基于序贯高斯条件模拟的GEDI数据联合Landsat8反演森林地上生物量

摘要：[目的]单一遥感技术估测森林生物量存在较大局限性，本研究旨在利用多源遥感协同技术互补激光雷达和光学遥感的优势，提高生物量估测精度。[方法]以星载激光雷达GEDI和光学遥感Landsat8数据为主要信息源，采用序贯高斯条件模拟方法实现GEDI光斑数据由“点”到“面”的空间扩展，结合地面1 38块生物量调查样地，利用随机森林回归方法估测云南省香格里拉云冷杉林的地上生物量。[结果]（1）采用序贯高斯条件模拟方法对GEDI光斑点进行空间插值，模拟的12个生物物理指标在空间分布上呈现出随机性、破碎化的特征，这与森林的空间分布聚集性非常相似，参与建模的9个指标OEC均大于0.90；（2）利用单一Landsat8光学遥感数据和地形因子构建的随机森林模型精度为：R2=0.82，RMSE=35.51 t·hm-2，P-0.77；Landsat8数据协同星载激光雷达GEDI数据构建的随机森林模型精度为：R2=0.86，RMSE=32.11 t·hm-2，P=0.80，模型精度明显提升；（3）利用多源遥感技术估测的香格里拉2019年云冷杉林地上的生物量总量为37 042 605.68 t，平均生物量为123.28 t·hm-2。[结论]基于地统计学的序贯高斯条件模拟方法考虑到研究对象的空间异质性、能克服一定的平滑效应，用于实现激光点由“点”到“面”的空间扩展是可行的。星载激光雷达GEDI与光学遥感Landsat8协同的多源遥感数据可有效填补单一遥感数据源的缺陷，提高森林生物量的估测精度，能为激光雷达联合光学遥感估测大范围、全覆盖的森林生物量提供参考。

关键词：GEDI;Landsat 8；序贯高斯条件模拟；随机森林；生物量

中图分类号：X87; S757.2 文献标识码：A 文章编号：1001-1498（2024）03-0049-12

森林生物量是森林生态系统长期物质循环与能量流动的结果，它是反映森林生产力和森林生态系统功能的重要参数，在森林的经营、监测与评价中起到重要作用。森林生物量的估测对合理保护、管理和利用森林资源有重要意义。传统的森林生物量估测主要采取实地测量法，虽然精度较高，但受到时间和空间差异的限制，需要耗费大量的人力和物力，并且对森林有一定的破坏作用。随着3S技术的迅速发展，结合遥感数据估测生物量成为一种有效的技术手段。目前，森林生物量估测大多是利用野外实测的生物量样地数据，通过分析样地数据与光学遥感因子之间的相关关系建立估测方程推广到区域尺度的生物量分布。然而，光学遥感易受天气影响，对植被的穿透性较弱，只能获取森林冠层表面的信息，而且信号容易饱和、数据存在混合像元等问题，导致模型估算精度不高。

激光雷达（LiDAR）能快速获取目标的空间三维信息、对植被的穿透力强、可实现对森林冠层垂直结构参数的精准监测，在森林参数反演、生物量估测方面具有很大优势。由于激光雷达是离散采样，很难提供全覆盖的数据，尽管它突破了光学遥感在森林应用领域的诸多缺陷，但依然存在数据获取和处理较难等问题。此外，任何单一的遥感技术都存在局限性，激光雷达技术也不例外，多传感器协同应用能够填补激光雷达在林业应用领域存在的不足，因此，激光雷达协同光学遥感以获取更高精度的森林参数迅速得到诸多研究者的应用与推广。池弘等利用GLAS星载激光雷达结合Landsat/ETM+数据估测了长白山地区森林地上生物量；边瑞利用无人机激光雷达结合Landsat数据估测了祁连山国家公园的森林生物量；巨一琳等联合机载LiDAR数据和多光谱数据估算了根河市大兴安岭生态观测站寒温带天然林的生物量。

全球生态系统动力学调查（GEDI）多波束激光雷达相比于GLAS，采样密度有较大提高，实现对森林冠层高度、生物量等参数的测量具有重要意义。本研究以典型高寒山区香格里拉市为研究区，在对GEDI波形数据处理的基础上提取大量的森林参数，利用基于地统计学的序贯高斯条件模拟方法对GEDI光斑点森林参数变量实现由“点”到“面”的空间扩展，获取全覆盖的GEDI数据，结合Landsat8-OLI数据和2016年调查的138个野外样地数据，利用随机森林方法估算香格里拉2019年云冷杉林地上生物量，评估多源遥感数据协同集成应用的优势。

1 材料与方法

1.1 研究区概况

香格里拉市位于云南省西北部“三江并流”区域，地理坐标为26°52＇～28°52＇N，99°20＇～100°19＇E，总面积11 613 km2，最高海拔为5545 m，最低海拔为1 503 m，属于山地寒温带季风气候，全年气温偏低，年平均气温仅为5.4 ℃，四季交替不明显，夏秋时节雨水较多，冬春时节较为干旱，年均降水量为268～945 mm。境内生态系统完好，森林资源丰富，森林覆盖率达74.99%，优势树种有云冷杉（Picea asperataMast. and Abies fabri （Mast.） Craib）、高山松（Pinus densata Mast.）、云南松 （Pinusyunnanensis Franch.）、高山栎（Quercussemicarpifolia Smith）等。研究区地理位置见图1。

1.2 数据来源与处理

1.2.1 星载激光雷达GEDI数据

全球生态系统动力学调查（GEDI）雷达由美国国家航天航空局（NASA）2018年12月5日发射成功，自2019年4月开始收集全球51.6°N～51.6°S之间基于足印的测量数据，轨道高度约400 km，光斑直径为25m，沿轨光斑点间隔距离为60 m，相邻轨道间距约600 m，提供了较高的采样密度和精度。GEDI包含的数据产品种类丰富，根据数据处理阶段的不同分为L1、L2、L3、L4四个级。本研究所使用的GEDI数据为2019年4-12月获取的L2B产品，它是光斑尺度的生物物理指标，包含陆地表面上每个激光足印点精确的经纬度坐标、高程、海拔、覆盖率和垂直剖面指标等，数据从NASA官方平台获取（https：//search.earthdata.nasa.gov）。

提取37个GEDI光斑数据指标，包括quality_flag、degrade_flag、sensitivity 3个质量指标，各指标的描述见表1。GEDI为波形数据，采集过程中会受到云层和其他环境因素的影响，存在部分地理定位精度差以及信号质量差的光斑，因此根据3个质量指标对初始的3 858个GEDI光斑点进行筛选。quality_flag有0和1两个值，值为O时表示光斑信号质量差，需剔除；degrade_flag有0和1两个值，值为1时表示卫星处于降轨状态，所获取的光斑点无效，需剔除；Sensitivity是灵敏度参数，阈值为0到1，在陆地上使用的GEDI光斑点灵敏度阈值一般在0.95以上，剔除阈值小于0.95的光斑点。筛选后共有1 302个光斑点，其中80%作为训练数据集，用于获取各变量的空间分布结果，剩余的20%作为验证数据集，用于检验空间插值的精度。筛选前后的研究区GEDI光斑点分布见图2。

1.2.2 Landsat8遥感影像数据

使用分辨率为30m的Landsat8-OLI数据，共需3景影像覆盖整个研究区，数据获取时间为2019年12月。遥感影像广泛被用于森林生物量的估测研究，在前期对影像预处理的基础上，提取NDVI、 RVI、 DVI、PVI、EVI、ARVI、RDVI、V13、SARVI、MSAVI、SLAVI、VIS234、B53、B345、ND43、ND563等16个植被指数用于云冷杉林地上生物量的估测。

1.2.3 野外调查样地数据

野外调查数据为2016年香格里拉市森林资源二类调查的角规控制样地数据，共有138块云冷杉样地，每块样地面积为1 hm2，包含样地点坐标、小班地类、优势树种、平均胸径、树高等参数，利用二元材积表计算样地的林分蓄积量。由于Landsat8-OLI和GEDI的数据采集时间均为2019年，所以根据《云南省主要树种材积生长率表》的材积净生长率（云杉0.85%，冷杉0.88%）推算2019年的林分蓄积量，然后利用林木生物量扩展因子法计算2019年的样地生物量。蓄积量一生物量转换模型见式（1），云冷杉生物量扩展因子及木材密度见表2，138块样地的2019年生物量统计参数见表3。为与遥感影像像元尺度相匹配，将每hm2样地生物量数据转换为900m2的样地生物量数据，样地点的空间分布见图1。

B= V·SVD·BEF（1）

式中，B为森林地上生物量，V为蓄积量，SVD为某一树种或树种组的木材密度，BEF为生物量扩展因子。

1.2.4 数字高程模型数据

数字高程模型（DEM）数据为美国国家航空航天局（NASA）和日本宇宙航空研究开发机构（JAXA）共同开发的ASTER GDEMV2数字高程数据产品，分辨率为30 m。该数据用于Landsat8影像的地形辐射校正，以及基于栅格表面算法提取坡度、坡向和高程3个地形因子（图3）。根据图3可知，香格里拉海拔较高，总体呈现出西北高东南低的趋势，海拔落差较大，境内地势陡峭，多为斜坡、陡坡和峭破，地形结构复杂。

1.3 研究方法

1.3.1 变异函数

变异函数的理论模型可通过球状模型、指数模型和高斯模型进行拟合，用于描述区域化变量的变化规律，变异函数的模型由模型类型、基台值（Sill）、变程（Range）和块金值（Nugget）唯一确定。本研究利用GS+9软件对GEDI光斑变量进行变异函数分析，采用3个模型进行预测，以模型的决定系数（R2），残差平方和（RSS）对模型进行评价。在分析前，对所有变量进行归一化处理以提高模型的收敛速度，即将所有变量的数值调至[0，1]或[-1，1]，计算方法见式（2），并对所有变量的数据结构进行正态分布性检验，将不符合正态分布的变量进行开立方变换，符合条件后方可进行变异函数分析。在前期对GEDI光斑变量进行变异函数分析和筛选的基础上，选择最优变异函数模型用于实现序贯高斯条件模拟。

Y=y-ymin/ymax-ymin（2）

式中，y为归一化后的值，y为初始值，ymin为初始数据的最小值，ymax为初始数据的最大值。

1.3.2 序贯高斯条件模拟

序贯高斯条件模拟（SGCS）以变异函数理论为基础，结合区域化变量，用随机模拟方法估计研究变量。根据已知数据构造高斯函数，将区域化随机变量Z（x）的每一个取值看作是符合高斯函数（即正态分布函数）F（x）的一次随机实现。在每一个模拟位置X处，F（x）是以n个已知数据Z （xi）（i=1，2，…，n）和此前的m-1个模拟值z （xj）（j=1，2…，m-1）为条件的累积条件概率密度函数，在模拟过程通过构建局部累计条件概率分布来实现最终结果。该方法对所预测的空间数据可能的取值结果及其概率进行度量，能在最大程度上正确反映区域化变量的空间波动性，再现真实资源特性空间变异曲线的波动。在大多数情况下，原始数据并不满足模拟过程的要求，需通过正态变换后再进行模拟，最后通过逆变换将模拟结果转换为初始值。

1.3.3 随机森林回归

随机森林（RF）是以决策树为基础的机器学习集成算法，最早由Breiman等提出。它利用Bootsrap重抽样方法从训练数据中随机抽取一部分样本，对每个Bootsrap样本进行决策树建模，然后组合多棵决策树的预测，通过投票决策的方式得出最终预测结果。在随机森林中，每棵决策树都是独立并在随机选择的子样本上进行训练的，从而有效降低过拟合的风险。该方法能够进行特征重要性分析，相比于神经网络和支持向量机等其它算法，它在分析变量关系上具有显著优势。

1.3.4 精度评估方法

（1）空间插值精度评估

序贯高斯条件模拟是基于Kriging估计的一种不确定性空间插值方法，会产生一定的估测误差。本研究采用总体估计值一致性（OEC）对GEDI各变量的序贯高斯条件模拟结果进行精度评估，即以GEDI光斑点模拟估计的总体总值（TGs）与光斑点估计的总体总值置信区间中值（TGp）的比值来衡量，OEC越接近于1，表明空间插值精度越高。各指标的计算公式如下：

式中，A为研究区总面积，a为光斑面积，n为光斑点观测数量，gci为光斑点实测值，m为像元总数，scj为第/个像元的模拟值。OEC代表模拟总值与光斑点估计总值的接近程度，二者相等时，OEC=1。

（2）生物量估测模型精度评估

模型精度的高低直接影响到生物量的估测结果，本研究基于138个云冷杉样地生物量数据，采用留一交叉验证法验证云冷杉林生物量非参数估测模型的精度，以模型的决定系数（R2）、均方根误差（RMSE）、以及总体估测精度（P）作为精度评价指标。

2 结果与分析

2.1 GEDI光斑点变量筛选及序贯高斯条件模拟

2.1.1 GEDI变量筛选

保留光斑点的quality_flag和degrade_flag两个指标只有1和0单一固定值，参与构建生物量估测模型的意义不大，因此仅用于光斑点的筛选，最终采用剩余的35个变量进行变异函数分析，各变量的模型拟合精度见图4。根据变异函数的分析结果，剔除R2小于0.5的变量，从而减小序贯高斯条件模拟方法对GEDI光斑点进行空间插值的误差，获取精度较高的空间分布结果。筛选后保留12个GEDI变量，其变异函数的拟合参数见表4。

2.1.2 GEDI变量序贯高斯条件模拟

基于GS+9软件拟合的12个GEDI变量的最优变异函数，根据模型参数（变程、块金值、偏基台值）进行简单Kriging插值，最后采用ArcGIS 10.5中“高斯地统计模拟”工具完成序贯高斯条件模拟。通过反复实验，依次设定模拟次数为10次、25次、50次、75次、100次、125次，对比分析后，发现50次模拟后像元均值的方差变化趋于稳定，因此确定模拟的次数为50次。为与GEDI光斑点的面积尺度相匹配，输出GEDI各变量模拟的栅格数据分辨率为22.15 m，并对经过归一化处理和正态变换后的变量进行反归一化处理和逆变换，最终得到真实的模拟结果。序贯高斯条件模拟的流程见图5，12个GEDI变量模拟的结果见图6。利用预留20%的实测光斑数据检验模拟结果的精度，其总体估计值一致性的检验结果见表5。

根据总体估计值一致性检验结果，pgap_theta_ error、rg_a2、rg_a6 3个变量的OEC较低，且经过检验总体估计值未在置信区间内，因此剔除OEC最低的3个变量，最终保留满足精度要求的9个变量参与构建云冷杉林生物量估测模型。

2.2 随机森林建模

2.2.1 利用Landsat8光学数据构建估测模型

利用Python 3.7软件，基于Landsat8影像的16个植被指数和3个地形因子，采用随机森林特征重要性分析方法对参与建模的19个特征进行特征重要性贡献度分析，其中对建模贡献度最大的变量是V13，贡献率为12.97%；贡献率最小的变量是VIS234，贡献率为1.52%，各特征的建模贡献度占比见图7。采用所有特征构建随机森林回归模型，模型的R2=0.82、RMSE=35.51 t·hm-2、P=0.77，模型精度拟合的散点图见图8。

2.2.2 联合GEDI和Landsat8数据构建估测模型

考虑到与样地面积和Landsat8数据的栅格大小相匹配，将基于2.1节GEDI光斑数据9个变量的序贯高斯条件模拟结果重采样成分辨率为30 m的栅格数据。针对保留的9个GEDI光斑变量、Landsat8-OLI影像的16个植被指数和3个地形因子，采用随机森林特征重要性分析方法对参与建模的所有特征进行特征重要性贡献度分析，其中对建模贡献率最大的变量是V13，贡献率为12.05%；贡献率最小的变量是RVI，贡献率为0.96%，各特征的建模贡献度占比见图9。

利用9个GEDI光斑变量联合Landsat8影像的16个植被指数特征、以及3个地形因子建立随机森林回归模型，模型的R2=0.86、RMSE=32.11t·hm-2、P=0.80，两种遥感数据联合后构建的模型精度明显提升。模型精度拟合的散点图见图10。

2.2.3 生物量估测模型地形因子贡献度分析

根据图7和图9可知，两个随机森林模型中3个地形因子的建模贡献度占比均较大，尤其以DEM较为显著，前后两个模型中的贡献度分别达到11 .18%和7.21%。香格里拉地处高寒山区，海拔落差较大，云冷杉为喜阴树种，主要分布在海拔较高的阴坡和半阴坡地带，生长趋势严重受到海拔高度的影响，垂直地带差异性强，因此DEM对模型的贡献度较大。

2.3 多源遥感数据协同的生物量估测

根据两个随机森林模型精度的对比分析结果，利用星载激光雷达GEDI数据协同Landast8光学遥感数据构建模型估测香格里拉2019年云冷杉林的地上生物量。估测的2019年云冷杉林地上生物量的总量为37 042 605.68 t，平均生物量为123.28 t·hm-2，生物量空间分布结果见图11。香格里拉地处高寒山区，平均海拔约为3 500 m，总体呈现出西北高东南低的趋势。云冷杉林主要集中分布于南部地区，根据图11可知，高生物量值并非全部集中分布于南部地区，而是分布在海拔3 600 m以上的高寒地区，这是因为冷杉林主要分布于海拔3 500～3 800 m的阴坡、半阴坡地带，云杉林主要分布于海拔3 200～3 700 m的阴坡、半阴坡以及沟谷地带。此外，香格里拉人口较多的乡镇主要分布在南部的建塘镇、虎跳峡镇、金江镇和三坝乡，受人文因素和自然因素的影响，云冷杉林地上生物量总体分布不均，区域差异较大。

3 讨论

3.1 序贯高斯条件模拟方法对激光点空间扩展的应用分析

地统计学是研究空间变异及格局的有效方法，最初被应用于地质研究领域，随着统计学的发展，以及森林资源变化监测的重要性日益显著，该方法也被广泛应用于林业研究领域。星载激光雷达GEDI只能提供离散的采样点，不具有成像性，所以研究采用基于地统计学的序贯高斯条件模拟方法实现了GEDI光斑由“点”到“面”的空间扩展，获取与地面调查样地对应的变量指标数据。相比于Kriging插值，该方法克服了一定的平滑效应，能够在最大程度上真实地展现研究对象的空间分布格局及变化。根据1 2个GEDI变量的模拟结果（图6）可知，每个生物物理指标都呈现出破碎化、随机性的特征，这与香格里拉森林空间分布聚集性非常相似，且参与建模的9个变量OEC均在0.9以上，精度满足要求，因此，序贯高斯条件模拟方法可用于激光点的空间扩展，为离散的“点”数据转向大区域尺度的“面”数据的应用提供了可靠的技术支持。

3.2 云冷杉材积生长率不确定性及驱动因素分析

本研究参考《云南省主要树种材积生长率表》云杉和冷杉的材积净生长率推算2019年云冷杉林地上生物量，该表按15个龄级分别计算各龄级的材积净生长率，其中云杉I龄级的材积净生长率为8.84%，XV龄级的净生长率为0.35%；冷杉I龄级的材积净生长率为14.64%，XV龄级的净生长率为0.12%，两个树种不同龄级的材积净生长率差距较大，因减少成本和受其他因素限制，统一采用总平均净生长率，相关研究表明，利用该表给出的净生长率计算的材积生长量结果偏小。根据2016年二类调查的林木材积推算2019年的林木材积，期间有3年的时间跨度，近年来气候变化多样，时间跨度导致的气候因素会影响林木的生长趋势，进而影响生物量的生长量，贾勃等的研究发现，生物量的生长量受到年均温度的正向影响；何潇等的研究也表明气候因子对林分生物量有显著影响。此外，研究利用生物量扩展因子法计算2019年云冷杉林的地上生物量，三重不确定性因素会导致估测结果存在微弱的偏差。胥丽等利用该方法获取了香格里拉2019年的栎类森林地上生物量精度较高的估测结果，但后续研究中，若条件允许建议采用实地验证法对估测结果进行检验，并对林木材积生长率变化气候驱动因素进行分析，降低生长率不确定性因素导致的生物量估测偏差。

3.3 生物量估测精度分析及模型优化建议

针对生物量的估测结果，本研究与前人的结果做了一些对比。谢福明等基于KNN模型估测了香格里拉2016年云冷杉林地上生物量，RMSE为41.64 t·hm-2．而本研究基于随机森林模型的RMSE为32.11 t·hm-2，估测精度明显提升。根据云冷杉材积净生长率，以及生物量扩展因子法，基于2016年香格里拉森林资源二类调查的小班数据计算了2019年云冷杉林地上生物量，其总生物量为39 961 196.39 t，138块样地的2019年生物量平均值为125.77 t·hm-2。相比于二者的结果，本研究的估测结果具有一定可靠性，但在精度方面依然有可提升的空间。非参数模型是一种黑箱操作原理，它是直接或间接地从实际系统的实验分析中得到响应，机械式的系统化难以提高模型的预测精度，为提高模型的预测精度，获取更加精准的生物量空间分布制图，可对模型进行优化。宋涵玥等采用随机搜索法（Random Searching）对模型参数进行优化，优化后的生物量估测精度提高2.01%。此外，还可采用遗传优化算法、粒子群优化算法和贝叶斯优化方法对模型进行优化。

4 结论

本研究采用星载激光雷达GEDI数据协同Landsat8-OLI影像估测了滇西北生态系统较为完善的高寒山区香格里拉云冷杉林的地上生物量，评估多源遥感数据协同估测森林生物量的优势。主要结论如下：

（1）利用地统计学GS+9软件对GEDI观测点做变异函数分析，经第一次筛选的12个GEDI变量变异函数的最优模型均为指数模型，模型决定系数均大于0.5。基于变异函数的分析结果，获取12个变量的序贯高斯条件模拟空间分布，采用OEC检验其模拟精度，并进行第二次筛选，筛选后9个GEDI变量的OEC均高于0.9，较接近于1，且模拟估计的总体总值均在置信区间内，满足精度要求，基于地统计学的序贯高斯条件模拟方法可实现GEDI光斑点由“点”到“面”的空间扩展。

（2）利用单一光学遥感Landsat8数据的1 6个植被指数和3个地形因子构建生物量估测模型，模型的R2=0.82，RMSE=35.51 t·hm-2，P-0.77；星载激光雷达GEDI数据协同光学遥感Landsat8数据构建的生物量估测模型R2=0.86，RMSE=32.11t·hm-2，P=0.80，精度明显提升。多源遥感数据协同能够互补不同的数据优势，提高森林生物量的估测精度。

（3）利用2019年的两种遥感数据，结合2016年野外调查数据，根据材积净生长率和生物量扩展因子法获取138块云冷杉样地2019年的生物量，利用随机森林回归方法估测的香格里拉2019年云冷杉林的地上生物量总量为37 042 605.68 t，平均生物量为123.28 t·hm-2，主要分布在38.54～282.19 t·hm-2之间。

（责任编辑：彭南轩）

基金项目：云南省农业联合专项一重点项目（20230IBD070001-002）;云南省教育厅科学研究基金项目（2023Y0728），