GEE实战：从零构建京津冀地区土地利用随机森林分类模型

1. 京津冀土地利用分类实战入门第一次接触Google Earth EngineGEE做土地利用分类时我被它强大的云端计算能力震撼到了。不用下载海量遥感数据直接在浏览器里就能完成从数据处理到模型训练的全流程。这次我们就用京津冀地区作为案例手把手带你完成一个完整的随机森林分类项目。京津冀作为我国重要的城市群土地利用类型复杂多样。我们选用2023年的Landsat 9数据因为它的空间分辨率30米和光谱特性非常适合中尺度土地利用分类。整个流程可以分为五个关键步骤数据准备、样本采集、特征工程、模型训练和精度验证。最棒的是所有操作都在GEE的云端完成你的笔记本电脑再也不用担心跑不动大数据量了。2. 数据准备与预处理2.1 研究区与时间窗口设置首先需要明确研究区域和时间范围。京津冀地区包含北京、天津和河北三个省级行政区我们可以用GEE提供的行政区划数据来划定范围// 加载中国省级行政区划数据 var province ee.FeatureCollection(users/cduthes1991/boundry/China_province_2019); // 筛选京津冀地区 var beijing province.filter(ee.Filter.eq(provinces,beijing)); var tianjin province.filter(ee.Filter.eq(provinces,tianjin)); var hebei province.filter(ee.Filter.eq(provinces,hebei)); // 合并三个区域 var roi beijing.merge(tianjin).merge(hebei); Map.centerObject(roi, 7); Map.addLayer(roi.style({color:red,fillColor:00000000}),{},京津冀边界);时间窗口选择很关键建议避开冬季12月-2月和雨季7-8月。我测试发现3-6月和9-10月的影像质量最好var startDate 2023-03-01; var endDate 2023-10-01;2.2 Landsat 9数据处理Landsat 9数据需要做三个预处理辐射定标、去云和计算光谱指数。这三个步骤可以封装成函数链式调用// 辐射定标函数 function applyScaleFactors(image) { var opticalBands image.select(SR_B.).multiply(0.0000275).add(-0.2); var thermalBands image.select(ST_B.*).multiply(0.00341802).add(149.0); return image.addBands(opticalBands, null, true) .addBands(thermalBands, null, true); } // 去云函数 function rmL9Cloud(image) { var cloudBitMask (1 3); var cloudShadowBitMask (1 4); var qa image.select(QA_PIXEL); var mask qa.bitwiseAnd(cloudBitMask).eq(0) .and(qa.bitwiseAnd(cloudShadowBitMask).eq(0)); return image.updateMask(mask); } // 光谱指数计算函数 function addIndices(image) { var ndvi image.normalizedDifference([SR_B5,SR_B4]).rename(ndvi); var ndbi image.normalizedDifference([SR_B6,SR_B5]).rename(ndbi); var mndwi image.normalizedDifference([SR_B3,SR_B6]).rename(mndwi); return image.addBands(ndvi).addBands(ndbi).addBands(mndwi); }3. 样本点采集与处理3.1 样本采集策略样本质量直接影响分类效果。建议每类至少采集100个样本点均匀分布在整个研究区。在GEE中可以通过手动绘制或导入已有样本// 示例创建水体样本 var water ee.FeatureCollection([ ee.Feature(ee.Geometry.Point([116.3, 39.9]), {landcover: 1}), ee.Feature(ee.Geometry.Point([117.2, 38.8]), {landcover: 1}), // 更多样本点... ]); // 建筑用地样本 var urban ee.FeatureCollection([ ee.Feature(ee.Geometry.Point([116.4, 39.9]), {landcover: 2}), // 更多样本点... ]); // 林地样本 var forest ee.FeatureCollection([ ee.Feature(ee.Geometry.Point([115.5, 40.2]), {landcover: 3}), // 更多样本点... ]);3.2 样本划分与平衡样本需要随机划分为训练集70%和验证集30%。我建议使用stratifiedSplit保持各类样本比例// 合并所有样本 var samples water.merge(urban).merge(forest); // 添加随机数列 var samples samples.randomColumn(random); // 分层抽样 var training samples.filter(ee.Filter.lte(random, 0.7)); var validation samples.filter(ee.Filter.gt(random, 0.7)); print(训练样本数, training.size()); print(验证样本数, validation.size());4. 随机森林模型构建4.1 特征选择与模型训练除了原始波段我们加入NDVI、NDBI和MNDWI三个指数作为特征。随机森林的关键参数包括树的数量300、每棵树考虑的特征数默认sqrt// 准备特征图像 var features img.select([SR_B2,SR_B3,SR_B4,SR_B5,SR_B6,SR_B7, ndvi,ndbi,mndwi]); // 训练随机森林分类器 var classifier ee.Classifier.smileRandomForest(300) .train({ features: training, classProperty: landcover, inputProperties: features.bandNames() }); // 执行分类 var classified features.classify(classifier);4.2 分类结果可视化为不同地类设置合适的颜色// 定义可视化参数 var visParams { min: 1, max: 3, palette: [blue, red, green] }; Map.addLayer(classified, visParams, 分类结果); Map.addLayer(foods, {color: yellow}, 验证样本);5. 精度评估与结果导出5.1 精度验证指标使用混淆矩阵计算各项精度指标// 验证集预测 var validated validation.classify(classifier); // 计算混淆矩阵 var confusionMatrix validated.errorMatrix(landcover, classification); // 输出精度指标 print(总体精度, confusionMatrix.accuracy()); print(Kappa系数, confusionMatrix.kappa()); print(用户精度, confusionMatrix.consumersAccuracy()); print(生产者精度, confusionMatrix.producersAccuracy());5.2 结果导出分类结果可以导出到Google Drive或GEE Assets// 导出到Google Drive Export.image.toDrive({ image: classified, description: 京津冀土地利用分类, fileNamePrefix: LUCC_京津冀_2023, region: roi, scale: 30, maxPixels: 1e13, fileFormat: GeoTIFF, formatOptions: { cloudOptimized: true } });6. 常见问题与优化建议在实际项目中有几个容易踩的坑需要注意样本不平衡问题如果某类样本过少可以设置classWeights参数平衡各类权重特征相关性高相关性的特征会影响模型表现可以用ee.Reducer.correlation()检查过拟合问题增加树的数量不一定能提高精度建议通过交叉验证找到最优参数季节性影响可以尝试加入多时相数据提高分类精度一个实用的技巧是在模型训练前对特征进行标准化// 特征标准化 var meanStd features.reduceRegion({ reducer: ee.Reducer.mean().combine({ reducer2: ee.Reducer.stdDev(), sharedInputs: true }), geometry: roi, scale: 30, maxPixels: 1e13 }); var featuresScaled features.map(function(image) { return image.subtract(meanStd.get(mean)).divide(meanStd.get(stdDev)); });7. 进阶扩展思路完成基础分类后还可以进一步后处理使用ee.Kernel.square()进行多数滤波消除小斑块变化检测比较不同年份的分类结果分析土地利用变化精度提升加入高程数据、夜间灯光数据等辅助特征自动化将整个流程封装为函数方便应用到其他区域比如加入DEM数据var dem ee.Image(CGIAR/SRTM90_V4).clip(roi); var slope ee.Terrain.slope(dem); var featuresWithDEM features.addBands(dem).addBands(slope);