当前以机载LiDAR点云数据求取林分和单木参数的过程中，DEM(数字高程模型)滤波求取是最基础的工作之一，此时需尽力得到纯地面激光点数据而去除地上植被和建筑等非地面数据的干扰。因后续森林参数的提取多是基于由DEM得到的nDSM(规格化表面模型)进行，所以非地面数据一旦进入DEM则会直接影响最终林木参数提取精度。

随着机载LiDAR硬件水平的提高, 以机载LiDAR数据为基础，生成理想的DSM(数字表面模型)已不是难题，而生成DEM绝大部分算法都是基于激光数据点的高程突变信息进行的，几乎每种方法都有其自身的缺陷，有待于进一步改进(张小红，2007)。尽管机载LiDAR量测技术在国外已相当成熟，但是量测数据处理的算法仍处于前期研究发展阶段，其中关键的问题之一就是机载LiDAR数据的滤波和分类(Axelsson，1999)。目前, 国内外比较典型的基于点云本身的DEM滤波算法主要包括总体算法和局部算法2种。从点云总体入手的算法主要有迭代线性最小二乘内插法(Kraus et al., 1998)等；从局部入手的算法主要有数学形态学方法(Lindenberger，1993；Kilian et al., 1996)、移动窗口法(Petzold et al., 2000)、基于地形坡度变化的滤波法(Vosselman，2000；Mass，1999)、移动曲面拟合滤波法(张小红等，2004)等。局部算法应用较广，这些算法从总体上看，起点大都是基于先认定局部(如某窗口尺度内)点云中高程最低数据点为地面点，并由此先确定一些地面数据，然后根据某些原则(例如高程、角度等阈值)，认定距离这些地面数据较近的其他数据点作为地面点，从而达到加密地面点的作用。采取不同加密原则产生了不同的算法, 2类算法通常都需要人工参与调整参数或阈值, 所以，如何减少人工干预、获取较高精度的DEM模型，一直是DEM滤波研究的热点之一。

目前已有不少LiDAR量测系统可记录多次回波信号，本文正是基于此类数据，提出了在点云数据中，依据各类回波数据的强度和高程属性分析数据性质，并自动进行数据粗略划分，然后使用判别分析来完成DEM提取工作，避免了人工确定阈值参数。本文的方法称为“多次回波判别分析法”，计算简便且运算速度快，取得了较为理想的结果。

1 研究区域及试验数据

试验地位于甘肃省张掖市大野口观测站(100°15′ E，38°32′ N)，属祁连山高山深谷地貌，树种分布为青海云杉(Picea crassifolia)天然纯林。本文所用数据涉及机载LiDAR数据和地面数据2种，地面数据又包含建模样地和验证样地数据。

机载LiDAR数据：试验区域共重叠飞行6次(6个航带)，获得近似的重叠航带数据。在航带内设置建模样地和验证样地。对验证样地应用样地内涉及的全部6次飞行数据(平均点间隔约为0.54 m，每平方米约3.47个点)，获取验证样地激光数据点145 404个(多次飞行数据叠加); 建模样地则只取其第5次飞行数据，获取激光数据点20 946个(单次飞行)。这些已经足够用于分析激光数据分布特点和建立判别模式的要求。

建模样地数据：用于分析激光数据分布特点及建立判别模式。在航带中选取一块既包含成片林木、又包含林中空地试验样地作为建模的训练样地，样地中间部分为空地，两端为大片林地，大小约为100 m×130 m。其机载LiDAR数据结构包括：每个数据点的三维坐标、回波强度、回波次序号和回波总数6个属性。通过此样地的LiDAR数据，得出单次和多次回波数据的空地和林内的分布特点。

验证样地数据：用于验证DEM判别精度及后续林地林木参数研究。实测了航带中1块郁闭林地，面积约为100 m×100 m，林木密集分布。使用DGPS和全站仪获取样地地面定位点41个，二者观测结果坐标差值的平均值为0.100 m(横坐标)、0.204 m(纵坐标)、-0.361 m(高程坐标); 并选定样地中单株林木基部地面点作为其位置点，进一步结合全站仪测定了验证样地中全部1 439株单木地面点的三维坐标，总计获取地面点三维位置实测数据1 480个。由实测数据计算得验证样地为平均海拔约2 800 m、坡度约9°、坡向约325°的阴坡，该验证样地内因有出露煤层，存在开采后遗留的塌陷坑，故地形表面也存在局部波动。

两样地数据点云的地面形状如图 1，图 1A为验证样地激光数据水平分布情况，其中网格顶点及网格内中心点为样地地面定位点，内部林木参数及单木位置三维坐标均为外业实测; 图 1B为建模样地激光数据水平分布，红色实心圆表示样地中心，连续的红色三角形画出贯穿样地的一条完整扫描线，并示范地给出样地右下角一个数据点的空间位置(WGS-84坐标)。

2 研究方法概述 2.1 激光点云数据分类

本文划分整体激光点云数据为单次回波点集和多次回波点集，进一步在多次回波点集中再分割出首次回波点集、中间回波点集和末次回波点集3种数据类型，分类方式见图 2。在此分类的基础上，分别提取出研究样地所涉及的各类型数据子点集加以分析利用。

2.2 强度信息分类

影响每个激光数据点强度的主要因素有探测距离、反射方位、物质密度变动等。相同物质的回波信号强度较为接近，故可用来区分不同属性的物质。对同一种物质的表面，对LiDAR系统来讲，物体表现的光谱特性接近激光波长时，会具有较强的反射性。一般认为，在小区域、飞行条件近似时，可认为回波强度只与物质表面有关(张小红，2007)。本文在分别分析不同数据点集的反射强度特征基础上，再探讨相互间联系。

2.3 地面或植被激光点云的分类过程

把全部激光点云分成2类：植被π_P和地面π_G。分3步进行计算：1)通过激光点性质分析，认定部分地面点T_G和部分植被点T_P; 2)用已经认定的地面点和植被点建立判别模式P(H，I)，H是光点高程，I是光点强度; 3)对未被认定的其余激光点，使用判别模式P(H，I)来判别每个点的类别归属，得到其判别分类的结果。

3 试验数据分析和计算过程 3.1 单次回波和多次回波的空间分布

单次回波：图 3A为建模样地内全部单次回波显示，可看出建模样地中部空地上激光点排列密集且规则，而两端林地内激光点分布松散，这是因为地面只能生成单次回波。多次回波只产生在林地内(图 3C)，其使用能加大激光点描述林地情况的密度，在水平分布上可见林内存在大量多次回波激光脚点，而林中空地则由单次回波点表现。

对样地总体及选取其中林地内、空地上2个小区域(图 3A中a，b)的单次回波数据，统计得到建模样地的单次回波激光点平均面密度D_sur和平均采样密度D_ave，公式为: D_sur=N/S，。其中N为区域内激光点个数; S为区域多边形面积。此参数反映了林地和空地对单次回波响应的大致情况，即空地内的单次回波平均采样点距离更小，而林地内单次回波平均采样点距离大(表 1)。

图 3B为按扫描线绘制的建模样地前200个数据点示意图，求得各扫描线中相邻数据点间水平距离最大极值是9.55 m，位置处于林下，更可看出单次回波在林地内沿扫描方向形成远大于采样平均距离的空洞，对整个建模样地此数值则达到了15.63 m，在图 3A中标出了这2个点的位置。

多次回波：建模样地内包含最多回波次数为5的回波数据。在一个发射脉冲返回的多个激光点中，首次回波点高程最高，不涉及地面，更能体现林木上表层的信息; 末次回波点高程最低，更接近地表; 中间数据点在高程上位于前两者之间，也不涉及地表，可反映林内垂直结构。图 3C为建模样地内全部多次回波显示，只发生在林木范围，通过计算每个发射脉冲所得首次回波和末次回波之间的高程差，可体现这些回波点所代表林木的空间垂直分布情况，尤其高程差极大值在总体上大致反映了林木的高度情况，因为末回波点不一定就落在地表上; 而其高程差极小值一方面体现了传感器的分辨精度，另一方面也体现了植被枝叶高度分布的复杂性。建模样地的多次回波高程差极小值为0.16 m，极大值为18.84 m。

3.2 样地激光点反射强度和空间高程关系

图 4为建模样地(上排，第5次飞行)和整条航带(下排，第5次飞行)激光数据的强度直方图。二者的强度和频数范围由表 2给出，可见二者有相似的强度曲线分布。

单次回波的强度曲线呈明显双峰形态，表征为地面和植被两大范畴，其中与多次回波强度重叠区间(具多次回波各数据类型)必然是植被范围，另一区间则可能是表征土壤、岩石或其混合等非植被的地表范围。

多次回波中的末次回波点相对在该多次回波光束中高程最低，更接近地面或者就是地面点，其频数曲线呈现出轻微的双峰形态，在强度属性上既有接近植被，也有接近地表; 多次回波中非末回波激光为多次回波中首次回波和中间回波的集合，空间分布上高程较高，只可能是植被点，这些点分布在植被冠层的中上部位。

建模样地和整个航带的强度跨度范围从大到小均依次为单次回波、多次回波的末次回波、多次回波的首次回波、多次回波的中间回波。这启示我们可以从强度曲线差异中找出一些信息加以利用。

3.3 激光点云数据距地面远近的相对空间分布划分

结合前面的分析结果，在单次回波分析中得知空地由单次回波表征，林地内易形成多次回波(也含单波); 而通过多次回波的垂直分布得知首次回波和中间回波点位置上不属地面点，距离地面较高; 再通过数据的强度分布和空间高程的关系可得出：图 2中，多次回波所包含的3个点集属于林地，各点集在强度跨度上各不相同，且多回波的强度分布范围小于单回波。因此可依据反射强度并结合高程垂直特点划分激光点云总体数据为4种类型，本文中称之为远离地表、近地表、疑似地表和绝对地表。划分原则如下。

远离地表：为首回波和中间回波强度最小值至中间回波强度最大值区间。此部分数据一定属于植被范围(首波和中间波)，空间位置距地面较高。单次回波和末次回波数据强度处于此强度区间的，也应表征林木。

近地表：为中间回波强度最大值至首次回波强度最大值区间。此部分虽为首次回波，但强度上已经不处于中间回波范围，更趋向远离植被而接近地表，故空间位置较低。这部分数据是加密地表数据的关键，此强度区间内的单次回波和末次回波数据几乎全部都在空间上接近地面。

疑似地表：为多次回波首回波强度最大值至多次回波末回波强度最大值区间。此部分数据强度上已经脱离了表征植被的首回波和中间回波强度范围，但仍包含多次回波的末波和单次回波数据。从末次回波角度看，数据在水平位置分布上接近植被，垂直位置分布上处于植被基部或地面。

绝对地表：为多次回波强度最大(末波)至单次回波强度最大区间。在强度上看，这部分数据完全脱离了多次回波强度范围，只有单次回波存在，从前面的单波分析来看，其表征的应为地面。

图 5给出了建模样地(A)和验证样地(B)强度直方图的多项式拟合曲线，目的是使数据区间段内的曲线描述更加平滑，在验证样地中进一步按上述划分依据给出了验证样地总体激光数据距离地表情况的4种划分区间示意。

3.4 激光点所属类型的判别分析

通过上述分析，可以把全部激光点分成3类(图 6): 1)多次回波首回波和中间回波数据集的空间位置是植被，称为植被训练样本; 2)强度大于多次回波首回波最大强度的末回波和单次回波数据集的空间位置是地面，称为地面训练样本; 3)其他的单次回波和多次回波末回波数据集的空间位置是待定样本。将利用植被训练样本和地面训练样本建立判别公式P(H, I)。用判别分析的方法，判定每个数据点的归属。本文采用Fisher和Bayes 2种判别方法进行探讨。

3.4.1 Fisher判别式的建立

Fisher判别原理如下：设2类总体π_i的均值和协方差矩阵为(μ⁽ⁱ⁾, Σ_i)。线性函数y=l′x=(μ⁽¹⁾-μ⁽²⁾)′Σ^-1x称为Fisher线性判别函数，取m为2个总体均值在投影方向上的中点，即μ⁽²⁾)′Σ^-1(μ⁽¹⁾+μ⁽²⁾)，则其判别规则为; 当μ⁽¹⁾, μ⁽²⁾, Σ未知时，可由总体中分别抽出样品并计算相应的样本均值和协方差阵来作为μ⁽¹⁾, μ⁽²⁾, Σ的估计。

3.4.2 Bayes判别式的建立

Bayes判别法的基本思想是计算待判样品x属于各总体的条件概率P(l|x)，l为地面或植被; 比较这2个概率的大小，从而将样品归入来自概率最大的总体。

设C(g|p)和C(p|g)分别为2种错判所带来的损失，即将地面点g判为植被点p或者相反，考虑2个总体的先验概率，则由某种判别准则因误判所带来的平均损失为ECM(expected cost of misclassificatiom)(期望误判损失，也称误判风险)：

一个优秀的判别准则应该使ECM最小，使上式达到最小的判别区域划分存在下面的定理(张润楚，2006)：

若认为2种错判的损失相等，在二元正态的情形下π_i~N₂(μ⁽ⁱ⁾, Σ_i)，i=1, 2，本文假定协方差阵Σ>0相等，可直接推出判别函数为y(i/x)=lnq_i-, i=1, 2。通过比较每个数据点2类判别函数值的大小，确定样点的归属。

3.4.3 验证样地判别结果的得出

在依据强度和高程划分了点云数据距离地表分布的基础上，本文采用疑似地表和绝对地表 2部分激光点数据集作为描述地表的地面训练样本T_G; 采用多次回波的首回波和中间回波激光点数据集作为描述植被的植被训练样本T_P进入判别分析。判别对象(待判样本)Q为验证样地激光点云的单次回波和多次回波的末回波激光点数据集(图 6)。选入判别的变量(指标)为激光点数据的强度和高程2个属性值。

Fisher判别：依据强度和高程2个指标直接就可进行，进入判别式的变量有μ⁽¹⁾, μ⁽²⁾, Σ，m和待判样本x。其中μ⁽¹⁾, μ⁽²⁾为T_G和T_P中回波高程的平均数和回波强度的平均数，Σ是T_G和T_P中回波高程和强度的协方差矩阵，m可用μ和Σ求出，x为待判样本Q中的每个样品。本文Fisher判别计算使用ForStat 2.1版进行(唐守正等，2008)，最终判得地面数据23 848个。

Bayes判别：需确定总体的概率分布表达式，这样的条件根据数据分析还不够得到其全部真实分布，依据初始总体的分位数图是检查数据是否呈现正态分布的有效手段(黄华兵，2008)，本文假定其接近正态，并试用强度和高程二维正态总体分布函数来进行Bayes判别。进入判别式的变量有q_g，q_p，μ⁽¹⁾, μ⁽²⁾, Σ和x。其中q_g，q_p为T_G和T_P中每类数据各自发生的概率，μ⁽¹⁾, μ⁽²⁾, Σ和x含义同Fisher判别，这部分判别为作者编程实现，最终判得地面数据13 647个。

3.5 DEM滤波的判别分析和直接栅格/窗口DEM滤波模拟方法结果的比较

用判别分析所得地表类别数据点进行0.25 m规则网格的TIN内插，生成验证样地DEM模型；并用此DEM模型计算出验证样地内实测地面点的模拟高程值(Fisher内插Z，Bayes内插Z)。

为了比较判别分析方法的优良性，本文同时对验证样地地表进行10，15和20 m 3种尺寸窗口格网划分，提取格网中最低激光点作为地面点，并进行0.25 m规则网格的TIN内插直接进行DEM模拟，得到DEM模型并算出实测地面点的模拟高程值(网格10，15和20 m内插Z)。因为，若取更小于10 m的窗口划分，则进入树木冠幅或相邻冠幅相交叉的地上范围，故不再更小划分。

然后用实测地面点所在栅格的高程减去上述的模拟高程差(高差)，进行统计比较，比较结果如表 3所示。

验证样地实测地面点和Fisher判别结果进行DEM模拟的结果如图 7所示，其中投影线条为等高线，二者均为0.25 m规则网格线性内插的结果。

4 结论与讨论

1) 对山区林地林木复杂环境样地的DEM提取采用判别分析的方法是可行的。本文用验证样地实测的地面点作为DEM精度检验的依据，这些地面点覆盖整个验证样地。由表 3可以看出，判别分析结果与实测地面点同位置的高差比较，其平均高差精度达到0.24 m，方差0.05 m，误差处于-0.5~0.5 m之间的数据(表 3频数统计，也从误差的大小和点数分布2方面说明了2个DEM的相合程度)占总量的93.85%，分布更为合理。

2) 多次回波点云数据的分类提取并进行判别分析避免了人工确定参数或阈值，有助于提高数据处理的自动化。基于多次回波和强度属性划分，并以其高程作为辅助的DEM判别分析方法，区别于传统绝大多数以高程突变的DEM生成依据的方式。将激光脚点依据回波或其他属性分为若干类型并分别提取，针对不同研究目的组合相关的数据集，不仅减少计算数据量，而且增加了数据中针对具体信息的含量。

3) 此次飞行数据为青海云杉纯林，实测郁闭度约0.8。从试验结果上看，在不计激光点云地面点数据概率分布和错判损失的情况下，直接使用Fisher判别即可取得较好的结果。从建模样地、验证样地和航带数据的强度波形上看(图 5)，各类数据在波形特征和跨度上都呈现相似的特征，所以初步认为，在实际应用中，只需确定待研究林木分布范围(或进行范围分割)，然后如能在点云文件中划分出2类训练样本，且激光数据点密度能够满足建模精度，则本方法都是适用的。