简介这是一份面向遥感与机器学习研究者的Python实战指南聚焦利用GEDI L4A星载激光雷达足迹级生物量数据、Sentinel-2多光谱影像、光谱指数及SRTM高程/坡度数据结合随机森林算法完成地上生物量密度AGBD建模。教程以Mafungautsi森林保护区为测试区完整覆盖Earth Engine账户初始化与API认证、Sentinel-2合成影像构建、光谱指数计算、训练与测试数据准备、随机森林模型运行、性能评估以及AGBD空间预测与可视化。文档同时解释了GEDI L4A数据的基本原理如足迹级AGBD估算、基于波形相对高度RH指标与模拟波形建模以及按全球区域和植物功能类型分层建模等有助于读者理解数据来源与模型假设。特别讨论了过拟合现象及提升泛化能力的思路适合关注森林碳储量和生产力评估的工程师深入学习。压缩包共1个PDF文件大小仅141KB轻量便携。目前已有101人学习这份资料对正在开展多源遥感生物量估算的读者具有直接参考价值。1. 基于GEDI和Sentinel-2的地上生物量密度建模为什么随机森林是第一个值得复现的基线大范围碳储量估算里实地样方不可能铺满研究区比较靠谱的替代方案是让GEDI激光雷达脚点提供生物量真值Sentinel-2光学影像提供空间连续的特征面再用随机森林把两者打通。用Python走完这条链路是遥感机器学习里最经典、也是绝大多数生态项目默认的第一版基线。无论你做毕业设计、城市森林碳密度图还是区域生态区划都会遇到同一个问题怎么把激光雷达的离散脚点变成一张连续AGBD空间分布图。这篇笔记不绕开数据处理细节直接把GEDI脚点筛选、Sentinel-2特征栈构建、脚点-像元匹配、模型训练和空间验证中容易翻车的地方讲透。2. 数据获取与预处理GEDI脚点质量筛选和Sentinel-2特征栈搭建预处理决定建模上限。GEDI官方产品按轨道发布HDF5文件Sentinel-2在各大云平台按分幅存储这两种数据格式完全不同想要喂给同一个Python训练流程得先把它们清洗成干净的表结构和栅格栈。2.1 GEDI L4A脚点提取从HDF5到干净样本表我习惯直接用GEDI L4A的AGBD产品它给出的就是Mg/ha单位的地上生物量密度省去从L2A波形反演AGBD那一步。L4A以HDF5格式分发文件里同时含有大量波形模拟参数但建模真正需要的字段不多。拿到轨道文件先别急着跑完整脚本应当打印HDF5的顶层键确认当前版本字段名和旧教程是否一致。import h5py import numpy as np import pandas as pd h5_path GEDI_L4A_2020_xxx_yyy.h5 with h5py.File(h5_path, r) as f: # L4A常用路径是 /gediL4A/ 下的一级字段 lon f[/gediL4A/lon_lowestmode][:] lat f[/gediL4A/lat_lowestmode][:] agbd f[/gediL4A/agbd][:] agbd_se f[/gediL4A/agbd_se][:] qflag f[/gediL4A/quality_flag][:] degrade f[/gediL4A/degrade_flag][:]lon_lowestmode和lat_lowestmode表示最低模式高程对应的经纬度实际使用中通常把它当作脚点位置。agbd_se是标准不确定度做样本加权时可以用也可以在筛选时设定阈值来过滤传感器观测质量差的脚点。质量筛选是第一个容易出错的动作也是最关键的valid ((qflag 1) (degrade 0)) valid (agbd 0) (agbd 1500) (~np.isnan(agbd)) # 经验上把agbd_se 50的样本剔除避免传感器不确定性过大 valid (agbd_se 50) df pd.DataFrame({ lon: lon[valid], lat: lat[valid], agbd: agbd[valid], agbd_se: agbd_se[valid], }) df df.drop_duplicates(subset[lon, lat])qflag等于1代表官方质量位通过degrade_flag为0代表卫星姿态和指向正常。agbd上限1500是一个粗筛森林AGBD通常不会超过这个值如果你做的是稀树草原或灌丛建议把上限降得更低。agbd_se小于50这个阈值属于我个人的经验值并不是官方标准样本量紧张时可以放宽到80但模型残差会明显增大。还有一点容易被忽略GEDI的轨道文件是沿轨60米步长采样同一轨道内相邻脚点的空间自相关很强。下载样本时不要只拿一条轨道我一般会收集研究区内两三个季节、多轨道的数据让脚点空间分布尽量分散否则后续随机森林很容易学到轨道方向上的虚假模式。2.2 Sentinel-2 L2A波段选择与辅助植被指数Sentinel-2我直接使用L2A地表反射率产品绕开自己跑大气校正。选波段时不需要把全部13个波段都放进去常用组合是B02、B03、B04、B08加B8A、B11、B12再额外算两个指数NDVI和NDMI基本覆盖植被绿度和水分状态。import rioxarray import xarray as xr band_map { B02: blue, B03: green, B04: red, B08: nir, B8A: nir_narrow, B11: swir1, B12: swir2, } data_vars [] for band_code, band_name in band_map.items(): da rioxarray.open_rasterio(fS2_L2A_{band_code}.tif).squeeze() da da.rename(band_name) da da.astype(float64) data_vars.append(da) stack xr.merge(data_vars)这里要注意Sentinel-2 L2A下载下来是整数DN值要先除以10000换算成反射率再参与指数计算。如果你用的是Level-1C产品必须自己跑Sen2Cor或依赖其他大气校正服务我为了减少无关配置直接选用L2A。所有波段必须重采样到同一网格。我一般以B08的10米网格作为基准把B8A、B11、B12从20米重采样到10米采用最近邻法而不是双线性。最近邻不会平滑掉植被边界对后续随机森林训练更稳妥。2.3 云掩膜与时间合成别让一个季相的云残留拖累训练单期Sentinel-2影像经常有云和云影残留直接用于特征提取会把红光和近红外波段拉出异常值。常规做法是选择研究区生长季内的多期无云L2A影像在像元级取中位数合成。import glob import numpy as np tif_files glob.glob(S2_L2A_2023_summer/*.tif) band_stack [] for f in tif_files: da rioxarray.open_rasterio(f) band_stack.append(da) # 先按像元判断有效值再取中位合成 s2_median xr.concat(band_stack, dimtime).median(dimtime)反演之后必须确认合成影像中每个像素是否在所有波段上有效。云掩膜后的边界区域通常会出现一个或几个波段是空值如果不处理后续按坐标提取特征时会拿NaN参与训练。一般做法是生成一个有效像素掩膜特征提取时判断窗口是否满足阈值。valid_mask np.isfinite(s2_median.to_array()).all(axis0) print(有效像元占比:, valid_mask.mean())如果有效像元占比低于80%说明合成影像残留空洞过多建议增加影像期数或换一个时间窗口而不是靠插值强行填洞。3. 训练样本构建把GEDI脚点映射到Sentinel-2像元的三种策略GEDI脚点是离散点Sentinel-2是连续栅格建模前必须把两者对齐到同一个样本表中。这里需要决策的不是算法而是“用多大窗口去匹配脚点”。窗口太小光学噪声大窗口太大把相邻地物混合进特征同样污染标签。3.1 单像元与3×3邻域窗口定位误差如何影响特征质量GEDI脚点的定位精度不是固定的受轨道指向和地形起伏影响实际位置可能偏移数米到十几米。如果直接让脚点坐标落在Sentinel-2的单像元上很可能会取到激光足迹边缘以外的地表信号。更稳的做法是用脚点坐标为中心提取3×3邻域窗口取均值近似覆盖GEDI约25米足印对应的光学信号。import rasterio from rasterio.windows import Window src rasterio.open(s2_feature_stack.tif) def extract_window_mean(src, lon, lat, half_window1): row, col src.index(lon, lat) w Window(col - half_window, row - half_window, 2 * half_window 1, 2 * half_window 1) data src.read(windoww, boundlessTrue) # 剔除NodeData边界避免越界像元进入统计 data[data src.nodata] np.nan data data.reshape(data.shape[0], -1).astype(float32) with np.errstate(invalidignore): mean_vals np.nanmean(data, axis1) return mean_valshalf_window等于1时是3×3窗口覆盖约900平方米和GEDI足迹量级匹配。如果脚点稀疏、影像空洞多可以放大到half_window2但窗口太大也会把不同树种的边界混合掉样本量和特征纯度之间需要权衡。另一个经常被忽视的问题是边界处理。脚点靠近影像边缘时Window读取会拿到空值或越界数据。上面代码用boundlessTrue把窗口外当作NaN处理再用np.nanmean过滤这样至少不会把无效窗口直接算成异常特征。3.2 构建特征矩阵与分组ID为空间交叉验证保留轨道信息把每条脚点的窗口均值聚合起来就得到特征矩阵。这一步的目标是生成一个干净的DataFrame行为脚点列为波段特征和AGBD标签。# 假设已经循环得到 features: (n_footprints, n_bands) X pd.DataFrame(features, columnsstack_band_names) y df[agbd].values # 按轨道块或空间网格生成分组ID供后续空间交叉验证使用 df_model[group_id] df.index // 500 df_model pd.concat([X, pd.Series(y, nameagbd)], axis1) df_model df_model.dropna(subsetstack_band_names) print(f有效样本: {len(df_model)})dropna这一步必须在建模前完成。很多新手在这一步跳过检查直到sklearn报错才回头处理NaN。更关键的是保留group_id后面做空间交叉验证时按轨道或空间区块分组而不是随机打乱所有样本。如果研究区包含明显的植被类型差异比如针叶林与阔叶林的反射率特征完全不同可以用土地覆盖类型生成更粗的分组ID。分组越粗空间验证越严格但训练数据量也会相应减少需要根据项目目标取舍。4. 随机森林与超参数调节从基线模型到能上线的预测器特征矩阵准备好之后就可以进入随机森林回归算法的训练环节。我的建议是不要一上来就做超参数大搜索先用一组保守参数把训练闭环跑通确认数据管道没有问题再做调优和空间验证。4.1 基线模型R2、RMSE与一次完整训练流程from sklearn.model_selection import train_test_split from sklearn.ensemble import RandomForestRegressor from sklearn.metrics import r2_score, mean_squared_error import numpy as np feature_cols [c for c in df_model.columns if c ! agbd] X df_model[feature_cols].values y df_model[agbd].values X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, random_state42 ) rf RandomForestRegressor( n_estimators300, max_depth12, min_samples_leaf2, max_features0.5, n_jobs-1, random_state42, ) rf.fit(X_train, y_train) y_pred rf.predict(X_test) print(R2:, r2_score(y_test, y_pred)) print(RMSE:, np.sqrt(mean_squared_error(y_test, y_pred)))max_features0.5是我常用的起始值。sklearn默认的“sqrt”在特征数量七八个时约等于3在特征数量二十个以上时会偏少0.5则能让树之间保持更强的随机性降低相关性。max_depth限制为12是为了防止单棵树在几万样本上把噪声一起背下来min_samples_leaf设为2则可以对离群脚点做一定压制。基线模型跑完先看RMSE而不是只看R2。R2对数据范围很敏感如果样本集中在一个窄生物量区间R2会非常好看但实际预测误差可能仍然很高。RMSE才是评估AGBD结果实用性的核心指标。4.2 超参数随机搜索先粗后细控制时间成本随机森林的超参数主要是n_estimators、max_depth、min_samples_leaf和max_features。网格搜索要跑大量组合样本量大时非常耗时我一般优先用RandomizedSearchCV做粗搜再在最优参数附近做一次细搜。from sklearn.model_selection import RandomizedSearchCV param_dist { n_estimators: [200, 300, 400, 500], max_depth: [8, 10, 12, 15, None], min_samples_leaf: [1, 2, 5], max_features: [0.2, 0.4, 0.6, 0.8, sqrt], } search RandomizedSearchCV( RandomForestRegressor(random_state42), param_dist, n_iter30, cv5, scoringneg_root_mean_squared_error, n_jobs-1, random_state42, ) search.fit(X_train, y_train) print(best params:, search.best_params_) print(CV RMSE:, -search.best_score_)注这里的cv5是普通K折并未考虑空间自相关存在一定程度的高估但它只用于快速筛选参数区间最终的模型评估要交给下一阶段的空间交叉验证。如果样本量超过5万n_estimators从300继续增加带来的提升有限真正影响精度的是max_depth和min_samples_leaf对异常点的压制。4.3 模型保存与整图预测训练完成后把模型保存下来再对整个Sentinel-2覆盖区做预测才能得到连续的AGBD分布图。直接逐像元循环预测很慢我一般按分块预测再拼接。import joblib joblib.dump(rf, rf_agbd_model.joblib)import numpy as np import xarray as xr def predict_raster(rf_model, s2_da, feature_order, chunk_size2000): ds s2_da[feature_order] arr ds.to_array().transpose(y, x, band).values h, w, _ arr.shape out np.full((h, w), np.nan, dtypefloat32) for i in range(0, h, chunk_size): for j in range(0, w, chunk_size): block arr[i:ichunk_size, j:jchunk_size] valid_mask np.isfinite(block).all(axis-1) if not valid_mask.any(): continue pred rf_model.predict(block[valid_mask]) temp np.full((block.shape[0], block.shape[1]), np.nan, dtypefloat32) temp[valid_mask] pred out[i:ichunk_size, j:jchunk_size] temp return s2_da.isel(xslice(0, w), yslice(0, h)).copy(dataout)这段代码以2000×2000的块为单位预测避免影像过大时内存溢出。np.isfinite判断所有波段有效才预测云洞和边缘像元直接输出NaN最后在GIS里可以单独筛选。5. 避坑清单GEDI与Sentinel-2随机森林建模的五条踩坑记录这一章是从实际项目中攒下来的问题每一条都经历过从现象到定位原因再到解决的完整过程。新手可以照着排查熟手也能对照自己的流程看看有没有漏掉。5.1 质量筛选后样本量骤减训练集空间分布严重失衡现象下载了整条GEDI轨道原始脚点几万个经过quality_flag和degrade_flag筛选后可用样本居然不足几百模型训练后RMSE高到离谱。原因GEDI轨道数据里本来就有相当一部分是低质量波形加上研究区地形起伏和云量影响degrade_flag1的点不少。如果只取一条轨道文件高质量脚点比例可能只有两到三成。解决不要用单条轨道拼一个模型。收集研究区内不同季节、多条轨道的样本再做筛选如果点还是不够可以引入L4A相邻时相的数据扩展时间窗口或是用GEDI L2B的RH95等结构指标做补充。宁可样本量稍大也要保证空间分布足够分散。5.2 特征栅格出现大量NaN训练时样本被静默丢弃现象提取特征矩阵后直接dropna发现很大比例脚点被丢掉模型训练完生成预测图发现大量区域是空值。原因Sentinel-2 L2A合成影像在云掩膜后边缘和空洞处的像元值本来就是NaN。10米像素的残留空洞在高植被覆盖区很常见尤其是回波边缘。解决先做多期影像中位数合成再做有效像元占比过滤。如果空洞范围超过窗口面积的20%我一般直接丢弃这个脚点不推荐用插值填洞。插值补出来的特征只是制造虚假平滑对随机森林没有任何信息增益。5.3 随机划分的测试集R2很好换到邻区却崩掉现象随机切分训练测试集时R2有0.85但把模型应用到相邻县或相邻年份影像上预测值和实测样地数据偏离严重。原因同一个轨道上相邻脚点的空间距离只有几十米随机划分测试集时很多空间邻近的样本被分到两边模型实际上记住了局部空间模式。另一个常见原因是B11、B12等短波红外波段在不同年份间波动很大模型对时间变化非常敏感。解决改用GroupKFold按轨道或空间区块分组。评估指标除了RMSE还要画残差与预测值的散点图看是否存在随生物量升高而增大的系统性偏差。如果要迁移到邻区最好补充目标区少量实测点做校准不要直接拿原模型做无约束外推。5.4 特征之间共线性强特征重要性不稳定现象feature_importances_里NDVI排第一但红光和近红外的重要度几乎一样换一批样本排序就变了。原因Sentinel-2原始波段和植被指数之间本就存在强相关随机森林面对高度可替代的特征组合时会把重要性随机分配给其中某一个导致结果不稳定。解决建模前先算特征相关性矩阵把相关系数绝对值大于0.9的特征成对剔除建模后再用permutation importance做一次验证。不要把随机森林的feature_importance当作因果解释工具它只能告诉你哪些特征在节点分裂中最常用。5.5 预测图出现沿轨道方向的条带状条纹现象生成AGBD分布图后发现高值和低值区域沿GEDI轨道方向呈现条带像一道道划痕。原因GEDI沿轨采样步长只有60米采样密度本身就携带了轨道方向的信息。如果特征矩阵里缺少描述森林垂直结构的变量随机森林会不自觉地使用脚点分布密度来拟合预测面自然出现条纹。另一个诱因是只用单条轨道数据训练地形和季节效应被模型记住。解决先判断条纹是沿Sentinel-2轨道还是沿GEDI轨道。属于GEDI条带时加入地形因子并汇入多轨道数据交叉训练通常能明显缓解。更严格的验证是用方向半变异函数检查预测面如果存在强方向性结构说明特征工程还有缺口而不是继续调参能解决的。6. 空间交叉验证与排列重要性把模型拷问一遍再出图随机森林训练完成后最忌讳直接拿着测试集R2去汇报。我的习惯是强制走一遍空间交叉验证和排列重要性分析这两步能让模型结论扎实很多。空间交叉验证用GroupKFold实现关键在于分组ID必须能代表空间块或轨道。from sklearn.model_selection import GroupKFold from sklearn.inspection import permutation_importance groups df_model[group_id].values gkf GroupKFold(n_splits5) cv_r2 [] cv_rmse [] for train_idx, val_idx in gkf.split(X, y, groupsgroups): rf_cv RandomForestRegressor(**search.best_params_, random_state42) rf_cv.fit(X[train_idx], y[train_idx]) pred rf_cv.predict(X[val_idx]) cv_r2.append(r2_score(y[val_idx], pred)) cv_rmse.append(np.sqrt(mean_squared_error(y[val_idx], pred))) print(空间CV R2:, np.mean(cv_r2)) print(空间CV RMSE:, np.mean(cv_rmse))GroupKFold确保同一个轨道或空间区块的样本不会同时出现在训练集和验证集里这样得到的精度才接近真实应用场景。排列重要性用于检查特征贡献的稳定性perm permutation_importance( rf, X_val, y_val, n_repeats20, n_jobs-1, random_state42 ) for i, val in enumerate(perm.importances_mean): print(f{feature_cols[i]}: {val:.4f})排列重要性的含义是把某个特征的值随机打乱后模型误差增加多少。增加越多说明模型对该特征依赖越大。它与随机森林自带的feature_importances不同不会因为特征之间存在共线性就把重要性集中到某一个上。从那以后我每次建模都会强制走一遍GroupKFold和排列重要性分析空间验证能过滤掉靠空间自相关刷出来的虚高精度排列重要性也能防止被共线波段误导。希望这些习惯能帮你在GEDI与Sentinel-2的地上生物量密度建模过程中少走弯路。本文还有配套的精品资源点击获取