叶面积指数（Leaf Area Index，LAI）作为反映作物群体长势特征的关键指标，已非常广泛地应用于作物长势监测^[1]，它通常被定义为单位地表面积上的单面叶片的面积之和，与作物光合作用、呼吸作用和水分利用等过程密切相关^[2-4]。目前，获取LAI的方式主要分为两大类：田间实测和遥感估算。田间实测方法有方格法、称重法、扫描仪法等，测量相对准确，但由于只能获取有限的采样点数据且费时费力、有一定破坏性因而在实际应用中受到局限^[5]，而借助遥感技术可以快速、无损地估算大尺度范围的LAI，尤其是无人机遥感技术，相比于卫星遥感、航空遥感等高空遥感技术，它具有便捷、高时效性、空间分辨率高等诸多的优势而更适合农田尺度的作物长势监测^[6]。无人机飞行平台可以通过搭载丰富的传感器来获取地面作物信息，目前使用较为广泛的无人机载传感器有普通数码相机、可见光—近红外多光谱相机、高光谱相机等^[7]。其中，光谱相机可以获取地物在不同波段下的光谱反射率，继而通过数学运算可以得到各种植被指数以突出植被特征，因此可在LAI反演中作为重要的模型输入参数。国内外诸多学者通过构建植被指数与实测LAI的经验统计模型来实现LAI的反演^[8-10]，这种模型结构简单、适用性强，尤其是在不清楚地面实况或遥感信号产生机理过于复杂的情况下，是一种较合适的反演方法^[11]。然而，由于多光谱遥感手段只能获取作物冠层信息、无法得到作物垂直结构信息而使其估测能力有限^[12]。而且，当监测区域的作物覆盖度高时，光谱信号会饱和^[13]，植被指数表现出对冠层LAI变化不敏感的现象，这也使得估测精度受到一定程度的影响。

激光雷达（Light Laser Detection and Ranging，LiDAR）作为一种主动式遥感技术，可借助无人机遥感平台实现短时间内获取作物高精度三维信息^[12]，恰好弥补多光谱遥感对作物垂直信息提取不足的缺陷，为基于无人机遥感的LAI反演提供了新的技术手段。基于机载LiDAR的LAI监测技术率先在林业研究中得到推广应用，且表现出较好的监测效果^[14-16]。此外，已有研究证明采用LiDAR数据反演森林地上部生物量时并未出现饱和现象或者有极大值，因此认为通过该法可以解决光谱遥感存在植被指数饱和的问题^[16]。然而，由于农作物高度较低，冠层结构不明显等原因，无人机载LiDAR对农作物的监测还处于探索阶段^[17]：苏伟等^[18]的研究中使用机载LiDAR提取玉米冠层点云结构参数，将其与实测LAI构建回归模型，结果表明预测LAI与实测LAI之间的决定系数（Coefficient of Determination，R²）达0.79，均方根误差（Root Mean Square Error，RMSE）为0.47，反演效果较好，证实了LiDAR点云数据反演玉米等中等高度作物LAI的可行性。更进一步的研究，若能够综合多光谱和LiDAR遥感技术在作物冠层光谱信息及作物空间结构信息获取方面的优势，将无人机多光谱和LiDAR数据融合应用来反演LAI等作物生长参数或许将为作物长势监测提供更优的手段。前期在玉米等较大株型农作物中的研究表明，与单独使用多光谱数据相比，将LiDAR与多光谱数据融合，可有效提高其地上部生物量的估测精度（R²分别增加0.30，RMSE减少54.92 g·m^–2，NRMSE下降7.65%）^[19]。因此，本研究旨在探索无人机载LiDAR与多光谱数据融合是否能提高株型较小的冬小麦LAI反演精度，以为冬小麦LAI反演提供更优模型。

试验设置于河南省新乡市封丘县（34°53′~35°14′ N，114°14′~114°46′ E）中国科学院封丘农业生态试验站内。试验区地势平坦，土壤类型为潮土，年平均气温在13.5℃~14.5℃之间，年降水量为615.1 mm，属于半干旱、半湿润的暖温带季风气候。选用试验站内的冬小麦氮肥梯度处理长期试验地作为研究对象，该试验地为自2005年开始的冬小麦-夏玉米轮作长期试验。冬小麦的生育期为每年的10月至次年6月，根据作物生长需要，设置了F1、F2、F3、F4和F5共5个施氮处理，单季作物施氮量分别为0、150、190、230和270 kg·hm^–2，应用尿素作为单一氮肥。其中，60%作为基肥，40%作为追肥，磷肥和钾肥等肥料不进行差异处理。F1、F2、F4、F5处理重复9次，F3处理重复12次，试验区内共48个小区（图 1）。

图 1 试验区布局 Fig. 1 Arrangement of the experiment zone

试验采用大疆M210四旋翼无人机作为飞行平台，该飞行平台可实现全自动起飞降落和规划航线飞行，续航20~30 min，最大载重2 kg，适合作业面积相对较小的，搭载微单相机、轻便多光谱仪等仪器的情况。飞行平台上搭载有Prime RedEdge-MX五通道型多光谱相机（MicaSense，美国）和扫描鹰HS-300 LiDAR扫描仪（RIEGL，奥地利）。其中，多光谱相机的相机分辨率约为130万像素（1280×960），焦距为5.4 mm，光谱参数如表 1所示；LiDAR扫描仪，其激光扫描单元最大视场角360°可调节，测距精度为15 mm，最大测距范围（60%反射率目标）为250 m。

表 1 Table 1

表 1 Prime RedEdge-MX五通道型多光谱相机的光谱参数 Table 1 Spectral parameters of the prime rededge-MX five-channel multispectral camera 谱带序号 Serial number of spectrum 谱带名称 Name of spectrum 中心波长 Central wavelength/nm 半峰全宽 Full width at half peak/nm 1 蓝Blue 475 20 2 绿Green 560 20 3 红Red 668 10 4 近红外Near infrared 840 40 5 红边Red edge 717 10

表 1 Prime RedEdge-MX五通道型多光谱相机的光谱参数 Table 1 Spectral parameters of the prime rededge-MX five-channel multispectral camera

于2019年4月28日，即冬小麦孕穗期进行试验区低空及地面数据的采集。搭载五通道多光谱相机和LiDAR扫描仪的无人机飞行平台，根据地面控制站预先设计好的航线在当天晴朗无云的情况下于10:00~14:00进行全程自动飞行。飞行高度设定为50 m，多光谱图像的地面分辨率约为3.49 cm，每次飞行前后，在地面上对校准的反射率白板进行拍摄来获取多光谱相机的辐射校准图像。无人机飞行结束后可以获取到涵盖整个试验区的若干张单波段光谱图像以及原始的LiDAR点云数据，同期使用SunScan冠层分析仪（Delta-T，英国）进行试验区冬小麦叶面积指数的地面测定工作，每个小区取五个点，于无人机飞行的同期在10:00~14:00测定，记录各点读数，求其平均值作为该小区LAI地面实测值。由此获取到试验区内冬小麦关键生育时期（孕穗期）的无人机多光谱图像、LiDAR点云数据和地面实测LAI这三部分低空和地面数据。考虑到后续的建模过程，将各不同处理内试验数据均按照7︰3的比例随机划分建模集和验证集，共32个小区试验数据为建模集，16个小区试验数据为验证集。

由无人机遥感平台获得的多光谱图像和LiDAR点云数据需要进行前期预处理。为了得到覆盖整个试验区的多光谱正射图像，需要将若干张单波段光谱图像进行拼接，使用Pix4Dmapper 4.5.6软件通过已知反射率值的辐射校准图像来对多光谱图像进行辐射校准以及全程自动化拼图。同时，通过机载LiDAR获取的点云数据需要预先转化成Las数据集。

针对LiDAR点云数据，首先用CloudCompare 2.10.2软件对其进行去噪处理，去掉点云数据中明显的高点和低点等异常点，然后利用基于不规则三角网滤波（TIN）算法对点云进行多次滤波处理^[20]，TIN的原理是通过迭代建立三角网的方式来过滤地面点：假设在一定范围内必有一个初始点，先由这些初始地面点构成三角网，然后根据限制条件（地形角度、迭代角和迭代距离）不断向三角网中加入点进行迭代，逐步细化三角网，而最终得到的三角网顶点为地面点^[21]。由于该滤波算法只是粗略地将地面点和植被点分开，因此需要使用人工交互的方法将错分的点云进行正确归类。由于LiDAR点云是离散的数据，为了得到试验区小麦的连续表面形态及垂直结构，借助ArcGIS 10.1将分离出的地面点和植被点进行栅格化处理，根据王安^[22]的研究结果，选用反距离插值算法进行插值，空间分辨率设为0.02 m，地面点建立数字高程模型（Digital Elevation Model，DEM），它可以反映试验区地势的起伏变化，除噪点之外的所有点建立数字表面模型（Digital Surface Model，DSM），DSM与DEM进行栅格差减运算得到冠层高度模型（Canopy Height Model，CHM），在去除了地表起伏的影响后，CHM在一定程度上可以反映小麦表面形态和冠层垂直结构。

将得到的冠层高度模型在R3.5.2中进行参数提取，提取的参数包括：三维体积指标（BIOVP）、平均株高（Hmean）、75%的分位数高度（H75）、激光穿透力指数（Laser Penetration Index，LPI）这四种点云结构参数。其中BIOVP是根据Han等^[23]提出的估测生物量的一个体积指标，它是指小麦植株高度在相应面积范围内积分和，Hmean是指小麦的平均株高，H75是指小麦植株的75%分位数高度，LPI是指地面激光回波点数与总回波点数之比。

预处理时已经拼接好的多光谱图像，在R3.5.2中将其按照试验小区划分兴趣区并求每个小区的平均植被指数，参照各研究LAI反演时常用的6种植被指数，包括：归一化植被指数（Normalized Difference Vegetation Index，NDVI）^[24]、土壤调节植被指数（Soil Adjusted Vegetation Index，SAVI）^[25]、改进的叶绿素吸收指数（Modified Chlorophyll Absorption Ratio Index，MCARI）^[26]、三角植被指数（Triangular Vegetation Index，TVI）^[27]、归一化红边差值植被指数（Normalized Difference Red Edge Index，NDRE）^[28]、比值植被指数（Resistant Vegetation Index，RVI）^[29]。各植被指数的计算公式见表 2。

表 2 Table 2

表 2 植被指数计算公式 Table 2 Calculation formula of vegetation indices 植被指数 Vegetation index 计算公式 Calculation formula 归一化植被指数NDVI ${\text{NDVI}} = ({\rho _{{\text{nir}}}} - {\rho _{{\text{red}}}})/({\rho _{{\text{nir}}}} + {\rho _{{\text{red}}}})$ 土壤调节植被指数SAVI ${\text{SAVI}} = (1 + L) \times ({\rho _{{\text{nir}}}} - {\rho _{{\text{red}}}})/({\rho _{{\text{nir}}}} + {\rho _{{\text{red}}}} + L)(L = 0.5)$ 改进的叶绿素吸收指数MCARI $ {\text{MCARI}} = ({\rho _{{\text{rededge}}}} - {\rho _{{\text{red}}}}) - 0.2 \times ({\rho _{{\text{rededge}}}} - {\rho _{{\text{green}}}}) \times ({\rho _{{\text{rededge}}}} - {\rho _{{\text{red}}}}) $ 三角植被指数TVI $ {\text{TVI}} = 0.5 \times [120 \times ({\rho _{{\text{nir}}}} - {\rho _{{\text{green}}}}) - 200 \times ({\rho _{{\text{red}}}} - {\rho _{{\text{green}}}})] $ 归一化红边差值植被指数NDRE ${\text{NDRE}} = ({\rho _{{\text{nir}}}} - {\rho _{{\text{rededge}}}})/({\rho _{{\text{nir}}}} + {\rho _{{\text{rededge}}}})$ 比值植被指数RVI $RVI = {\rho _{{\text{nir}}}}/{\rho _{{\text{red}}}}$ 注：ρred ，ρgreen ，ρrededge ，ρnir分别代表小麦冠层在可见光的红光、绿光、红边、近红外波段的反射率值。Note：ρred，ρgreen ，ρrededge and ρnir stands respectively for reflectance in red，green，red edge，and near-infrared bands of wheat canopy.

表 2 植被指数计算公式 Table 2 Calculation formula of vegetation indices

在用遥感手段反演作物生长参数中常用的两类方法为物理模型法和经验统计模型法，相比原理复杂的物理模型，经验统计模型由于结构简单、模型适用性强而应用较为广泛^[11]，本研究使用经验统计模型对冬小麦LAI进行反演。因需要输入多个变量构建经验统计模型，因此建模方法采用多元线性回归（Multivariable Linear Regression，MLR）和偏最小二乘回归（Partial Least Squares Regression，PLSR）这两种多变量统计回归方法，使用MLR方法的原因是它对建模结果容易解释，可以较为直观地反映各自变量与因变量之间的关系^[23]。而与传统的MLR方法相比，PLSR模型在建模过程中结合了主成分分析法和多元线性回归法的特征，可以解决MLR方法中自变量之间存在的多重共线性问题^[30]。建模过程全部在R 3.5.2中进行。

在该研究中，使用决定系数R²评估预测LAI与实测LAI之间的拟合优度，R²的范围在0~1之间，数值越高表明拟合效果越好。同时，使用均方根误差RMSE衡量观测值与真实值之间的偏差。RMSE越小，模型估算越准确。

用反距离插值算法将LiDAR点云数据进行栅格化处理后，得到空间分辨率为0.02 m的数字高程模型（DEM）（图 2）、数字表面模型（DSM）（图 3）以及冠层高度模型（CHM）（图 4）。由DEM图可以看出，试验区地势整体较为平坦，而且在去除了地表起伏影响之后，CHM图可以看出小麦冠层高度均不高于63 cm，F1处理与其他处理的小麦冠层高度差异明显。

图 2 试验区数字高程模型（DEM） Fig. 2 Digital elevation model of the experiment zone

图 3 试验区数字表面模型（DSM） Fig. 3 Digital surface model of the experiment zone

图 4 试验区冠层高度模型（CHM） Fig. 4 Canopy height model of the experiment zone

为筛选合适的LiDAR点云结构参数及植被指数作为LAI估测模型的输入参数，将点云结构参数及植被指数与实测LAI做了相关性分析，相关分析矩阵图如下（图 5），由各变量之间的相关分析矩阵图可以看出：实测LAI与各建模参数的相关系数绝对值均达到了0.90以上，且均为极显著相关；各植被指数之间呈极显著的强自相关性；各点云结构参数之间相关性同样较强，且极为显著。

注：图中对角线区域代表变量自身的分布，对角线的左下部分为横纵轴两个变量的散点图和拟合曲线，对角线的右上部分代表横纵轴两个变量的相关系数，“***”代表变量在P < 0.01的水平下显著相关。  Note：The areas alongside the diagonal represent distribution of the variable per se；in the lower left part of the diagonal are the scatter plots and the curves fitting between the two variables；the upper-right part of the diagonal shows the relationship between the two variables；“***” indicates that variables are significantly correlated（P < 0.01）. 图 5 叶面积指数与各建模参数的相关分析矩阵图 Fig. 5 Matrix diagram of the correlation between leaf area index and various modeling parameters

基于相关性分析结果，由于各植被指数之间和点云参数之间存在不同程度的自相关性，因此采用MLR方法建模时，为了防止过度拟合、降低模型的复杂程度，选用与实测LAI相关性最强的NDVI和SAVI这两种植被指数以及同样与实测LAI相关性最强的H75和BIOVP这两种点云结构参数作为模型输入参数。PLSR法建模时，将涉及到的所有参数一同作为输入参数，然后根据交叉验证结果选择模型的主成分个数^[31]，根据分析结果在融合LiDAR点云结构参数与植被指数建模时选用3个主成分，在仅使用植被指数建模时选用1个主成分。

首先，在只有植被指数NDVI、SAVI以及融合了LiDAR点云结构参数H75、BIOVP和植被指数NDVI、SAVI作为模型输入参数这两种情形下，分别由建模集各建模参数与实测LAI进行多元线性回归建立LAI反演模型，同时为检验模型是否可靠，利用验证集数据对该模型进行验证。MLR回归方程如下（表 3）：

表 3 Table 3

表 3 利用不同建模参数反演LAI的多元线性回归方程 Table 3 Multiple linear regression equation for LAI inversion using different modeling parameters 建模参数Modeling parameters 回归方程Regression equation LAI~NDVI+SAVI LAI=11.807×NDVI–2.864×SAVI–1.271 LAI~NDVI+SAVI+BIOVP+H75 LAI=–2.920×NDVI+4.842×SAVI+0.183×BIOVP+3.352×H75

表 3 利用不同建模参数反演LAI的多元线性回归方程 Table 3 Multiple linear regression equation for LAI inversion using different modeling parameters

类似地，接下来用偏最小二乘法建模，首先将通过交叉验证结果中筛选的所有植被指数中1个主成分与实测LAI进行偏最小二乘回归，然后从所有点云结构参数与植被指数中筛选出3个主成分作为模型输入参数，进行偏最小二乘回归，PLSR回归方程如表 4所示。

表 4 Table 4

表 4 利用不同建模参数反演LAI的偏最小二乘回归方程 Table 4 Partial least squares regression equation for LAI inversion using different modeling parameters 建模参数Modeling parameters 回归方程Regression equation LAI~NDVI+SAVI+MCARI+TVI+NDRE+RVI LAI=0.165×NDVI+0.165×SAVI+0.164×MCARI+0.164×TVI+0.163× NDRE+0.165×RVI LAI~BIOVP+Hmean+H75+LPI+NDVI+SAVI+ MCARI+TVI+NDRE+RVI LAI=0.313×BIOVP–0.009×Hmean+0.170×H75–0.073×LPI+0.147× NDVI+0.067×SAVI+0.171×MCARI–0.069×TVI–0.045×NDRE+ 0.178×RVI

表 4 利用不同建模参数反演LAI的偏最小二乘回归方程 Table 4 Partial least squares regression equation for LAI inversion using different modeling parameters

最后，为比较分析模型效果，将建模集与验证集的估测LAI与实测LAI做散点图及拟合曲线。图 6可以看到，MLR建模方法下，在只有植被指数参与建模时（图 6 a、图 6 b），建模集和验证集估测LAI与实测LAI之间的拟合曲线R²分别为0.937和0.897，RMSE分别为0.382和0.492。而融合了LiDAR点云结构参数与植被指数后（图 6 c、图 6 d），R²分别为0.950和0.901，RMSE分别为0.341和0.480，相比单独使用植被指数建模时R²均有所升高，且估测LAI与实测LAI之间的偏差降低，总体而言提高了LAI的估算精度。

注：红色线条为1︰1参考线。以植被指数为输入参数的建模集a和验证集b，以LiDAR点云结构参数与植被指数为输入参数的建模集c和验证集d。  Note：The red line represents the 1︰1 reference line. Modeling Set a and Validation Set b with Vegetation Indices as Input Parameters，Modeling Set c and Validation Set d with LiDAR Point Cloud Structure Parameters and Vegetation Indices as Input Parameters. 图 6 基于MLR建模方法的LAI预测值与实测值散点图 Fig. 6 Scatter plot of estimated LAI and field-observed LAI based on the MLR modeling method

PLSR建模方法下，在只有植被指数参与建模时（图 7 a、图 7 b），建模集和验证集估测LAI与实测LAI之间的拟合曲线R²分别为0.957和0.892，RMSE分别为0.319和0.486。而融合了LiDAR点云结构参数与植被指数后（图 7 c、图 7 d），R²分别为0.968和0.909，RMSE分别为0.278和0.445，较单独使用植被指数建模时其拟合优度同样有所提高，并且估测LAI与实测LAI之间的偏差也有所降低，提高了LAI的估算精度。

注：红色线条为1︰1参考线。以植被指数为输入参数的建模集a和验证集b，以LiDAR点云结构参数与植被指数为输入参数的建模集c和验证集d。  Note：The red line represents the 1︰1 reference line. Modeling Set a and Validation Set b with Vegetation Indices as Input Parameters，Modeling Set c and Validation Set d with LiDAR Point Cloud Structure Parameters and Vegetation Indices as Input Parameters. 图 7 基于PLSR建模方法的LAI预测值与实测值散点图 Fig. 7 Scatter plot of estimated LAI and field-observed LAI based on the PLSR modeling method

总之，无论是多元线性回归法还是偏最小二乘回归法，在融合了机载LiDAR数据与植被指数后均较仅使用植被指数建模效果好，这同时也证明了LiDAR数据对于冠层垂直结构提取的优势能够弥补多光谱遥感在此方面的缺陷，但同时也可以看到模型评价指标之间的差别最多不到0.1，建模效果提升不算明显，通过分析认为这可能与建模和验证数据集过于离散有关，可以看出除F1处理的LAI实测值在1附近，F2~F5处理LAI大部分均介于4~6之间，因此数据集的扩充问题还有待解决。

为提高田块尺度冬小麦LAI的反演精度，弥补无人机多光谱遥感在LAI反演时对作物垂直信息提取不足且易产生植被指数饱和的问题以及挖掘LiDAR遥感在作物垂直信息提取方面的潜力，本研究将无人机载LiDAR提取的作物冠层点云结构参数与机载多光谱提取的植被指数融合，借助多元线性回归法和偏最小二乘回归法进行建模，实现了冬小麦关键生育时期——孕穗期LAI的反演。整体而言，LAI反演模型效果较好，说明借助经验统计模型的方法使用无人机载LiDAR和多光谱数据建模反演冬小麦LAI是可行的。无论是多元线性回归还是偏最小二乘回归的建模方法，融合了LiDAR点云结构参数与植被指数的模型要优于只有植被指数的模型，虽然模型对比差异不是很明显，但一定程度上验证了LiDAR数据可以弥补光谱遥感数据在作物垂直方向上信息提取不足的缺陷，二者的融合应用可以提高冬小麦叶面积指数的反演精度。然而，由于本研究中所用样本数据集中的LAI数据整体没有均匀分布在一定范围内，而是过于离散，导致在一定程度上影响了对模型效果提升程度的估计，所以未来仍需更丰富的研究来支持和证实此结论。