本發(fā)明涉及一種對激光雷達回波波形數(shù)據(jù)處理的方法，是一種通過LM算法實現(xiàn)激光雷達波形數(shù)據(jù)分解的方法，屬于機載激光雷達技術(shù)領(lǐng)域。

背景技術(shù)：

水深測量是海洋測繪的重要內(nèi)容。傳統(tǒng)的船載多波束回聲測深方法雖然具有較高的測量精度，但無法進入沿岸淺水區(qū)域和海礁密集區(qū)域。近年來，機載激光雷達的出現(xiàn)和發(fā)展填補了沿海淺水區(qū)域水深測量技術(shù)的空白，已成為一種快速高效的水深及海底地形探測方法。機載激光雷達探測系統(tǒng)可以廣泛用于海洋水文勘測(包括海水深度測量、海底地貌測繪和海水光學(xué)參數(shù)遙測等)、水雷和魚群的探測、海洋環(huán)境污染監(jiān)測等眾多領(lǐng)域，因而在海洋科學(xué)研究中有巨大的作用。

機載激光雷達系統(tǒng)在海洋探測領(lǐng)域有著探測精度高、探測效率高、不受海水溫度和鹽度影響等優(yōu)勢。機載激光雷達系統(tǒng)的工作原理是機載激光發(fā)射器向海面發(fā)射激光，通過計算海面回波和海底回波的時間差求解出海水深度。為了提高水深測量精度，需求解出激光的精確回波時間，但是由于海水的組成十分復(fù)雜，增加了對回波波形的數(shù)據(jù)處理難度，無法求出精確的回波時間。因此對波形數(shù)據(jù)的處理方法就變得尤為重要。

目前波形數(shù)據(jù)的處理算法主要有三種：一是閾值法，二是形心法，三是波形分解法。閾值法在面對較為復(fù)雜的回波時，由于未考慮光斑內(nèi)不同目標(biāo)物回波間的相互影響，使得處理結(jié)果的準(zhǔn)確性上存在一定問題。對于形心法，當(dāng)待測目標(biāo)的回波波形出現(xiàn)上升沿或者下降沿變緩、波形發(fā)生展寬和畸變時，可能會出現(xiàn)較大的偏差。波形分解法將不同高度處的目標(biāo)物回波從接收的波形信號中提取出來，采用高斯函數(shù)近似刻畫單個目標(biāo)的回波波形。波形分解法可以將由于水深較淺而重疊為一個波峰的海面與海底回波分解以判斷海面與海底回波的對應(yīng)位置。傳統(tǒng)方法使用EM(Expectation Maximization)算法進行波形分解，但是EM算法存在收斂速度較慢等問題。

技術(shù)實現(xiàn)要素：

本發(fā)明使用LM算法進行波形分解。LM(Levenberg-Marquardt)算法是一種利用標(biāo)準(zhǔn)數(shù)值優(yōu)化技術(shù)的快速算法，既有高斯-牛頓法的局部收斂性，又有梯度下降法的全局特性。在設(shè)定合理初值的情況下，能得到可靠的最優(yōu)參數(shù)。本發(fā)明中，當(dāng)?shù)玫郊す獾幕夭úㄐ螘r，有效地將海面回波與海底回波分解為兩個波形，利用兩個回波的時間差求出海水深度。

為了解決上述技術(shù)問題，本發(fā)明提出的一種基于LM算法分解激光雷達波形確定海水深度的方法，通過激光雷達接收獲得海面回波與海底回波構(gòu)成的總體回波，并包括以下步驟：

步驟一、波形預(yù)處理，包括波形的噪聲去除和波形平滑；其中，波形的噪聲去除中在進行噪聲估計時，將總體回波數(shù)據(jù)采集點的前2.5％的點和后2.5％的點幅值的平均值作為噪聲的均值；選擇高斯函數(shù)對上述總體回波的波形進行平滑處理；

步驟二、對于海面回波及海底回波波形參數(shù)的初始化，包括：

對于海面回波，采用式(1)所示的高斯函數(shù)與指數(shù)函數(shù)的卷積模擬海面回波波形，以增加時間常數(shù)τ以模擬后向散射對海面回波波形的影響；

式(1)中的初始參數(shù)包括：函數(shù)幅值h_G、時間常數(shù)τ、海面回波波峰的對應(yīng)時間t_G及高斯函數(shù)的標(biāo)準(zhǔn)差σ_G；自變量為t

總體回波最大峰值點對應(yīng)的幅值為A，總體回波最大峰值點對應(yīng)時間為t_p，求出幅值為A'＝α×A處所對應(yīng)的時間t_a和時間t_b，t_a＜t_b，幅值系數(shù)α＝0.1；

設(shè)：

W_α＝t_b-t_a (2)

A_α＝t_p-t_a (3)

B_α＝t_b-t_p (4)

式中：t_b與t_a的時間間隔為W_α、t_p與t_a的時間間隔為A_α、t_b與t_p的時間間隔為B_α及中間參數(shù)μ；

則

t_G＝t_p-σ_G[-0.193(B_α/A_α)²+1.162(B_α/A_α)-0.545] (8)

將式(6)、式(7)和式(8)代入式(1)中，求解出函數(shù)幅值h_G；

對于海底回波，使用式(9)所示高斯函數(shù)代表海底回波信號，

y_G(t)＝A_maxexp[-(t-t_max)²/2σ²] (9)

式(9)中的初始參數(shù)包括：高斯函數(shù)的最大幅值A(chǔ)_max，海底回波波峰的對應(yīng)時間t_max，σ標(biāo)準(zhǔn)差；

將總體回波波形減去初始的海面回波波形得到海底回波波形，海底回波波形的最大峰值點的波峰對應(yīng)的幅值為A_max，海底回波波形對應(yīng)的時間為t_max，若激光雷達發(fā)射脈沖的半寬為ω，則取

步驟三、對步驟二初始化后的參數(shù)進行最優(yōu)化：

利用LM算法通過迭代對初始波形參數(shù)進行優(yōu)化，設(shè)x為含有n個參數(shù)的向量；e_k(x)＝y(tǒng)(x，t_k)-y(t_k)為每個采樣點的殘差，t_k為在第k個采樣點的時間；目標(biāo)函數(shù)為

E(x)的二次模型為m(x+s)＝E(x)+g(x)^Ts+0.5s^TH(x)s,其中s為步長，g(x)為E(x)的梯度，H(x)為E(x)的Hessian矩陣；

設(shè)e(x)的雅克比矩陣為J(x)，e(x)＝[e₁(x)，…,e_k(x)，…，e_m(x)]^T，其中m為采樣點的數(shù)量；

則E(x)＝e(x)^Te(x)，g(x)＝2J(x)^Te(x)，H(x)＝2J(x)^TJ(x)；

對于LM算法，步長為s(λ)＝-[λI+J(x)^TJ(x)]-¹J(x)^Te(x)，λ為阻尼參數(shù)，I為n階單位矩陣，n為參數(shù)個數(shù)；

所述LM算法的具體過程如下：

3-1)設(shè)置參數(shù)組成的向量x₀及置信半徑初值δ，求E(x₀)；

3-2)設(shè)λ＝0，求s(λ)和||s(λ)||₂；

3-3)若||s(λ)||₂≤1.5δ，則設(shè)調(diào)整向量s＝s(λ)，δ＝min{δ，s(λ)}，轉(zhuǎn)到步驟3-5)；

3-4)若||s(λ)||₂＞1.5δ，尋找一個當(dāng)次迭代阻尼參數(shù)λ_k，使||s(λ_k)||₂∈[0.75δ，1.5δ]，設(shè)s＝s(λ_k)；

3-5)設(shè)調(diào)整后的參數(shù)向量x₊＝x₀+s，求E(x₊)，計算實際變化量ΔE＝E(x₊)-E(x₀)；

3-6)若ΔE≥0，則取δ∈[0.1δ，0.5δ]，并返回步驟3-4)；

3-7)計算m(x₊)，設(shè)預(yù)計變化量ΔE_pred＝m(x₊)-E(x₀)，實際變化量與預(yù)計變化量的比值R＝ΔE/ΔE_pred若R≤0.1，則δ＝0.5δ，若R≥0.75，則δ＝2δ，否則δ不變；

3-8)若|s|/(0.001+|x₀|)≥0.1，則設(shè)x₀＝x₊，并返回步驟3-2)，否則最優(yōu)化的參數(shù)x_f＝x₊，E_f＝E(x₊)，停止計算，得到最優(yōu)化的參數(shù)x_f，并將該最優(yōu)化的參數(shù)x_f分別帶入式(1)和式(9)，從而獲得海面回波波形與海底回波波形；

步驟四、求解海水深度D：通過步驟三得到的海面回波波峰的對應(yīng)時間t_G和海底回波波峰的對應(yīng)時間t_max求出海水深度其中c為光速，n_w為海水折射率。

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果是：

利用LM算法可以得到最優(yōu)參數(shù)，實現(xiàn)海面回波與海底回波的分解，大大減少算法復(fù)雜程度。在海面回波與海底回波十分接近的情況下(如圖2)，仍可以分解出海面波形和海底波形。

附圖說明

圖1是本發(fā)明中總體回波波形分解流程圖；

圖2是本發(fā)明中對海面回波進行初始參數(shù)的估計；

圖3是本發(fā)明中對參數(shù)進行最優(yōu)化后獲得的波形分解效果圖。

具體實施方式

下面結(jié)合附圖和具體實施例對本發(fā)明技術(shù)方案作進一步詳細描述，所描述的具體實施例僅對本發(fā)明進行解釋說明，并不用以限制本發(fā)明。

本發(fā)明提出的一種基于LM算法分解激光雷達波形確定海水深度的方法，其設(shè)計思路如圖1所示，在得到回波波形后，進行回波波形分解操作，即，首先對波形進行降噪和平滑處理；然后，根據(jù)回波波形設(shè)定初始波形參數(shù)；其次，通過LM算法擬合，在每次迭代調(diào)整參數(shù)，當(dāng)?shù)嬎惴辖Y(jié)束條件時，求得最優(yōu)參數(shù)；最后，將回波波形分解為海面波形與海底波形，進而求解出海水深度。具體步驟如下：

步驟一、波形預(yù)處理，包括波形的噪聲去除和波形平滑；其主要目的是減少波形中的噪聲對參數(shù)初值設(shè)定的干擾。波形的噪聲去除中在進行噪聲估計時，將總體回波數(shù)據(jù)采集點的前2.5％的點和后2.5％的點幅值的平均值作為噪聲的均值(若有2000個采樣點則取前50和后50個采樣點求解平均值)；需選擇高斯函數(shù)作為濾波器，因此，選擇高斯函數(shù)對上述總體回波的波形進行平滑處理；

步驟二、對于海面回波及海底回波波形參數(shù)的初始化，包括：

對于海面回波：

采用式(1)所示的高斯函數(shù)與指數(shù)函數(shù)的卷積模擬海面回波波形，以增加時間常數(shù)τ以模擬后向散射對海面回波波形的影響；

式(1)中的初始參數(shù)包括：函數(shù)幅值h_G、時間常數(shù)τ、海面回波波峰的對應(yīng)時間t_G及高斯函數(shù)的標(biāo)準(zhǔn)差σ_G；自變量為t

如圖2所示，求出的總體回波最大峰值點對應(yīng)的幅值為A，總體回波最大峰值點對應(yīng)時間為t_p，求出幅值為A'＝α×A處所對應(yīng)的時間t_a和時間t_b，t_a＜t_b，幅值系數(shù)α＝0.1；

設(shè)：

W_α＝t_b-t_a (2)

A_α＝t_p-t_a (3)

B_α＝t_b-t_p (4)

式中：t_b與t_a的時間間隔為W_α、t_p與t_a的時間間隔為A_α、t_b與t_p的時間間隔為B_α及中間參數(shù)μ；

則

t_G＝t_p-σ_G[-0.193(B_α/A_α)²+1.162(B_α/A_α)-0.545] (8)

將式(6)、式(7)和式(8)代入式(1)中，求解出函數(shù)幅值h_G。

對于海底回波：

使用式(9)所示高斯函數(shù)代表海底回波信號，可以減少算法復(fù)雜度。

y_G(t)＝A_maxexp[-(t-t_max)²/2σ²] (9)

式(9)中的初始參數(shù)包括：高斯函數(shù)的最大幅值A(chǔ)_max，海底回波波峰的對應(yīng)時間t_max，σ標(biāo)準(zhǔn)差；

將總體回波波形減去初始的海面回波波形得到海底回波波形，海底回波波形的最大峰值點的波峰對應(yīng)的幅值為A_max，海底回波波形對應(yīng)的時間為t_max，若激光雷達發(fā)射脈沖的半寬為ω，則取

步驟三、對步驟二初始化后的參數(shù)進行最優(yōu)化：

利用LM算法通過迭代對初始波形參數(shù)進行優(yōu)化，設(shè)x為含有n個參數(shù)的向量；e_k(x)＝y(tǒng)(x，t_k)-y(t_k)為每個采樣點的殘差，t_k為在第k個采樣點的時間；目標(biāo)函數(shù)為

E(x)的二次模型為m(x+s)＝E(x)+g(x)^Ts+0.5s^TH(x)s,其中s為步長，g(x)為E(x)的梯度，H(x)為E(x)的Hessian矩陣；

設(shè)e(x)的雅克比矩陣為J(x)，e(x)＝[e₁(x)，…,e_k(x)，…，e_m(x)]^T，其中m為采樣點的數(shù)量；

則E(x)＝e(x)^Te(x)，g(x)＝2J(x)^Te(x)，H(x)＝2J(x)^TJ(x)；

對于LM算法，步長為s(λ)＝-[λI+J(x)^TJ(x)]-¹J(x)^Te(x)，λ為阻尼參數(shù)，I為n階單位矩陣，n為參數(shù)個數(shù)；

所述LM算法的具體過程如下：

3-1)設(shè)置參數(shù)組成的向量x₀及置信半徑初值δ，求E(x₀)；

3-2)設(shè)λ＝0，求s(λ)和||s(λ)||₂；

3-3)若||s(λ)||₂≤1.5δ，則設(shè)調(diào)整向量s＝s(λ)，δ＝min{δ，s(λ)}，轉(zhuǎn)到步驟3-5)；

3-4)若||s(λ)||₂＞1.5δ，尋找一個當(dāng)次迭代阻尼參數(shù)λ_k，使||s(λ_k)||₂∈[0.75δ，1.5δ]，設(shè)s＝s(λ_k)；

3-5)設(shè)調(diào)整后的參數(shù)向量x₊＝x₀+s，求E(x₊)，計算實際變化量ΔE＝E(x₊)-E(x₀)；

3-6)若ΔE≥0，則取δ∈[0.1δ，0.5δ]，并返回步驟3-4)；

3-7)計算m(x₊)，設(shè)預(yù)計變化量ΔE_pred＝m(x₊)-E(x₀)，實際變化量與預(yù)計變化量的比值R＝ΔE/ΔE_pred若R≤0.1，則δ＝0.5δ，若R≥0.75，則δ＝2δ，否則δ不變；

3-8)若|s|/(0.001+|x₀|)≥0.1，則設(shè)x₀＝x₊，并返回步驟3-2)，否則最優(yōu)化的參數(shù)x_f＝x₊，E_f＝E(x₊)，停止計算，得到最優(yōu)化的參數(shù)x_f，并將該最優(yōu)化的參數(shù)x_f分別帶入式(1)和式(9)，從而獲得海面回波波形與海底回波波形，波形分解的效果如圖3所示。

步驟四、求解海水深度D：通過步驟三得到的海面回波波峰的對應(yīng)時間t_G和海底回波波峰的對應(yīng)時間t_max求出海水深度其中c為光速，n_w為海水折射率。

盡管上面結(jié)合附圖對本發(fā)明進行了描述，但是本發(fā)明并不局限于上述的具體實施方式，上述的具體實施方式僅僅是示意性的，而不是限制性的，本領(lǐng)域的普通技術(shù)人員在本發(fā)明的啟示下，在不脫離本發(fā)明宗旨的情況下，還可以做出很多變形，這些均屬于本發(fā)明的保護之內(nèi)。