基于地物散射的全波LiDAR波形分解方法

摘要

本发明提供一种基于地物散射的全波LiDAR波形分解方法，包括：激光雷达发射信号到地物目标上反馈回来的信号叠加形成全波LiDAR波形；采用有效散射单元来描述全波LiDAR波形：将全波LiDAR波形转化为各采样点对应的具有相同散射能量的有效散射单元；采用阶梯函数对全波LiDAR波形中的有效散射单元个数建模；利用RJMCMC算法对全波LiDAR波形进行分解，得到地物目标相对高程以及散射特性。本发明针对采用有效散射单元来描述全波LiDAR波形，将全波LiDAR波形转化为各采样点对应的具有相同散射能量的有效散射单元，能够将地物散射特性与波形分解相结合，并提出采用阶梯函数对有效散射单元个数建模，打破了使用传统高斯统计分布模型来拟合原始波形的思想。

说明书

技术领域

本发明所属全波LiDAR数据处理技术领域，具体是一种基于地物散射的全波LiDAR波 形分解方法。

背景技术

波形分解是全波LiDAR数据研究的主要方面，是全波LiDAR数据应用的前期工作，因 此波形分解成为全波LiDAR数据处理领域研究的热点及难点问题。波形分解是指用统计分布 模型拟合原始波形数据的过程。

LiDAR波形分解包括两个关键问题：(1)实现波形数据的最佳分解；(2)分解结果能够反 映被测目标的散射特性。目前波形分解方法中多数采用高斯分布来拟合原始波形数据，尽管 高斯分解模型具有模型简单，拟合容易等优点，但在实际激光反射中，当激光光束照到较为 复杂的地物目标时，波形会发生倾斜、拖尾等形变，在这种情况下，采用高斯函数拟合原始 波形数据将产生较大的拟合误差。广义高斯模型、Nakagami和Burr等统计模型也常用来拟 合原始波形数据。

但这些波形分解方法仅仅考虑波形拟合，并没有从地物散射机理上考虑波形的形成，不 能充分的反映被测目标的散射特性。并且，目前的波形分解方法仅仅针对波形数据本身，并 没有考虑研究区域内地物目标及其散射特征等特性，因此难于挖掘蕴含于波形数据中的地物 信息。

发明内容

针对现有技术存在的不足，本发明提供一种基于地物散射的全波LiDAR波形分解方法。

本发明的技术方案是：

基于地物散射的全波LiDAR波形分解方法，包括：

步骤1、激光雷达发射信号到地物目标上反馈回来的信号叠加形成全波LiDAR波形，包 括分布在有效高程区域内的采样点及各采样点处的散射能量；

步骤2、采用有效散射单元来描述全波LiDAR波形：将全波LiDAR波形转化为各采样 点对应的具有相同散射能量的有效散射单元；

步骤3、采用阶梯函数对全波LiDAR波形中的有效散射单元个数建模：用阶梯的高度表 示各个高程点对应的有效散射单元个数的均值，阶梯的位置表示地物目标的相对高程；

步骤4、利用RJMCMC(可逆跳转马尔科夫链蒙特卡罗)算法对全波LiDAR波形进行分解， 即求解阶梯函数，得到地物目标相对高程以及散射特性。

所述步骤2的具体步骤如下：

步骤2-1、设各采样点处的散射能量中的散射能量最小值c＝min(y_i)为一个有效散射单元；

步骤2-2、计算每个采样点对应的有效散射单元个数z_i＝y_i/c,i＝1,2,…,n；即采样点处的 散射能量y_i除以有效散射单元的散射能量c；每个采样点对应的有效散射单元个数最多为z_max；

步骤2-3、设有效高程区域内每两个采样点间有z_max个高程点，且z_i个有效散射单元随机 分布在z_max个高程点上，则每个高程点对应的有效散射单元个数w_j为0或1，其中，j＝1,2,…, t，j为高程点索引，t为高程点个数，即：t＝z_max×n。

所述步骤3的具体步骤如下：

步骤3-1、将全波LiDAR波形分解为m+1个阶梯，包括m个跳变点和m+1个阶梯高度；

步骤3-2、定义全波LiDAR波形或然率服从均值是阶梯函数的泊松分布，建立全波LiDAR 波形或然率模型，即参数集Φ＝{s,h,m}条件下高程点对应有效散射单元个数的条件概率，其 中，s为跳变点，h为阶梯高度，m为跳变点个数；

步骤3-3、根据贝叶斯定理构建已知高程点对应的有效散射单元个数的条件下参数集Φ＝ {s,h,m}的联合后验概率密度函数。

有益效果：

本发明针对采用有效散射单元来描述全波LiDAR(LightDetectionAndRanging)波形， 将全波LiDAR波形转化为各采样点对应的具有相同散射能量的有效散射单元，能够将地物散 射特性与波形分解相结合，并提出采用阶梯函数对有效散射单元个数建模，打破了使用传统 高斯统计分布模型来拟合原始波形的思想，为全波LiDAR波形建模提供了有效的新思路。

附图说明

图1为本发明具体实施方式的基于地物散射的全波LiDAR波形分解方法流程图；

图2为本发明具体实施方式的步骤2的具体流程图；

图3为本发明具体实施方式的步骤4-1的具体流程图；

图4为本发明具体实施方式的步骤4-2的具体流程图；

图5为本发明具体实施方式的增加阶梯操作实例；

图6为本发明具体实施方式的仿真实验所采用的数据集以及实验区域对应的遥感影像， 其中，图6(a1)至(d1)为原始波形数据，图6(a2)至(d2)为实验区域分别对应的遥感影像，其中 不同灰度代表不同高程的地物目标；

图7为本发明具体实施方式的由原始数据得到的每个采样点对应有效散射单元个数均值 以及由本发明得到的阶梯函数的叠加图，其中，图7(a)至(d)分别为图6(a1)至(d1)的分解结果；

图8为本发明具体实施方式的由数据计算得到的每个阶梯对应散射单元个数与本发明仿 真实验结果得到的每个阶梯对应散射单元个数的叠加图，(a)～(d)分别为不同实验区域的叠加 图；

图9为本发明具体实施方式的10000次迭代过程中后验概率取对数后的变化情况，(a)～(d) 分别为不同实验区域后验概率取对数后的变化情况。

具体实施方式

下面结合附图详细说明本发明的具体实施方式。

基于地物散射的全波LiDAR波形分解方法，如图1所示，包括：

步骤1、激光雷达发射信号到地物目标上反馈回来的信号叠加形成全波LiDAR波形，包 括分布在有效高程区域内的采样点及各采样点处的散射能量。实验数据采用ICESat-GLAS测 高数据，每个波形有544个采样点，前151个采样点间隔增加到4ns(60cm)，后393个采样 点间隔为1ns(15cm)，有效高程区间为[0,L]，L＝150m；

x＝{x_i；i＝1,2,…,n}，y＝{y_i；i＝1,2,…,n}，其中，i为采样点索引，x_i为分布在有效高 程区间[0,L]内的采样点，y_i为在采样点x_i处的散射能量，n＝544为采样点个数。

步骤2、采用有效散射单元来描述全波LiDAR波形：将全波LiDAR波形转化为各采样 点对应的具有相同散射能量的有效散射单元；

根据雷达方程原理将全波LiDAR波形用有效散射单元描述，如图2所示，具体如下：

步骤2-1、设各采样点处的散射能量中的散射能量最小值c＝min(y_i)为一个有效散射单元， c为一个有效散射单元的散射能量；

步骤2-2、计算每个采样点x_i对应的有效散射单元个数z_i＝y_i/c,i＝1,2,…,n；即采样点处 的散射能量y_i除以有效散射单元的散射能量c；每个采样点对应的有效散射单元个数最多为 z_max；

步骤2-3、设有效高程区域内每两个采样点间有z_max个高程点，且z_i个有效散射单元随机 分布在z_max个高程点上，则每个高程点对应的有效散射单元个数w_j为0或1，其中，j＝1,2,…, t，j为高程点索引，t为高程点个数，即：t＝z_max×n。

步骤3、采用阶梯函数对全波LiDAR波形中的有效散射单元个数建模：用阶梯的高度表 示各个高程点对应的有效散射单元个数的均值，阶梯的位置表示地物目标的相对高程；

步骤3-1、将全波LiDAR波形分解为m+1个阶梯，包括m个跳变点和m+1个阶梯高度；

步骤3-1-1、确定跳变点个数m。

假设跳变点个数m服从均值为λ的泊松分布，其概率密度函数为： 其中，m≤m_max，m_max为预先指定的最大阶梯数；

步骤3-1-2、生成m个跳变点s＝{s_k；k＝1,2,…,m}。

假设跳变点s_k(k＝1,2,…,m)服从均匀分布在有效高程区间[0,L]上的2m+1个随机点的偶 数有序统计分布，各跳变点的联合概率密度函数为： p(s)＝p(s₁,s₂,...,s_m)＝(2m+1)！L^-2m-1s₁(s₂-s₁)...(s_m-s_m-1)(L-s_m)，即：在有效高程区间[0,L] 上随机生成2m+1个随机点，把其中偶数点的位置作为跳变点s，跳变点设为0(s₀)<s₁<s₂<… <s_m<L(s_m+1)，其中，k为跳变点索引；

步骤3-1-3、生成m+1个阶梯高度h＝{h_l；l＝0,1,…,m}。

跳变点子区间[s₀,s₁)对应的阶梯高度为h₀，跳变点子区间[s_k,s_k+1)，k＝1,2,…,m对应的 阶梯高度为h_k，假设阶梯高度h_l(l＝0,1,…,m)服从形状参数和尺度参数分别为α和β的 Gamma分布且相互独立，各阶梯高度的联合概率密度函数为： $p\left(h\right)=p(h_0,h_1,...,h_m)=\underset{l=0}{\overset{m}{\Pi }}\frac{h_l^{\alpha -1}}{\Gamma \left(\alpha \right){\beta }^{\alpha }}\exp (-\frac{h_l}{\beta }).$

步骤3-2、建立全波LiDAR波形或然率模型，即参数集Φ＝{s,h,m}条件下高程点对应有 效散射单元个数的条件概率；

定义或然率p(w_j|Φ)(j＝1,2,…,t)为服从均值为g(·)的泊松分布，g(·)为定义在区间[0,L] 上的阶梯函数，在[0,L]区间内有m+1个阶梯，包括m个跳变点s以及m+1个阶梯高度h， 其中，参数集Φ＝{s,h,m}，s＝{s₁,s₂,…,s_m}，h＝{h₀,h₁,…,h_m}。或然率可表示为： $p\left(w\right|\Phi )=\underset{j=1}{\overset{t}{\Pi }}p\left(w_j\right|\Phi )=\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j\right)}^{w_j}}{\left(w_j\right)!}\exp (-g_j).$

步骤3-3、根据贝叶斯定理构建已知高程点对应的有效散射单元个数w的条件下参数集Φ ＝{s,h,m}的联合后验概率密度函数。

根据Bayes定理，已知高程点对应的散射单元个数w条件下参数集Φ＝{s,h,m}的后验概 率密度函数为，

$\left(\begin{array}{c}p\left(\Phi |w\right)\propto \left(\log \left(p\left(w|\Phi \right)\right)\right)p\left(s,h,m\right)=\left(\log \left(p\left(w|\Phi \right)\right)\right)p\left(s|m\right)p\left(h|m\right)p\left(m\right)\\ =\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j\right)}^{w_j}}{w_j}\exp \left(-g_j\right)\times \left(2m+1\right)!L^{-2m-1}{s}_1\left(s_2-s_1\right)\mathrm{...}\left(s_m-s_{m-1}\right)\left(L-s_m\right)\\ \times \underset{l=0}{\overset{m}{\Pi }}\frac{h_l^{\alpha -1}}{\Gamma \left(\alpha \right){\beta }^{\alpha }}\exp \left(-\frac{h_l}{\beta }\right)\times \frac{{\lambda }^m}{m!}\exp \left(-\lambda \right)\end{array}\right)$

步骤4、利用RJMCMC(可逆跳转马尔科夫链蒙特卡罗)算法对全波LiDAR波形进行分解， 即求解阶梯函数，得到地物目标的相对高程以及散射特性。

步骤4-1，更新参数集Φ＝{s,h,m}，包括更新跳变点s和更新阶梯高度h，如图3所示， 具体步骤如下：

步骤4-1-1，更新跳变点s。

从当前m个跳变点中随机抽取一个跳变点s_k作为候选跳变点s_k^*，k∈{1,2,…,m}，该操 作只改变该跳变点的位置，而保持其它参数不变，候选跳变点s_k^*服从区间[s_k-1,s_k+1]上的均匀 分布，h_l和h_l+1对应的候选跳变点区间分别为[s_k-1,s_k^*]和[s_k^*,s_k+1]；

步骤4-1-2，计算更新跳变点s的接受率其中， $R_s=\frac{\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j^{*}\right)}^{w_j}}{w_j}\exp \left(-g_j^{*}\right)}{\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j\right)}^{w_j}}{w_j}\exp \left(-g_j\right)}\times \frac{\left(s_{k+1}-s_k^{*}\right)\left(s_k^{*}-s_{k-1}\right)}{\left(s_{k+1}-s_k\right)\left(s_k-s_{k-1}\right)}.$在[0,1]之间生成一个随机数e，如果a(s_k,s_k^*) ＞e，则接受更新跳变点结果，s_k＝s_k^*；反之，拒绝更新跳变点，且保持原有参数不变；

步骤4-1-3，更新阶梯高度h。

从当前m+1个阶梯高度中随机抽取一个阶梯高度h_l，k∈{0,1,…,m}，该操作只改变该 阶梯高度，而保持其它参数不变，候选阶梯高度为h_l^*，log(h_l^*/h_l)服从区间[-1/2,1/2]上的均匀 分布；

步骤4-1-4，计算阶梯高度h的接受率a(h_l,h_l^*)＝min{1,R_h}，其中， $R_h=\frac{\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j^{*}\right)}^{w_j}}{w_j}\exp (-g_j^{*})}{\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j\right)}^{w_j}}{w_j}\exp (-g_j)}\times {\left(\frac{h_l^{*}}{h_l}\right)}^{\alpha -1}\exp \{-\frac{h_l^{*}-h_l}{\beta }\}.$在[0,1]之间生成一个随机数e，如果a(h_l,h_l^*) ＞e，接受更新阶梯高度结果，h_l＝h_l^*；反之，拒绝更新阶梯高度，且保持原有参数不变。

步骤4-2，增加或删除阶梯，如图4所示，具体步骤如下：

步骤4-2-1，设置阈值参数，包括：设置跳变点个数最大值T_m，最大迭代次数T；

步骤4-2-2，判断执行增加阶梯操作或执行删除阶梯操作：设选择增加阶梯操作的概率和 选择删除阶梯的概率分别为b_m和d_m，当m＝1时，d_m＝0；当m＝T_m时，b_m＝0；当m≠0且 m≠T_m时，b_m＝d_m；在[0,1]之间生成一个随机数u，如果u＞0.5，则执行增加阶梯操作，转 至步骤4-2-3；否则执行删除阶梯操作，转至步骤4-2-4；

步骤4-2-3，执行增加阶梯操作：

首先选取候选跳变点s^*，s^*均匀分布在区间[0,L]上，并以概率1位于区间(s_k,s_k+1)上。如 果该候选跳变点被接受，则将其标记为s_k+1^*，s_k+1,s_k+2,…,s_m重新标号为s_k+2^*,s_k+3^*,…,s_m+1^*， 与之对应阶梯高度的标号做同样的操作。对于子区间(s_k,s_k+1^*)和(s_k+1^*,s_k+2^*)生成候选阶梯高度 h_l^*和h_l+1^*。

设跳变前后变化的区间及阶梯高度满足加权几何平均关系： $(s^{*}-s_k)log\left(h_l^{*}\right)+(s_{k+1}-s^{*})log\left(h_{l+1}^{*}\right)=(s_{k+1}-s_k)log\left(h_l\right),$并假定h_l^*和h_l+1^*间存在微小的扰动，使 得h_l+1^*/h_l^*＝(1-v)/v，其中，v为均匀分布于[0,1]区间的随机数。生成阶梯的接受概率为： a_b＝min{1,R_b}，其中，

$\left(\begin{array}{c}{R}_b=\frac{p\left(w|{\Phi }^{*}\right)p\left({\Phi }^{*}\right)r_b\left({\Phi }^{*}\right)}{p\left(w|\Phi \right)p\left(\Phi \right)r_d\left(\Phi \right)q\left(u\right)}\left|\frac{\partial {\Phi }^{*}}{\partial (\Phi ,u)}\right|\\ =\frac{\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j^{*}\right)}^{w_j}}{w_j}\exp \left(-g_j^{*}\right)}{\underset{j=1}{\overset{t}{\Pi }}\frac{{\left(g_j\right)}^{w_j}}{w_j}\exp \left(-g_j\right)}\times \frac{\lambda }{m+1}\times \frac{2\left(m+1\right)\left(2m+3\right)}{L^2}\frac{\left(s^{*}-s_k\right)\left(s_{k+1}-s^{*}\right)}{s_{k+1}-s_k}\\ \times \frac{1}{\Gamma \left(\alpha \right){\beta }^{\alpha }}{\left(\frac{h_l^{*}{h}_{l+1}^{*}}{h_l}\right)}^{\alpha -1}\exp \left\{-\frac{h_l^{*}+h_{l+1}^{*}-h_l}{\beta }\right\}\times \frac{d_{m+1}L}{b_m\left(m+1\right)}\times \frac{{\left(h_l^{*}+h_{l+1}^{*}\right)}^2}{h_l}\end{array}\right)$

在[0,1]之间生成一个随机数e，如果a_b＞e，接受增加阶梯操作结果，则s^*＝{s₁,…,s_k,s_k+1^*, s_k+2^*,s_k+3^*,…,s_m+1^*}，h^*＝{h₀,h₁,…,h_l^*,h_l+1^*,…,h_m^*}。

步骤4-2-4，执行删除阶梯操作：

在m个跳变点中随机抽取一个跳变点s_k+1，如果删除该跳变点被接受，顺序改变后跳变 点的标号为s_k+1^*,s_k+2^*,…,s_m^*，使得(s_k^*,s_k+1^*)＝(s_k,s_k+2)，其概率为1/m+1。子区间(s_k^*,s_k+1^*) 对应的新高度为h_l^*，并满足如下关系： (s_k+1-s_k)log(h_l)+(s_k+2-s_k+1)log(h_l+1)＝(s_k+2-s_k)log(h_l^*)。

删除阶梯操作是生成阶梯操作的对偶操作，因此其接受概率为：在[0,1] 之间生成一个随机数e，如果a_d＞e，接受删除阶梯操作结果，则s^*＝{s₁,…,s_k,s_k+1^*,…,s_m^*}， h^*＝{h₀,h₁,…,h_l-1,h_l^*,…,h_m^*}。

图5为增加阶梯操作的实例，其中s^*、h_l^*和h_l+1^*分别为新增加的跳变点以及增加跳变点 后新生成的两个阶梯高度。

步骤4-3，当迭代指次数达到最大迭代次数T时，停止迭代，反之，重复步骤4-1，继续 迭代。

步骤4-4，构成阶梯的两个跳变点之间的中点为该阶梯对应的地物目标的相对高程，阶梯 的高度为高程点对应的有效散射单元个数的均值，有效散射单元个数的均值与高程点个数的 乘积为有效散射单元个数，有效散射单元个数与根据遥感影像得到的实际地物目标的面积之 比为地物的散射系数，即散射特性。

全波LiDAR波形的最后一个返回波为地面的返回波，设其相对高程为0，得到其它地物 目标的相对高程，当前阶梯高度即有效散射单元个数的均值与高程点个数的乘积为有效散射 单元个数。根据仿真实验得到的有效散射单元个数与根据遥感影像得到的实际地物目标的面 积之比为地物的散射系数。若两个地物目标的散射特性相近则认为两个地物目标相同(材质 或类型相同)。

通过实验，进一步说明本发明的有效性：

1.实验环境：在CPU为Corei5-34703.20GHz、内存4GB、Windows7旗舰版系统上使 用MATLAB7.12.0软件编程实现仿真。

2.选择4个实验区域进行波形分解的仿真实验，并将分解结果与原始数据进行叠加，以 及将实验得到地物目标的相对高程和散射系数与遥感影像及实地测量得到的对应地物相比 较，对实验结果进行评价。

3.如图6所示，图6(a1)至(d1)为4个实验区域((a1)至(d1)分别区域Ⅰ、区域Ⅱ、区域 Ⅲ、区域Ⅳ)的原始波形数据，图6(a2)至(d2)为4个实验区域((a1)至(d1)分别区域Ⅰ、区域 Ⅱ、区域Ⅲ、区域Ⅳ)对应的遥感影像，其中不同灰度代表不同目标地物。

4.如图7所示，为每个实验区域由原始数据得到每个采样点对应有效散射单元个数与由 本发明得到阶梯函数的叠加图，其中，图7(a)至(d)分别为图6(a1)至(d1)的分解结果。可以看 出，由本发明得到的阶梯函数跳变点的位置能够和原始数据重合，阶梯函数高度与数据得到 的每个采样点对应的散射单元个数趋势相同；根据图6(a2)至(d2)的遥感影像可以看出本发明 的仿真实验得到阶梯函数的每个阶梯对应一个地物目标，但是地物有高低起伏变化时会被分 成多个阶梯，如图7(a)和(b)的最后三个阶梯对应的是有坡度的地面。

5.如图8(a)～(d)所示，为每个实验区域由数据计算得到的每个阶梯对应散射单元个数与 本发明实验结果得到的每个阶梯对应散射单元个数的叠加图。可以看出，本发明得到的每个 阶梯对应散射单元的个数与数据每个回波对应有效散射单元个数基本相同。

6.如图9所示，为每个实验区域10000次迭代过程中后验概率取对数后的变化情况，(a)～ (d)分别为不同实验区域后验概率取对数后的变化情况。可以看出，本发明方法在迭代过程中 后验概率很快收敛，并趋于稳定。

7.为了精确验证本发明方法有效性，根据实验结果计算每类地物目标对应的相对高程以 及散射系数。计算结果如表1所示，

表1每类地物目标对应的相对高程以及散射系数

可以看出，本发明方法得到的相对高程与经过实地勘测得到的实际地物高程基本一致； 并且相同材质的目标地物散射系数基本相同，区域Ⅰ和区域Ⅱ白色地物的沥青屋顶以及黑色 地物石灰地面，区域Ⅲ浅灰地物和区域Ⅳ白色地物的水泥屋顶，如图6所示。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何本 技术领域人士在本发明所揭露的技术范围内，可理解到的替换，都应涵盖在本发明的包含范 围之内，例如基于有效散射单元的波形分解、基于散射系数的波形分解、基于阶梯函数的波 形分解等。因此，本发明的保护范围应该以权利要求书的保护范围为准。