使用卫星嵌入数据集进行回归

在 GitHub 上修改
报告问题
网页历史记录
本教程是有关卫星嵌入数据集的一系列教程中的一篇，另请参阅简介、非监督式分类、监督式分类和相似性搜索。

嵌入字段可像用于分类一样，用作回归的特征输入/预测变量。

在本教程中，我们将学习如何使用 64D 嵌入字段层作为输入，进行预测地上生物量 (AGB) 的多元回归分析。

NASA 的 Global Ecosystem Dynamics Investigation (GEDI) 任务沿地面样线以 30 米的空间分辨率每隔 60 米收集一次 LIDAR 测量数据。我们将使用 GEDI L4A 地上生物量密度栅格数据集，其中包含地上生物量密度 (AGBD) 的点估计值，这些值将用作回归模型中的预测变量。

选择区域

我们先来定义一个感兴趣的区域。在本教程中，我们将选择印度西高止山脉中的一个区域，并将一个多边形定义为几何变量。或者，您也可以使用代码编辑器中的绘图工具在感兴趣的区域周围绘制一个多边形，该多边形将作为几何变量保存在导入中。我们还使用卫星底图，以便轻松找到植被覆盖区域。

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');


图：选择地上生物量预测的感兴趣区域

选择时间段

选择要运行回归分析的年份。请注意，卫星嵌入是按年间隔进行汇总的，因此我们将整个年份定义为时间段。

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

准备卫星嵌入数据集

64 波段卫星嵌入图像将用作回归的预测变量。我们加载卫星嵌入数据集，并过滤出所选年份和区域的图片。

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

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

卫星嵌入图片会网格化为图块，并以图块的 UTM 区域投影形式提供。这样一来，我们便可获得覆盖感兴趣区域的多个卫星嵌入图块。为了获得单张图片，我们需要将它们拼接起来。在 Earth Engine 中，输入图像的镶嵌图会被分配默认投影，即比例为 1 度的 WGS84。由于我们将在本教程的后续部分中聚合并重新投影此镶嵌图，因此保留原始投影很有帮助。我们可以从其中一个图块中提取投影信息，并使用 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);

准备 GEDI L4A 镶嵌

由于 GEDI 生物量估计值将用于训练我们的回归模型，因此在使用之前过滤掉无效或不可靠的 GEDI 数据至关重要。我们应用了多个掩码来移除可能存在错误的测量结果。

  • 移除所有不符合质量要求的衡量结果（l4_quality_flag = 0 且 degrade_flag > 0）
  • 移除所有相对误差较高（“agbd_se”/“agbd”> 50%）的测量结果
  • 根据 Copernicus GLO-30 数字高程模式 (DEM) 移除坡度 > 30% 的所有测量结果

最后，我们选择感兴趣的时间段和区域的所有剩余测量结果，并创建一个镶嵌图。

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);

对输入进行重采样和聚合

在对像素进行抽样以训练回归模型之前，我们会对输入进行重采样和重新投影，使其位于同一像素网格中。GEDI 测量结果的水平精度为 +/- 9 米。在将 GEDI AGB 值与卫星嵌入像素匹配时，这会造成问题。为了克服这一问题，我们会对所有输入图片进行重采样和聚合，以获得更大的像素网格，其中包含原始像素的平均值。这也有助于去除数据中的噪声，并有助于构建更好的机器学习模型。

// 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));

重新投影和聚合像素是一项开销很大的操作，因此最好将生成的堆叠图像导出为资源，并在后续步骤中使用预先计算的图像。这有助于解决处理大型区域时出现的计算超时或用户内存超出错误。

// 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
});

启动导出任务，并等待其完成。完成后，我们导入资源并继续构建模型。

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

提取训练特征

我们已准备好输入数据，可以提取训练特征了。我们在回归模型中使用卫星嵌入频段作为因变量（预测变量），并使用 GEDI AGBD 值作为自变量（被预测变量）。我们可以提取每个像素处的重合值，并准备训练数据集。我们的 GEDI 图像大部分被遮盖，仅在少量像素处包含值。如果我们使用 sample()，则会返回大部分为空的值。为了克服这一问题，我们从 GEDI 掩码创建了一个类波段，并使用 stratifiedSample() 来确保我们从非掩码像素进行采样。

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());

训练回归模型

现在，我们可以训练模型了。Earth Engine 中的许多分类器都可用于分类和回归任务。由于我们想要预测一个数值（而非类别），因此可以将分类器设置为以 REGRESSION 模式运行，并使用训练数据进行训练。模型训练完成后，我们可以将模型的预测值与输入值进行比较，并计算均方根误差 (rmse) 和相关系数 r^2，以检查模型的性能。

// 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);


图：观测到的 AGBD 值与模型预测的 AGBD 值

针对未知值生成预测

对模型感到满意后，我们可以使用训练好的模型，根据包含预测变量波段的图片，在未知位置生成预测结果。

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

包含每个像素的预测 AGBD 值的图片现已准备好导出。我们将在下一部分中使用此数据来直观呈现结果。

// 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
});

启动导出任务，并等待其完成。完成后，我们导入资产并直观呈现结果。

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');


图：预测的地上生物量密度 (AGBD)。深绿色表示预测的生物质密度较高

估计总生物量

现在，我们获得了图像中每个像素的预测 AGBD 值，可用于估计该区域的地上总生物量 (AGB) 存量。但我们必须先移除所有属于非植被区域的像素。我们可以使用 ESA WorldCover 地表覆盖数据集并选择植被像素。

// 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)');

GEDI AGBD 值的单位为每公顷兆克 (Mg/ha)。如需获得总 AGB，我们将每个像素乘以其面积（以公顷为单位），然后将这些值相加。

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);

在 Earth Engine 代码编辑器中试用本教程的完整脚本。