一种改进的全波形LiDAR数据贝叶斯分解算法

一种改进的全波形lidar数据贝叶斯分解算法
技术领域
1.本发明涉及激光探测技术领域，特别涉及机载激光雷达波形分解算法。


背景技术：

2.目前，激光雷达系统(light detected and ranging，简称lidar，)主要由脉冲发射和接收装置、信号处理部分组成，而机载激光雷达系统包括激光雷达系统、全球定位系统、惯性导航系统还有数码相机等光谱成像设备等，简称alb(airborne lidar bathymetry)。机载激光雷达系统的历史始于hickman和hogg在1969年首次测试激光脉冲在近岸水深测量中的能力。机载激光雷达技术成本低廉、高效，备受如林业、城市三维建模，水利电网建设等领域所喜爱，是研究者对大面积地表进行研究的理想工具，在陆地上的运用已十分广泛。机载激光雷达系统通常采用1064nm的红外激光系统进行测量，但随着450nm～580nm“海洋窗口”的发现及海洋开发的需求，其中532nm的蓝绿激光波段能很好的穿透水体，一些系统便开始使用该波段。如此，能同时满足水深和陆地测量需要的双频机载激光雷达系统应时而生，如vq和czmil系列。专业的激光雷达系统软件包括realm和terrascan等，它们在水深小于50米的浅水区具有很大的优势，可以在近海岸、暗礁、岩礁或是测量船无法行驶的水域完成高效而快速的测量工作。在我国，水深50米以内的海域面积有近50万平方公里，如此宽广的海域面积都适合采用双频机载激光雷达系统进行测量，其中蕴含着巨大的经济价值。同时，双频机载激光雷达系统可以在陆地和水域的测量中无缝衔接，是一种极具诱惑力的新技术。1976年，nabrit等人的实验表明，对同样大小的海底地形进行测量，传统船载声呐系统所需费用是机载激光雷达系统的五倍。但是激光雷达在水深测量和水下目标探测等方面的应用目前还处在起步阶段，双频激光雷达回波数据处理的研究较为缺乏，相关的研究主要集中在欧美等早期就开始研究激光雷达的发达国家。
3.早前，机载激光雷达系统只能记录离散的回波信号，通常为两次，并且用户不能使用回波的强度信息。riegl公司在2004年率先发布了一台全波形激光雷达系统，该系统一经发布便能投入到商业化使用。全波形激光雷达系统的物理原理与传统的激光雷达系统类似，但全波形激光雷达系统能够通过数字采样，记录不同时间分辨率(如1/2/4纳秒)被照射物体散射的全部回波。因此，这些波形不仅包括地面的响应，而且还包括沿激光线被照射表面的多次散射响应。更为重要的是系统记录了激光回波波形的特征信息，它不仅仅包含了反射地物的空间位置，还包含了一些特征值如波形的振幅、脉宽等信息。通过分析这些特征信息可以计算出目标的空间分布或者几何关系如波宽越大反映着树木越高，树木冠层可以根据距离进行计算，还可以用于分析地物地貌的形状、反射率及相互关联关系等，为激光雷达分割或分类算法提供目标的特征信息。在海洋测量中，全波形数据处理后能获取水深数值，还能获取脉冲传输中的各种媒介信息。这些信息可以服务于海洋环境的保护、海浪潮汐的检测、海岸线监测等相关海洋研究领域。
4.由此可见，与传统的激光雷达系统相比，全波形激光雷达系统能够对回波信号进行完整而密集的存储，数据加工后能提供更为精确、高密度和信息丰富的点云。通过对全波
形回波进行波形分解后可以获取更加精确的点云坐标以及更加丰富目标特征信息，尤其是在森林地区，它能提供更多的地物信息。虽然双频机载全波形激光雷达技术的开发异常艰难，但其优越的性能以及科研和社会价值极具竞争力，国家和科研人员应竭尽所能去共克这些难关。全波形激光雷达系统的发展也让波形处理算法受到更多人的关注。早在2000年，guenther等人就阐述了精度对于机载激光雷达系统的重要性，而数据的后处理算法的研究是其突破点。mallet等人在2010年也指出激光雷达信号重建是一个非常重要的研究课题。经过数十年的发展，已经发表了不少波形分解算法，但这些算法仍很难提取完整的全波形信息，从而导致信息丢失。
5.在这一背景下，本发明的主要目标是基于实测激光雷达数据对已有的全波形数据分解算法进行改进，通过定性和定量的对比分析，结合贝叶斯算法优化波形分解模型，提出基于优化马尔可夫链蒙特卡罗(markov chain monte carlo，简称mcmc)算法的全波形lidar数据贝叶斯分解方法，为促进算法的实践应用提供参考价值。


技术实现要素：

6.本发明公开了机载激光雷达波形分解算法改进方法。主要包括以下步骤。
7.第一，针对具有动态范围大、非稳态特点的lidar信号，基于传统阈值法构建新的一阶可导阈值函数，开发一种适用于双频激光雷达系统降噪的自适应阈值方法，为后续的波形分解进行有效去噪。
8.第二，基于mcmc抽样的全波形lidar数据贝叶斯分解算法优化。其中，先进行mcmc采样。具体步骤如下：
9.1)输入：设定状态转移次数阈值n1，样本个数n2，平稳分布π(x)，任意选定的马尔科夫链状态转移矩阵q。
10.2)通过峰值识别算法得到初始状态值x0(参数a/u/σ/λ等同理)。
11.3)for t＝0到n1 n2-1:
12.a)从条件概率分布q(x|x
t
)中采样得到样本x
*
；
13.b)从均匀分布采样u～uniform[0,1]；
[0014]
c)如果u＜α(x
t
,x
*
)＝π(x
*
)q(x
*
,x
t
)，则接受转移x
t
→
x
*
，即x
t
1＝x
*
；
[0015]
d)否则不接受转移，即x
t
1＝x
t
；
[0016]
那么，样本集(x
n1
,x
n1 1
,...,x
n1 n2-1
)就是所求的平稳分布对应的样本集。
[0017]
优化收敛条件。当接受率过小如0.1，这导致数十万次迭代后可能还未收敛。针对激光雷达通常存在海量数据需要处理的情况，这样的采样效率很低，因此进行方法改进，假设a(i,j)为0.1，a(j,i)
[0018]
为0.2，代入前面的条件公式中得到：
[0019]
π(i)q(i,j)
×
0.1＝π(j)q(j,i)
×
0.2
ꢀꢀꢀ
(1)
[0020]
两边同时扩大2倍后，等式同样成立：
[0021]
π(i)q(i,j)
×
0.2＝π(j)q(j,i)
×4ꢀꢀꢀ
(2)
[0022]
那么，针对a(i,j)可以做如下改进：
[0023]
[0024]
假设马尔科夫链状态转移矩阵q对称，即q(i,j)＝q(j,i)，公式简化为：
[0025][0026]
进行延迟拒绝方法改进。具体过程如下：
[0027]
1)输入：设定状态转移次数阈值n1，需要的样本个数n2，平稳分布π(x)，任意选定的马尔科夫链状态转移矩阵q；
[0028]
2)通过峰值识别算法得到初始状态值x0(参数a、u、σ和λ等参数同理)；
[0029]
3)for t＝0到n1 n2-1:
[0030]
a)从条件概率分布q(x|x
t
)中采样得到样本x
*
；
[0031]
b)从均匀分布采样u～uniform[0,1]；
[0032]
c)如果则接受转移x
t
→
x
*
，即x
t
1＝x
*
；
[0033]
d)否则不接受转移，即x
t
1＝x
t
；
[0034]
e)延迟拒绝。
[0035]
通过上面系列步骤的计算，最终得到符合平稳分布的样本集
[0036]
(x
n1
,x
n1 1
,...,x
n1 n2-1
)。
[0037]
第三，基于mcmc算法的全波形lidar数据贝叶斯分解模型构建。
[0038]
其中，进行先验分布。通过mcmc采样可以得到未知参数如峰值a、峰值位置μ和波宽σ等的样本集。以参数μ为例，将分解得到的各个波形分量的位置按时间先后排序可得0＜μ1＜μ2＜...＜μm＜τ，如果μk(k＝1,2,...,m)在区间[0,τ]上服从均匀分布，而且彼此之间是独立的。那么，可以通过下式求得概率密度函数：
[0039][0040]
如果波宽σk(k＝1,2,...,m)服从(μ
σ
，σ
σ2
)正态分布，并且彼此之间独立，那么可以通过下式求得概率密度函数：
[0041][0042]
如果是在高斯混合模型中，会在不同的高斯类函数前增加权重π，分解得到的波形分量的权重πk(k＝1,2,...,m)在区间[0,1]上如果服从均匀分布，并且彼此之间独立，那么可以通过下式求得概率密度函数：
[0043][0044]
如果波形分量的数目服从方差和期望都为λ的泊松分布(poisson distribution)，那么可以通过下式求得概率密度函数：
[0045][0046]
式中，2≤m≤m
max
，m
max
为波形分量数目的最大值。
[0047]
至此，通过上述方法可以获取未知参数a、u和σ等的先验分布，通过贝叶斯算法进
行求解。根据贝叶斯原理，非线性的激光雷达模型可以用下式表示：
[0048]
yi＝f(xi,β)*εiꢀꢀꢀꢀꢀꢀꢀ
(9)
[0049]
式中，yi是观测数据，f(xi,β)是含有未知参数β和预测因子xi的非线性函数，εi为服从(0，τ2)分布的独立误差，τ是logεi的标准差，xi是lidar波形的第i个采样点。
[0050]
在贝叶斯框架下，通过线性或非线性等数学方程来指定确定性模型，未知的模型参数则由各种概率分布随机处理。以高斯模型为例来进行贝叶斯分解在lidar数据处理中的应用。高斯模型可以用下式表达：
[0051][0052]
由于要求解未知的模型参数，在此对公式中的yi进行对数变换，yi的可能性函数的对数形式如下式：
[0053]
logεi＝logy
i-logf(xi,β)
ꢀꢀꢀꢀ
(11)
[0054][0055]
通过对波形原始数据的统计，可以缩小这些参数的合理范围：a[10,150]，u[15,100]，σ[4,15]。为比较先验信息对模型性能造成的影响，可以为每个参数分配一个均匀分布，aj～u(10,100)、uj～u(10,150)、σj～u(4,15)。本发明采取峰值检测算法获取参数的先验分布，其中由于σj难以取决，采用连续峰间差值的三分之一来表示σj的先验信息。综合来看，将aj、uj、σj的先验分布设定为aj～n(aj,102)、uj～n(uj,52)、σj～n(σj/3,32)和τ～n(0,0.52)的正态分布。那么，通过前面求得的所有兴趣参数(a、u和σ)的先验分布和lidar波形数据分布p(y|x,β,τ)可以求得后验分布的对数形式，后验分布即为波形数据的最可能函数，如下式：
[0056][0057][0058]
式中，m代表每个波形的采样数，n代表波形分量的数量，1/τm是波形分量的概率密度函数。
[0059]
本发明的有益之处是为量化激光雷达波形数据获取过程中产生的误差以及波形分解导致的模型不确定性，提出了基于mcmc算法的贝叶斯分解改进方法，并结合在激光雷达数据处理中发现的问题对算法进行了优化，包括收敛条件的优化和延迟拒绝。利用mcmc模拟通过峰值检测法获得的参数，代入贝叶斯推断中即可得到参数的估计值。
附图说明
[0060]
图1是本发明的流程图
[0061]
图2是本发明的基于mcmc算法的全波形激光雷达数据贝叶斯分解优化方法采样示意图。
[0062]
图3是本发明的五种算法对不同地区的数据进行波形分解所耗费的时间。
具体实施方式
[0063]
为使本发明的目的、技术方案及优点更加清楚明白，以下举出优选实施例，结合附图对本发明具体实施作进一步详细说明。
[0064]
实施例：
[0065]
结合图2，说明本发明机载激光雷达波形分解算法改进方法实例结果：
[0066]
第一，针对具有动态范围大、非稳态特点的lidar信号，基于传统阈值法构建新的一阶可导阈值函数，开发一种适用于双频激光雷达系统降噪的自适应阈值方法，为后续的波形分解进行有效去噪。
[0067]
应用500组波形样本数据进行自适应阈值法处理。结果表明(见表1)，本发明降噪方法在激光雷达模拟数据中比传统硬阈值法提高了近2db的信噪比(signal-to-noise ratio，snr)，比软阈值法提高了近4db的信噪比，提高了信号特征的可分离性。
[0068]
表1三种阈值函数的去噪指标
[0069][0070]
第二，基于mcmc抽样的全波形lidar数据贝叶斯分解算法优化。
[0071]
基于mcmc抽样，使用贝叶斯分解方法对500组的样本波形处理的实验结果表明(见表2)。只有在高斯模型中个别波形出现发散的情况(的收敛率为99.4％)，整体上三种模型在一致的实验参数下均能很好的收敛。其中，高斯模型参数估计值的标准差稍大于自适应高斯模型，但其整体的waic值和误差的均值最小，分别为777.6和11.9。虽然威布尔模型也能完全收敛，但其处理时间为高斯算法的4倍之多，达29168秒。此外，当使用自适应高斯或威布尔模型时，参数的物理意义也受到限制。gold和rl两种反卷积模型分解结果的waic值都比另外三种数学模拟法大，分别为1180.6和1004.4，其原因可能是反卷积后波形的振幅值发生了较大的变化。其中，gold模型的处理时间与高斯模型相似，具有较快的处理速度，而rl模型的迭代时间接近高斯模型的4倍，达26435秒。
[0072]
表2贝叶斯框架下高斯、自适应高斯和威布尔模型的分解结果统计
[0073][0074]
分解的实验结果表明，威布尔和自适应高斯模型都需要花费更多的时间来寻找最
优参数以及使额外的参数取得收敛。从精度、不确定性、物理意义和处理效率等方面分析，高斯模型是全波形激光雷达数据分解的最佳模型。
[0075]
第三，基于mcmc算法的全波形lidar数据贝叶斯分解模型构建。
[0076]
选取了描述波形分量的参数：峰值a、峰值的位置μ和波宽σ，形状参数λ和k，并选取了均方根误差(root mean squared error，rmse)、相关系数c和拟合优度r2作为评价指标来定量评估高斯、自适应高斯、威布尔、richardson-lucy(rl)和gold等五种算法分解结果的精度，结果见表3。
[0077]
表3五种算法的分解精度比较(数据：奥果韦)
[0078][0079]
结果表明：
[0080]
1)高斯算法的拟合误差最大(9.96)，gold算法的误差最小，仅为4.802。这结果表明gold算法具有很强的鲁棒性。rl算法的误差为7.424，相比gold算法偏高的原因可能是该算法对噪声更加敏感。
[0081]
2)weibull算法的拟合信号与原始信号的相关系数在五种算法中最低，仅为0.71，自适应高斯算法具有最高的相关系数(0.851)。这结果进一步表明自适应高斯算法分解波形的准确性。
[0082]
3)五种算法的拟合优度都达到了0.9以上(越接近1效果越好)。这结果表明5种算法都具有很好的拟合效果。
[0083]
基于激光雷达数据运行后，五种算法的分解速度对比结果如图3所示。结果表明：
[0084]
1)在处理一万条以下的数据时，五种算法的处理速度的差异可以忽略不计；而在波形数量达到上万条时开始有了明显差异，这个阈值因不同的计算平台或硬件水平而变化。
[0085]
2)rl在运算数据处理10万条数据在4个地区所花费的时间分别为1.8万秒、3.5万秒、4万秒和3.5万秒，所耗费的时间几乎是最多的；weibull算法随着数据量的增加，其处理速度迅速下降，在奥果韦地区甚至超过了rl算法，达4.4万秒；这些结果表明weibull和rl算法并不适用于实时计算的lidar系统
[0086]
3)高斯算法在处理路易斯安那州和奥果韦地区10万条数据时所花费的时间最少，分别为0.4万秒和0.7万秒；而在拉森港和耶洛奈夫地区高斯算法花费时间最少的是gold算法(0.1万秒和0.2万秒)。两种算法都具有很快的信号分解速度。
[0087]
以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变形。因此，所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。
[0088]
本发明未详尽描述的技术内容均为公知技术。