基于單波束聲吶的沉積層參數序貫估計

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

(1.中國科學院 水聲環境特性重點實驗室，北京 100190； 2.中國科學院聲學研究所，北京 100190； 3.中國科學院大學，北京 100049)

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

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

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

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

1 高頻海底反向散射包絡模型

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

(1)

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

(2)

式中：ρw和vw分別為水層的密度和聲速；Mz為海底底質平均粒徑，其他參數根據經驗公式[14]計算。

2 沉積層參數估計

2.1 全局優化方法

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

利用包絡反演海底聲學和散射特性參數通常采用最小二乘均方誤差函數作為優化問題中的代價函數[5]：

(3)

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

2.2 無跡卡爾曼濾波

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

基于隱馬爾可夫模型[15]，卡爾曼濾波首先建立狀態空間方程。假設非平穩隨機過程中Xk為系統狀態，狀態空間方程表示為：

Xk=fk(Xk-1,Vk)

(4)

Yk=hk(Xk-1,Wk)

(5)

式中：Xk和Yk的方程表達式又分別稱為系統的狀態方程和觀測方程；Vk和Wk分別為模型中狀態噪聲與觀測噪聲，兩者之間互不相關。利用無跡變換選取2n+1個采樣點作為Sigma點集：

(6)

采樣點的權值ω(i)為：

(7)

UKF在此求出一組Sigma點的預測，并對他們加權求均值，得到系統狀態和協方差的預測值：

(8)

(9)

(10)

根據預測值，再次利用UT變換產生新的Sigma點集，并代入觀測方程得到預測的觀測量：

(11)

(12)

接著，通過加權求和得到系統預測的均值及協方差，并計算卡爾曼增益和系統狀態的協方差更新：

(13)

(14)

(15)

(16)

(17)

最后對系統狀態量的預測值進行校正：

(18)

3 海上實驗

3.1 實驗概況

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

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

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

圖2 換能器姿態監視Fig.2 Pitch and roll angle of SBS along the track

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

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

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

3.2 數據預處理

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

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

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

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

圖7 整個測線海底散射回波包絡Fig.7 The experimental backscatter envelope on the track

3.3 序貫估計和全局優化反演

圖8給出了海底底質平均粒徑值Mz分別取3.56、4.70、7.40時高頻海底反向散射包絡模型預報的散射包絡。當Mz增大時，單波束聲吶接收到的聲強度增大，散射包絡的面積增大。在該模型的基礎上可利用預處理后的80組散射信號包絡分別進行海底聲學參數全局反演和序貫估計。

圖8 不同海底底質平均粒徑值的預報包絡Fig.8 Prediction envelope of different mean grain sizes

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

圖9(a)和圖9(b)分別給出了SA和UKF預測包絡與實測包絡的最小二乘均方誤差隨距離變化的針狀圖。由于SA在誤差收斂到一定閾值時終止搜索，其預測精度較為穩定。UKF根據預測與實測數據之間的誤差反饋，自適應地調整待估參數的激勵，由圖可見在拐點附近UKF可及時調整狀態參數的更新量以減小預測誤差，但在海底底質劇變區測量誤差增大，這是因為狀態轉移方程假定動態模型為緩變過程。SA和UKF的最小二乘均方誤差均控制在0.05以下，表明預測結果可信度尚好。

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

利用序貫估計的海底聲學參數預測海底散射包絡隨距離如圖10所示。通過對比圖7可知預測的海底散射包絡與實測結果對應較好。

圖10 利用序貫估計的海底聲學參數預測海底散射包絡Fig.10 The seafloor backscattering envelope predicted by sequential estimation

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

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

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

圖12 海底底質類型分布Fig.12 Core sampler type distribution

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

表1 SA與UKF的計算時間代價對比Table 1 Calculating time cost comparison

4 結論

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

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

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