1 引 言

与传统激光雷达(light detection and ranging,LiDAR)系统记录点云数据不同，全波形LiDAR系统能够以波形的形式记录一定高程范围内不同高程点上的后向散射能量^{[1, 2]}。记录的波形信息可以提供一定范围(光斑)内的地表高程信息，以及不同高程点上反射能量。根据不同高程点反射能量的大小，能够检测出激光光束传播过程中目标垂直方向(高程方向)的分布情况^{[3, 4]}。根据获取数据时激光束照射在地面光斑的大小，全波形LiDAR系统可分为大光斑系统和小光斑系统。大光斑LiDAR系统的光斑直径一般为8~70 m，脉冲发射频率为40~500 Hz，光斑之间间距在170 m左右，光斑分布稀疏；小光斑LiDAR系统的光斑直径一般为0.3~1 m，脉冲发射频率为40~500 kHz，光斑间的间距小于1 m，光斑分布密集^[5]。目前，常用的配有全波形数字化仪的LiDAR系统包括NASA机载的LVIS(laser vegetation imaging sensor)系统^{[6, 7]}、SLICER(scanning LiDAR imager of canopies by echo recover)系统^{[8, 9]}和星载的GLAS(geoscience laser altimeter system)系统^{[10, 11]}等大光斑LiDAR系统；Optech ALTM3100、Leica ALS-70、TopEye Mark Ⅱ、Riegl LMS-Q560以及Falcon Ⅲ等小光斑LiDAR系统^{[12, 13, 14]}。

全波形LiDAR数据记录的是地物在不同高程点上后向散射能量，因此全波形LiDAR数据不仅蕴含地物的高程信息，而且能够反映地物的散射特性。波形分解是获取这些信息的重要手段。波形分解是指用有限一致分布统计模型(Gaussian、Nakagami、Burr等)的组合来拟合原始波形数据。

目前，已经有许多学者对全波形LiDAR数据的分解进行了研究。文献[15]通过雷达方程证明高斯函数可以拟合全波形LiDAR数据，为分解波形数据提供了理论依据。文献[16]提出了基于高斯模型的波形分解方法。该方法虽然能够将原始波形数据分解为若干个高斯分量，但容易陷于局部最优解，无法得到高精度的拟合结果。文献[17]针对波形数据的高斯特性，对原始波形数据进行了分解，该方法虽然提高了拟合精度，但是容易遗漏部分重叠的返回波，分解结果不能充分反映地物沿高程方向的信息。

通过高斯分解，分析高斯分量数、均值和标准差等参数，能够得到地物高程信息和后向散射特性。上述文献中高斯分量数均为给定值，然而由于受地形起伏复杂、回波信号噪声以及地物多次漫反射等诸多因素影响，很难人为确定波形分解数，并且容易遗漏部分重叠的返回波，降低波形拟合精度。因此，为了自动确定波形分解数，获得精确的参数估计值，本文提出了一种基于可变分量高斯混合模型和可逆跳转马尔科夫链蒙特卡洛(reversible jump Markov chain Monte Carlo,RJMCMC)的全波形LiDAR数据分解方法。

给定的全波形LiDAR数据集x={x_i;i=1,2,…,n}，y={y_i;i=1,2,…,n}，其中，x_i为分布在有效高程L内的采样点，y_i为在采样点x_i处的反射能量，n为采样点个数。由此，全波形LiDAR数据可进一步表示为：z={(x_i,y_i);i=1,2,…,n}。一束激光的波形数据为由采样点x作为横坐标，能量值y作为纵坐标形成的点集，如图 1所示。

图 1 波形数据 Fig. 1 Waveform data

图选项

在波形分解过程中，首先假设理想波形为高斯混合模型，用理想波形来拟合实际波形数据；定义能量函数以刻画二者的拟合程度；通过能量函数定义吉布斯分布，表征全波形LiDAR数据的或然率；根据贝叶斯定理构建刻画波形分解的后验概率模型；设计RJMCMC算法模拟该后验概率模型，从而确定波形分解个数及各高斯分量参数，完成波形分解。图 2所示为该算法总体流程，其中，实线框为完成的功能，虚线框为完成功能需要的模型和函数。

图 2 波形分解流程 Fig. 2 Waveform

图选项

设理想波形为高斯混合模型，表达为m个高斯函数分量加权和

点击浏览公式

为了控制理想模型与波形数据拟合程度定义能量函数。设参数集Φ={μ,σ,w,m}，其中，μ={μ₁,μ₂,…,μ_m}；σ={σ₁,σ₂,…,σ_m}；w={w₁,w₂,…,w_m}。定义控制y_i和 之间差异程度的能量函数为U(z|Φ)，由试验分析能量函数用绝对误差和残差二次型的平方根的效果相同，但采用绝对误差可以减少计算量，因此，用绝对误差代替残差二次型的平方根，即

点击浏览公式

由能量函数U(z|Φ)定义表征波形数据或然率的吉布斯分布h(z|Φ)为

点击浏览公式

为了分解波形，需要构建已知波形数据z条件下高斯分量的位置μ、标准差σ、权重w、个数m的联合后验概率密度函数。根据Bayes定理，已知波形z条件下的后验概率为

点击浏览公式

式中的先验概率讨论如下：

(1) 高斯分量位置μ。高斯分量的均值μ_k为对应第k个分解波形最大反射能量的高程，称为该高斯分量的位置。将m个高斯分量的位置按照增序排列为0<μ₁<μ₂<…<μ_m<L。假设μ_k(k=1,2,…,m)服从区间[0,L]上的均匀分布且相互独立，则其联合概率密度函数为

(2) 高斯分量标准差σ。假设各高斯分量的标准差σ_k(k=1,2,…,m)服从均值和标准差分别为μ_σ和σ_σ的正态分布且相互独立，则其联合概率密度函数为

(3) 高斯分量权重w。假设各高斯分量的权重w_k(k=1,2,…,m)服从[0,1]上均匀分布并相互独立，则其概率密度函数为

(4) 高斯分量个数m。假设高斯分量的个数m服从均值为λ的泊松分布，则其概率密度函数为 ，其中，2≤m≤k_max，k_max为预先指定的最大高斯分量数。

综上，式(4)的后验概率可重写为

点击浏览公式

波形分解是通过对构建的后验概率模型进行模拟求解，本文采用RJMCMC算法进行后验概率模拟。

RJMCMC算法实现参数空间维数改变时，对参数集Φ的相关采样满足概率分布(式(5))。提取候选参数矢量Φ^*，根据文献[18]，其接受率表示为

点击浏览公式

当参数空间维数不变时，称为M-H(Metropolis-Hastings)算法，式(6)接受率简化为

点击浏览公式

采用RJMCMC算法对式(5)后验概率模型进行模拟，为了实现完备采样，设计4种移动操作类型：①更新高斯分量参数；②更新高斯分量权重w；③合并或分裂高斯分量；④生成或删除高斯分量。

更新高斯分量参数包括更新高斯分量位置μ和更新高斯分量标准差σ，该操作中参数空间维数没有发生改变，其过程包括：在m个高斯分量中随机抽取一个高斯分量k∈{1,2,…,m}，该操作改变该高斯分量位置μ_k或标准差σ_k，而保持其他参数不变。候选μ_k^*和σ_k^*分别服从以当前μ_k和σ_k为均值，σ_t和σ_s为标准差的正态分布。其中σ_t和σ_s分别为μ_k和σ_k的迭代步长，是预先定义的常数。当σ_t过大时可以很快遍历整个μ_k取值范围，但是得到μ_k的精度较低；当σ_t过小时能够取到较高精度的μ_k，但是需要耗费较长时间。同理，σ_s具有相同的作用。通过试验，当σ_t和σ_s分别取50和1时得到的结果比较理想。则改变μ和σ的接受概率分别为：a(μ_k,μ_k^*)=min{1,R_μ}，a(σ_k,σ_k^*)=min{1,R_σ}，其中

点击浏览公式

点击浏览公式

在更新高斯分量权重w的操作中，参数空间维数没有发生改变，其过程包括：在m个高斯分量中随机抽取两个高斯分量，假设选取第k和k+1个高斯分量，该操作改变w_k和w_k+1，保持其他参数不变。随机生成s∈(0,1)，则改变高斯分量权重前后需满足平衡条件w_k^*=(w_k+w_k+1)s和w_k+1^*=(w_k+w_k+1)(1-s)。则改变w的接受概率为：a(w_k,w_k^*)=min{1,R_w}，其中

点击浏览公式

在合并或分裂高斯分量操作中，高斯分量的个数m发生变化，同时参数空间维数也发生改变。在分裂处理中，首先以等概率选择高斯分量k，将其分裂成两个高斯分量，分别记为k^*和(k+1)^*，它们的位置、标准差以及权重分别为(μ_k^*,μ_k+1^*)、(σ_k^*,σ_k+1^*)、(w_k^*,w_k+1^*)。将μ_k+1、μ_k+2、…、μ_m重新标号为μ_k+2^*、μ_k+3^*、…、μ_m+1^*，σ、w作同样的操作，则参数μ={μ₁,μ₂,…,μ_m}、σ={σ₁,σ₂,…,σ_m}和w={w₁,w₂,…,w_m}分别变为μ^*={μ₁,μ₂,…,μ_k^*,μ_k+1^*,…,μ_m+1^*}、σ^*={σ₁,σ₂,…,σ_k^*,σ_k+1^*,…,σ_m+1^*}和w^*={w₁,w₂,…,w_k^*,w_k+1^*,…,w_m+1^*}。分裂后，高斯分量的位置μ_k^*和μ_k+1^*均服从以当前μ_k为均值，σ_t为标准差的正态分布，其中σ_t为预先定义的常数；标准差σ_k^*和σ_k+1^*均服从以当前σ_k为均值，σ_s为标准差的正态分布，其中σ_s为预先定义的常数；权重w_k^*和w_k+1^*满足平衡条件w_k^*=w_ku₁和w_k+1^*=w_k(1-u₁)，u₁为保证维度平衡而设置的随机数，u₁∈(0,1)。则分裂高斯分量的接受概率为a_s=min{1,R_s}，其中

点击浏览公式

合并高斯分量是分裂高斯分量的对偶操作。其过程包括：在m+1状态下选择合并高斯分量操作，其概率为m_m；然后在原有的m+1个高斯分量中以等概率随机抽取一个高斯分量，其概率为1/m+1，然后在与其相邻的两个高斯分量中随机抽取一个与其合并，其概率为1/2。由于合并高斯分量操作是分裂高斯分量操作的对偶操作，因此其接受概率为：a_c=min{1,1/R_s}。

在生成或删除高斯分量操作中，高斯分量数m发生变化，同时参数空间维数也发生改变。对于生成高斯分量操作，首先以等概率选择一个高斯分量k，设新生成的高斯分量的位置、标准差和权重分别为μ_k+1^*、σ_k+1^*和w_k+1^*。μ_k+1^*服从(μ_k,μ_k+1)之间的均匀分布，σ_k+1^*由其服从的先验分布中抽取，权重w_k^*和w_k+1^*满足平衡条件w_k^*=w_ku₂和w_k+1^*=w_k(1-u₂)，u₂为保证维度平衡而设置的随机数，u₂∈(0,1)，则生成一个高斯分量的接受概率为a_b=min{1,R_b}，其中

点击浏览公式

删除高斯分量是生成高斯分量的对偶操作,其过程包括：在m+1状态下选择删除操作，其概率为d_m；在原有m+1个高斯分量中以等概率随机抽取一个高斯分量k，并删除该高斯分量，然后将位置μ_k+1、μ_k+2、…、μ_m+1重新标号为μ_k^*、μ_k+1^*、…、μ_m^*，其概率为1/m+1。由于删除高斯分量操作是生成高斯分量操作的对偶操作，因此其接受概率为a_d=min{1,1/R_b}。

试验数据采用ICESat-GLAS测高数据，来源于美国冰雪中心(national snow and ice data center，NSIDC)(http：//nsidc.org/data/icesat/)获取，数据采集时间为2003年10月05日至2005年10月23日。GLAS系统为首个极地轨道星载大光斑LiDAR系统，每秒发射40次激光脉冲，每个激光脉冲在地面上形成直径约为70 m的光斑，光斑之间的间隔为170 m。每个光斑对应一个波形，每个波形有544个采样点。每个数据采集时间段对应一个操作周期。L1A和L2A操作周期的波形数据采样点间隔为1 ns(15 cm)，有效高程范围为81.6 m。从L2B操作周期以后，对波形进行了压缩处理，前151个采样点间隔增加到4 ns(60 cm)，有效高程范围增加到150 m^{[19, 20]}。ICESat-GLAS共提供了15种标准数据产品(GLA01-GLA15)，用户可以根据不同需求和目的选择相应的标准数据产品。试验数据采用33版本的GLA01和GLA14两种标准数据产品。其中GLA01中记录了光斑号、时间、脉冲强度、经纬度等信息，GLA14中记录高程、大地经纬度、回波振幅等地表高度信息^[21]。

本文的研究区针对城市区域，图 3所示为研究区的遥感影像和ICESat-GLAS光斑分布，图 3(a)中由上至下4条线分别为2004年5月20日、2005年5月22日、2005年10月23日和2003年10月05日的卫星轨迹。由遥感图像可以看出，该城区主要地物为道路、建筑物、广场、公园、矿区等。本文选择了小区、广场、公园、郊区、矿区等几个具有代表性的区域进行研究，图 3(b)-(f)分别为选择5个典型区域的遥感影像，圆圈内为光斑所覆盖的范围。图 4为5个研究区域对应的ICESat-GLAS原始波形数据，表 1为5个研究区域的数据信息。

图 3 研究区遥感影像及ICESat-GLAS光斑的分布 Fig. 3 Remote sensing image of study area and distribution for ICESat-GLAS shots

图选项

图 4 GLAS原始波形 Fig. 4 GLAS raw waveform

图选项

为便于对不同返回波能量进行比较和加权，增强不同波形之间的可比性，需对波形数据进行归一化处理，使能量值落入到[0,1]区间。

为了提高波形拟合的精度，需要对波形数据进行去噪处理。定义噪声阈值为单个回波的4倍标准差，即：ε=4σ_r，其中，σ_r为前100个高程点能量值的标准差。当脉冲能量小于阈值时认为是噪声，将其能量视为0。

首先，采用本文提出的算法对5个区域的波形数据进行波形分解试验。图 5为最终分解结果。

图 5 分解结果 Fig. 5 Decomposition results

图选项

图 5(a)为区域1波形分解结果，分解成5个高斯分量。根据该区域遥感影像(图 3(b))可以看出该区域主要地物有高低不齐的平顶建筑物，车棚以及树木。为了验证地物覆盖，对实地进行了勘测，并测得其平均高程分别为20、14.5、3和1 m。因此该地区分解的5个高斯分量分别代表高建筑物返回回波、低建筑物返回回波、车棚返回回波、树木返回回波以及地面返回回波。试验得到的地物高程分别为20.3、14.4、2.7和1.2 m。图 5(b)为区域2波形分解结果，分解成两个高斯分量，根据该区域遥感影像(图 3(c))可以看出该区域主要地物为低矮植被，实测高程为0.3 m。因此该区域分解的两个高斯分量分别为低矮植被返回波和地面返回回波。试验得到的植被高程为0.3 m。图 5(c)为区域3波形分解结果，分解成7个高斯分量，根据该区域遥感影像(图 3(d))可以看出该区域垂直方向上主要地物有亭子和较高的植被，实测高程分别为8和4.5 m。因此该地区分解的7个高斯分量分别代表亭子返回回波，中间的几个高斯分量为植被返回回波，最后一个高斯分量为地面返回回波。试验得到亭子和植被的高程分别为8和4.4 m。图 5(d)为区域4波形分解结果，分解成两个高斯分量，根据该区域遥感影像(图 3(e))可以看出该区域主要地物为尖顶的低建筑物，经实地勘测，其高程为3 m。因此该区域分解的两个高斯分量分别为建筑物顶部返回回波和地面返回回波，试验得到建筑物的高程为2.7 m。图 5(e)为区域5波形分解结果，分解成7个高斯分量，根据该区域遥感影像(图 3(f))可以看出该区域垂直方向上主要地物是地面倾斜的矿区，实测高差为8.5 m。因此该地区分解的7个高斯分量分别代表矿区不同高程地面返回回波。试验得到地面倾斜矿区的高差为8.6 m。

通过分解后波形可以看出，低矮植被返回信号能量比较集中，振幅较高，脉冲较窄(图 5(b))；较高植被返回信号能量较分散，振幅较低，脉冲较宽(图 5(c))；平顶建筑物返回信号能量较大，振幅较高，脉冲较窄(图 5(a))；尖顶建筑返回信号能量较小，振幅较低，脉冲较宽(图 5(d)；地面倾斜的矿区根据地面倾斜程度的不同返回信号能量、振幅、脉冲的宽度都有所不同(图 5(e)。

从图 5的分解结果可以看出，对于不同地物类型区域的波形数据，由本文算法得到的理想波形都能很好地拟合原始波形数据，并且分解的波形能够与实际地物相对应，能够得到地物沿高程方向的散射特性，说明本文算法具有可行性和有效性。

试验进行了10 000次迭代，在迭代循环过程中，误差很快减小。图 6为5个研究区域10 000次迭代误差变化曲线。由图 6可以发现，前4000次迭代误差有一定波动，在后6000次迭代中误差达到稳态，为了使试验结果更加准确，忽略前4000次的迭代结果，比较后6000次高斯分量个数，并计算其分别出现的概率(表 2为区域1不同高斯分量个数出现的概率)，将出现概率最大的高斯分量个数作为最终分解的波形个数。

图 6 10 000次迭代误差变化曲线 Fig. 6 Change in error during 10 000 iterations

图选项

该算法基于Matlab 7.12.0 开发，运行环境为Intel Core i5-3470-3.20 GHz-4 G，Window7操作系统。以区域1为例，该区域有544个采样点，一次迭代完成5种移动操作的时间约为0.02 s。因此，完成10 000次迭代总共的操作时间约为4 min。更新高斯分量位置、更新高斯分量标准差、更新高斯分量权重、合并或分裂高斯分量以及生成或删除高斯分量分别占用时间的比例约为：10%、6%、4%、45%以及35%。其中，合并或分裂高斯分量和生成或删除高斯分量两种操作中不仅需要判断是否满足进行该操作的条件，而且还要判断操作后高斯分量的均值、标准差以及权重是否满足条件，因此，需要的操作时间略长。

为了验证本文算法的可靠性和准确性，将本文方法与Matlab提供的CFTOOL高斯拟合方法进行比较。由于CFTOOL高斯拟合方法不能自动确定分解成高斯分量数，所以将本文算法得到的高斯分量数作为该方法高斯分量数，然后将原始波形数据分解为固定数的高斯分量之和。图 7为CFTOOL高斯拟合方法的波形分解结果。表 3分别列出了两种方法得到的5个试验区域高斯分布的均值估计值和标准差估计值。对比两种分解方法结果，本文方法不仅能自动获得高斯分量数，还可以较好地将波形分解，分解的波形能够与不同高程地物相对应，同时得到的高斯曲线能够较好地拟合原始波形。CFTOOL高斯拟合方法不能自动获得高斯分量数，虽然分解波形能较好地拟合原始波形，但未能分解出与不同高程地物对应的波形分量，如区域1的第3个高斯分量标准差估计和区域2的第2个高斯分量标准差估计值过大，振幅接近于平值，并且对部分波形未能拟合，特别是有细微能量变化的波形，如区域3的第1个脉冲和区域4的第1个脉冲，区域5的第2个脉冲。因此，CFTOOL高斯拟合方法分解的波形不能与地物相对应，无法反映地物沿高程方向的散射特性，其原因在于CFTOOL高斯拟合方法没有对分解波形的参数进行约束，因此出现标准差估计值较大等情况。

图 7 CFTOOL方法分解结果 Fig. 7 Decomposition results from CFTOOL

图选项

为了对波形分解结果进行定量评价，分别计算每个试验区域原始波形与拟合波形之间的相关系数ρ和KS距离^[22]，计算公式分别为

点击浏览公式

点击浏览公式

表 4为计算结果，根据计算结果可以看出，得到的相关系数均在0.98以上，KS距离均在0.2以下，说明本文算法具有有效性和准确性。

本文提出了一种基于可变分量混合高斯模型和RJMCMC的全波形LiDAR数据分解方法，首先假定波形数据服从高斯分布，然后对LiDAR波形数据实现可变个数的分解。本文方法不仅能自动确定波形分解数，同时还可以有效地获取最优的模型参数估计。利用本文算法分别对不同地物类型的区域进行了波形分解试验，并和Matlab中CFTOOL的高斯拟合方法进行了对比，定性和定量的测试结果验证了本文方法的有效性、可靠性和准确性。虽然本文算法能够准确地确定波形分解的个数，实现波形的分解，并且分解的结果能够与地物相对应，同时具有很高的拟合精度，但是对于一些非对称、倾斜或者拖尾的波形，本文的方法还不能非常准确地拟合。在未来的工作中，对如何进一步提高波形拟合精度，使算法更普遍地适用于各种形状的波形，还需要继续从理论上进行研究，如选择其他的随机模型来模拟原始波形数据，使分解的波形更准确地拟合原始波形数据，从而提高波形拟合精度。