Phân loại không giám sát bằng tập dữ liệu nhúng vệ tinh

Chỉnh sửa trên GitHub
Báo cáo vấn đề
Nhật ký trang
Tác giả: spatialthoughts
Hướng dẫn này nằm trong một loạt hướng dẫn về tập dữ liệu Nhúng vệ tinh. Bạn cũng có thể xem Giới thiệu, Phân loại có giám sát, Hồi quy và Tìm kiếm tương tự.

Trong hướng dẫn trước (Giới thiệu), chúng ta đã thấy cách Satellite Embeddings ghi lại quỹ đạo hằng năm của các quan sát từ vệ tinh và các biến khí hậu. Điều này giúp bạn dễ dàng sử dụng tập dữ liệu để lập bản đồ cây trồng mà không cần mô hình hoá kiểu hình cây trồng. Việc lập bản đồ loại cây trồng là một nhiệm vụ khó khăn, thường đòi hỏi phải mô hình hoá kiểu hình cây trồng và thu thập mẫu thực địa cho tất cả các loại cây trồng được trồng trong khu vực.

Trong hướng dẫn này, chúng ta sẽ sử dụng phương pháp phân loại không giám sát để lập bản đồ cây trồng, cho phép chúng ta thực hiện nhiệm vụ phức tạp này mà không cần dựa vào nhãn trường. Phương pháp này tận dụng kiến thức địa phương về khu vực cùng với số liệu thống kê tổng hợp về cây trồng, vốn có sẵn ở nhiều nơi trên thế giới.

Chọn một khu vực

Trong hướng dẫn này, chúng ta sẽ tạo một bản đồ loại cây trồng cho Quận Cerro Gordo, Iowa. Quận này nằm trong vùng trồng ngô của Hoa Kỳ, nơi có 2 loại cây trồng chính là ngô và đậu nành. Kiến thức địa phương này rất quan trọng và sẽ giúp chúng tôi quyết định các thông số chính cho mô hình của mình.

Hãy bắt đầu bằng cách lấy ranh giới cho quận đã chọn.

// Select the region
// Cerro Gordo County, Iowa
var counties = ee.FeatureCollection('TIGER/2018/Counties');

var selected = counties
  .filter(ee.Filter.eq('GEOID', '19033'));
var geometry = selected.geometry();
Map.centerObject(geometry, 12);
Map.addLayer(geometry, {color: 'red'}, 'Selected Region', false);


Hình: Khu vực đã chọn

Chuẩn bị tập dữ liệu Nhúng vệ tinh

Tiếp theo, chúng ta sẽ tải tập dữ liệu Satellite Embedding (Nhúng vệ tinh), lọc hình ảnh cho năm đã chọn và tạo một hình ảnh ghép.

var embeddings = ee.ImageCollection('GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL');

var year = 2022;
var startDate = ee.Date.fromYMD(year, 1, 1);
var endDate = startDate.advance(1, 'year');

var filteredembeddings = embeddings
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

var embeddingsImage = filteredembeddings.mosaic();

Tạo mặt nạ đất trồng trọt

Để mô hình hoá, chúng ta cần loại trừ các khu vực không phải là đất trồng trọt. Có nhiều tập dữ liệu toàn cầu và khu vực mà bạn có thể dùng để tạo mặt nạ cắt. ESA WorldCover hoặc GFSAD Global Cropland Extent Product là những lựa chọn phù hợp cho các tập dữ liệu toàn cầu về đất trồng. Một sản phẩm mới bổ sung gần đây là ESA WorldCereal Active Cropland (Đất trồng trọt đang hoạt động của ESA WorldCereal) có bản đồ theo mùa về đất trồng trọt đang hoạt động. Vì khu vực của chúng tôi nằm ở Hoa Kỳ, nên chúng tôi có thể sử dụng một tập dữ liệu khu vực chính xác hơn là Lớp dữ liệu về đất trồng trọt của USDA NASS (CDL) để lấy mặt nạ cây trồng.

// Use Cropland Data Layers (CDL) to obtain cultivated cropland
var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var croplandMask = cdl.select('cultivated').eq(2).rename('cropmask');

// Visualize the crop mask
var croplandMaskVis = {min: 0, max: 1, palette: ['white', 'green']};
Map.addLayer(croplandMask.clip(geometry), croplandMaskVis, 'Crop Mask');


Hình: Khu vực được chọn có mặt nạ đất trồng

Trích xuất mẫu huấn luyện

Chúng tôi áp dụng mặt nạ đất trồng trọt cho ảnh ghép nhúng. Giờ đây, chúng ta chỉ còn lại tất cả các điểm ảnh đại diện cho đất trồng trọt trong quận.

// Mask all non-cropland pixels
var clusterImage = embeddingsImage.updateMask(croplandMask);

Chúng ta cần lấy hình ảnh Nhúng vệ tinh và thu thập các mẫu ngẫu nhiên để huấn luyện một mô hình phân cụm. Vì khu vực mà chúng ta quan tâm có nhiều pixel bị che, nên việc lấy mẫu ngẫu nhiên đơn giản có thể dẫn đến các mẫu có giá trị rỗng. Để đảm bảo có thể trích xuất số lượng mẫu không có giá trị rỗng mong muốn, chúng tôi sử dụng phương pháp lấy mẫu phân tầng để thu được số lượng mẫu mong muốn ở những khu vực không bị che khuất.

// Stratified random sampling
var training = clusterImage.addBands(croplandMask).stratifiedSample({
  numPoints: 1000,
  classBand: 'cropmask',
  region: geometry,
  scale: 10,
  tileScale: 16,
  seed: 100,
  dropNulls: true,
  geometries: true
});

Xuất mẫu sang một thành phần (không bắt buộc)

Việc trích xuất mẫu là một thao tác tốn nhiều tài nguyên tính toán và bạn nên xuất các mẫu huấn luyện đã trích xuất dưới dạng Tài sản rồi sử dụng các tài sản đã xuất trong các bước tiếp theo. Điều này sẽ giúp khắc phục lỗi tính toán hết thời gian chờ hoặc lỗi vượt quá bộ nhớ người dùng khi làm việc với các khu vực rộng lớn.

Bắt đầu tác vụ xuất và đợi tác vụ này hoàn tất rồi mới tiếp tục.

// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';

var samplesExportFc = 'cluster_training_samples';
var samplesExportFcPath = exportFolder + samplesExportFc;

Export.table.toAsset({
  collection: training,
  description: 'Cluster_Training_Samples',
  assetId: samplesExportFcPath
});

Sau khi hoàn tất tác vụ xuất, chúng ta có thể đọc các mẫu đã trích xuất trở lại mã của mình dưới dạng một tập hợp đối tượng.

// Use the exported asset
var training = ee.FeatureCollection(samplesExportFcPath);

Hiển thị các mẫu một cách trực quan

Cho dù chạy mẫu một cách tương tác hay xuất sang một tập hợp đối tượng, giờ đây, bạn sẽ có một biến huấn luyện với các điểm mẫu. Hãy in mẫu đầu tiên để kiểm tra và thêm các điểm huấn luyện vào Map.

print('Extracted sample', training.first());
Map.addLayer(training, {color: 'blue'}, 'Extracted Samples', false);


Hình: Các mẫu ngẫu nhiên được trích xuất để phân cụm

Thực hiện phân cụm không giám sát

Giờ đây, chúng ta có thể huấn luyện một trình phân cụm và nhóm các vectơ nhúng 64D thành một số lượng cụm riêng biệt đã chọn. Dựa trên kiến thức địa phương, có 2 loại cây trồng chính chiếm phần lớn diện tích, cùng với một số loại cây trồng khác chiếm phần còn lại. Chúng tôi có thể thực hiện phân cụm không giám sát trên Satellite Embedding để thu được các cụm pixel có quỹ đạo và mẫu thời gian tương tự. Các pixel có đặc điểm quang phổ và không gian tương tự cùng với kiểu hình tương tự sẽ được nhóm vào cùng một cụm.

ee.Clusterer.wekaCascadeKMeans() cho phép chúng ta chỉ định số lượng cụm tối thiểu và tối đa, đồng thời tìm ra số lượng cụm tối ưu dựa trên dữ liệu huấn luyện. Lúc này, kiến thức địa phương của chúng ta sẽ hữu ích trong việc quyết định số lượng tối thiểu và tối đa của các cụm. Vì chúng tôi dự kiến sẽ có một số loại cụm riêng biệt (ngô, đậu nành và một số loại cây trồng khác), nên chúng tôi có thể sử dụng 4 làm số lượng cụm tối thiểu và 5 làm số lượng cụm tối đa. Bạn có thể phải thử nghiệm với những con số này để tìm ra con số phù hợp nhất với khu vực của mình.

// Cluster the Satellite Embedding Image
var minClusters = 4;
var maxClusters = 5;

var clusterer = ee.Clusterer.wekaCascadeKMeans({
  minClusters: minClusters, maxClusters: maxClusters}).train({
  features: training,
  inputProperties: clusterImage.bandNames()
});

var clustered = clusterImage.cluster(clusterer);
Map.addLayer(clustered.randomVisualizer().clip(geometry), {}, 'Clusters');


Hình: Các cụm thu được từ quá trình phân loại không giám sát

Chỉ định nhãn cho các cụm

Sau khi kiểm tra bằng mắt, các cụm thu được ở các bước trước đó gần như trùng khớp với ranh giới trang trại trong hình ảnh có độ phân giải cao. Dựa trên kiến thức địa phương, chúng tôi biết rằng hai cụm lớn nhất sẽ là ngô và đậu nành. Hãy tính diện tích của từng cụm trong hình ảnh.

// Calculate Cluster Areas
// 1 Acre = 4046.86 Sq. Meters
var areaImage = ee.Image.pixelArea().divide(4046.86).addBands(clustered);

var areas = areaImage.reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'cluster',
    }),
    geometry: geometry,
    scale: 10,
    maxPixels: 1e10
    });

var clusterAreas = ee.List(areas.get('groups'));

// Process results to extract the areas and create a FeatureCollection

var clusterAreas = clusterAreas.map(function(item) {
  var areaDict = ee.Dictionary(item);
  var clusterNumber = areaDict.getNumber('cluster').format();
  var area = areaDict.getNumber('sum')
  return ee.Feature(null, {cluster: clusterNumber, area: area})
})

var clusterAreaFc = ee.FeatureCollection(clusterAreas);
print('Cluster Areas', clusterAreaFc);

Chúng tôi chọn 2 cụm có diện tích lớn nhất.

var selectedFc = clusterAreaFc.sort('area', false).limit(2);
print('Top 2 Clusters by Area', selectedFc);

Nhưng chúng ta vẫn chưa biết cụm nào là loại cây trồng nào. Nếu có một vài mẫu ruộng ngô hoặc đậu nành, bạn có thể phủ chúng lên các cụm để xác định nhãn tương ứng. Nếu không có mẫu thực địa, chúng ta có thể tận dụng số liệu thống kê tổng hợp về cây trồng. Ở nhiều nơi trên thế giới, số liệu thống kê tổng hợp về cây trồng được thu thập và xuất bản thường xuyên. Đối với Hoa Kỳ, Cục Thống kê Nông nghiệp Quốc gia (NASS) có số liệu thống kê chi tiết về cây trồng cho từng quận và từng loại cây trồng chính. Năm 2022, Quận Cerro Gordo, Iowa có Diện tích trồng ngô là 161.500 mẫu Anh và Diện tích trồng đậu nành là 110.500 mẫu Anh.

Dựa vào thông tin này, giờ đây chúng ta biết rằng trong số 2 cụm hàng đầu, cụm có diện tích lớn nhất có nhiều khả năng là ngô và cụm còn lại là đậu nành. Hãy chỉ định các nhãn này và so sánh các vùng được tính toán với số liệu thống kê đã xuất bản.

var cornFeature = selectedFc.sort('area', false).first();
var soybeanFeature = selectedFc.sort('area').first();
var cornCluster = cornFeature.get('cluster');
var soybeanCluster = soybeanFeature.get('cluster');

print('Corn Area (Detected)', cornFeature.getNumber('area').round());
print('Corn Area (From Crop Statistics)', 163500);

print('Soybean Area (Detected)', soybeanFeature.getNumber('area').round());
print('Soybean Area (From Crop Statistics)', 110500);

Tạo bản đồ cắt

Giờ đây, chúng ta đã biết nhãn cho từng cụm và có thể trích xuất các pixel cho từng loại cây trồng rồi hợp nhất chúng để tạo bản đồ cây trồng cuối cùng.

// Select the clusters to create the crop map
var corn = clustered.eq(ee.Number.parse(cornCluster));
var soybean = clustered.eq(ee.Number.parse(soybeanCluster));

var merged = corn.add(soybean.multiply(2));
var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(merged.clip(geometry), cropVis, 'Crop Map (Detected)');

Để giúp diễn giải kết quả, chúng ta cũng có thể sử dụng các phần tử trên giao diện người dùng để tạo và thêm chú thích vào bản đồ.

// Add a Legend
var legend = ui.Panel({
  layout: ui.Panel.Layout.Flow('horizontal'),
  style: {position: 'bottom-center', padding: '8px 15px'}});

var addItem = function(color, name) {
  var colorBox = ui.Label({
    style: {color: '#ffffff',
      backgroundColor: color,
      padding: '10px',
      margin: '0 4px 4px 0',
    }
  });
  var description = ui.Label({
    value: name,
    style: {
      margin: '0px 10px 0px 2px',
    }
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')}
)};

var title = ui.Label({
  value: 'Legend',
  style: {fontWeight: 'bold',
    fontSize: '16px',
    margin: '0px 10px 0px 4px'}});

legend.add(title);
legend.add(addItem('#ffd400', 'Corn'));
legend.add(addItem('#267300', 'Soybean'));
legend.add(addItem('#bdbdbd', 'Other Crops'));


Hình 1: Bản đồ cây trồng được phát hiện có cây ngô và đậu nành

Xác thực kết quả

Chúng tôi có thể thu thập bản đồ loại cây trồng bằng tập dữ liệu Satellite Embedding mà không cần nhãn trường chỉ bằng cách sử dụng số liệu thống kê tổng hợp và kiến thức địa phương về khu vực. Hãy so sánh kết quả của chúng ta với bản đồ chính thức về loại cây trồng của Lớp dữ liệu về đất trồng (CDL) của NASS thuộc Bộ Nông nghiệp Hoa Kỳ (USDA).

var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var cropMap = cropLandcover.updateMask(croplandMask).rename('crops');

// Original data has unique values for each crop ranging from 0 to 254
var cropClasses = ee.List.sequence(0, 254);
// We remap all values as following
// Crop     | Source Value | Target Value
// Corn     | 1            | 1
// Soybean  | 5            | 2
// All other| 0-255        | 0
var targetClasses = ee.List.repeat(0, 255).set(1, 1).set(5, 2);
var cropMapReclass = cropMap.remap(cropClasses, targetClasses).rename('crops');

var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(cropMapReclass.clip(geometry), cropVis, 'Crop Landcover (CDL)');

Mặc dù có sự khác biệt giữa kết quả của chúng tôi và bản đồ chính thức, nhưng bạn sẽ nhận thấy rằng chúng tôi có thể thu được kết quả khá tốt mà không tốn nhiều công sức. Bằng cách áp dụng các bước xử lý hậu kỳ cho kết quả, chúng ta có thể loại bỏ một số vùng nhiễu và lấp đầy các khoảng trống trong đầu ra.

Hãy thử toàn bộ tập lệnh cho hướng dẫn này trong Trình chỉnh sửa mã Earth Engine.