Combined GEDI Data and Landsat 8 for Inversion of Forest Aboveground Biomass Based on Sequential Gaussian Condition Simulation

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

Figure 1.  Geographical location of the study area

全球生态系统动力学调查(GEDI)雷达由美国国家航天航空局(NASA)2018年12月5日发射成功，自2019年4月开始收集全球51.6° N～51.6° S之间基于足印的测量数据，轨道高度约400 km，光斑直径为25 m，沿轨光斑点间隔距离为60 m，相邻轨道间距约600 m^[23]，提供了较高的采样密度和精度。GEDI包含的数据产品种类丰富，根据数据处理阶段的不同分为L1、L2、L3、L4四个级^[24]。本研究所使用的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两个值，值为0时表示光斑信号质量差，需剔除；degrade_flag有0和1两个值，值为1时表示卫星处于降轨状态，所获取的光斑点无效，需剔除；Sensitivity是灵敏度参数，阈值为0到1，在陆地上使用的GEDI光斑点灵敏度阈值一般在0.95以上，剔除阈值小于0.95的光斑点^[21]。筛选后共有1 302个光斑点，其中80%作为训练数据集，用于获取各变量的空间分布结果，剩余的20%作为验证数据集，用于检验空间插值的精度。筛选前后的研究区GEDI光斑点分布见图2。

Figure 2.  GEDI light spot distribution maps before (a) and after (b) screening in the study area

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

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

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

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

Figure 3.  Elevation (a), slope (b) and aspect (c) maps of the study area

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

式中，Y为归一化后的值，y为初始值，$ {y}_{min} $为初始数据的最小值，$ {y}_{max} $为初始数据的最大值。

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

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

（1）空间插值精度评估

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

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

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

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

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

Figure 4.  Fitting accuracy of variogram of GEDI variables

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

Figure 5.  Flowchart of Sequential Gaussian Condition Simulation with GEDI spot variable

Figure 6.  The results of Sequential Gaussian Condition Simulation for GEDI variables

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

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

Figure 7.  Contribution proportion of modeling features’ importance

Figure 8.  Scatterplot of random forest model fit for using Landsat8 data

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

Figure 9.  Contribution proportion of modeling features’ importance

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

Figure 10.  Scatterplot of random forest model fit for combining two types of remote sensing data

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

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

Figure 11.  Spatial distribution of aboveground biomass estimation results of spruce-fir forest

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

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

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

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

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

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

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