距離依賴的聲速場反演與運動聲源的跟蹤定位

戴 淼， 李亞安

(西北工業大學 航海學院，西安 710072)

水下聲傳播環境是一種地區性、參數值多變的復雜環境，現場儀器測量難以大范圍的精確獲取所需要的海洋環境參數。聲學反演技術將聲傳播模型與觀測數據相結合，能夠快速簡便的獲取所需要的參數。聲速剖面 (Sound Speed Profile, SSP) 是影響海洋聲傳播的重要因素之一，其結構特征決定了水下聲傳播路徑和傳播損失[1]。經驗正交函數 (Empirical Orthogonal Function, EOF) 是表示十分有效的基函數，通常只需要前幾階便可以精確的反演出SSP，大大減少了待反演參數[2]。已有研究成果表明：在水平變化環境下的平均SSP反演[3]和SSP特征量統計分布規律的獲取[4]與同步觀測數據相比，具有很好的一致性。在這里反演得到的是平均意義上的估計，SSP所具有的時空變化特性，尤其是聲速快變區域的存在將影響反演精度。

序貫濾波方法基于系統的觀測值對系統狀態進行估計和更新，能夠根據相關的物理統計模型預測未知海洋聲學參數[5]。卡爾曼濾波器用于處理跟蹤問題[6]，要求狀態方程和觀測方程均為線性。擴展卡爾曼濾波(Extended Kalman Filter, EKF)把非線性問題線性化，在海洋聲學參數反演方面取得了一定效果[7]。Yardim等[8]通過EKF、無跡卡爾曼濾波(Unscented Kalman Filter, UKF)和粒子濾波(Particle Filter, PF)將地聲反演轉換為跟蹤問題，對具有時空變化特性的環境參數進行估計。Carrière等[9]將聲學測量數據同化問題建模為狀態-空間模型，把沿海環境距離依賴聲速場參數作為狀態變量，利用EOF進行參數低維化，并首次應用集合卡爾曼濾波(Ensemble Kalman Filter, EnKF)估計了SSP隨時間的演化[10]。Li等[11]對比多種序貫濾波算法分析了距離獨立環境下地聲參數反演跟蹤情況。需要指出的是上述研究過程均采用固定收發裝置，即聲源與接收陣之間是相對靜止的。

傳統的固定收發裝置得到的反演結果為聲源和接收陣之間的平均SSP，空間分辨率存在局限性；在實際情況中，受季節、溫度以及動態海洋(存在渦旋、內波、冷熱水團等)的影響，收發裝置之間存在的SSP會隨著時間/距離的變化而變化，尤其是聲速快變區域的存在，這就要求更多的待反演參數，故而反演算法的速度和精度都將受到制約。針對上述問題，由于SSP的時空演化特性可在一定條件下近似的建模為EOFs系數隨時間/距離的演化，在已有序貫濾波算法基礎上，本文提出借助運動聲源的時變性，實現距離依賴快變聲速場反演與運動聲源位置的聯合估計。通過引入運動聲源在不同空間位置處發射信號，利用垂直陣進行接收，分別采用EKF和EnKF同時處理多個待反演參數。該方法相比固定收發布設裝置增加了聲場的空間信息，同時聲場中包含的信息又可以逆推聲源的位置，因而能夠在重構真實聲速場的同時對運動聲源進行跟蹤定位，并且根據估計結果的均方根誤差對算法的估計性能進行分析、比較。研究表明：EnKF表現出較優越的估計性能，即使在陣元數目減少的局限條件下，該算法仍然能夠很好的跟蹤估計出運動聲源參數和聲速場在時間/距離維度上的演化軌跡。

1 理論與算法

1.1 距離依賴條件下的聲傳播模型

海洋環境復雜多變導致無法對其精確建模，一般采用量化水下各種物理過程對聲傳播影響的方法，運用合理的聲傳播模型最大程度的重構水下物理場。簡正波模型通過求解與深度有關的特征方程，得到一組振動模式，振動的頻率給出與模式傳播相關聯的水平波數，對加權的每個模式的貢獻進行疊加來構成總聲場。簡正波聲傳播模型本質是利用簡正波模式求解亥姆霍茲方程[12],頻率為ω的空間任意點源s(r,z)在柱坐標下聲壓可表示為

(1)

式中：c(r,z)表示聲速，在距離獨立條下c(r,z)=c(z)，模態方程為

(2)

根據邊界條件，聲壓正比于各階模態的和與漢克爾函數的乘積

(3)

式中：φm(z)和km分別表示為深度方向上的模態和水平波數。

將距離獨立環境推廣至距離依賴簡正波聲傳播理論：水平距離分成段，各段均近似視為距離獨立環境，并對每一段以標準簡正波求解聲場，最后根據邊界條件將各段的解連接在一起[10]：

(4)

(5)

(6)

聲壓在第l個界面上的連續性條件為

(7)

耦合簡正波的數值求解過程非常復雜[13]，在實際計算中，通常采用絕熱簡正波理論[14]忽略模態之間的能量耦合項，即各階模態從當前距離至下一距離是絕熱耦合的。

1.2 非線性濾波算法

濾波算法將不同性質的測量數據有效結合起來，并且能夠給出相關尺度范圍內與觀測數據吻合的動力過程的確定性描述，目的在于更好的估計所需要的參數，改進動力模型的性能。考慮到距離依賴SSP的動態演化特性，采用序貫濾波的思想對其進行預測估計，下面給出兩種非線性濾波算法。

1.2.1 擴展卡爾曼濾波

擴展卡爾曼濾波 (Extended Kalman Filter, EKF)方法實質是將非線性模型在其狀態變量均值的范圍內進行泰勒級數展開，并在當前估計值處線性化，即對過程方程或測量方程求偏導，從而將非線性問題線性化。EKF算法可表示為以下遞推形式：

首先給出兩個非線性預設方程：

xk=f(xk-1,uk,wk)

(8)

yk=h(xk,vk)

(9)

函數f用k-1時刻的狀態和k時刻的輸入uk預測k時刻的狀態，相應的，函數h以k時刻的狀態預測k時刻的測量值，wk和vk為噪聲向量。實際中并不知道每一時刻狀態噪聲和觀測噪聲的值，可以將其假設為零。

k時刻的預測值：

xk|k-1=f(xk-1|k-1,uk, 9)

(10)

yk|k-1=h(xk|k-1, 0)

(11)

狀態協方差矩陣的預測值

(12)

k時刻的更新：

新息

rk=yk-yk|k-1

(13)

新息的協方差矩陣

(14)

卡爾曼增益

(15)

狀態協方差矩陣

Pk|k=(I-GkFk)Pk|k-1

(16)

狀態向量

xk|k=xk|k-1+Gkrk

(17)

上述函數用到了狀態轉換函數和觀測函數的Jacobian:

(18)

(19)

式中：下標k|k-1表示由時刻的值預測得到k時刻的值，下標k|k表示使用k時刻新息后的更新值，下標k則表示k時刻的值。

1.2.2 集合卡爾曼濾波

集合卡爾曼濾波(Ensemble Kalman Filter, EnKF)基本思想:初始化產生一組系統狀態的采樣作為背景集合，根據集合的預測結果估計狀態向量和觀測向量的協方差，然后利用觀測信息和協方差信息計算分析集合，通過系統模型傳遞分析集合并得到下一時刻的背景數據集合，遞推向前進行預測。EnKF算法可表示為：

(20)

(21)

(22)

卡爾曼增益

(23)

式中：Rk為k時刻的觀測誤差協方差矩陣。

根據k時刻的觀測信息，更新背景集合，得到k時刻的分析集合：

(24)

(25)

對比兩種濾波算法可以看出：EnKF使用非線性集合估計真實統計值，省略了對模型算子和觀測算子的線性化過程，在提高估計精度的同時降低了計算量。

1.3 聲速剖面的經驗正交函數表示

SSP反演作為一種多維優化問題，反演算法的復雜度和精度通常取決于待反演參數的個數。EOF用于表征海水的SSP，大大減少了描述聲速垂直結構所需要的參數，其原理如下：

EOF將SSP隨深度變化的關系參數化。通常可以由觀測海域的實測SSP計算得到，不同時刻測得的SSP可表示為：

(26)

(27)

式中：fn為矩陣R的特征向量，即經驗正交函數；λn為特征向量對應的特征值。協方差矩陣R表征了該海域聲速起伏的不確定量。將上述的N個特征值λn按從大到小的順序排列，選取前K個特征值對應的特征向量來表示SSP，則海域內任意一點的聲速值可表示為：

(28)

式中：αk(x,y)為第k階EOF系數。

由式(28)可知，EOF是SSP在深度維上的離散化。根據距離依賴的聲傳播理論，現將SSP在距離維進行離散化處理，即把聲源和接收陣之間的水平距離分成段，相應的聲速可表示為

(29)

ul(r)為門函數，定義為：

(30)

在第l個區域，有

(31)

式中：αk, l為相應距離處的EOFs系數。

在最小均方意義下，EOF用于表征SSP在保證反演精度的前提下，可降低算法復雜度及計算時間，從而提高算法性能。

1.4 運動聲源的狀態方程

(32)

式中：s為狀態變量，Δt為狀態更新的時間增量 (即兩次測量的時間間隔) ，vz和va分別表示速度在深度方向的擾動和徑向的加速度擾動，H為狀態轉移矩陣。

狀態方程的第二部分用于描述SSP的狀態，將EOFs系數作為狀態參量，SSP的狀態方程可由前n階EOFs系數的狀態方程代替：

(33)

式中：wk為EOFs系數在空間上演化的過程噪聲。

在這里，將觀測方程用于表征垂直線列陣 (VLA) 接收到的聲壓信號，根據簡正波理論，合成聲壓可以表示為水平距離R，深度Z，頻率ω，聲速c以及海底邊界條件參數(BCs)的非線性方程：

yk=h(Z,R,ω,c1, …，ck,BCs)+vk

(34)

式中：h為簡正波聲傳播模型，ck為第k個距離上的聲速剖面，vk為測量噪聲。

運動聲源的引入可增加水下聲場的空間信息，同時聲場中所包含的信息又能夠逆推聲源的位置，通過觀測方程將狀態向量 (待反演參數) 和觀測向量 (VLA接收到的聲壓) 聯系起來。

2 結果與分析

2.1 實測聲速剖面的EOF反演

為證明EOFs系數表示SSP的有效性，將給出EOF反演SSP與實測聲速剖面的對比。中美聯合水聲實驗組于2001年6月在東中國海進行寬帶信號傳播實驗(ASIAEX,ECS,2001)[15]，圖1為現場CTD實測的聲速樣本和平均聲速剖面，其中細線為實測的54組聲速剖面，粗線為平均聲速剖面。可見，實驗期間溫躍層 (30～60 m) 范圍內聲速隨時間的起伏較大，變化最大峰值達到了11 m/s以上。

圖1 東海不同時刻聲速剖面樣本Fig.1 Sound speed profiles in different time of the East China Sea

圖2給出對該組聲速剖面所提取到的EOFs：其中圖2(a)為圖(1)所示的聲速起伏，并給出如圖2(b)所示的前六階EOFs值，由圖2(c)可見，前三階EOFs可表征近90%的聲速起伏。利用圖2提取的6階EOFs對54條聲速剖面中任意時刻CTD所測得的其中一條聲速剖面進行擬合，并將反演得到的SSP與實測SSP對比于圖3中。

驗證結果表明：只需利用前幾階EOFs即較少EOFs系數便能夠較為精確的表征該海域海水SSP的演變過程，此方法相比于直接將聲速剖面表征為深度的函數大大減少了待反演參數的個數。在實際應用過程中，通常希望在保證一定精度的前提下用盡量少的參數來表示復雜的聲場環境，EOF階數的選取還應考慮計算量的大小。一般的，只需選取3～6階EOFs便可較為準確的表征聲速場的起伏情況。

(a) 聲速剖面樣本差異

(b) 前六階EOFs

(c) 前10階EOFs所占聲速剖面起伏的能量比圖2 提取的海水聲速剖面經驗正交函數Fig.2 EOFs of SSPs in the water volume

圖3 某一時刻CTD實測聲速剖面(實線)與EOF反演聲速剖面(虛線)Fig.3 CTD measured SSP (solid line) and EOF invert SSP (dotted line)

2.2 距離依賴聲速場重構和運動聲源的跟蹤定位

由于沒有獲取到同步的聲壓測量數據，因此本文中聲壓根據環境模型擬合得到，波導環境參考Benchmarking距離依賴淺海波導模型，海底參數根據東中國海(ASIAEX,ECS,2001)實測數據設置。假設仿真波導環境如圖4所示：在平整海面條件下，海面邊界為簡單的壓力釋放邊界。海深106 m,16元VLA等間距間隔4 m布設于距離海面15 m處，陣長60 m。運動聲源初始深度為40 m，從水平方向初始位置1 km處以2.5 m/s的徑向速度遠離接收陣，并每隔20s發射一次信號，聲源輻射頻率為400 Hz，20 min內共移動60次。

圖4 距離依賴的淺海環境參數模型Fig.4 Experimental configuration of the range dependent shallow water environment

為了減少待反演樣本數目，擬采用前三階EOFs系數作為聲速參數來實現聲速場的重構。假定過程噪聲和測量噪聲服從零均值的高斯分布，聲場模型為絕熱簡正波模型[17]，聲壓作為觀測值由聲學軟件KRAKENC根據環境參數合成。分別用EKF和EnKF濾波算法對運動聲源參數(深度、距離和速度)以及EOFs前三階系數同時進行跟蹤，并將跟蹤估計結果與真實值進行比較。其中，EnKF樣本大小為50。

圖5為EKF和EnKF的參數聯合估計結果，其中，(a)～(c) 分別為運動聲源深度、距離和速度隨時間的演化，(d)～(f) 分別為EOFs前三階系數隨距離的演化。實線代表運動聲源參數和前三階EOFs系數隨時間/距離的真實演化過程，點線和虛線分別是EKF和EnKF估計結果。可以看出，EKF除了運動聲源的速度以外，其余參數均能夠跟蹤上，產生這一現象的主要原因是由于擬合的聲壓中并不包含聲源的速度信息。EnKF對于所有反演參數均能夠跟蹤上，且估計結果相對穩定，波動較小，相比于EKF更能逼近待反演參數真實狀態隨時間/距離的演化軌跡。

圖6給出上述兩種濾波算法估計的運動聲源參數在時間維的均方根誤差，可以看出，EnKF在聲源的深度、距離和速度隨時間演化過程中均表現出較EKF優越的估計性能，幾乎與聲源的真實位置 (深度、距離) 相吻合，均方根誤差在實時觀測過程中接近于零。

圖5 EKF和EnKF參數聯合估計結果,圖中實線為參數真實演化軌跡Fig.5 Simultaneous estimation results of EKF (dash-dotted line) and EnKF (dashed line). The true values trajectories are in solid line

圖6 運動聲源參數Fig.6 Root-mean-square-errors of the moving source parameters estimation in time dimension

由EKF和EnKF估計得到的聲速參數重構的聲速場和其相對誤差如圖7所示。由(a)，(c)可知，兩種濾波算法均能夠重構出真實的聲速場，同時可以看到聲速場個別區域存在較大的聲速擾動現象。在濾波估計相對誤差(b)，(d)中顯示：EnKF估計效果整體優于EKF，相對誤差均較小，尤其在聲速變化較為劇烈區域，EnKF表現出更為理想的估計效果。

圖8給出重構的聲速場在距離維度上的均方根誤差，兩種濾波算法均能夠大致跟蹤上EOFs系數的真實演化趨勢。其中，EnKF表現出更優越的估計性能，均方根誤差波動非常小。經計算，EKF和EnKF在距離維度上的平均均方根誤差分別為0.573 m/s和0.011 m/s。

圖7 EKF和EnKF重構的距離依賴聲速場Fig.7 Inversion of range dependent sound speed field

圖8 濾波估計重構的聲速場在距離維的均方根誤差Fig.8. Root-mean-square-errors in range dimension of the invert sound speed fields from different filters estimation

需要指出的是：EnKF的估計性能整體優于EKF，主要是由于測量方程的非線性導致。EKF只適用于弱非線性系統，因其保留了非線性函數的一階近似項，忽略了高階項從而引起更大的誤差；而EnKF不需要對模型算子和觀測算子線性化，具有較高的處理效率和處理高維非線性系統的能力，方差和均值的估計精度至少可達2階以上。

上述仿真均為16陣元的VLA得到的估計結果。但在實際應用中通常不能滿足大量陣元數目的要求，因此需要驗證EnKF濾波在減少陣元數目的局限條件下算法的寬容性。保留原有VLA的有效陣長60 m，將陣元個數減少為8個，陣間距增大至8.57 m, 由圖9、10可知，在相同的時間/距離維度，EnKF仍然能夠保持較為穩健的估計結果。

圖9 不同陣元個數VLA的EnKF估計運動聲源參數Fig.9 Root-mean-square-errors in time dimension of the moving source parameters EnKF estimation

圖10 不同陣元個數VLA的EnKF估計重構聲速場在距離維的均方根誤差Fig.10 Root-mean-square-errors in range dimension of the EnKF estimate sound speed fields from different elements VLA

與16陣元VLA估計結果相比，8陣元VLA在相同時間/距離維度內的EnKF對運動聲源參數/聲速場估計結果的均方根誤差均略有增大，經計算，在相同距離維度上重構的聲速場平均均方根誤差為0.024 m/s。陣元數目的減少降低了空間采樣率，使估計結果的誤差略有增大。但即使在陣元數目減少的局限條件下，EnKF仍然能夠很好的跟蹤估計出運動聲源參數和EOFs系數即距離依賴的快變聲速場的演化軌跡，體現了該算法的魯棒性和寬容性。

3 結 論

受海洋內部復雜動力學的影響，在大范圍海域或者聲速擾動劇烈區域，收發裝置之間存在的聲速場會隨著距離的變化而變化，因而要求更多的待反演參數，反演算法的速度和精度都將受到制約。本文針對距離依賴聲速快變區域，對ASIAEX, ECS, 2001實驗現場CTD測量的聲速樣本數據進行EOFs提取，明確只需前幾階EOFs便可較為準確的表示SSP，驗證了具有時空變化特性的SSP在一定條件下可近似建模為經驗正交函數系數隨時間/距離演化的可行性，同時根據實驗中觀測的水文資料合成聲壓作為觀測值。提出通過將運動聲源參數(深度、距離和速度)和聲速參數 (EOFs前三階系數) 作為測量值放入狀態方程，利用EKF和EnKF濾波算法實現了運動聲源位置與聲速場的聯合估計。通過分析對比，表明本文提出的反演方法在重構真實聲速場的同時可對運動聲源進行跟蹤定位。其中，EnKF表現出更優越的估計性能，并驗證了算法的魯棒性和寬容性。運動聲源的引入，可改善利用傳統的固定收發裝置反演只能夠得到聲源和接收陣之間的平均SSP的局限性，使反演結果更接近于具有時空變化特性SSP的真實演化狀態。同時，運動聲源的靈活性可增加局部區域聲壓的采樣密度，豐富聲學和環境的觀察信息，能夠改善傳統聲層析法固定網絡采樣分辨率和觀測范圍的局限性，提高預報精度。

[ 1 ] CHIU L Y S, LIN Y T, CHEN C F, et al. Focused sound from three-dimensional sound propagation effects over a submarine canyon [J]. J. Acoust. Soc. Am., 2011, 129(6): 260-266.

[ 2 ] LEBLANC L R, MIDDLETON F H. An underwater acoustic sound velocity data model [J]. J. Acoust. Soc. Am., 1980, 67(6): 2055-2062.

[ 3 ] 何利，李整林，彭朝暉，等. 南海北部海水聲速剖面聲學反演[J]. 中國科學：物理學 力學 天文學, 2011, 41(1): 49-57.

HE Li, LI Zhenglin, PENG Chaohui, et al. Inversion for sound speed profiles in the northern of South China Sea [J]. Scientia Sinica Pysica,Mechanica & Astronomica, 2011, 41(1): 49-57.

[ 4 ] 李佳，楊坤德，雷波，等. 印度洋中北部聲速剖面結構的時空變化及其物理機理研究[J]. 物理學報, 2012, 61(8): 282-264.

LI Jia, YANG Kunde, LEI Bo, et al. Research on the temporal-spatial distributions and the physical mechanisms for the sound speed profiles in north-central Indian Ocean [J]. Acta Physica Sinica, 2012, 61(8): 282-264.

[ 5 ] YARDIM C, MICHALOPOULOU Z H, GERSTOFT P. An overview of sequential Bayesian filtering in ocean acoustics [J]. IEEE J. Ocean Eng., 2011, 36(1): 71-89.

[ 6 ] CULVER R L, HODGKISS W S. Comparison of Kalman and least squares filters for locating autonomous very low frequency acoustic sensors [J]. IEEE J. Ocean. Eng.,1988, 13(4) :282-290.

[ 7 ] CANDY J V, SULLIVAN E J. Model-based environment inversion: A shallow water ocean application [J]. J. Acoust. Soc. Am.,1995, 98(3): 1446-1454.

[ 8 ] YARDIM C, GERSTOFT P, HODGKISS W S. Tracking of geoacoustic parameters using Kalman and particle filters [J]. J. Acoust. Soc. Am.,2009, 125(2): 746-760.

[11] LI J L, ZHOU H. Tracking of time-evolving sound speed profiles in shallow water using an ensemble Kalman-particle filter [J]. J. Acoust. Soc. Am., 2013, 133(3): 1377-1386.

[12] PEKERIS C L. Theory of propagation of explosive sound in shallow water [J]. Geologi11cal Society of America Memoris, 1948, 27: 1-116.

[13] EVANS R B. A coupled mode solution for acoustic propagation in a waveguide with stepwise depth variations of a penetrable bottom [J]. J. Acoust. Soc. Amer.,1983, 74(1): 188-195.

[14] PORTER A D, REISS E L. Extension of the method of normal modes to sound propagation in an almost-stratified medium [J]. J. Acoust. Soc. Amer., 1965, 37(1): 19-27.

[15] DAHL P, ZHANG R H, MILLER J, et al. Overview of results from the Asian Seas International Acoustics Experiment in the East China Sea [J]. IEEE J. Oceanic Eng.,2004, 29(4): 920-928.

[16] Becker K M, Frisk G V. The impact of water column variability on horizontal wave number estimation and mode based geoacoustic inversion results [J]. J. Acoust. Soc. Amer., 2008, 123(2) : 658-666.

[17] PORTER M B. The KRAKEN normal mode program[C]. Technical Report, SACLANT Undersea Research Center, La Spezia, Italy, 1991.

[18] LUO X, MOROZ I M. Ensemble Kalman filter with the unscented transform [J]. Physica D: Nonlinear Phenomena, 2009, 238(5): 549-562.