杨泽至，舒清态，邱 霜，席 磊

(西南林业大学 林学院，云南 昆明 650224)

叶面积指数(leaf area index,LAI)作为影响植被光合作用、呼吸作用、蒸腾作用和碳氮水循环等生理过程的关键因素以及反映植被净初级生产力和物种多样性水平的重要因子，是利用大气物理模型模拟植被与气候相互作用关系的关键输入参数[1-4].大尺度森林LAI 分布信息对于深入揭示地表植被的时间变化对气候变化的响应机制，进而深化对植被与气候相互作用关系的认识至关重要[1].因此，如何大面积科学、高效地获取森林LAI 空间分布信息是地学领域长期以来的研究热点.传统的LAI地面直调查方法虽然能够取得较高的测量精度，但存在耗时长、成本高、仅限于“点状”尺度等难以克服的固有缺陷，难以获取LAI 的空间格局信息，阻碍了LAI 的具体应用[5].当前，遥感技术凭借其快速获取大范围空间连续数据的技术优势已经成为快速提取森林LAI 空间分布信息的主流手段.包括光学遥感、机载激光雷达和合成孔径雷达等在内的多种遥感技术在植被LAI 定量反演中得到了不同程度的应用，尤其以基于被动光学数据的LAI估测应用最为广泛[5].然而，基于以上遥感数据的森林LAI 的估算方法均存在各自的局限性.光学遥感虽然能获取大尺度的空间连续数据，被广泛应用于多种生境条件下的LAI 估测，但植被指数等光谱变量常出现光饱和效应，进而容易导致LAI 的高值低估[6-7].机载激光雷达作为当前快速提取森林结构参数的先进遥感手段，尽管具有高精度获取不同高度的植被结构信息的独特优势，可实现LAI 的高精度反演，但其硬件昂贵、数据获取成本高，且通常只能获取范围较小的带状数据，使其技术推广受到制约[8-10].合成孔径雷达虽然不受气象条件限制，能够全天时进行植被LAI 建模，但模型精度易受到传感器特性、冠层结构特征及地表性质等诸多因素的影响，因而不同程度限制了该技术的广泛应用[5].随着遥感技术进一步发展，星载激光雷达为森林结构参数大尺度高效提取提供了新的技术途径[10].然而，以星载激光雷达作为主要数据源进行森林LAI 估测的研究尚不多见，缺乏对星载激光雷达估算森林LAI 潜力和效果的充分探讨.

新一代激光雷达卫星ICESat-2（Ice，Cloud，and land Elevation Satellite-2）搭载的传感器ATLAS（Advanced Terrain Laser Altimeter System）是全球首个光子计数激光雷达.作为新型星载激光雷达系统，与全波形激光系统GLAS 相比，ATLAS 光子数据具有更小的光斑直径（17 m）、更高的脉冲频率（10 kHz）和光子密度，在森林资源监测中具有重要的应用潜力[11].目前，ICESat-2 数据已逐渐在极地冰盖与海冰高程检测、地表高程测量、湖泊水位动态监测以及森林结构参数提取等多个领域发挥着关键作用[12].在森林结构参数提取方面，ICESat-2 数据的应用研究主要集中在区域尺度的林分树高分布制图和森林地上生物量反演两个方向[12-15]，而应用于森林LAI 估测的研究非常有限，鲜见报道.因此，基于ICESat-2 数据估测森林LAI 的有效性和潜力仍有待深入探讨.此外，与光学影像和机载LiDAR 点云数据不同，ICESat-2 光子数据缺乏空间连续性，无法直接形成连续的面状或带状数据，而必须借助一定数据外推方法才能使目标反演参数由点状分布拓展到面状分布，进而形成其预测性空间分布.目前，实现遥感变量或反演结果外推的方法主要包括两大类，即建立反演结果与卫星光学影像特征因子的统计关系[16-17]和克里格插值[18].其中，克里格插值在遥感变量及其反演结果外推中的应用较少，作为对未抽样点的线性无偏最优估计[19]，其不仅可以有效克服光学影像所固有的易饱和缺陷，而且能够实现ICESat-2 光子数据的充分利用.

综上，利用ICESat-2/ATLAS 数据进行森林结构参数反演具有重要的应用潜力，其在大尺度森林LAI 估测中的有效性需待深入探讨.对此，本文提出了基于ICESat-2 数据并结合克里格插值，实现森林LAI 连续估测的方法.其基本思路为：结合LAI 样地的对应ATLAS 光斑特征因子与其47 个LAI 实测数据，分别利用逐步回归和随机森林回归算法构建森林LAI 估测模型，随后通过克里格插值实现光斑特征因子属性的尺度外推，进而实现基于ICESat-2/ATLAS 数据的研究区森林LAI 的整体预测.基本研究过程如图1 所示.

图1 主要研究思路示意图Fig.1 Flow chart of main research ideas

1 数据与方法

1.1 研究区概况本研究以地处滇西北的香格里拉市为试验区.香格里拉市位于云南省西北部、迪庆藏族自治州中东部，西邻德钦县和维西傈僳族自治县，南接丽江市，东部与四川省稻城、木里二县毗邻，跨越99°8′2″E～100°21′15″E 和26°49′7″N～28°54′38″N，是青藏高原向云贵高原的过渡地带和横断山区的腹地（图2）.香格里拉市境内地形复杂多样，以山地地形为主，地势中间低，四周高，平均海拔3 000 m 以上，多种气候带镶嵌交错，土壤类型多样.该地区森林覆盖率高，地处亚热带常绿阔叶林植被带与青藏高原高寒植被带的交错地带，低纬高原的地理特征造就了其显著的植被垂直分异特征，主要树种为高山松Pinus densata、冷杉Abies fabri、云杉Piceaasperata、和落叶松Larix gmelinii等，并伴有杨树Populus tremula等伴生树种，尤其以高山松分布最为广泛，是我国重要的生态单元之一，具有重要的研究意义.

图2 研究区位置及LAI 样地分布Fig.2 Location of the study area and distribution of LAI sample sites

1.2 数据获取及处理研究数据主要包括ICESat-2/ATLAS 的2 种数据产品：全球地理定位产品ALTA03 数据陆地与植被冠层高度产品ATLA 08数据；试验区ATLAS 光斑尺度下的森林LAI 样地数据.

1.2.1 ICESat-2/ATLAS 数据获取 本研究采用ATLAS 采集的ATLA 03 和ATLA 08 数据作为研究区森林LAI 估测的遥感数据源，以探讨其在区域尺度的森林LAI 估测中的潜力.目前，ATLAS 按照4 个级别（Level0～Level3）已生成了共计21 种产品数据集，依次编号为ATL00～ATL21[20-21].研究所采用的覆盖整个试验区的ATL03、ATL08 产品数据集均由美国国家冰雪中心官网（https://nsidc.org/data/ ATL03/versions/1，https://nsidc.org/data/ATL08/versions/1https://nsidc.org/data/icesat-2/datasets）下载得到，并与试验区LAI 样地实测数据保持时间同步.为覆盖研究区整个林地区域，本试验获取了2020 年1 月至2021 年6 月之间香格里拉市境内的所有ATL03 和ATL08 产品数据集.其中，ATL03 数据共计118 条记录，含354 条轨道、708条光子轨道波束；ATL08 数据共计118 条记录，含354 条轨道、708 条光子轨道波束.

1.2.2 研究区ATLO8 光子数据去噪和分类处理

（1）光子点云去噪 ATLAS 光子数据在采集过程中易受太阳辐射、大气散射、地形条件和传感器特性等多种因素影响，包含大量非植被光子的背景噪声.因此，需进行强背景噪声的过滤，才能确保光子分类和地表信息提取的有效性[22-23].研究表明，基于改进OPTICS 的光子点云去噪算法具有较强的自适应性，且对聚类参数不敏感，尤其适用于存在强背景噪声和地形条件复杂的山地，并在分离信号光子与噪声光子、估算地面高程及树高方面表现均优于改进的DBSCAN 算法[12,24].本研究区地形复杂，山地广布，严重影响信号光子信息的高精度提取.因此，本文采用基于改进OPTICS 的光子点云去噪算法对研究区的ATL08 数据进行噪声滤除，以有效提取森林信号光子.该算法很好利用了噪声光子点和信号光子点在空间分布规律上的显著差异，具有较强机理性.其去噪过程如图3所示.

图3 基于改进OPTICS 的光子点云去噪过程[13]Fig.3 Photonic point cloud denoising process based on improved OPTICS[13]

（2）光子点云分类 经过噪声滤除后的ATLAS 光子数据包含了植被冠层光子和地面光子等信号光子.因此，需对信号光子进行分类处理，以有效提取森林冠层光子.研究显示，经改进的渐进三角不规则网络加密[progressive triangular irregular network (TIN) densification，PTD]光子点云分类算法在去除低层物体上的非地面点和在复杂环境中保留地面特征点方面比经典PTD 方法表现得更好.与经典的PTD 算法相比，改进型PTD 方法具有更低的分类误差和更高的光子分类精度，对地形条件复杂的区域具有较强的适用性[25-26].因此，本研究基于改进的PTD 光子分类算法对研究区去噪后的ATLAS 光子数据进行分类处理，进而有效提取植被冠层光子.该算法包括3 个核心步骤：参数指定、种子点选择和初始TIN 构造以及TIN 的迭代加密.

通过光子去噪和分类处理后，获得覆盖香格里拉市的上千万个ATL03 产品的有效光子，并基于ATL08 产品的数据组织结构，对所获取的ATL03产品光子进行100 m 的抽稀采样，进而提取得到研究区内有效光斑94 086 个.随后根据研究区最新林地范围Shp 文件提取得到覆盖研究区内有效林地光斑74 230 个，非林地有效光斑19 856 个.最终的研究区ATLAS 光斑分布如图4 所示.

图4 研究区林地ATLAS 光斑分布Fig.4 ATLAS spot distribution in the study area

1.2.3 特征因子提取与优选 本文采用Python语言，基于Jupyter 集成开发环境建立了ICESat-2/ATLAS 特征参数提取模块程序，以提取模型特征因子，并采用皮尔逊相关分析法对ICESat-2/ATLAS-LAI 模型因子进行优选.

1.2.4 样地数据采集及处理

（1）LAI 的样地设计 为了与ATLAS 光斑位置和大小完全相一致，本研究采用光斑脚印在地面上的投影作为调查样地，即LAI 样地是半径为8.5 m 的样圆，并以ATLAS光斑的圆心作为样地的几何中心.在调查过程中利用千寻星矩SR3 差分定位仪对光斑中心点坐标进行精确采集，以确保LAI 样地与ATLAS 光斑在空间上的重合.本次样地调查采用分层随机抽样法进行样地布设，共调查了覆盖研究区主要森林类型的LAI 样地54 块，整个样本数据能够有效反映研究区森林类型的分布情况.

（2）样地实测LAI 的计算 采用半球摄影法实现样地LAI 实测数据的获取[10].首先确保光照条件稳定，天色光较为均匀，且无阳光直射，随后对每个样地均匀选取13 个点位，并逐一在每个点位上将鱼眼镜头呈竖直向上各取1 张林冠鱼眼半球照片，形成13 张鱼眼照片，即每个LAI 样地均对应13 张林冠半球照片，以此作为计算样地LAI 值的基础.对照片质量合格的47 块样地的鱼眼照片进行预处理，分别统一照片行列号数和格式为3 000×4 000 和JPG 格式.最后，利用基于MATLAB 程序开发的 CANEYE V6.1 的半球图像分析软件对47块有效样地的鱼眼照片进行LAI 计算，进而得到47 块样地的LAI 实测值，其描述统计特征见表1.

表1 有效样地实测LAI 描述统计Tab.1 Descriptive statistics of measured LAI in effective sample plots

1.3 研究方法

1.3.1 建模方法 研究分别利用成熟、稳健的逐步回归和随机森林回归[27-30]（random forest regression，RFR）模型构建森林LAI 估测模型.

1.3.2 变异函数 变异函数又称半方差函数，其定量描述了地理数据的空间自相关性，可以表征区域化变量的空间分布格局，是地质统计学的核心工具和克里格插值的理论基础[19].在地学变量满足本征假设和二阶平稳的假设下，变异函数可定义为区域化变量在对应采样间隔下样本方差的数学期望，其函数表达式为：

式中:γ(h)为变异函数，h为成对样点的采样间距，n(h)是当样点间距为h时的样点总对数，z(xi)和z(xi+h)分别为区域化变量在空间位置xi和xi+h处的属性值.

经变异函数计算得到的实验变异函数需通过球状模型、指数模型线性模型和高斯模型等方法进行拟合，建立其理论变异函数模型，从而为地学变量的结构分析及其空间局部估计即克里格插值提供理论基础[19].变异函数理论模型的结构特征由函数类型、块金方差（Nugget variance，记作C0）、基台值（Sill，记作C0+C，C为偏基台值）和变程（Range，记作a）等4 个参数来表征.其中，块均方差指样本间距h=0 时的变异函数值，来自测量误差或局部微小变异；基台值是当采样间距增加到某一临界值时变异函数趋于稳定而所对应的常数；而变程则是变异函数达到基台值时的样本间隔距离，是地学变量空间自相关性质的作用距离.

1.3.3 克里格插值 克里格插值是建立在地学变量数据分布及其变异函数结构分析基础上，对未知点属性值的无偏最优估计.其插值过程中既充分考虑未知点与所有已知样本点之间的空间位置关系，又充分利用各样本点数据的空间结构特征，对每个已知点赋予一定的插值权重，并采用加权平均的方法对未知点的值进行预测[19,31].本研究采用克里格插值法，以光斑采样尺度（17 m×17 m）为插值尺度，对拟参与LAI 建模的ICESat-2 光子特征参数进行线性无偏最优估计，进而形成对应的光子特征参数的栅格图层.其插值表达式如下：

式中：Z0为拟参与建模的ICESat-2 光子特征参数在未知点的属性值，λi为参与插值的样地点对估测点属性的贡献权重，Xi为样地点位置，Z(Xi)为ICESat-2 光子特征参数在位置Xi处的观测值.在确定已知点的贡献权重时不仅要参照估测点和已知点之间的距离，还需考虑到已知点属性值及其空间分布格局.

1.3.4 模型评价方法 以可决系数R2和均方根误差(root mean square error,RSME)作为模型精度评价指标，并采用留一交叉验证法（leave-one-out method of cross-validation,LOOCV）对模型性能进行验证评估，以有效避免建模数据集和测试数据集划分过程中LAI 建模样本选取的随机性和不确定性[32]，提高小样本建模结果评价的可靠性.相关统计指标的表达式如下所示：

式中:yi为因变量实测值，为模型的预测值，为因变量实测值的均值，n为样本数.

2 结果与分析

研究结果包括模型因子相关性分析、森林LAI 估测模型的建立、LAI 模型特征因子克里格插值和研究区森林LAI 空间分布预测等4 个部分.

2.1 ATLAS08 光斑特征参数与样地实测LAI 相关性分析在提取47 个LAI 样地对应的ICESat-2 光斑特征参数的基础上，通过SPSS 25.0 版统计软件，采用皮尔逊相关分析法（Pearson correlation analysis）对47 个样地中的LAI 实测值与其对应的56 个光斑特征因子进行线性相关分析，以从中优选出与样地实测LAI 相关性较为显著的光斑特征因子作为LAI 的建模因子，其优选结果见表2.Pearson 相关性分析结果显示,与样地实测LAI 相关性最高光斑特征因子为Landsat percent，相关系数为0.346 (P＜0.05)；其次为冠层光子比例（Canopy photon rate），相关系数为0.327（P＜0.05）；相关性最低光斑特征因子为信噪比SNR，相关系数为0.01.共有4 个光斑特征因子与样实测LAI 的相关性较为显著，按其相关系数由大至小依次为排列为：Landsat percent ＞ Canopy photon rate ＞ Minimum canopy＞ Number of canopy photons.因此，本文将上述4个ICESat-2 光斑特征因子作为LAI 的模型因子.

表2 ATLAS 特征因子优选结果Tab.2 ATLAS feature factor preference results

2.2 LAI 模型预测结果与分析根据2.1 节的光斑特征因子优选结果，以Landsat percent、Canopy photon rate、Minimum canopy height 和Number of canopy photons 等4 个与实测LAI 相关性相对显著ICESat-2/ATLAS 特征变量为建模因子，以样地实测LAI 为模型因变量，分别采用逐步线性回归和随机森林回归算法，建立基于ICESat-2 光子数据的森林LAI 遥感估测模型.其建模结果见表3～表7.

2.2.1 基于逐步回归的LAI 估测模型 利用Minitab 统计分析软件，以经优选的4 个光斑特征因子为模型自变量，以47 个样地实测LAI 真值为响应变量，建立森林LAI 估测模型，以检测ICESat-2 光子数据与森林LAI 的线性关系水平.其建模结果表明（表3），逐步回归分析一共产生了3 个LAI线性估测模型.最高LAI 拟合精度为0.389，表明入选建模的ICESat-2 光斑特征因子能够解释样地LAI 的38.9%的方差变化.由其他因素引起LAI 方差变化占到了61.1%.由表4 并结合模型拟合优度及显著性可知，LAI 的最优线性模型为：LAI=-0.074+0.006×Landsat percent+0.242×Canopy photon rate+0.11×Minimum canopy height，其LOOCV 评价见表5.

表4 LAI 逐步回归模型参数Tab.4 Parameters of LAI stepwise regression model

表5 LAI 逐步回归模型留一交叉评价Tab.5 LAI stepwise regression model leave-one-out crossevaluation

2.2.2 基于随机森林回归的LAI 估测模型 为探究森林LAI 与ICESat-2/ATLAS 光斑特征因子的非线性映射关系，研究借助RFR 算法构建了基于ICESat-2 光斑特征因子的森林LAI 反演模型.基于2.1 节的光斑特征因子优选结果，结合47 份LAI-ICESat-2 样本数据，利用Python 语言，借助Jupyter book 平台，通过Scikit-learn 机器学习库（https://scikit-learn.org.cn）调用RandomForestRegressor 模块实现了LAI 的随机森林建模.其模型训练的关键超参数设置见表6.结果表明，留一交叉验证下的LAI 随机森林模型的R2和RMSE 分别为0.557 和0.175（图5、图6）.可见，RFR 模型在整体上较好地拟合了ICESat-2 光斑特征因子与LAI 的关系.

表6 LAI 随机森林回归模型参数取值Tab.6 Parameter values of LAI random forest regression model

图5 RFR 预测LAI 的效果点线图Fig.5 Point line diagram of the effect of RFR to predict LAI

图6 基于RFR 模型的LAI 预测精度Fig.6 Accuracy of LAI prediction based on RFR model

表7 显示，留一交叉验证下，LAI 的逐步线性回归模型和随机森林回归模型的可决系数R2分别为0.187 和0.557；RMSE 分别为0.231 和0.175；随机森林回归模型对LAI 的整体预测效果优于逐步回归模型.因此，本研究以随机森林回归模型作为估测整个研究区森林LAI 分布格局的建模方法.

表7 LAI 估测模型精度对比Tab.7 Comparison LAI estimation model accuracy

2.3 LAI 模型光斑特征因子克里格插值结果获取研究区ICESat-2/ALTAS 光子特征参数的连续分布数据是实现研究区森林LAI 空间分布制图的基础条件.为此，本文基于光斑特征因子优选结果，采用变异函数对经优选并参与LAI 建模的4 个光斑特征因子（表8）进行空间结构分析，并通过克里格插值对4 个光斑特征因子的属性值进行连续内插，继而得到4 个光斑特征因子的栅格图层.本节的变异函数分析通过GS+9.0 软件实现，该软件为用户提供了地学变量的变异函数值计算及其理论模型的拟合工具.4 模型因子的变异函数结构分析及其克里格插值结果分别见表9 和图7.

表8 待插值的4 个光斑特征因子Tab.8 The 4 spot feature factors to be interpolated

表9 LAI 模型光斑特征因子变异函数结构分析Tab.9 Structural analysis of the variance function of the LAI model light spot feature factor

图7 研究区LAI 光斑特征因子属性空间分布Fig.7 Spatial distribution of LAI spot feature factor attributes in the study area

变异函数结构分析表明，上述4 个ICESat-2 光子特征属性具有高度的空间异质性,且表现出不同的变异特征.研究区4 个光斑特征因子中除了冠层光子数（Number of canopy photons）以外，其余3 个光斑特征因子的变异函数最优拟合模型均为指数模型，其可决系数R2均在0.75 以上，且残差平方RSS 和在4 个变异函数理论模型中均为最小值.对于Number of canopy photons，其球状模型和高斯模型的决定系数R2均为0.68，其高斯模型的残差略高于球状模型，因而球状模型是其变异特征的最佳表征模型.此外，块金方差与基台值之比[C0/(C0+C)]是描述地学变量总空间变异性基本构成的度量，直接反映了随机因素对总空间异质性的影响强度，亦称块金效应.块金效应可分3 个水平：[75%，100%]；[25%，75%）；[0，25%），依次分别代表地学变量呈强空间自相关、中等空间自相关和弱空间自相关.此效应越强，表明随机因素对总空间异质性的影响越大，反之，表明总空间异质性中由变量空间自相关性引起的部分占比越高.块金效应表明，4 个光斑特征因子在研究区内由其空间自相关性引起的异质性在总空间异质性的占比均达80%以上，而由随机因素引起的空间变异性仅占到了总空间异质性的约6%～12%.4 个光斑特征因子在研究区内的空间自相关作用距离，即变程a处于5 700～7 500 m.可见，研究区4 个光斑特征因子的空间分布特征主要由ICESat-2 光子数据本身的结构性因素（如地、植被覆盖和冠层结构特征等）控制，即空间结构特性是形成其异质性的决定性因素，而受随机因素影响较小.

4 个光斑特征因子属性的变异特征分析显示，LAI 模型的4 个光斑特征因子属性均具有显著的空间自相关性，具备对其进行克里格插值的充分条件.因此，本文在变异结函数结果分析的基础上，借助ArcGIS 10.7 中的地统计分析（Geostatistical Analyst）模块分别对Landsat percent、Canopy photon rate、Minimum canopy height 和Number of canopy photons 等4 个光斑特征因子的属性进行克里格插值，进而形成4 个特征因子的栅格属性图层，如图7所示.

2.4 研究区森林LAI 空间分布预测基于2.2 节LAI 建模结果和2.3 节的ICESat-2 光斑特征因子属性的克里格内插结果进行香格里拉市的森林LAI 的空间分布预测.以研究区林地范围内的光斑特征因子连续分布数据为自变量，利用RFR 模型逐像元地对研究区森林LAI 进行预测，并以17 m输出分辨率，以确保LAI 的制图尺度与ICESat-2的光斑采样尺度一致，进而实现研究区森林LAI的连续制图（图8）.

图8 研究区森林LAI 预测分布图Fig.8 Forest LAI prediction distribution map in the study area

图8 所示的LAI 反演结果较好地反映了研究区内森林LAI 的空间分布，直观、有效地展示了森林LAI 的空间分布差异，与ICESat-2 光斑特征因子一样，具有高度的空间异质性.总体而言，研究区的LAI 整体分布不均衡，整体呈现“中间高、四周低”的分布格局.LAI 高值区集中分布在香格里拉市的中部偏北区域和南部区域，且多为森林腹地，而LAI 低值区主要分布在其中部、与城镇毗邻的林地区域，符合地理学第一定律（空间自相关性）和地理学第二定律（空间异质性）.

3 讨论与结论

3.1 讨论

（1）基于ICESat-2/ATLAS 光子数据估测森林LAI 具有重要的应用潜力.本研究基于新一代星载激光雷达ICESat-2 光子数据进行森林LAI 估测，以探讨星载光子数据在森林LAI 建模中的潜力，是区域尺度森林LAI 反演方法的新尝试.在横向对比上，本研究得到的研究区森林LAI 空间分布格局与谭德宏等[33]和赵洪莹等[34]等利用卫星光学数据对同一区域森林LAI 的预测分布高度一致，尤其在LAI 低值区更为吻合.这从侧面反映了ICESat-2 光子数据在森林LAI 估测中的可行性.此外，与Zhang 等[35]同样利用ICESat-2 数据并结合Beer-Lambert 定律进行大尺度森林LAI 反演不同，本文以林地ICESat-2 光斑内的森林LAI 的实测样本与ICESat-2 光子数据的统计关系进行建模，估测结果能够进行更为科学的地面验证，因而具有显著的生态学意义.在纵向对比上，ATLAS 光子计数激光雷达的数据采集性能得到显著提升，能获得直径更小且密度远大于第一代星载激光雷达GLAS全波形数据的光斑数据，为大尺度森林结构参数的快速提取提供了更为有效的数据支撑[11].早期研究，如骆社周等[16]利用直径为70 m 的GLAS 光斑数据，基于LPI 激光穿透指数建立LAI 估测模型，并结合TM 光学影像进行了大尺度森林LAI 反演，是基于星载激光雷达数据反演区域森林LAI 的一次重要尝试.然而，该研究受制GLAS 全波形数据自身特性，存在地面采样成本较高，反演精度地域性强等问题，其方法技术的可移植性不强.

（2）随机森林回归模型对LAI 的估测效果优于逐步回归模型.随机森林和逐步回归都是性能稳健、并被广泛应用的建模方法.本文分别采用逐步回归和随机森林模型进行研究区森林LAI 建模，并对比两种模型对LAI 的预测效果.留一交叉验证显示，逐步回归和随机森林回归模型对研究区森林LAI的预测精度R2分别为0.187 和0.557，RMSE 分别为2.34 和1.75，可见随机森林回归模型的预测效果优于逐步回归.然而，本文应用的两种模型对研究区森林LAI 的估测精度均不高，其主要原因在于LAI 实测样本容量有限，限制了随机森林模型预测性能的充分发挥，不易深入发现ICESat-2 光子数据与实测LAI 的内在关系.然而，较大的实测样本量，必然以大量财力物力的消耗为代价，降低了研究的可行性与遥感估测模型的实用性.因此，如何在有限条件下，充分发掘遥感数据与LAI 实测值的潜在关系、积极探讨新兴人工智能方法（如梯度提升GBTD[36]、极端随机森林Extree[37]或深度学习算法[38]）在基于ICESat-2 数据估测森林LAI中的潜力将是今后的研究重点.

（3）克里格插值法能够充分利用ICESat-2 数据对LAI 估测结果进行外推.在利用星载激光雷达进行森林LAI 反演的早期研究中，有关学者借助卫星光学影像实现基于星载激光雷达的反演结果的尺度上推，以使反演结果呈连续分布.如骆社周等[16]利用TM 影像将基于GLAS 波形数据反演得到的森林LAI 外推至试验区，从而形成了LAI 的预测分布图.Narine 等[15]借助Landsat-8 影像将基于ICESat-2 数据的AGB 反演结果尺度上推至整个研究区.然而，此类方法对估测模型的预测精度要求较高，且会造成估测结果外推过程中预测误差的二次传递，并带来光学数据的饱和问题，适用范围有限.针对上述问题，本研究基于变异函数结构分析，对经优选的4 个呈离散分布的ATLAS 光斑特征因子进行克里格插值，进而形成其各自属性覆盖研究区的连续数据，并通过RFR 实现了研究区森林LAI 的连续估测.本文方法不仅充分利用了ICESat-2 光子数据，而且有效避免了借助卫星光学影像对LAI 估测结果进行尺度外推过程中新的误差引入以及其带来的光饱和效应.此外，研究仅采用了普通克里格法对4 个参与LAI 建模的光斑特征因子进行插值，而未对多种克里格法的插值效果进行比较，其它克里格插值法，如协同克里格、泛克里格法在光斑特征因子插值中的效果需进一步探讨.

此外，ATLAS 的强、弱波束光子在森林LAI预测中的效果需分别探究；RFR 模型参数的优化与否对LAI 模型精度是否有显著提升也需进一步探究，是今后研究的重点之一.

3.2 结论针对传统光学遥感易受光饱和影响、机载激光雷达与微波数据获取成本高等问题，本文提出了一种基于ICESat-2/ATLAS 光子数据进行森林LAI 估测的新方法.该方法在研究区林地光子点云数据去噪、分类等预处理的基础上，提取了各个林地光斑的特征因子，并结合47 份森林LAI 样本数据，建立了基于样地光斑尺度的森林LAI 随机森林估测模型，其留一交叉验证下的R2和RMSE分别为0.557 和1.75.由此表明，ICESat-2 光子数据在区域尺度的森林LAI 预测中具有重要潜力.为实现研究区森林LAI 的连续估测，本文基于变异函数结构分析，对经优选的4 个参与LAI 建模的光斑特征因子进行克里格插值，形成其各自属性的对应栅格图层，并利用已训练的RFR 模型实现了研究区森林LAI 的整体估测，继而形成了17 m 空间分辨率的研究区森林LAI 空间分布图.本研究估测得到的森林LAI 在空间梯度上与同类研究以及研究区森林覆盖度基本一致，结论可信，是以较小成本实现大尺度森林LAI 快速反演的一个新思路.