基于單波束聲吶的沉積層參數(shù)序貫估計(jì)

肖旭， 任群言， 蘇林， 馬力， 鹿力成

(1.中國科學(xué)院 水聲環(huán)境特性重點(diǎn)實(shí)驗(yàn)室，北京 100190； 2.中國科學(xué)院聲學(xué)研究所，北京 100190； 3.中國科學(xué)院大學(xué)，北京 100049)

淺海水聲環(huán)境往往呈現(xiàn)水平非均勻性，其中海底底質(zhì)聲學(xué)特性的水平變化是其主要表現(xiàn)之一。由于直接采樣測量海底底質(zhì)聲學(xué)特性的方法過程繁瑣、成本高昂，許多研究聚焦于利用聲學(xué)信號反演沉積層聲學(xué)參數(shù)的間接方法。一種典型間接方法是基于單波束測深儀的沉積層地聲參數(shù)反演方法[1]。它基于海底高頻反向散射的復(fù)合粗糙度模型[2]，給出小入射角高頻聲波海底反向散射聲信號包絡(luò)時(shí)域建模方法[3]，并通過海底反向散射信號包絡(luò)進(jìn)行地聲參數(shù)反演[4]。

地聲參數(shù)反演問題通?？闯梢粋€(gè)非線性優(yōu)化問題，將各種全局優(yōu)化算法(模擬退火、遺傳算法等)應(yīng)用于求代價(jià)函數(shù)最優(yōu)解問題[5]。此類方法在尋優(yōu)過程中需反復(fù)調(diào)用正演模型計(jì)算，普遍存在計(jì)算開銷大、收斂速度慢，以及易陷入局部最優(yōu)、超參數(shù)調(diào)優(yōu)困難等問題?？紤]到在實(shí)際的海洋測試中，海底底質(zhì)特性在局部區(qū)域隨距離變化的過程較為緩慢，這為地聲參數(shù)的序貫估計(jì)提供了可行性。

卡爾曼濾波是一種典型序貫估計(jì)算法[7]，可對線性或非線性系統(tǒng)中狀態(tài)變量和參數(shù)進(jìn)行動(dòng)態(tài)估計(jì)[8]，利用卡爾曼濾波框架進(jìn)行序貫估計(jì)的方法已被有效應(yīng)用于水聲目標(biāo)測距[9]、聲速剖面反演[10]等領(lǐng)域?？柭鼮V波過程分為預(yù)測和校正2個(gè)步驟，可根據(jù)正演模型預(yù)測與實(shí)測數(shù)據(jù)之間的誤差反饋，自適應(yīng)地更新對狀態(tài)變量的激勵(lì)，從而校正模型參數(shù)的預(yù)測值。

本文提出了一種利用卡爾曼濾波對海底反向散射回波進(jìn)行沉積層參數(shù)序貫估計(jì)的方法，將模型應(yīng)用于2018年8月采集的海試數(shù)據(jù)，獲得測線上海底底質(zhì)類型平均粒徑分布，并與全局優(yōu)化算法中性能較優(yōu)的模擬退火算法進(jìn)行了計(jì)算效率和反演精度的對比。

1 高頻海底反向散射包絡(luò)模型

I(t)=Ii(t)+Iv(t)

(1)

式中：Ii為海底粗糙界面散射包絡(luò)；Iv為海底介質(zhì)體積散射包絡(luò)，假設(shè)兩者相互獨(dú)立。Ii是一個(gè)亥姆霍茲衍射積分方程的Kirchhoff近似解，考慮了單波束聲吶系統(tǒng)的安放位置、指向性、信號特征、水體擴(kuò)展和吸收等，可以很好地描述沉積層界面反向散射的物理過程，特別是在小入射角條件下。體積散射信號包絡(luò)Iv考慮了界面粗糙度對體積散射的影響。根據(jù)Sternlicht等[3]的模型，散射回波包絡(luò)聲壓最終表示為：

(2)

式中：ρw和vw分別為水層的密度和聲速；Mz為海底底質(zhì)平均粒徑，其他參數(shù)根據(jù)經(jīng)驗(yàn)公式[14]計(jì)算。

2 沉積層參數(shù)估計(jì)

2.1 全局優(yōu)化方法

地聲參數(shù)反演問題通?？醋魇欠蔷€性數(shù)值優(yōu)化問題，以待反演參數(shù)為決策變量，確立代價(jià)函數(shù)和優(yōu)化的可行域，即決策變量的搜索區(qū)間，并給定其須滿足的約束條件，通過全局優(yōu)化方法進(jìn)行迭代計(jì)算，最終得到模型參數(shù)的全局最優(yōu)解。

利用包絡(luò)反演海底聲學(xué)和散射特性參數(shù)通常采用最小二乘均方誤差函數(shù)作為優(yōu)化問題中的代價(jià)函數(shù)[5]：

(3)

常用的非線性全局優(yōu)化算法包括模擬退火算法[11]、遺傳算法[12]、粒子群算法[13]等。利用此類算法對海底底質(zhì)平均粒徑Mz進(jìn)行搜索，通過多次迭代使式(3)中的代價(jià)函數(shù)最小化。在與序貫估計(jì)進(jìn)行對比分析時(shí)，本文采用模擬退火作為反演模型的優(yōu)化算法。

2.2 無跡卡爾曼濾波

無跡卡爾曼濾波(unscented Kalman filter，UKF)是通過無跡變換(unscented transform，UT)產(chǎn)生一系列確定樣本來逼近狀態(tài)的后驗(yàn)概率密度，既沒有線性化忽略高階項(xiàng)的過程，也不用求解雅可比矩陣[7]。綜合考慮濾波性能和運(yùn)算效率，UKF是3種方法中的首選。

基于隱馬爾可夫模型[15]，卡爾曼濾波首先建立狀態(tài)空間方程。假設(shè)非平穩(wěn)隨機(jī)過程中Xk為系統(tǒng)狀態(tài)，狀態(tài)空間方程表示為：

Xk=fk(Xk-1,Vk)

(4)

Yk=hk(Xk-1,Wk)

(5)

式中：Xk和Yk的方程表達(dá)式又分別稱為系統(tǒng)的狀態(tài)方程和觀測方程；Vk和Wk分別為模型中狀態(tài)噪聲與觀測噪聲，兩者之間互不相關(guān)。利用無跡變換選取2n+1個(gè)采樣點(diǎn)作為Sigma點(diǎn)集：

(6)

采樣點(diǎn)的權(quán)值ω(i)為：

(7)

UKF在此求出一組Sigma點(diǎn)的預(yù)測，并對他們加權(quán)求均值，得到系統(tǒng)狀態(tài)和協(xié)方差的預(yù)測值：

(8)

(9)

(10)

根據(jù)預(yù)測值，再次利用UT變換產(chǎn)生新的Sigma點(diǎn)集，并代入觀測方程得到預(yù)測的觀測量：

(11)

(12)

接著，通過加權(quán)求和得到系統(tǒng)預(yù)測的均值及協(xié)方差，并計(jì)算卡爾曼增益和系統(tǒng)狀態(tài)的協(xié)方差更新：

(13)

(14)

(15)

(16)

(17)

最后對系統(tǒng)狀態(tài)量的預(yù)測值進(jìn)行校正：

(18)

3 海上實(shí)驗(yàn)

3.1 實(shí)驗(yàn)概況

本試驗(yàn)于2018年8月在黃海進(jìn)行，通過走航式測量海底散射回波，主要目的是對沉積物聲學(xué)特征遙測的可行性進(jìn)行研究。走航試驗(yàn)開始于試驗(yàn)海區(qū)E1站位點(diǎn)(36°00.180 4′N 121°36.804 0′E)，利用收發(fā)合置換能器發(fā)射 25 kHz單頻信號和線性調(diào)頻信號，于次日到達(dá)E2站位點(diǎn)，全程約為130 km。

本次試驗(yàn)采用的換能器聲源級和接收靈敏度較高，25 kHz聲源級為212.72 dB，接收靈敏度為-186.1 dB，收發(fā)合置換能器具有較好的指向性，如圖1所示。圖2給出了收發(fā)合置換能器上羅經(jīng)監(jiān)測的換能器姿態(tài)，在實(shí)驗(yàn)的后半部分隨海況惡化其姿態(tài)有所偏離，在反向散射包絡(luò)的建模中引入了對指向性的糾正。

圖1 換能器指向性Fig.1 Directionality pattern of SBS

圖2 換能器姿態(tài)監(jiān)視Fig.2 Pitch and roll angle of SBS along the track

通過對發(fā)射的線性調(diào)頻信號進(jìn)行處理得到換能器距離海底深度，如圖3所示。水深在80 km處有相對較小的波動(dòng)，其他較為平穩(wěn)。圖4給出了E1和E2站位測量的水文聲速剖面，在15～20 m海深處存在明顯的躍變層，在此聲速剖面基礎(chǔ)上計(jì)算信號的傳播時(shí)間。

圖3 測線海水深度Fig.3 Water depth along the track

圖4 聲速剖面曲線Fig.4 Speed of sound profile

3.2 數(shù)據(jù)預(yù)處理

實(shí)驗(yàn)中發(fā)射了頻率為25 kHz，時(shí)間脈沖長度為2 ms的單頻信號。信號采樣頻率為500 kHz，同時(shí)記錄了功率放大器監(jiān)視信號和收發(fā)合置換能器信號。散射信號包絡(luò)的獲取是通過對濾波后的原始信號進(jìn)行Hilbert變換并取絕對值獲得。圖5給出了第40組100個(gè)信號的海底散射回波包絡(luò)圖，可以看出由于海底的深度變化和換能器的上下浮動(dòng)等因素，接收信號的包絡(luò)到達(dá)時(shí)刻變化很大。由于收發(fā)合置換能器具有較高的源級，接收信號信噪比較高，故對信號包絡(luò)可采用最小門限法[5]進(jìn)行對齊。圖6給出了信號的第40組包絡(luò)對齊結(jié)果。在走航過程中海底的散射不同獲得的信號包絡(luò)差別較大，對獲取的散射回波包絡(luò)每100個(gè)樣本進(jìn)行平均，共得到80組信號包絡(luò)。

圖5 第40組信號的海底散射回波包絡(luò)Fig.5 Original backscatter envelopes of color plot for the 40th dataset

圖6 第40組海底散射回波包絡(luò)對齊結(jié)果Fig.6 Aligned backscatter envelopes of color plot for the 40th dataset

基于最小門限法對齊的信號，對獲取的信號回波包絡(luò)進(jìn)行80%的重疊平均處理，獲得整個(gè)E1—E2測線海底散射回波包絡(luò)，如圖7所示。利用獲得的信號的80組包絡(luò)分別進(jìn)行海底聲學(xué)參數(shù)全局優(yōu)化反演和序貫估計(jì)。

圖7 整個(gè)測線海底散射回波包絡(luò)Fig.7 The experimental backscatter envelope on the track

3.3 序貫估計(jì)和全局優(yōu)化反演

圖8給出了海底底質(zhì)平均粒徑值Mz分別取3.56、4.70、7.40時(shí)高頻海底反向散射包絡(luò)模型預(yù)報(bào)的散射包絡(luò)。當(dāng)Mz增大時(shí)，單波束聲吶接收到的聲強(qiáng)度增大，散射包絡(luò)的面積增大。在該模型的基礎(chǔ)上可利用預(yù)處理后的80組散射信號包絡(luò)分別進(jìn)行海底聲學(xué)參數(shù)全局反演和序貫估計(jì)。

圖8 不同海底底質(zhì)平均粒徑值的預(yù)報(bào)包絡(luò)Fig.8 Prediction envelope of different mean grain sizes

首先，利用全局優(yōu)化法分別對80組散射包絡(luò)進(jìn)行地聲參數(shù)反演，所用算法為SA，搜索區(qū)間為 [0,10]，初始溫度T0=100，溫度變化率ΔT=0.98，設(shè)多個(gè)終止條件，當(dāng)10次迭代內(nèi)目標(biāo)函數(shù)平均變化率δ<10-3或代價(jià)函數(shù)值小于0.02時(shí)終止迭代，且最大迭代次數(shù)Nmax=50。其次，利用序貫估計(jì)法對80組散射包絡(luò)依序進(jìn)行地聲參數(shù)估計(jì)，所用算法為UKF，為獲取相對準(zhǔn)確的初始狀態(tài)值，將第1組包絡(luò)的全局優(yōu)化反演結(jié)果作為UKF初始狀態(tài)，UKF初始狀態(tài)協(xié)方差設(shè)為單位矩陣I。序列長度設(shè)為80，對應(yīng)信號的80組散射包絡(luò)。由于地聲參數(shù)在短距離內(nèi)的狀態(tài)變化可視為緩變過程，令式(4)中狀態(tài)方程為Xk=Xk-1+Vk。其中Vk表示k時(shí)刻過程中的隨機(jī)擾動(dòng)，服從均值為零、協(xié)方差矩陣為Q=qI的高斯分布；式(5)測量方程設(shè)為Yk=hk(Xk-1)+Wk，hk為式(3)給出的函數(shù)，Wk表示測量噪聲，服從均值為零、方差為r的高斯分布。根據(jù)狀態(tài)轉(zhuǎn)移過程和測量特性將q設(shè)為0.01，r試設(shè)為1和2。

圖9(a)和圖9(b)分別給出了SA和UKF預(yù)測包絡(luò)與實(shí)測包絡(luò)的最小二乘均方誤差隨距離變化的針狀圖。由于SA在誤差收斂到一定閾值時(shí)終止搜索，其預(yù)測精度較為穩(wěn)定。UKF根據(jù)預(yù)測與實(shí)測數(shù)據(jù)之間的誤差反饋，自適應(yīng)地調(diào)整待估參數(shù)的激勵(lì)，由圖可見在拐點(diǎn)附近UKF可及時(shí)調(diào)整狀態(tài)參數(shù)的更新量以減小預(yù)測誤差，但在海底底質(zhì)劇變區(qū)測量誤差增大，這是因?yàn)闋顟B(tài)轉(zhuǎn)移方程假定動(dòng)態(tài)模型為緩變過程。SA和UKF的最小二乘均方誤差均控制在0.05以下，表明預(yù)測結(jié)果可信度尚好。

圖9 2種方法預(yù)測的最小二乘均方誤差隨距離變化針狀圖Fig.9 Stem figure of least square mean square error at different distance predicted by two methods

利用序貫估計(jì)的海底聲學(xué)參數(shù)預(yù)測海底散射包絡(luò)隨距離如圖10所示。通過對比圖7可知預(yù)測的海底散射包絡(luò)與實(shí)測結(jié)果對應(yīng)較好。

圖10 利用序貫估計(jì)的海底聲學(xué)參數(shù)預(yù)測海底散射包絡(luò)Fig.10 The seafloor backscattering envelope predicted by sequential estimation

圖11給出了全局優(yōu)化法和序貫估計(jì)法對海底底質(zhì)平均粒徑Mz的預(yù)測與海底底質(zhì)分類結(jié)果。從圖中可以看出UKF與SA的底質(zhì)分類結(jié)果基本一致，從大約60 km處海底底質(zhì)平均粒徑變小，表明海底底質(zhì)由軟變硬。測量噪聲方差r=2時(shí)相比r=1時(shí)測量結(jié)果曲線更加平滑，說明r的增大可使測量噪聲對模型的影響增大，使卡爾曼增益Kk減小，從而減弱系統(tǒng)的校正權(quán)值，降低了系統(tǒng)的瞬態(tài)響應(yīng)和穩(wěn)態(tài)值。根據(jù)Wentworth分類標(biāo)準(zhǔn)，獲得海底底質(zhì)平均粒徑為極細(xì)粉砂到中粉砂范圍。

圖12給出了試驗(yàn)海區(qū)位置的海底底質(zhì)類型分布，E1—E2測線為試驗(yàn)走航測線，在試驗(yàn)起始端底質(zhì)類型為黏土質(zhì)粉砂，與反演得到的平均粒徑一致，如圖11中左側(cè)方框所示；并很快過渡到砂質(zhì)粉砂，后半段測線為粉砂質(zhì)砂，與圖11中的右側(cè)方框一致。序貫估計(jì)結(jié)果與試驗(yàn)預(yù)報(bào)結(jié)果符合很好。

圖11 海底底質(zhì)平均粒徑Fig.11 Estimation results of mean grain size

圖12 海底底質(zhì)類型分布Fig.12 Core sampler type distribution

利用UKF進(jìn)行一步估計(jì)時(shí)，由于UKF無需同EKF一樣計(jì)算非線性函數(shù)的雅可比矩陣，其計(jì)算復(fù)雜度僅來源于UT變換中產(chǎn)生Sigma點(diǎn)集時(shí)和預(yù)測校正時(shí)所調(diào)用的正演函數(shù)計(jì)算。因此，序貫估計(jì)法所需的時(shí)間代價(jià)同全局優(yōu)化法一樣由正演程序的計(jì)算次數(shù)決定。對于序貫估計(jì)法，每組包絡(luò)固定執(zhí)行2L+1次正演模型計(jì)算，L為待估參數(shù)個(gè)數(shù)；對于全局優(yōu)化法，每組包絡(luò)需執(zhí)行的正演次數(shù)由尋優(yōu)過程的迭代次數(shù)N決定，且根據(jù)預(yù)定的終止條件有所不同，每組反演所需的迭代次數(shù)根據(jù)具體收斂情況、偽峰和鞍點(diǎn)個(gè)數(shù)而不定，可達(dá)數(shù)十次到上百次不等。在Intel Core i5-9400F 2.90Ghz處理器的PC機(jī)環(huán)境下應(yīng)用Matlab計(jì)算，二者所用計(jì)算時(shí)間對比如表1所示。UKF在單次估計(jì)中平均計(jì)算時(shí)間相比SA法縮短83.4%，總體時(shí)間開銷縮短79.6%。

表1 SA與UKF的計(jì)算時(shí)間代價(jià)對比Table 1 Calculating time cost comparison

4 結(jié)論

1)序貫估計(jì)預(yù)測的海底散射回波包絡(luò)與實(shí)測包絡(luò)符合較好，全局優(yōu)化反演和序貫估計(jì)法預(yù)測的最小二乘均方誤差均可控制在0.05以下，預(yù)測精度尚好。

2)海底底質(zhì)平均粒徑的預(yù)測值與海底底質(zhì)分類分布對比具有較高的可信度。序貫估計(jì)的結(jié)果顯示，海底底質(zhì)類型在試驗(yàn)起始端為黏土質(zhì)粉砂，并很快過渡到砂質(zhì)粉砂，后半段測線為粉砂質(zhì)砂，與試驗(yàn)預(yù)報(bào)結(jié)果符合很好。

3)序貫估計(jì)法在保證和全局優(yōu)化法精度相似的情況下有效提升算法效率，參數(shù)估計(jì)消耗的平均時(shí)間顯著縮短，具有良好的研究和應(yīng)用價(jià)值。