简介面向遥感与机器学习交叉领域研究者的Python实战指南以津巴布韦Mafungautsi森林保护区为测试区系统讲解如何整合GEDI L4A星载激光雷达脚印、Sentinel-2多光谱影像、光谱指数、SRTM海拔与坡度数据并应用随机森林算法完成地上生物量密度AGBD的建模与空间预测。资源仅含1个PDF文件141KB内容为完整的GEEPython操作流程覆盖Earth Engine账户初始化与认证、哨兵二号合成影像创建、特征变量构建、训练/测试集划分、模型性能评估以及过拟合现象的诊断与缓解策略。全文注重代码与原理结合既给出可复现的实现步骤也说明GEDI L4A产品中RH指标与AGBD参数化模型的背景适合希望掌握多源遥感数据融合与机器学习建模的林业、生态遥感研究人员。现已有101人学习作为轻量级入门材料可在短时间内快速建立从数据获取到结果可视化的完整认知。1. 用GEDI、Sentinel-2和随机森林估算森林地上生物量从像元到碳储量的捷径森林地上生物量密度AGBD一直是碳汇监测、林业调查和生态研究里最让人头疼的指标。传统做法是跑野外样地测胸径、树高、密度再用异速生长方程换算成每公顷生物量一亩一亩地推。样地布设密度稍微薄一点插值结果就没法看。2019年之后GEDI激光雷达给了另一条路它把全球森林切成一串直径约25米的足迹波形数据直接反演出每个足迹的地上生物量密度。把这些足迹当标签配上Sentinel-2的多光谱特征用随机森林回归做外推就能把点状观测铺成一张连续的生物量密度图。整套流程用Python就能跑通不需要深度学习环境一台普通笔记本加一个下午的数据清洗就能出第一版结果。真正花时间的不是建模而是把GEDI的HDF5文件、Sentinel-2的投影体系和随机森林的调参逻辑理顺。2. 数据准备GEDI足迹与Sentinel-2影像的坐标对齐和文件组织2.1 为什么用GEDI当标签、Sentinel-2当特征先说我为什么这么搭配。GEDI全称是全球生态系统动力学调查搭载在国际空间站上的激光雷达。它发射的激光脉冲穿过林冠、反射回波形波形里包含了冠层高度、垂直结构信息这些信息与生物量的统计关系比任何光学指数都直接。NASA的L2B产品直接给出了每个足迹的agbd字段单位是Mg/ha也就是每公顷干生物量省去了自己反演波形的功夫。它的脚印是一个约25米直径的圆这个尺度恰好和Sentinel-2的10米/20米混合分辨率对得上比MODIS那种500米像元精细一个数量级。Sentinel-2这边我常用的是L2A表面反射率产品。它的波段配置对植被特别友好B8近红外和B4红波段能算NDVIB5到B7的红边波段对冠层叶绿素敏感B11/B12短波红外对含水量和结构敏感。为什么不用Landsat也不是不行但30米分辨率对GEDI的25米足迹来讲混合像元更严重而且Sentinel-2五天一次的重访更容易凑出无云合成影像。随机森林是这三件套里最省心的部分。它不需要对特征做标准化能处理特征之间的共线性对几千到几万个样本的训练集来说几分钟就能跑完。更重要的是它能输出特征重要性这对后续向领导解释“为什么这片区域生物量高”很有用。遥感随机森林这条路线在很多碳汇项目里已经是被验证过的常见做法了。2.2 文件组织与数据下载后的目录规划下载数据这一步没有太多技巧GEDI从NASA的LP DAAC节点下载Sentinel-2从Copernicus Data Space Ecosystem下载或者直接在GEE里导出研究区合成影像。我一般会在项目根目录下建这样的结构保持可复现agbd_project/ ├── gedi/ │ ├── GEDI02_B_2019108041955_O03015_02_T01917_02_002_002.h5 │ └── GEDI02_B_202001... ├── s2/ │ ├── L2A_T50xxx_B02_10m.tif │ ├── L2A_T50xxx_B03_10m.tif │ ├── L2A_T50xxx_B04_10m.tif │ ├── L2A_T50xxx_B05_20m.tif │ └── ... ├── output/ │ ├── gedi_footprints.geojson │ └── s2_feature_stack.tif └── scripts/ ├── 01_extract_gedi.py ├── 02_build_features.py └── 03_train_rf.pySentinel-2原始产品是一片一片的SAFE文件夹里面分包着不同分辨率。我建议先把B5、B6、B7、B11、B12这些20米波段重采样到10米和B2、B3、B4、B8保持同一网格。这一步用gdal warp或rasterio都能做核心是统一投影和像元尺寸。GEDI的坐标是WGS84经纬度而Sentinel-2 L2A产品是按UTM分带的坐标对齐是后面所有步骤的前提。2.3 坐标基准检查别让投影不一致毁掉匹配我见过不止一次有人拿着WGS84的足迹经纬度直接去切UTM的Sentinel-2影像切出来偏了几百米还在找原因。GEDI给的是椭球坐标Sentinel-2的GeoTIFF头里写的投影通常是UTM。先把GEDI足迹读进来转成和目标影像一致的投影再去做缓冲区。代码如下import geopandas as gpd from shapely.geometry import Point import rasterio # 读取或构建GEDI足迹点 gdf gpd.GeoDataFrame( {agbd: [100.3, 85.4]}, geometry[Point(112.34, 27.61), Point(112.35, 27.62)], crsEPSG:4326 ) # 读一张Sentinel-2影像获取目标投影 with rasterio.open(s2/L2A_T50xxx_B02_10m.tif) as src: target_crs src.crs transform src.transform gdf_utm gdf.to_crs(target_crs) # 检查坐标是否落在影像范围 with rasterio.open(s2/L2A_T50xxx_B02_10m.tif) as src: bounds src.bounds x, y gdf_utm.geometry.x.values, gdf_utm.geometry.y.values valid ( (x bounds.left) (x bounds.right) (y bounds.bottom) (y bounds.top) ) print(f落在影像内的足迹比例: {valid.mean():.2f})这段代码的逻辑是把经纬度点统一到Sentinel-2所在的UTM投影然后检查有多少足迹落在影像范围里。如果这个比例低于90%先别急着建模大概率是研究区跨了两个UTM分带或者GEDI文件里有日期与研究区不匹配的轨道数据。比例异常时优先检查数据时间范围而不是调整模型。另外要注意这里的transform在计算窗口时有用后面切像元时会反复用到。3. 从HDF5提取GEDI生物量标签质量筛选与标准化处理3.1 GEDI L2B的HDF5结构和关键字段GEDI的HDF5文件打开后顶层是一组BEAM开头的光束目录从BEAM0000到BEAM0019。每个光束内部是并行的数组字段所有字段长度相同对应一次激光发射的shot。L2B产品里最核心的是agbd和agbd_se前者是每个足迹的地上生物量密度单位Mg/ha后者是它的标准误。除了生物量本身还有三个质量字段必须处理quality_flags、degrade_flag、sensitivity。quality_flags是算法质量标志等于1才说明这束激光的波形反演通过质检。degrade_flag是卫星姿态和轨道控制的状态标志等于1说明该时段卫星处于降级模式footprint位置精度会明显变差这类数据要剔除。sensitivity则是波形对地表的灵敏度反映激光脉冲穿透林冠后还能不能探测到地面数值在0到1之间越接近1说明越可靠。高生物量密度的林分sensitivity很容易掉到0.95以下直接把阈值卡死会丢掉大量高值样本这个矛盾留到第5章细说。3.2 用 h5py 遍历光束并构建样本表读取GEDI的标准姿势是h5py遍历光束再把筛选后的记录拼成DataFrame。我写的提取脚本长这样import h5py import numpy as np import pandas as pd h5_path gedi/GEDI02_B_2019108041955_O03015_02_T01917_02_002_002.h5 recs [] with h5py.File(h5_path, r) as f: # 过滤顶层目录里非BEAM开头的元数据组 beam_names [k for k in f.keys() if k.startswith(BEAM)] for beam in beam_names: g f[beam] qf g[quality_flags][:].astype(np.uint8) dg g[degrade_flag][:].astype(np.uint8) sen g[sensitivity][:] agbd g[agbd][:] agbd_se g[agbd_se][:] lat g[lat_lowestmode][:] lon g[lon_lowestmode][:] shot g[shot_number][:] for i in range(len(agbd)): recs.append( { beam: beam, shot_number: shot[i], lon: lon[i], lat: lat[i], agbd: agbd[i], agbd_se: agbd_se[i], quality_flags: qf[i], degrade_flag: dg[i], sensitivity: sen[i], } ) df pd.DataFrame(recs) print(f原始记录数: {len(df)})这段的逻辑是逐光束读取把每束光的全部shot记录收集起来最后统一过滤。注意每个数组都要加[:]才会真正读入内存否则拿到的是HDF5的数据集对象后面没法按元素判断。GEDI一个文件动辄几十万条记录个人经验是一次只处理一个文件处理完保存中间结果不要把所有文件同时读进内存不然16G内存的机器会很吃力。光束之间样本量差异很大这是轨道扫描几何造成的正常现象不需要刻意均衡。3.3 质量筛选与样本清洗拿到原始表之后清洗逻辑比读数据更重要。我按下面的顺序做每一步都打印删除数量df df[df[quality_flags] 1].copy() df df[df[degrade_flag] 0].copy() df df[df[agbd].notna()].copy() df df[(df[agbd] 0) (df[agbd] 800)].copy() # sensitivity 阈值0.90 起步0.95 太严时回退 df df[df[sensitivity] 0.90].copy() # 相对误差过滤标准误超过生物量一半的样本不可信 df df[df[agbd_se] / df[agbd] 0.5].copy() # 去除明显离群的经纬度跑不在研究区附近的 df df[(df[lon] 110) (df[lon] 114) (df[lat] 25) (df[lat] 30)].copy() print(f筛选后样本量: {len(df)})参数说明agbd上限设800 Mg/ha是给极端高密度林留的余量如果你研究区是温带林可以收紧到400sensitivity阈值0.90是我在热带和亚热带数据上的折中方案0.95会丢掉太多高值样本0.88以下又会混入波形没穿透冠层的数据相对误差0.5这条很多教程不写但实际效果很好它能把agbd_se异常大的低质量回归样本过滤掉。清洗完把结果存成GeoJSON后面特征提取直接读这个文件import geopandas as gpd from shapely.geometry import Point gdf gpd.GeoDataFrame( df, geometry[Point(x, y) for x, y in zip(df[lon], df[lat])], crsEPSG:4326, ) gdf.to_file(output/gedi_footprints.geojson, driverGeoJSON) print(f已保存 {len(gdf)} 个足迹到 GeoJSON)这里有个细节GeoJSON只支持WGS84坐标为了后续和Sentinel-2对齐我通常会再存一份转成UTM的版本gdf_utm gdf.to_crs(EPSG:32650) gdf_utm.to_file(output/gedi_footprints_utm50.gpkg)EPSG:32650是UTM 50N对应中国中东部的大多数研究区。如果你在研究区不知道带号用elevation数据先查一下或者直接用Sentinel-2影像自带的crs这是最不易出错的方案。4. Sentinel-2特征工程与随机森林把像元值变成训练特征4.1 特征设计波段、指数和纹理怎么取舍Sentinel-2进入模型的特征要分三组。第一组是核心波段我常用B2、B3、B4、B5、B6、B7、B8、B11、B12这里故意没放B1和B9它们气溶胶和水汽吸收波段对生物量建模贡献极低。第二组是植被指数NDVI、EVI、NDMI、Cire这四项足够覆盖绿度、冠层结构和水分特征。第三组是纹理特征比如红边波段的窗口标准差能反映林分破碎度和冠层粗糙度对区分同龄林和异龄林有帮助。特征不是越多越好随机森林虽然在特征多时不容易马上过拟合但无关特征会稀释特征重要性的可解释性。我的经验是先用下面这张表做最小集合跑完第一版再加纹理不要一步到位堆到三四十个特征。特征组变量名来源波段/计算方式用途可见光B2 B3 B410m波段基础光谱近红外B810m波段植被高反射区红边B5 B6 B720m重采样10m叶绿素和冠层结构短波红外B11 B1220m重采样10m含水量和结构绿度指数NDVI(B8-B4)/(B8B4)植被覆盖度增强指数EVI2.5*(B8-B4)/(B86B4-7.5B21)抗土壤背景水分指数NDMI(B8-B11)/(B8B11)冠层含水量红边指数CIreB7/B5-1叶绿素含量4.2 构建特征堆栈并提取足迹内像元先做一步预处理把S2的20米波段重采样到10米再和10米波段一起合成一个多波段GeoTIFF。重采样用rasterio就够方法选bilinear别选nearest因为20米到10米是加密像元最近邻会出现块状锯齿。构建堆栈的代码import rasterio import numpy as np from rasterio.warp import reproject, Resampling def resample_to_10m(src_path, dst_path): with rasterio.open(src_path) as src: profile src.profile.copy() scale src.width / (src.bounds.right - src.bounds.left) * 10 new_height int(src.height / (src.res[0] / 10)) new_width int(src.width / (src.res[1] / 10)) profile.update(heightnew_height, widthnew_width, transformsrc.transform * src.transform.scale( (src.width / new_width), (src.height / new_height) )) # 实际用 reproject 更稳妥这里简写 data src.read( out_shape(1, new_height, new_width), resamplingResampling.bilinear ) with rasterio.open(dst_path, w, **profile) as dst: dst.write(data)这里逐步说明参数含义new_height和new_width由原始分辨率与目标分辨率10米的比例算出重采样方法Resampling.bilinear适合连续型反射率数据profile继承原文件的坐标系和CRS避免投影出问题。各波段重采样后用gdal_merge或rasterio的merge把九波段四个指数合成一个12波段的feature_stack.tif合成后必须检查波段顺序后面特征矩阵的列名要和它一一对应。提取足迹内像元值我推荐数组切片而不是逐个调用rasterio.mask.mask因为足迹数量几万条时mask函数逐次读盘会慢到怀疑人生。先整幅读入内存按中心行列号切片import rasterio import numpy as np from rasterio.transform import rowcol src rasterio.open(output/s2_feature_stack.tif) stack src.read() # 形状: (12, height, width) stack[stack -9999] np.nan gdf gpd.read_file(output/gedi_footprints_utm50.gpkg) # 影像像元尺寸按10米算GEDI足迹半径12.5米约1.25个像元 radius_px 2 # 取3x3窗口 rows, cols rowcol(src.transform, gdf.geometry.x.values, gdf.geometry.y.values) feat_matrix [] valid_idx [] for i, (r, c) in enumerate(zip(rows, cols)): if r - radius_px 0 or c - radius_px 0: continue if r radius_px 1 stack.shape[1] or c radius_px 1 stack.shape[2]: continue # 3x3窗口内各波段的均值 patch stack[:, r - radius_px:r radius_px 1, c - radius_px:c radius_px 1] patch patch.reshape(patch.shape[0], -1) if np.isnan(patch).all(axis1).any(): continue feat_matrix.append(np.nanmean(patch, axis1)) valid_idx.append(i) feat_matrix np.array(feat_matrix) print(f提取成功特征样本数: {feat_matrix.shape})参数radius_px取2是经验值——实际上GEDI足迹直径25米在10米影像上是2.5个像元取2或3差别不大取3会引入更多邻域混合取1则只覆盖足迹中心16%。我用2也就是5x5像元窗口算下来覆盖面积2500平方米比足迹略大但有缓冲余量。窗口太小时足迹中心落在像元边缘会导致提取值和实际地物偏差大这是做生物量建模很容易踩的坑。4.3 随机森林回归训练、验证与特征重要性特征矩阵准备好标签从GeoJSON里按valid_idx对应取然后开训。随机森林回归算法在sklearn里实现得很成熟不需要自己写。我的参数设置如下from sklearn.ensemble import RandomForestRegressor from sklearn.model_selection import train_test_split from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error import numpy as np feature_names [B2,B3,B4,B5,B6,B7,B8,B11,B12, NDVI,EVI,NDMI,CIre] X feat_matrix y gdf.iloc[valid_idx][agbd].values.astype(float) # 过滤非有限值 mask np.isfinite(X).all(axis1) np.isfinite(y) X, y X[mask], y[mask] X_train, X_test, y_train, y_test train_test_split( X, y, test_size0.2, random_state42 ) rf RandomForestRegressor( n_estimators300, max_depth15, min_samples_leaf3, min_samples_split5, max_featuressqrt, n_jobs-1, random_state42, ) rf.fit(X_train, y_train) y_pred rf.predict(X_test) rmse np.sqrt(mean_squared_error(y_test, y_pred)) mae mean_absolute_error(y_test, y_pred) r2 r2_score(y_test, y_pred) print(fR2: {r2:.3f}) print(fRMSE: {rmse:.2f} Mg/ha) print(fMAE: {mae:.2f} Mg/ha)参数默认值为什么调这几个n_estimators在300棵之后RMSE下降速度已经很慢再往上只是线性增加训练时间max_depth限制在15可以避免树对单条样本死记硬背尤其当样本量只有几千时这个约束很关键min_samples_leaf3和min_samples_split5是防止树在低样本区再分裂回归树在样本稀疏区域容易拟合出尖刺。max_featuressqrt是随机森林的经典选择降低树间相关性比增加每棵树的能力更重要。n_jobs-1是让所有CPU核心参与。看模型输出不要只看R2要结合RMSE和MAE一起判断。RMSE远大于MAE说明预测在高值区有大幅度偏差这大概率是随机森林对极端高生物量区域系统性低估——这是随机森林回归的老毛病高值样本被平均值拉回来。特征重要性此刻能帮你判断哪些变量在主导imp pd.Series(rf.feature_importances_, indexfeature_names).sort_values(ascendingFalse) print(imp.head(6))通常NDVI、B11和红边波段会排在前列。如果B2、B3这些可见光排第一先检查是不是有云阴影或裸地样本混进来了这种样本会让模型学会区分“有没有植被”而不是“有多少生物量”。5. 常见翻车现场与排查别让三个看似正确的步骤毁掉模型5.1 quality_flags 忘了筛选导致模型学到噪声现象训练集R2能到0.8验证集R2掉到0.3特征重要性分布完全不符合常识。原因GEDI L2B的quality_flags等于0的记录是反演失败或低质量的波形但这些记录里agbd字段不是空值而是算法硬算出来的填充值。把它们当标签训练等于把大量随机噪声教给模型。这个问题在第一版建模里最常见因为HDF5读取时agbd字段看起来有数不会报错。解决在提取环节就必须按quality_flags1过滤且每个文件都检查过滤前后样本量差异。如果过滤掉的样本超过一半说明这期数据质量不佳考虑更换时间段的GEDI数据而不是靠模型去吸收噪声。5.2 提取Sentinel-2像元时只取中心点尺度全碎现象模型在本地验证集上RMSE表现尚可但整幅外推成图后出现大量椒盐状噪点空间连续性很差。原因GEDI一个足迹是25米直径圆你在Sentinel-2上只取中心那个10米像元的反射率值相当于用四分之一的信息代表整个足迹。地物在25米范围内不均匀时光谱值和生物量的对应关系被严重削弱。外推到面上时模型对相邻像元的光谱微变反应过度所以出现椒盐效果。解决提取特征时必须做窗口聚合半径至少取2个像元5x5窗口。我做项目时的血泪经验是宁可窗口取大一点引入邻域混合也不要取小导致代表性不足。窗口大小对结果影响比树的数量大得多这个超参数值得用交叉验证专门搜一轮。5.3 sensitivity 阈值越高越好现象阈值设0.95之后样本量从一万掉到两千且样本的agbd分布明显偏向低值区模型在高生物量区间的预测能力崩溃。原因sensitivity反映的是激光脉冲从冠层顶部到地面的衰减程度。高生物量的浓密林冠本来就会让地面回波变弱sensitivity天然偏低。你把阈值卡得越高保留下来的样本越是低矮稀疏植被训练数据丢失了高值段模型自然学不到高生物量长什么样。解决把阈值定在0.90到0.95之间同时看筛选后的agbd直方图保证高值区还有足够样本。如果研究区确实以高密度林为主可以考虑放宽到0.88观察交叉验证RMSE的变化。这个参数没有普适最优和你研究区的冠层密度分布强相关。5.4 时间没对齐GEDI和Sentinel-2差了好几个季节现象模型在训练集表现很好但分季节验证时叶落期的样本预测偏差明显增大。原因GEDI的采集时间是轨道过境瞬间Sentinel-2反射率则是某个时相或某段时间的合成影像。落叶林的生物量全年相对稳定但NDVI、红边指数受物候影响波动非常大。如果标签是夏天的GEDI特征用的是冬天影像模型学到的是“这个光谱形态对应某生物量”而不是真正的生物量与结构关系。解决要么把Sentinel-2影像按月份匹配到GEDI过境时间附近要么给特征里加一列“影像月份”作为辅助变量。更省事的做法是用整个生长季的合成影像做特征让回归关系更能容忍物候漂移。5.5 用随机分布的验证集骗自己现象整体随机切分验证的R2很高但按区域分块验证时R2骤降0.3以上。原因遥感像元存在强烈的空间自相关相邻足迹的光谱和生物量高度相似。随机切分时同一个林分的样本同时出现在训练和验证集里模型相当于开了“天眼”。这在所有遥感机器学习任务里都成立不是随机森林特有的问题而是验证策略的黑匣子。解决做空间交叉验证按足迹聚成格网例如把研究区切城3km乘以3km的块按块划分训练和验证集。如果按块验证后模型表现依然可接受说明泛化能力是真实的。外推的全图预测精度永远应该以空间交叉验证的结果为准而不是随机切分的那个数字。6. 从足迹到一张图外推预测与精度验证的收尾技巧模型训好后最后一步是把全研究区的Sentinel-2特征堆栈喂给模型逐像元预测生物量密度并输出GeoTIFF。这里有个操作细节大影像不能一次性成块predict内存会爆。按行分块是通用做法import rasterio import numpy as np rf trained_rf # 已训练好的模型 block_size 512 with rasterio.open(output/s2_feature_stack.tif) as src: profile src.profile.copy() profile.update(dtypenp.float32, count1, compresslzw) with rasterio.open(output/agbd_predicted.tif, w, **profile) as dst: for row in range(0, src.height, block_size): rows min(block_size, src.height - row) window rasterio.windows.Window(0, row, src.width, rows) arr src.read(windowwindow, boundlessTrue, fill_value-9999) # 转换成模型输入的二维数组 (n_pixels, n_features) orig_shape arr.shape[1:] flat arr.reshape(12, -1).T # 无效像元跳过 valid np.isfinite(flat).all(axis1) pred np.full(flat.shape[0], np.nan, dtypenp.float32) pred[valid] rf.predict(flat[valid]) dst.write(pred.reshape(orig_shape).astype(np.float32), 1, windowwindow)分块后要注意每个窗口读取时如果影像边缘不足一块用boundlessTrue配合fill_value填充保证窗口形状一致。预测后的整幅GeoTIFF直接在GIS软件里打开检查有没有异常的高值斑块比如道路或水体被预测出了生物量——这说明训练样本里混了非林地足迹需要回到第3章加一步土地利用掩膜过滤。最后一件事把预测图与GEDI验证足迹的散点图画出来。横轴是GEDI观测的agbd纵轴是模型预测的agbd按2D密度着色我习惯给散点图叠加1:1线。如果高值区的点明显整体偏向1:1线下方不要怀疑数据有问题这是随机森林回归的固有权衡——想在高值区拉回偏差可以尝试用极端梯度提升对高值样本加权或者在数据层面补充更多高生物量样本。我自己做第一个项目时就是因为没做空间交叉验证高估了模型后来按块重验才看清真实精度。这个流程里的每一步都值得你亲手跑一遍原代码再出门希望帮到你。本文还有配套的精品资源点击获取