前言
单纯依赖样地估测法或机载激光雷达测量方法进行某个县、市的森林植被碳储量估算，仍面临时间过长的问题。
在区域尺度上进行森林植被碳储量的反演，许多研究已经开始利用多光谱卫星影像提取植被的光谱特征，并基于一些植被指数与实测样方的生物量进行建模，以估算区域尺度的碳储量。随着GEDI和GLAS激光雷达数据的公开，激光雷达数据能够提取冠层高度这一重要的结构性指标，但由于其分辨率较低，通常需要结合卫星影像才能进行区域尺度上森林冠层高度的连续分布反演，进而实现森林植被碳储量或生物量估算的反演。
Aboveground-carbon-stock-distribution-map-for-secondary-Andean-forests-in-the-Tabio-test_Q320.jpg
本文将基于Google Earth Engine（GEE）平台，利用Sentinel数据与GEDI数据，介绍讨如何进行森林碳储量的反演。
反演思路：
- 利用GEDI数据结合Landsat或Sentinel光学影像数据，生成区域尺度上连续分布的森林冠层高度模型；
- 获取目标县市内树种分布的矢量数据；若没有此类数据，可使用已公开的土地覆盖类型数据；
- 对不同森林类型的冠层高度与相应的实测碳储量或生物量数据进行建模，最终通过训练好的模型反演整个区域的森林植被碳储量。
实际操作：
1、森林冠层高度反演
代码思路：
- 提取光学影像中的重要指数和波段特征，常见的处理方法可参考相关文献。
- 对GEDI数据进行质量筛选，并考虑坡度等因素进行进一步的筛选。
- 因为GEDI光斑存在位置偏移，需对其进行重采样到更粗的分辨率。有一篇文章为了提高分辨率，筛选光斑中心位置与影像像元中心位置很近的光斑，以减小位置偏移的影响，并将其仅重采样至30m分辨率。并将GEDI提取的冠层高度指数与影像计算的指数进行图层叠加
- 根据筛选后的光斑位置进行样本提取，即提取改光斑位置对应的树高指数和光谱指数，用来构建森林冠层高度的反演模型。
Map.centerObject(boundary,9)
var startDate=ee.Date.fromYMD(2021,4,1)
var endDate=ee.Date.fromYMD(2022,10,30)
// 加载Sentinel-2光谱反射率数据
//这一部分是为了合成光学影像提取植被指数和光谱属性
// sentinal_2 data mask clouds
var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED");
var s2Filter=s2.filterDate(startDate,endDate)
.filterBounds(boundary)
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
var s2projection=s2Filter.first().select("B3").projection()
function maskS2clouds(image) {
var qa = image.select('QA60');
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
.and(qa.bitwiseAnd(cirrusBitMask).eq(0));
return image.updateMask(mask).divide(10000)
;}
// Function to compute spectral indices
var addIndices = function(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('ndvi');
var evi = image.expression(
'2.5 * ((NIR - RED)/(NIR + 6*RED - 7.5*BLUE + 1))', {
'NIR': image.select('B8'),
'RED': image.select('B4'),
'BLUE': image.select('B2')}).rename('evi');
return image
.addBands(ndvi)
.addBands(evi);};
// remove cloudy and then add indices
var S2Processed = s2Filter
.map(maskS2clouds)
.map(addIndices);
// 坐标转换
var S2_composite = S2Processed.median().reproject({crs: 'EPSG:4326', scale: 100})
print('Projection, crs, and crs_transform:', S2_composite.projection());
Map.addLayer(S2_composite.clip(boundary), {bands: ['B11', 'B8', 'B3'], min: 0, max: 0.3});
// 这一部分加载DEM数据，并计算高程，坡度和坡向
// 加载 SRTM
var SRTM = ee.Image("USGS/SRTMGL1_003");
//裁剪研究区
var elevation = SRTM.clip(boundary);
//坐标转换
var elevation = elevation.reproject({crs: 'EPSG:4326',scale: 100});
print('Projection, crs, and crs_transform:', elevation.projection());
// 提取坡度
var slope = ee.Terrain.slope(SRTM).clip(boundary);
var slope = slope.reproject({crs: 'EPSG:4326',scale: 100});
print('Projection, crs, and crs_transform:', slope.projection());
// 加载土地覆盖类型数据，筛选出森林类别
森林
// 加载ESA World Cover数据
var dataset = ee.ImageCollection("ESA/WorldCover/v100").first();
// 裁剪研究区
var ESA_LC_2020 = dataset.clip(boundary);
// 从土地覆盖中提取森林区域
var forest_mask = ESA_LC_2020.updateMask(ESA_LC_2020.eq(10) );
// 显示森林
var visualization = {bands: ['Map'],};
Map.addLayer(forest_mask, visualization, "Trees");
// 合并预测变量，并导入GEDI中的RH95变量。
// 合并预测变量
var mergedCollection = S2_composite
.addBands(elevation
.addBands(slope));
// 准备训练数据集
// 定义mask
var qualityMask = function(im) {
return im.updateMask(im.select('quality_flag').eq(1))
.updateMask(im.select('degrade_flag').eq(0));
};
// 导入数据。筛选gedi光斑数据
var dataset = ee.ImageCollection('LARSE/GEDI/GEDI02_A_002_MONTHLY')
.filterDate(startDate, endDate)
.map(qualityMask)
.select('rh95').filterBounds(boundary);
var gediMosaic = dataset.mosaic().clip(boundary)
光斑
// var clippedmergedCollection=mergedCollection.addBands(gediMosaic)
// .clip(boundary)
var stacked=mergedCollection.addBands(gediMosaic).reproject({crs: 'EPSG:4326',scale: 100})
var predictors = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8',
'B11', 'B12','ndvi','evi', 'elevation', 'slope'];
var predicted =['rh95']
var predictorImage = stacked.select(predictors);
var predictedImage = stacked.select(predicted);
var classMask = predictedImage.mask().toInt().gt(0).rename('class').clip(boundary);
//有值为1无值为0,加上.gt(0)才会有问题
print(classMask);
// print(predictedImage);
// print(predictedImage.mask().toInt());
Map.addLayer(classMask, {min: 0, max: 1}, 'Class Mask');
var numSamples = 1000
//分层抽样，class必须为整数
var training = stacked.addBands(classMask)
.stratifiedSample({
numPoints: numSamples,
classBand: 'class',
region: boundary,
scale: 100,
classValues: [0,1],
classPoints: [0, numSamples],
dropNulls: false,//这里我选择false，否则失败，采样为0
tileScale: 16
});
// output mode to REGRESSION
var model = ee.Classifier.smileRandomForest(50)
.setOutputMode('REGRESSION')
.train({
features: training,
classProperty: 'rh95',
inputProperties: predictors
});
// Get model predictions for training samples
var predicted = training.classify({
classifier: model,
outputName: 'height_predicted'
});
// Calculate RMSE
var calculateRmse = function(input) {
var observed = ee.Array(
input.aggregate_array("rh95"));
var predicted = ee.Array(
input.aggregate_array('height_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(['rh95', 'height_predicted']),
xProperty: 'rh95',
yProperties: ['height_predicted'],
}).setChartType('ScatterChart')
.setOptions({
title: 'Tree height (m)',
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);
var predictedImage = stacked.classify({
classifier: model,
outputName: 'height'
});
var exportPath = 'users/foundation/gedi_map/';
Export.image.toAsset({
image: predictedImage.updateMask(forest_mask).clip(boundary),
description: 'Predicted_Image_Export',
assetId: exportPath + 'predicted_height',
region: boundary,
scale: 100,
maxPixels: 1e10
});
// // 导出图像，指定比例和区域
// Export.image.toDrive({
// image: predictedImage.updateMask(forest_mask).clip(geometry),
// description: 'height',
// scale: 100,
// crs: 'EPSG:4326',
// maxPixels: 6756353855,
// region: boundary
// });
结果：我这里将结果导入到asset里，方便计算生物量时调用
2、估测森林生物量
导入生物量的实测数据，建模估测，这部分代码与上面估测冠层高度类似。
var col = ee.FeatureCollection(table); // 实测生物量数据
var chm = image; // 冠层高度模型数据
Map.centerObject(col, 9);
// 1.提取CHM值与生物量数据
// 使用`sample`方法从CHM图像中提取每个样地中心点的CHM值，并与生物量数据匹配
var samples = chm.sampleRegions({
collection: col,
properties: ['carbon'], // 假设生物量数据存储在'biomass'属性中
scale: 100 // CHM的分辨率
});
// 2.拆分训练和验证数据
var split = 0.8; // 80%数据用于训练，20%用于验证
var training = samples.filter(ee.Filter.lt('system:index', split));
var validation = samples.filter(ee.Filter.gte('system:index', split));
// 3. 构建回归模型
var model = ee.Classifier.smileLinear().train({
features: training,
classProperty: 'carbon',
inputProperties: ['CHM']
});
// 4.用回归模型进行预测，用训练好的模型预测生物量
var carbonPrediction = chm.classify(model);
Map.addLayer(carbonPrediction, {min: 0, max: 100}, 'Predicted carbon');
// 5. 使用验证数据评估模型
var validationPredictions = validation.classify(model);
var rmse = validationPredictions
.errorMatrix('carbon', 'classification')
.accuracy();
print('RMSE:', rmse);
这里我没有生物量实测数据，也就没必要展示结果了
参考：
《全球陆地碳汇的遥感和优化计算方法》陈镜明
《GEDI L4A Footprint Level Aboveground Biomass Density, Version2.1》
https://spatialthoughts.com/2024/02/07/agb-regression-gee/
《Assessing Protected Area's Carbon Stocks and Ecological Structure at Regional-Scale Using Gedi Lidar》Mengyu Liang
