做森林高度这件事我一直有个执念遥感出来的高度如果不是经过严格清洗的那宁可别用。GEEGoogle Earth Engine上手简单拉数据也快但“拉到数据”和“能用的数据”之间差着一整套数据处理和分析功夫。这篇文章就把我在GEE里处理森林高度数据以激光雷达脚印数据为主的完整思路和实操过程捋一遍包括数据怎么选、怎么清洗、怎么聚合、怎么跟辅助数据做交叉分析以及那些文档里很少写、但实测才会踩到的坑。如果你是刚接触遥感数据处理或者已经会用GEE但总觉得聚出来的高度图“心里没底”这篇文章可以直接当操作手册看。1. 项目思路与数据方案高度数据到底该怎么用1.1 先想清楚一个问题你要的高度是什么高度很多人一上来就问“GEE里有没有森林高度数据”这个问题其实没法直接回答因为“高度”至少有三种理解林冠高度从地面到树冠顶的高度即CHLCanopy Height Layer常用于森林生物量、碳储量估算。地面海拔DEM比如SRTM、COP30、NASADEM这类数据不是“树高”而是地形高度。体元高度/相对高度比如RH100、RH95这类从激光雷达波形导出的相对高度百分位数。我处理的项目最终目的是做区域林分高度制图回答“这个区域林分平均冠层高度是多少、空间上怎么分布”所以锚定的是第一种和第三种。如果你只是想看地形那完全不需要走激光雷达的路子直接加载DEM再计算山体阴影就行。这个前置判断很重要少走很多弯路。另外GEE里可用的森林高度相关数据源有好几个我做了一个对比方便你按需选数据源类型空间密度优点主要限制GEDI激光雷达脚印点稀疏但全球覆盖直接测量冠层结构精度高不是连续面域脚印之间有空隙ICESat-2光子计数激光雷达更稀疏的沿轨数据覆盖纬度极高地段夜间也能测适用于大尺度趋势区域制图难度大TanDEM-X雷达干涉生成连续面域全球连续高度模型主要反映地表少量冠层树高反演需额外建模光学影像哨兵/兰德间接反演连续面域空间连续、可多期对比需要训练样本反演不确定性高结论GEDI精度最高但输出是离散脚印面域连续成图一定要靠插值和辅助建模。项目里我选择“GEDI为真实高度样本 其他连续栅格作为辅助特征”的路线。1.2 GEDI做主角但为什么要“借助”它GEDI是美国国家航空航天局在国际空间站上搭载的激光雷达系统发射的激光脉冲打到地面后返回通过波形反演可以得到不同相对高度上的冠层能量分布。它的脚印大概直径25米左右采样密度取决于轨道分布。听起来很美但它有三个天然限制沿轨采样GEDI数据沿着卫星轨迹分布轨迹之间往往有大片空白直接加载出来就是一条条“光带”不是整片林子。云和地形影响云遮挡、坡度大、地表粗糙都会削弱信号导致部分脚印质量不合格。不是逐年均匀覆盖探测周期跟空间站轨道有关不同年份同一地区的脚印数量差异很大。所以GEDI数据在GEE里本质上是一个稀疏的、高质量的高度“实测样本集”咱们要做的是通过清洗、聚合和辅助数据建模把这些点“铺满”到整个研究区得到一份连续、可用的森林高度栅格。2. GEDI高度数据清洗比聚合更关键的一步2.1 在GEE里把GEDI脚印拉出来我先说加载。GEDI数据在GEE中通常以FeatureCollection形式存在每个Feature就是一个脚印属性里带高度值和一堆质量参数。我这里用典型的RH100数据来做演示代码大致是var region ee.FeatureCollection(你的研究区矢量); var gedi ee.FeatureCollection(GEDI/001/RH100); Map.centerObject(region, 9); Map.addLayer(gedi.limit(200), {color: red}, GEDI原始脚印);跑完这段你会发现地图上看脚印是有的但零零散散而且很多看起来高度极高或为负——这就是没做清洗的原始数据。不要指望它直接能用。注意不同GEDI产品版本或数据分发的路径可能不同加载前先print一下看看数据结构确认属性名比如RH100、sensitivity、beam别直接写死属性名。2.2 质量过滤一个脚印要不要看这四个指标这是整个数据处理里我认为价值最大的一步。GEDI脚印从发射到返回会经历大气、地形、地表覆盖等多重干扰GEE里能用的质量属性通常至少有这几个rh100冠层相对高度树高这是要用的核心值。sensitivity灵敏度表示该波形对冠层的探测能力。低于0.9的基本上无法分辨冠层结构这类脚印必须剔除。beam光束编号。GEDI有不同强度的光束弱光束power beam信噪比相对低有时需要单独考虑或者干脆用强光束。quality_flag / degrade_flag数据质量标志和降级标志。凡是标记为质量差或已降级的脚印直接排除。我自己的过滤器写下来是长这样的var gediFilt gedi .filterBounds(region) .filterDate(2020-01-01, 2022-12-31) .filter(ee.Filter.and( ee.Filter.gte(rh100, 0), // 高度非负 ee.Filter.lte(rh100, 60), // 树高不可能超过60米按区域情况调整 ee.Filter.gte(sensitivity, 0.9), ee.Filter.eq(degrade_flag, 0), ee.Filter.stringStartsWith(beam, BEAM) // 或者单独选强光束 )); print(过滤后的脚印数量, gediFilt.size());这句话看着简单但每个字段都有讲究。比如rh100上限不能拍脑袋定如果研究区是热带雨林上限可以放到70米要是研究区在比较干旱的区域上限设到40米都算宽。盲目保留超高值聚合的时候非常容易被少数异常点带偏。2.3 地形因子被忽略的“隐形杀手”除了数据本身的质量字段还有一个我强烈建议做的过滤地形坡度。原因很简单GEDI激光雷达在高坡度地区返回波形会拉伸变形反演高度误差会急剧增大。哪怕它的quality_flag显示正常面对陡坡依然要打问号。处理方式是结合SRTM或COP30做坡度计算再过滤掉坡度大的脚印var slope ee.Terrain.slope(ee.Image(NASA/SRTMGL1_003)); var gediFiltered gediFilt.map(function(f) { var s slope.sample(f.geometry(), 30).first(); return f.set(slope, s.get(slope)); }).filter(ee.Filter.lt(slope, 15)); // 坡度小于15度的保留这里30是一个采样半径参考了GEDI脚印大小。坡度阈值到底取多少建议做了敏感性测试再说比如分别用10度、15度、20度过滤看高度统计量变化大不大。如果变化剧烈说明数据对地形太敏感需要更严格。2.4 一个容易被漏掉的细节脚印点偏移GEDI脚印的坐标是波束在地面的投影中心但实际覆盖是一个大概直径25米的圆形区域。如果你拿点位置去跟高分辨率影像比较有时候会觉得“点怎么落在了路或河流上”这不是坐标错误而是脚印本身有空间不确定度和地表覆盖混合的问题。在这个项目里我采取了一个笨但有效的办法把脚印做缓冲再取其中心与土地覆盖数据判别。比如用ESA WorldCover或GlobeLand30判断脚印中心如果落到了非森林区域就把该脚印剔除。var landcover ee.Image(ESA/WorldCover/v200); var gediForest gediFiltered.map(function(f) { var lc landcover.sample(f.geometry().centroid(), 10).first(); return f.set(lc, lc.get(Map)); }).filter(ee.Filter.oneOf(lc, [10, 20])); // 10树覆盖, 20灌木等, 按分类设置这样清洗完之后剩下的脚印才算真正“能用”。3. 森林高度成图从离散脚印到连续栅格3.1 聚合方式选择均值不是唯一解数据清洗完了下一步就是把这些离散点转成影像。GEE里最直接的方法是reduceToImage就是按照像素取点的统计值。但怎么取有讲究均值mean常规做法适合脚印密度比较均匀、异常值已经被清洗掉的场景。最大值max如果想反映林冠顶部的极值可能用得上但对异常值非常敏感。百分位数percentile比如取P90表示区域内较高的冠层水平能抑制较低灌丛的影响。我实测下来的体会是如果研究区森林结构复杂直接用mean会偏低因为脚印可能打到林间的空隙用P90或P95会更贴近真实林冠高度。所以我的标准流程是同时生成mean和P90后续跟实测样地对比看哪个更适合。代码参考var heightImage gediClean .reduceToImage({ properties: [rh100], reducer: ee.Reducer.percentile([90]) }) .rename(height_p90) .float(); Map.addLayer(heightImage, {min: 0, max: 40, palette: [#2c7bb6, #ffffbf, #d7191c]}, 林冠高度P90);这里有个概念要明确reduceToImage的输出分辨率取决于输入集默认投影和后续计算如果不指定scale可能在投影转换时出现分辨率不一致的问题。我一般会在聚合后显式重投影到目标栅格的坐标系比如UTM。3.2 空间插值的另一种思路用影像作为载体GEDI点变成栅格后往往会有空值因为脚印与脚印之间本来就有间隙。很多人会想着插值但GEE里没有传统意义上的IDW或克里金“开箱即用”的全局工具不过有几种替代方案领域均值对空值像素做邻域统计填充。焦点均值用reduceNeighborhood对小范围空洞做填补。辅助数据回归填充用连续的辅助数据地形、气候、光学影像波段与GEDI高度建立机器学习回归模型预测每个像素的高度。这是我认为最严谨的办法后续我会单独写一节。如果你的需求只是做一个概览图用焦点填充就够了。代码大致这样var filled heightImage.unmask(-999) .reduceNeighborhood({ reducer: ee.Reducer.mean(), kernel: ee.Kernel.square(90, meters) }) .updateMask(heightImage.mask().not()) // 只在空洞区覆盖 .where(heightImage.mask().not(), ee.Image(-999)); var finalHeight heightImage.where(finalHeight.unmask(-999).eq(-999), filled);注意焦点填充只是“图面平滑”不是“科学插值”边界值会有一定失真别在精度要求高的项目里直接用。3.3 数据导出别等跑完再想坐标系导出连续森林高度栅格时我踩过最大的坑是坐标系选择。如果你用经纬度坐标直接导出面积量算和后续坡度分析会出问题。建议统一在区域对应的UTM投影下导出。var exportRegion region.geometry(); Export.image.toDrive({ image: finalHeight, description: forest_height_p90, folder: GEE输出, region: exportRegion, scale: 30, crs: EPSG:32650, // 根据研究区纬度带换带号 maxPixels: 1e13 });项目里我用的是30米分辨率因为GEDI脚印本身就相当于25米量级强行输出10米只会造出“看起来精细、实际上空洞百出”的假高分辨率产品。输出分辨率应匹配数据源固有尺度这是处理遥感数据最基本的原则。4. 交叉分析与辅助数据把高度数据用出价值4.1 高度与地形要素的联合分析森林高度图出来后我习惯先做一轮“高度 vs 地形”的交叉分析这能快速发现数据异常。具体做法是按坡度、坡向、海拔分层统计高度均值。比如如果某个区域内林冠高度P90在中海拔突然变成0那就很可能是GEDI脚印缺失而非真的没有树。var elev ee.Image(NASA/SRTMGL1_003); var demByHeight finalHeight.addBands(elev.rename(elevation)); var zones ee.Image(你的土地覆盖类型或高程带栅格).rename(zone); var stats finalHeight.addBands(zones).reduceRegions({ collection: region, reducer: ee.Reducer.mean().group(1), // 按zone字段分组 scale: 30 }); print(stats);这类分析的价值在于让你对数据在空间上的“行为”有直觉。比如我知道阔叶林在陡坡处常表现为高度偏低如果出图显示陡坡区平均高度反而异常高那八成是有异常脚印残留果断回头改清洗参数。4.2 结合NDWI做林分状态分析GEDI高度确实跟林分蓄积、生物量直接相关但单独高度不够。为什么因为树高跟生理状态关系不紧密一片健康的成熟林和一片受胁迫的稀疏林可能高度相同。所以我在分析阶段会叠加NDWI归一化水体差异指数这类水分敏感指标。项目里我用的哨兵2影像计算NDWIvar s2 ee.ImageCollection(COPERNICUS/S2_SR) .filterBounds(region) .filterDate(2021-01-01, 2021-12-31) .filter(ee.Filter.lt(CLOUDY_PIXEL_PERCENTAGE, 10)) .median(); var ndwi s2.normalizedDifference([B3, B8]).rename(NDWI); var heightNdwi finalHeight.addBands(ndwi); print(heightNdwi.sample({scale: 30}).first());把高度变量和NDWI放进同一个样本集后我可以用散点图或简单的统计分析来判断哪些区域是“高树但水分不足”哪些是“中等高度但水分条件好”。虽然只是皮毛级的生境描述但加在一起能为后续的碳汇分析或者森林健康评估打下基础。4.3 与机器学习数据处理衔接别让特征“打架”关于热词里总提到的“机器学习中的数据处理”我在森林高度项目里的体会是如果下一步你要用随机森林或深度网络反演高度那么特征之间的尺度一致性和空缺值处理比算法本身更决定成败。常见的做法是把GEDI高度样本作为标签影像波段、地形、气候作为特征逐像素采样后做回归。关键处理步骤包括所有特征要统一到同一个投影和分辨率。特征缺失值不能简单填0因为树高模型对0很敏感。我习惯填整个特征分布的中位数或者用掩膜剔除该像素。标签GEDI高度不要只用RH100可以同时输出RH60、RH80、RH95等在模型里做多输出反而能提高模型稳定性和可解释性。这些内容展开又是一篇文章但我先把方向指出来数据清洗是机器学习建模的前置条件模型的精度天花板从清洗时就定下了。5. 常见问题与排查技巧5.1 为什么GEDI脚印在GEE里越过滤越少很多人过滤完发现脚印数量从几万降到几百直接怀疑自己代码写错了。其实很正常。GEDI在对全球采样时本来就有多种无效脚印再加上研究区小、云雨多、时间跨度短最后剩下的高质量脚印确实不多。对策是扩大时间窗比如从2年扩到3年或者放宽非核心条件比如坡度从10度放宽到15度但要记得先做敏感性分析。我遇到过有项目区脚印太少、统计分析失去意义的情况最后方案是并入邻近区域数据或者在结论里明确说明“研究区内直接观测覆盖有限结果以趋势参考为主”。5.2 聚合结果大面积空洞怎么办聚合后的空洞有两种一种是数据处理中掩膜造成的一种是脚印本身覆盖不到。处理方式前面讲过用焦点填充或辅助回归。这里要补一句填充图一定要在元数据里留个属性标记哪些像素是实测、哪些是填充。否则后面拿到成图的人没法判断置信度工程上这是大忌。5.3 导出影像在ArcGIS/QGIS里看是黑的这个纯粹是可视化拉伸问题。GEE显示的色带是动态拉伸的导出的GeoTIFF默认按实际像元值存储。用桌面软件打开如果没设置拉伸低值区域看起来就会是一团黑。别急着改数据先在软件里做“百分比裁剪”或“线性拉伸”基本都能解决。5.4 时间轴不一致导致趋势分析出错做森林高度时序趋势分析时时间一致性非常头疼。不是每年同一地区都有足够的GEDI脚印。我在处理过程中会把多年数据按季节或年聚合但前提是每年的数据量要基本均衡。如果某年脚印极少宁愿把该年剔除也不要拿一个很不均衡的样本量去算趋势。这种统计陷阱会直接让看似平滑的趋势线失真。6. 实操心得一定要补齐的几件事第一件保留中间过程输出。我在第一次跑全流程时只保留了最终成图结果客户追问“你这个高度图哪些地方可信度高、哪些是插值出来的”我一时语塞。后来调整为导出清洗后脚印、聚合前高度点、聚合后空值掩膜三个中间图层整个分析链条都清晰了。第二件没有实测样地不要硬说精度。GEDI虽然是光达数据但它不是绝对真值。地形起伏、冠层内多次散射、密集林冠下的地面反射弱都会带来系统偏差。如果经费和野外条件允许至少拉几条样带跟GEDI脚印位置对一对算一下RMSE比任何R方和图像指标都让人安心。第三件基于实测后你会发现GEDI的RH100在某些区域往往偏高。我一开始不信后来仔细了解发现波形并未完全触底、被落叶层或灌木层拖住导致地面回波不明显算法估出来的地面位置低于真实地面位置高度自然偏大。所以如果你后续做生物量模型建议把RH100乘以一个区域系数或者直接用RH80、RH90作为因变量系统误差会小很多。至于下一步扩展我个人比较看好的方向是把GEDI高度跟长时间序列雷达后向散射比如哨兵1做联合时间分析因为雷达对冠层结构和水分变化敏感可以弥补光达脚印稀疏、不能高频监测的短板。GEE在这条路线上的潜力很大数据都在同一个平台里不需要来回下载搬运。最后分享一个小技巧每次写完GEE脚本我会抽时间把代码块整理成一个带参数的函数数据路径、时间范围、过滤阈值全部做成变量。这样下次换研究区只需要改一行参数而不是重新写一遍整个处理流程。这个习惯帮我节省了大量重复劳动也让我回看旧项目时更快定位当时筛选的细节。