Hồi quy 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 không giám sát, Phân loại có giám sát và Tìm kiếm tương tự.

Bạn có thể dùng các trường nhúng làm đầu vào/yếu tố dự đoán cho hồi quy theo cách tương tự như cách bạn dùng cho phân loại.

Trong hướng dẫn này, chúng ta sẽ tìm hiểu cách sử dụng 64 lớp trường nhúng 3D làm dữ liệu đầu vào cho một phân tích hồi quy đa biến để dự đoán sinh khối trên mặt đất (AGB).

Sứ mệnh Điều tra động lực học hệ sinh thái toàn cầu (GEDI) của NASA thu thập các phép đo LIDAR dọc theo các đường cắt ngang trên mặt đất với độ phân giải không gian 30 m ở khoảng thời gian 60 m. Chúng ta sẽ sử dụng tập dữ liệu GEDI L4A Raster Aboveground Biomass Density chứa các điểm ước tính về mật độ sinh khối trên mặt đất (AGBD) sẽ được dùng làm biến dự đoán trong mô hình hồi quy.

Chọn một khu vực

Hãy bắt đầu bằng cách xác định một khu vực quan tâm. Trong hướng dẫn này, chúng ta sẽ chọn một khu vực ở Western Ghats của Ấn Độ và xác định một đa giác làm biến hình học. Ngoài ra, bạn có thể sử dụng Công cụ vẽ trong Trình chỉnh sửa mã để vẽ một đa giác xung quanh khu vực mà bạn quan tâm. Đa giác này sẽ được lưu dưới dạng biến hình học trong phần Nhập. Chúng tôi cũng sử dụng lớp cơ sở Vệ tinh để dễ dàng xác định vị trí của các khu vực có cây cối.

var geometry = ee.Geometry.Polygon([[
  [74.322, 14.981],
  [74.322, 14.765],
  [74.648, 14.765],
  [74.648, 14.980]
]]);

// Use the satellite basemap
Map.setOptions('SATELLITE');


Hình: Chọn khu vực quan tâm để dự đoán sinh khối trên mặt đất

Chọn một khoảng thời gian

Chọn một năm mà chúng ta muốn chạy hồi quy. Xin lưu ý rằng Satellite Embeddings được tổng hợp theo khoảng thời gian hằng năm, vì vậy, chúng ta sẽ xác định khoảng thời gian cho cả năm.

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

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

Hình ảnh nhúng vệ tinh 64 băng tần sẽ được dùng làm yếu tố dự đoán cho hồi quy. Chúng tôi tải tập dữ liệu Nhúng vệ tinh, lọc hình ảnh cho năm và khu vực đã chọn.

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

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

Hình ảnh Nhúng vệ tinh được chia thành các ô và phân phát trong phép chiếu cho các vùng UTM của ô. Do đó, chúng ta sẽ có nhiều ô Nhúng vệ tinh bao phủ khu vực mà chúng ta quan tâm. Để có được một hình ảnh duy nhất, chúng ta cần ghép các hình ảnh đó lại với nhau. Trong Earth Engine, một tổ hợp hình ảnh đầu vào được chỉ định phép chiếu mặc định, đó là WGS84 với tỷ lệ 1 độ. Vì chúng ta sẽ tổng hợp và chiếu lại bức ảnh ghép này ở phần sau của hướng dẫn, nên việc giữ lại phép chiếu ban đầu là điều hữu ích. Chúng ta có thể trích xuất thông tin phép chiếu từ một trong các ô và đặt thông tin đó trên ảnh ghép bằng hàm setDefaultProjection().

// Extract the projection of the first band of the first image
var embeddingsProjection = ee.Image(embeddingsFiltered.first()).select(0).projection();

// Set the projection of the mosaic to the extracted projection
var embeddingsImage = embeddingsFiltered.mosaic()
  .setDefaultProjection(embeddingsProjection);

Chuẩn bị ảnh ghép GEDI L4A

Vì các số liệu ước tính về sinh khối của GEDI sẽ được dùng để huấn luyện mô hình hồi quy của chúng tôi, nên việc lọc ra dữ liệu không hợp lệ hoặc không đáng tin cậy của GEDI trước khi sử dụng là rất quan trọng. Chúng tôi áp dụng một số mặt nạ để loại bỏ những phép đo có thể bị lỗi.

  • Xoá tất cả các phép đo không đáp ứng yêu cầu về chất lượng (l4_quality_flag = 0 và degrade_flag > 0)
  • Xoá tất cả các phép đo có sai số tương đối cao ("agbd_se" / "agbd" > 50%)
  • Xoá tất cả các phép đo trên sườn dốc > 30% dựa trên Mô hình độ cao kỹ thuật số (DEM) GLO-30 của Copernicus

Cuối cùng, chúng tôi chọn tất cả các phép đo còn lại cho khoảng thời gian và khu vực quan tâm rồi tạo ra một bức ảnh ghép.

var gedi = ee.ImageCollection('LARSE/GEDI/GEDI04_A_002_MONTHLY');
// Function to select the highest quality GEDI data
var qualityMask = function(image) {
  return image.updateMask(image.select('l4_quality_flag').eq(1))
      .updateMask(image.select('degrade_flag').eq(0));
};

// Function to mask unreliable GEDI measurements
// with a relative standard error > 50%
// agbd_se / agbd > 0.5
var errorMask = function(image) {
  var relative_se = image.select('agbd_se')
    .divide(image.select('agbd'));
  return image.updateMask(relative_se.lte(0.5));
};

// Function to mask GEDI measurements on slopes > 30%

var slopeMask = function(image) {
  // Use Copernicus GLO-30 DEM for calculating slope
  var glo30 = ee.ImageCollection('COPERNICUS/DEM/GLO30');

  var glo30Filtered = glo30
    .filter(ee.Filter.bounds(geometry))
    .select('DEM');

  // Extract the projection
  var demProj = glo30Filtered.first().select(0).projection();

  // The dataset consists of individual images
  // Create a mosaic and set the projection
  var elevation = glo30Filtered.mosaic().rename('dem')
    .setDefaultProjection(demProj);

  // Compute the slope
  var slope = ee.Terrain.slope(elevation);

  return image.updateMask(slope.lt(30));
};

var gediFiltered = gedi
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

var gediProjection = ee.Image(gediFiltered.first())
  .select('agbd').projection();

var gediProcessed = gediFiltered
  .map(qualityMask)
  .map(errorMask)
  .map(slopeMask);

var gediMosaic = gediProcessed.mosaic()
  .select('agbd').setDefaultProjection(gediProjection);

// Visualize the GEDI Mosaic
var gediVis = {
  min: 0,
  max: 200,
  palette: ['#edf8fb','#b2e2e2','#66c2a4','#2ca25f','#006d2c'],
  bands: ['agbd']
};

Map.addLayer(gediMosaic, gediVis, 'GEDI L4A (Filtered)', false);

Lấy mẫu lại và tổng hợp dữ liệu đầu vào

Trước khi lấy mẫu các pixel để huấn luyện mô hình hồi quy, chúng tôi lấy mẫu lại và chiếu lại các đầu vào vào cùng một lưới pixel. Các phép đo GEDI có độ chính xác theo chiều ngang là +/- 9 m. Điều này gây ra vấn đề khi so khớp các giá trị GEDI AGB với các pixel Nhúng vệ tinh. Để khắc phục điều này, chúng tôi lấy mẫu lại và tổng hợp tất cả hình ảnh đầu vào thành một lưới pixel lớn hơn với các giá trị trung bình từ các pixel ban đầu. Điều này cũng giúp loại bỏ nhiễu khỏi dữ liệu và giúp xây dựng một mô hình học máy hiệu quả hơn.

// Choose the grid size and projection
var gridScale = 100;
var gridProjection = ee.Projection('EPSG:3857')
  .atScale(gridScale);

// Create a stacked image with predictor and predicted variables
var stacked = embeddingsImage.addBands(gediMosaic);

//  Set the resampling mode
var stacked = stacked.resample('bilinear');

// Aggregate pixels with 'mean' statistics
var stackedResampled = stacked
  .reduceResolution({
    reducer: ee.Reducer.mean(),
    maxPixels: 1024
  })
  .reproject({
    crs: gridProjection
});

// As larger GEDI pixels contain masked original
// pixels, it has a transparency mask.
// We update the mask to remove the transparency
var stackedResampled = stackedResampled
  .updateMask(stackedResampled.mask().gt(0));

Việc chiếu lại và tổng hợp các pixel là một thao tác tốn kém, vì vậy, bạn nên xuất hình ảnh xếp chồng kết quả dưới dạng một Tài sản và sử dụng hình ảnh được tính toán trước trong các bước tiếp theo. Điều này sẽ giúp khắc phục lỗi computation timed out (hết thời gian tính toán) hoặc user memory exceeded (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.

// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';
var mosaicExportImage = 'gedi_mosaic';
var mosaicExportImagePath = exportFolder + mosaicExportImage;
Export.image.toAsset({
  image: stackedResampled.clip(geometry),
  description: 'GEDI_Mosaic_Export',
  assetId: mosaicExportImagePath,
  region: geometry,
  scale: gridScale,
  maxPixels: 1e10
});

Bắt đầu tác vụ xuất và đợi cho đến khi tác vụ này hoàn tất. Sau khi hoàn tất, chúng ta sẽ nhập Tài sản và tiếp tục xây dựng mô hình.

// Use the exported asset
var stackedResampled = ee.Image(mosaicExportImagePath);

Trích xuất các đặc điểm huấn luyện

Chúng ta đã chuẩn bị sẵn dữ liệu đầu vào để trích xuất các đặc điểm huấn luyện. Chúng tôi sử dụng các dải nhúng vệ tinh làm biến phụ thuộc (biến dự đoán) và các giá trị AGBD của GEDI làm biến độc lập (biến được dự đoán) trong mô hình hồi quy. Chúng ta có thể trích xuất các giá trị trùng khớp tại mỗi pixel và chuẩn bị tập dữ liệu huấn luyện. Hình ảnh GEDI của chúng tôi hầu hết được che phủ và chỉ chứa các giá trị ở một số ít pixel. Nếu chúng ta sử dụng sample(), hàm này sẽ trả về hầu hết các giá trị trống. Để khắc phục vấn đề này, chúng ta tạo một dải lớp từ mặt nạ GEDI và sử dụng stratifiedSample() để đảm bảo chúng ta lấy mẫu từ các pixel không được che.

var predictors = embeddingsImage.bandNames();
var predicted = gediMosaic.bandNames().get(0);
print('predictors', predictors);
print('predicted', predicted);

var predictorImage = stackedResampled.select(predictors);
var predictedImage = stackedResampled.select([predicted]);

var classMask = predictedImage.mask().toInt().rename('class');

var numSamples = 1000;

// We set classPoints to [0, numSamples]
// This will give us 0 points for class 0 (masked areas)
// and numSample points for class 1 (non-masked areas)
var training = stackedResampled.addBands(classMask)
  .stratifiedSample({
    numPoints: numSamples,
    classBand: 'class',
    region: geometry,
    scale: gridScale,
    classValues: [0, 1],
    classPoints: [0, numSamples],
    dropNulls: true,
    tileScale: 16,
});

print('Number of Features Extracted', training.size());
print('Sample Training Feature', training.first());

Đào tạo mô hình hồi quy

Giờ đây, chúng ta đã sẵn sàng huấn luyện mô hình. Bạn có thể sử dụng nhiều trình phân loại trong Earth Engine cho cả nhiệm vụ phân loại và hồi quy. Vì muốn dự đoán một giá trị bằng số (thay vì một lớp) – chúng ta có thể đặt trình phân loại để chạy ở chế độ REGRESSION và huấn luyện bằng dữ liệu huấn luyện. Sau khi huấn luyện mô hình, chúng ta có thể so sánh kết quả dự đoán của mô hình với các giá trị đầu vào và tính toán sai số bình phương căn trung bình (rmse) và hệ số tương quan r^2 để kiểm tra hiệu suất của mô hình.

// Use the RandomForest classifier and set the
// output mode to REGRESSION
var model = ee.Classifier.smileRandomForest(50)
  .setOutputMode('REGRESSION')
  .train({
    features: training,
    classProperty: predicted,
    inputProperties: predictors
  });

// Get model's predictions for training samples
var predicted = training.classify({
  classifier: model,
  outputName: 'agbd_predicted'
});

// Calculate RMSE
var calculateRmse = function(input) {
    var observed = ee.Array(
      input.aggregate_array('agbd'));
    var predicted = ee.Array(
      input.aggregate_array('agbd_predicted'));
    var rmse = observed.subtract(predicted).pow(2)
      .reduce('mean', [0]).sqrt().get([0]);
    return rmse;
};
var rmse = calculateRmse(predicted);
print('RMSE', rmse);

// Create a plot of observed vs. predicted values
var chart = ui.Chart.feature.byFeature({
  features: predicted.select(['agbd', 'agbd_predicted']),
  xProperty: 'agbd',
  yProperties: ['agbd_predicted'],
}).setChartType('ScatterChart')
  .setOptions({
    title: 'Aboveground Biomass Density (Mg/Ha)',
    dataOpacity: 0.8,
    hAxis: {'title': 'Observed'},
    vAxis: {'title': 'Predicted'},
    legend: {position: 'right'},
    series: {
      0: {
        visibleInLegend: false,
        color: '#525252',
        pointSize: 3,
        pointShape: 'triangle',
      },
    },
    trendlines: {
      0: {
        type: 'linear',
        color: 'black',
        lineWidth: 1,
        pointSize: 0,
        labelInLegend: 'Linear Fit',
        visibleInLegend: true,
        showR2: true
      }
    },
    chartArea: {left: 100, bottom: 100, width: '50%'},

});
print(chart);


Hình: Giá trị AGBD được ghi nhận so với giá trị AGBD dự đoán của mô hình

Tạo giá trị dự đoán cho các giá trị không xác định

Sau khi hài lòng với mô hình này, chúng ta có thể sử dụng mô hình đã huấn luyện để tạo ra các dự đoán tại những vị trí chưa biết từ hình ảnh có chứa các dải dự đoán.

// We set the band name of the output image as 'agbd'
var predictedImage = stackedResampled.classify({
  classifier: model,
  outputName: 'agbd'
});

Giờ đây, hình ảnh chứa các giá trị AGBD được dự đoán tại mỗi pixel đã sẵn sàng để xuất. Chúng ta sẽ sử dụng thông tin này trong phần tiếp theo để hình dung kết quả.

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

Export.image.toAsset({
  image: predictedImage.clip(geometry),
  description: 'Predicted_Image_Export',
  assetId: predictedExportImagePath,
  region: geometry,
  scale: gridScale,
  maxPixels: 1e10
});

Bắt đầu tác vụ xuất và đợi cho đến khi tác vụ này hoàn tất. Sau khi hoàn tất, chúng ta sẽ nhập Tài sản và trực quan hoá kết quả.

var predictedImage = ee.Image(predictedExportImagePath);

// Visualize the image
var gediVis = {
  min: 0,
  max: 200,
  palette: ['#edf8fb','#b2e2e2','#66c2a4','#2ca25f','#006d2c'],
  bands: ['agbd']
};

Map.addLayer(predictedImage, gediVis, 'Predicted AGBD');


Hình: AGBD dự đoán. Màu xanh lục đậm hơn cho biết mật độ sinh khối dự kiến cao hơn

Ước tính tổng sinh khối

Giờ đây, chúng ta đã có các giá trị AGBD dự đoán cho từng pixel của hình ảnh và có thể dùng các giá trị này để ước tính tổng lượng sinh khối trên mặt đất (AGB) trong khu vực. Nhưng trước tiên, chúng ta phải xoá tất cả các pixel thuộc những khu vực không có thảm thực vật. Chúng ta có thể sử dụng tập dữ liệu về độ che phủ mặt đất ESA WorldCover và chọn các pixel có thực vật.

// GEDI data is processed only for certain landcovers
// from Plant Functional Types (PFT) classification
// https://doi.org/10.1029/2022EA002516

// Here we use ESA WorldCover v200 product to
// select landcovers representing vegetated areas
var worldcover = ee.ImageCollection('ESA/WorldCover/v200').first();

// Aggregate pixels to the same grid as other dataset
// with 'mode' value.
// i.e. The landcover with highest occurrence within the grid
var worldcoverResampled = worldcover
  .reduceResolution({
    reducer: ee.Reducer.mode(),
    maxPixels: 1024
  })
  .reproject({
    crs: gridProjection
});

// Select grids for the following classes
// | Class Name | Value |
// | Forests    | 10    |
// | Shrubland  | 20    |
// | Grassland  | 30    |
// | Cropland   | 40    |
// | Mangroves  | 95    |
var landCoverMask = worldcoverResampled.eq(10)
    .or(worldcoverResampled.eq(20))
    .or(worldcoverResampled.eq(30))
    .or(worldcoverResampled.eq(40))
    .or(worldcoverResampled.eq(95));

var predictedImageMasked = predictedImage
  .updateMask(landCoverMask);
Map.addLayer(predictedImageMasked, gediVis, 'Predicted AGBD (Masked)');

Đơn vị của giá trị AGBD GEDI là megagram trên hécta (Mg/ha). Để tính tổng AGB, chúng tôi nhân từng pixel với diện tích của pixel đó tính bằng héc ta rồi cộng các giá trị lại với nhau.

var pixelAreaHa = ee.Image.pixelArea().divide(10000);
var predictedAgb = predictedImageMasked.multiply(pixelAreaHa);

var stats = predictedAgb.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: geometry,
  scale: gridScale,
  maxPixels: 1e10,
  tileScale: 16
});

// Result is a dictionary with key for each band
var totalAgb = stats.getNumber('agbd');

print('Total AGB (Mg)', totalAgb);

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.