森林冠层结构形态是森林的重要几何参数，与森林资源调查的重要测树因子——树高、冠幅直接相关，准确地描述和表达森林冠层结构形态，有利于树高、冠幅、生物量等参数的提取和估算(庞勇等，2008； Chen，2007； Falkowski et al.，2006； Leeuwen et al.，2010)。机载激光雷达(light detection and ranging，LiDAR)系统发射的脉冲信号与森林冠层的枝、叶及地面发生作用，通过接收来自森林冠层枝、叶及地面的回波信息，可获取森林冠层及地面的空间位置信息。这些回波信号形成离散点云数据，可反映森林冠层结构形态。冠层高度是指森林冠层枝、叶离地面的垂直距离，是反映森林冠层结构形态的主要参数。对于离散点云数据，通常通过对点云进行归一化得到冠层高度，即冠层点高程减去地面点高程。归一化的冠层高度点云是树高、冠幅、生物量等森林参数提取和估算的基础数据。在树高、冠幅、生物量等森林参数的提取和估算时，通常需要对归一化的冠层高度点云进行空间插值，生成规则格网数据，即森林冠层高度模型(canopy height model，CHM)，然后再利用CHM进行森林参数提取和估算(Anderson et al.，2005; Gatziolis et al.，2010)。

CHM是LiDAR林业应用中单木参数提取的主要数据源，通常采用图像处理手段对CHM进行图像分割，得到分割后的CHM，然后根据分割结果提取单木冠幅、位置及树高等参数。CHM的质量影响冠幅分割，利用未准确表达森林冠层结构形态的CHM进行冠幅分割，通常会出现过分割或欠分割现象，从而直接影响单木冠幅、位置和树高等提取的参数准确性(Chen et al.，2006； Koch et al.，2006)。CHM是否准确表达森林冠层结构形态不仅与原始离散点云数据的质量有关，而且与空间插值方法有关。国内外研究者对有关离散点云生成数字地面模型(digital terrain model，DTM)空间插值方法进行了大量研究，如B样条插值(B-Spline)(朱长青等，2006； Kraus et al.，2001； Pirotti et al.，2013)、普通克里金插值法(ordinary Kriging，OK)(Rasib et al.，2013； Kim et al.，2013； Guo et al.，2010)、线性插值三角网法(triangulation with linear interpolation，TLI)(Bater et al.，2009)和反距离权重插值法(inverse distance weighted，IDW)(Mikita et al.，2013)等。DTM采用的点云数据为地面点，而CHM采用的点云数据为森林冠层点云。目前对森林冠层离散点云生成CHM插值方法研究较少。一些研究者直接采用生成DTM的插值方法生成CHM(Vazirabad et al.，2015)，事实上这并不一定合适。空间插值方法有很多，对于具有不同特点的数据而言，没有绝对最优的空间内插方法，应针对数据空间分布特点选择合适的插值方法。基于激光雷达离散点云数据具有点云分布、密度均匀的特点。同时，在森林冠层边缘局部区域包括属于森林冠层的点云和属于地面的点云，属于冠层的点云高度值较大，而属于地面的点云高度值较小；并且单株木之间存在林冠空隙，即分属不同单株木的点云之间存在林冠空隙，即森林冠层点云在冠层边缘处存在高度值突变。

由离散点云数据插值生成的CHM应能准确地反映和表达森林冠层表面形态及冠幅之间的林冠空隙，不同插值方法对由离散点云数据生成CHM是否准确表达森林冠层形态有重要影响。为了研究不同插值方法的插值效果，本文基于森林区域离散点云的特点，利用B-Spline，OK，TLI和IDW 4种插值方法构建冠层高度模型，从插值影像像元统计量、平面视图、三维视图、剖面图等方面进行比较、分析和评价，为森林冠层高度模型插值方法选择提供参考。

1 试验数据

试验数据为位于湖南省平江县芦头林场的机载LiDAR数据，采用的机载LiDAR设备为IGI公司(德国，www.igi.eu)制造的LiteMapper 6800。分别选取针叶林、针阔混交林、阔叶林3个森林类型面积为30 m×30 m的样地作为试验数据，试验样地基本信息如表 1所示。野外实地测量了样地内树高及平面位置等参数。

2 试验方法 2.1 数据预处理

利用TerraScan软件(芬兰TerraSolid公司，http://terrasolid.fi/)对原始点云数据进行噪声去除，如飞鸟点、低点，然后进行点云分类。由于森林区域目标地类为地面和植被，需将建筑物、电力线等点云去除，保留地面点和植被点。将分类后的地面点云利用不规则三角网法(triangulated irregular network，TIN)插值生成DTM。对植被点及地面点进行归一化处理，即利用点云与DTM对应平面位置的高程之差森林高度代替点云高程，得到归一化植被高度点云。然后对归一化植被高度点云按回波次数提取所有的归一化首次回波植被高度点云子集，首次回波包括只有1次回波的点云和多次回波中第1次回波点云。前者往往属于地面和森林冠层，而后者往往属于森林冠层。

2.2 插值方法

空间插值法是研究区域变量空间分布的基本方法。本文选择B-Spline，TLI，OK和IDW 4种插值方法进行比较分析。

B-Spline插值是采用样条函数插值，对一些限定的点通过控制估计方差，对采样点进行多项式拟合的方法来进行插值(Childs，2004)。由于 B 样条插值函数具有很好的连续性，因此被广泛应用于地学模型的构建。B-Spline插值采用正规化(Regularized)样条函数，主要参数有权重和搜索点数。权重是影响表面插值特征的参数，采用Regularized样条函数时，用来指定最小化期间附加到三阶导数项的权重，权重越大，输出表面越平滑。搜索点数是每个插值像元时所使用的点数，搜索点数越多，输出表面越平滑。本研究B-Spline插值权重设为 0.1，搜索点数为12个。

TLI是将采样点连接成Delaunay三角形，形成一个覆盖采样区域的三角网，然后利用三角网进行线性插值的方法(Popescu et al.，2002)。

克里金(Kriging)插值法是地学统计中最常用的插值法之一，是以空间自相关为基础，考虑了区域化变量的相关性和变异性，使用变异函数度量样点间的空间相关性，对插值点的区域化变量进行线性无偏最优估计的一种空间插值方法。本质上，Kriging插值法是对插值点周围的采样点进行加权以得出插值点的估计值，权重不仅取决于采样点之间的距离、插值点的位置，还取决于采样点的整体空间排列。用于拟合采样点空间关系的拟合模型包括线性模型、高斯模型、球状模型、指数模型、圆形模型等，地学插值通常选择球状模型(Kravchenko et al.，1999)。Kriging插值有很多种，其中OK插值是一种假设数据变化呈正态分布、变量期望值未知时的Kriging插值法。OK插值参数为：采用球面半变异模型，最大搜索距离0.8 m，最小搜索点数4，最大搜索点数20。

IDW(段祝庚等，2014； 王勇等，2008)是在一定距离范围内，以插值点与样本点间的距离的幂为权重进行加权平均，离插值点越近的样本点赋予的权重越大，见式(1)：

(1)

(2)

式中：(X，Y，Z)为插值点坐标；(x_i，y_i，z_i)搜索样本点坐标； p为权重； q为幂次； i为样本点编号； n为以插值点为中心、R为半径的搜索区域内样本点总数，搜索半径。

IDW插值法分为全局搜索和局部搜索，全局搜索为所有样本点参与插值； 局部搜索为仅插值点一定范围内的样本点参与插值，参与插值的点数可由参数搜索半径或插值点数决定。本研究IDW插值采用局部搜索法，沿任意方向搜索，最大搜索半径为0.8 m，距离幂次为2。

CHM主要用于描述冠层高度及其顶层结构形态，通常采用首次回波点归一化植被高度点云生成CHM(段祝庚等，2014)。针对试验数据的首次回波归一化植被高度点云，利用开源地理信息系统软件SAGA(system for automated geoscientificanalyses，www.saga-gis.org)分别采用B-Spline，TLI，OK和IDW 4种插值方法进行插值生成CHM，并三维显示。试验区机载LIDAR的点云密度平均为8～10 m^-2，即原始点云平均间距为0.3～0.35 m，CHM插值像元大小与原始点云平均间距相当，设为0.3 m。为了评价以上插值方法针对离散点云数据生成CHM的效果，通过CHM栅格影像平面视图、三维视图、剖面图及其像元统计量进行比较、分析及评价。

3 结果与讨论 3.1 插值结果及评价

利用去除噪声的原始点云与地面点插值生成的DEM求差可得到归一化点云。样地1原始点云和归一化点云三维视图如图 1a，b所示。然后利用归一化点云分别采用B-Spline，TLI，OK和IDW 4种插值方法生成CHM。样地1不同插值方法生成CHM的三维视图分别如图 1c-f所示。

不同森林类型样地1，2，3各种不同插值法CHM插值结果如图 2，3，4所示，其像元值统计量如表 2所示。评价插值算法是否适合离散点云生成CHM的原则是生成的CHM与森林冠层的自然形态的一致性，二者的一致性可以通过CHM栅格图及其剖面图目视评价及像元值统计量反映。

CHM与森林冠层自然形态的一致性要求CHM冠层形态可以适当填充，对原始点云适当进行平滑，而冠层边缘点云不被过度平滑，保留高度突变，同时林冠空隙仍然保留也不被过分填充。从三维视图图 1和平面视图图 2，3，4可以看出，B-Spline和TLI插值法生成的CHM冠层形态保留过多细节，未适度平滑，而林冠空隙进行了一定程度填充。特别是B-Spline插值的CHM显得比较破碎，未构成完整冠层顶部。同时将所有无首次回波点云数据的空值(no data)区域都进行了填充，林冠空隙(gap)也被过分填充。OK插值法生成的CHM影像模糊，出现过度平滑现象。而IDW插值法对冠层顶部进行了适当填充和平滑，但冠层边缘不被过度平滑，保留高度突变，同时林冠空隙仍然保留也不被过分填充，相对其他3种插值方法而言，较好地保留了森林冠层的真实自然形态。空值像元数量可以一定程度反映林冠空隙是否被过分填充。由表 1也可以看出，B-Spline和TLI插值法的空值像元数均为0，空值区域被填充，与目视判读结果一致。

CHM剖面图可以详细反映沿某剖面线插值结果细节。图 5，6，7分别为不同森林类型的样地1，2，3由B-Spline，TLI，OK和IDW插值法生成的CHM剖面图[样地1 CHM第70行(Line=70)、样地2 CHM第75行(Line=75)、样地3 CHM第40行(Line=40)]，从剖面图上可以更清楚地得出上述相同的结论。

另一方面，栅格数据的像元统计量像元最大值(高度值)与归一化点云高度最大值、实测最大树高的一致性可评价插值方法的效果。由表 2可以看出，样地1，2，3由B-Spline插值生成的CHM像元最大值分别为20.022，20.670，39.871 m，比归一化点云高度最大值15.04，16.02，30.58 m大4.982，4.650，9.291 m，经过插值后像元最大值明显增加。B-Spline插值的CHM像元最大值明显偏离归一化点云高度最大值和实测最大树高； 而IDW插值的CHM像元高度最大值与归一化点云高度最大值和实测最大树高相差较小，基本一致。

另外，无论哪一种插值方法，插值结果均受参与插值的原始点云数量及密度的影响。本研究采用首次回波进行插值生成CHM。由表 1可以看出，样地1(针叶林)、2(针阔混交林)和3(阔叶林)的首次回波数与点云总数之比分别为79.5%，77.5%和86.1%，首次回波的点云密度分别为8.6，8.0和7.3 m^-2，不同森林类型的点云密度差别不大。森林结构具有复杂性，不同森林类型具有不同的森林结构。而有关森林类型对插值方法的影响本文未予研究与讨论，这将是下一步研究的内容。

综上所述，对于森林区域空间分布均匀且存在高度突变的点云数据，IDW插值法优于B-Spline，OK，TLI插值法，生成的CHM能较准确反映森林冠层的真实自然形态。

3.2 IDW插值搜索半径选择

由式(1)可知搜索半径R和幂次q是IDW插值的关键参数，二者共同决定参与插值的点的权重，其中幂次q通常设为2。R决定参与插值点的数量，即搜索点数。确定搜索半径R有2种方式：一是直接设定搜索半径R的大小； 二是通过设定参与插值点的总数，间接确定搜索半径。当点云密度均匀时宜采用直接设定搜索半径R的方式； 当点云密度不均匀时，通过设定参与插值点的总数可以自适应地确定搜索半径。一般来说，机载LiDAR点云密度较为均匀，可通过设定搜索半径R的大小的方式进行IDW插值。由于搜索半径R的大小决定了参与插值的点的数量，下面讨论搜索半径R大小的选择。

以样地1为例，点云总数为9 698点，首次回波点数为7 704点，平均点云密度为10.77 m^-2，平均点间距为0.3 mm。图 8为样地1搜索半径R设为0.3～1.5 m IDW插值CHM。表 3为对应的IDW插值不同搜索半径冠层高度统计。不同插值半径的CHM冠层高度最小值、最大值、算数平均值、标准差均大致相同，无太大变化；但随着搜索半径的减小，空值像元(no data cell)逐渐增多。搜索半径为1.5，1.2和1.0 m时均无空值像元，而搜索半径为0.8，0.5和0.3 m时，空值像元数分别为221，674和1 686个。

图 9为搜索半径R为0.3～1.5 m时IDW插值生成的CHM中Line=70处剖面及其对应点云，参与插值的点云宽度(width)为搜索半径R的2倍，即点云宽度为0.6～3.0 m。由图 9可以发现，随着IDW插值搜索半径的增加，参与插值的点云迅速增加。尽管剖面宽度内的点云都参与了插值，但IDW是按距离平方加权插值，离插值点较远的样本点起作用较小，插值点高度值主要由较近样本点的高度值决定。特别地，当样本点与插值点重合时，样本点的权重为1，即插值点的高度与该样本点的高度完全一致，故IDW插值是无误差插值法。

搜索半径为0.3～1.5 m时，CHM剖面存在稍许差别，主要表现在搜索半径较小如0.3，0.5，0.8 m时 CHM剖面存在空值，但当搜索半径增加到1.0 m及以上时，CHM剖面便不再出现空值情况。空值像元的形成是由于在搜索半径范围内无首次回波点存在。当搜索半径增加到一定程度时，空值像元完全被填充，CHM剖面不再变化，CHM也就趋于稳定； 但稳定的CHM并不是最佳的CHM，原因是林冠空隙被过分填充，且冠幅边缘均被向外扩张了，向外扩张的最大值为搜索半径R。

搜索半径影响参与插值IDW插值的样本点数主要表现在2方面：一方面是空值像元的产生； 另一方面是林冠空隙被过分填充及冠幅边缘的向外扩张。过小的搜索半径会导致CHM存在过多的空值像元，产生冠层空洞，不能准确表达冠层的结构形态，对森林参数提取不利； 搜索半径过大，冠层会过度平滑，冠层信息丢失较多，且冠幅边缘向外扩张较大，从而导致较小林冠空隙完全填充，改变冠幅实际大小，不能真实表达冠层结构形态，也不利于冠幅、树高等森林参数提取。搜索半径与点云密度有较大关系，通常点云密度大，需要选择较小的搜索半径； 点云密度小，需要选择较大的搜索半径。对比搜索半径为0.3～1.5 m生成的CHM(图 8)可以发现，搜索半径为0.3 m时产生了过多空值像元，搜索半径为1.2～1.5 m时冠幅明显向外扩张了。合理选择IDW插值时搜索半径对CHM质量有重要影响。搜索半径为原始点云平均间距的1.5～2.5倍，即搜索半径设为0.5～0.8 m较为合适。

4 结论

由于森林冠层离散点云存在空间分布均匀且高度突变的特点，由离散点云数据插值生成CHM时，插值算法对CHM空间结构形态有较大影响，影响CHM准确表达森林冠层自然形态。通过对B-Spline，TLI，OK和IDW等插值算法的分析和比较，得出以下结论：

1) 对于森林区域空间分布均匀且存在高度突变的点云数据，B-Spline插值对空值(no data)区域都进行了填充，林冠空隙(gap)也被过分填充，且CHM像元最大值明显偏离原始插值数据；TLI插值的CHM显得比较破碎；OK插值法对影像过度平滑，生成的CHM影像模糊；而IDW插值法对冠层顶部进行了适当填充和平滑，但冠层边缘不被过度平滑，保留高度突变，同时林冠空隙仍然保留也不被过分填充。IDW插值法优于B-Spline，OK，TLI插值法，生成的CHM能较准确反映森林冠层的真实自然形态。

2) 搜索半径影响参与IDW插值的点云点数量，搜索半径过小，可能导致CHM存在大量空值像元，产生冠层空洞； 搜索半径过大，冠层会过度平滑，冠层信息丢失较多，且冠幅边缘向外扩张较大，不能真实表达冠层结构形态，也不利于冠幅、树高等森林参数提取。IDW插值应选择合适的搜索半径，搜索半径为原始点云间隔的1.5～2.5倍较为合适。