一種基于序貫估計的直達聲區水面艦船被動測距方法?

張雪冬 牛海強 吳立新

(1 中國科學院聲學研究所 聲場聲信息國家重點實驗室 北京 100190)

(2 中國科學院大學 北京 100049)

0 引言

直達聲區水面聲源的被動測距一直是研究熱點問題，在實際應用中具有重要意義[1?4]。傳統的匹配場測距方法通過拷貝場和測量場之間的匹配對距離進行估計，但該方法依賴于對海洋環境參數、尤其是海底參數的準確估計。近年來，基于射線的盲解卷積方法(Ray-based blind deconvolution,RBD)[5?6]被廣泛應用。它利用常規寬帶波束成形，通過將垂直陣(Vertical line array,VLA)的導向向量指向特定的射線路徑(例如直達波)來估計隨機聲源的相位，進而得到信道響應(Channel impulse response,CIR)。當聲源為可控通信聲源或船舶噪聲時，RBD方法可以只利用接收陣型和陣列位置處的聲速剖面(Sound speed profile,SSP)來估計CIR。已有的研究已經證明了利用RBD方法估計CIR的可靠性[7?12]。目前RBD方法被成功應用于聲源定位[7?9]、接收陣陣型估計[10?11]和海底參數反演[12?13]中。在利用RBD方法進行被動測距方面，已有的研究主要是通過基于匹配場[7]、基于聲線多途時延差[8]和陣列不變量[9]等方法，對RBD估計到的CIR(簡稱為RBD-CIR)進行處理。但海底參數未知時，匹配場方法和基于聲線多途時延差的方法不再適用；在復雜波導中，如聲速剖面存在躍層或海底為多層分布時，對陣列不變量的估計也變得復雜。本文利用RBD-CIR的直達波相對參考陣元的到達時間差和聲速剖面信息，基于序貫方法，對聲源-接收距離進行估計，避免了對海底參數和陣列不變量的估計，并利用2016年美國圣巴巴拉海峽實驗數據驗證了該方法的有效性。由于受實驗環境海深的限制，直達聲區僅在<3.5 km的范圍內，但此方法仍適用于深海環境的直達聲區，測距范圍可達幾十千米[14?15]。

序貫方法利用狀態-空間方程來估計和更新系統的狀態值，利用前一狀態的估計值和更新值來對下一個狀態進行預測，其估計值和更新值由關聯了聲學測量量與待同化量的聲傳播模型來得到[16]。常用的序貫方法有卡爾曼濾波[16?17]、集合卡爾曼濾波(Ensemble Kalman filter,EnKF)[18?19]和粒子濾波[17,20]等。這幾種方法均在海洋目標跟蹤定位方面進行了應用并取得了很好的效果。本文選擇了EnKF作為序貫估計濾波器來對聲源-接收距離進行估計。EnKF 基于貝葉斯原理[21]，利用集合成員分布來模擬待估計參數，并假設過程中的先驗概率分布、似然函數和后驗概率分布均服從高斯分布，在非線性系統中有較好的表現。

本文主要分為以下幾個部分：首先介紹基于射線的盲解卷積方法和基于狀態-空間模型的EnKF方法；第二部分對2016年美國圣巴巴拉海峽實驗數據進行處理，并將序貫估計結果與測量值和幾何方法的結果進行了比較，并對結果進行了分析；第三部分比較了本文方法和傳統匹配場方法的測距結果；最后對本文的工作進行了總結。

1 理論

1.1 RBD方法[5?6]

假設位于海洋信道中點源發射的時域信號為s(t)，其頻域表達式為S(ω)=|S(ω)|eiΦs(ω)。則第j(j=1,2,···,N，N為陣元個數)個VLA陣元所接收到的信號可表示為

其中，zs為聲源深度，zj為第j個VLA陣元深度，G(zj,zs,ω)是聲源和第j個VLA陣元之間的格林函數。對于適用于射線方法的高頻信號，各個射線路徑的到達角θ可以通過N個陣元的接收信號做常規波束形成得到

其中，到達第j個VLA陣元的第k條射線與到達參考陣元的時延差τj,k可以通過常規波束形成得到，T(θk)與頻率幾乎無關，主要由聲源和VLA之間的射線傳播時間決定[22]。RBD方法利用接收信號波束形成B(ω;θk;N)后的相位

來對接收信號的相位進行旋轉，從而抵消未知的聲源相位譜Φs(ω)，因此有

在大多數淺海信道中，信道響應(公式(4))中的分母在頻率上足夠恒定，所以可以通過平方根項的歸一化消除聲源信號幅度對估計結果的影響。最后，將RBD-CIR的直達波的相對到達時間差作為測量值，利用EnKF，可對聲源-接收距離進行估計。

1.2 EnKF[21]

序貫估計提供了一個狀態-空間模型，當出現新的觀測值時可對狀態值進行預測和更新，它模擬了狀態向量與觀測向量之間的關系，該模型由狀態方程(5a)和測量方程(5b)構成：

在本文中，狀態向量xk為聲源深度和聲源-接收距離在第k個時間步長時的值。狀態方程預測算子f(·)模擬了狀態向量隨時間步長的變化過程，由于本文不考慮聲源的運動模型，所以該算子取為單位矩陣I。狀態噪聲vk表示了狀態方程預測結果和真實結果的誤差，通常情況下，認為其服從高斯分布。本文中觀測向量yk為直達波在不同接收深度上的相對到達時間差(以中間陣元為參考陣元)，該時間差從RBD-CIR中提取得到。測量算子h(·)代表狀態向量xk與觀測向量yk之間的關系，在本文中利用射線數值模型Bellhop[23]來得到。測量噪聲wk描述了測量值yk與狀態方程預測值h(xk)之間的誤差。

EnKF 本質上利用采樣集合的概率密度分布來描述狀態向量xk和測量向量yk，將集合的均值作為估計結果。假設狀態向量的初始集合為=1,2,···,M}并服從正態分布，其中i表示集合成員序號，M表示集合成員個數。當vk和wk服從均值為0，方差為Vk和Wk的高斯分布且互不相關時，公式(5a)和公式(5b)可改寫為

其中，和表示由狀態方程得到的預測狀態向量及其協方差矩陣，表示由測量方程更新后得到的狀態向量，H是h(·)的線性化算子，K表示Kalman增益，上標T表示矩陣轉置。

由于EnKF 用集合中成員的概率分布來表達概率密度函數，則預測狀態向量集合在時間步長k的采樣均值和協方差矩陣可表示為

可以看出，當M趨于無窮時，由于收斂于真實的狀態向量值，因此趨于真實的誤差協方差矩陣。而且狀態-空間模型中包含算子H的項都可由以下公式代替

由此可看出，每個由狀態方程產生的狀態向量成員可以被更新為

則更新后的狀態向量集合會作為初值來遞歸地產生下一步預測狀態向量集合的初值最終的估計結果為狀態向量在收斂階段的后驗概率密度(Probably density function,PDF)的平均值。由于EnKF 無需得到狀態預測算子f(·)和測量算子h(·)的線性化表達式，在非線性系統中具有很大的優勢。

2 實驗數據處理結果

為了驗證上述方法的有效性，本文基于序貫估計利用RBD-CIR提取的直達波相對到達時間差來估計聲源-接收距離，并對實驗數據結果進行了分析。實驗數據來自于2016年9月美國圣巴巴拉海峽實驗(SBCEx-16)，實驗區域為以圣巴巴拉海峽為中心、由34.3?N、120.2?W和34.2?N、120.0?W圈出的方形區域，如圖1(a)所示。該區域位于圣羅莎島北部、加利福尼亞海岸線17 km以外，南北長約11 km、東西長約18 km。實驗區域的海深從西北部(500 m)向東南部逐漸增加(590 m)，而本文所研究的區域位于東南部分，如圖1(b)所示。圖1(b)中虛線為航道邊界，航道寬1 n mile，航道間的間隔也為1 n mile。4個32元的VLA 均勻錨系在航道中心線處的海底。本文僅對其中一個VLA(圖1(b)紅點)的接收信號進行了分析。該VLA錨系在深580 m的海底上方約5.8 m處，由兩個子陣構成，上方的子陣由16個間隔1 m的陣元構成，其中第一個陣元在實驗期間沒有收集數據；下方的子陣由16個間隔3.75 m的陣元構成；兩個子陣之間間隔5.6 m，總有效陣長為75.85 m。接收信號采集系統的采樣頻率為25 kHz。實驗中溫鹽深儀(Conductivity,temperature,depth,CTD)測量到的16個剖面的平均聲速剖面見圖1(d)?？梢钥闯?，聲速在水面附近為1510 m/s，在海深50 m附近迅速降為1491 m/s，并逐漸過渡到海底580 m處的1485 m/s，因此聲速剖面在50 m以上存在躍層。

圖1(b)中的藍色圓點表示商船Anna Maersk在9月16日12:14:30–12:19:30期間每間隔10 s的航跡，由自動測量系統(Automatic identification system,AIS)測得，圖中藍色箭頭表示船行進方向。在此期間，Anna Maersk 距離VLA最遠的位置(點P)有3440 m，而最近有1634 m。受航道條件的限制，實驗過程中不存在距離VLA 小于1.6 km的航船；受該區域海深的限制，當聲源-接收距離大于3.5 km時，航船超出了直達聲區的范圍。所以僅對有限距離范圍內的航船噪聲進行處理。圖1(c)為圖1(b)中綠色圓點所示航跡的海底回聲定位圖像?？梢钥闯觯摥h境下海底呈現出多層分布，且海深隨距離的變化而變化。

以聲源-接收距離r=3440 m為例(即圖1(b)中的點P)，圖2介紹了利用RBD方法估計CIR的詳細步驟。圖2(a)為VLA上第一個陣元(498.35 m)處接收到的5 s 長的艦船輻射噪聲信號。對VLA所有陣元的接收信號做常規寬帶波束形成，處理頻帶為150～1000 Hz，信號時間長度為5 s，得到波束圖如圖2(b)所示，其中直達波和海底反射波的到達角分別為θD=?10?和θB=12?(以水平到達方向為0?，垂直到達方向為90?)。對VLA所有陣元的接收信號在直達波到達方向θD上做波束形成(公式(2))，將波束形成后的信號與接收信號做匹配濾波，然后對幅度進行歸一化(公式(3)～(4))，得到RBD-CIR如圖2(c)所示。從圖中可以觀察到直達波和兩組到達角相同的海底反射波，這也印證了海底多層分布的情況。

圖1 2016年美國圣巴巴拉海峽實驗Fig.1 2016 Santa Barbara Channel Experiment

圖2 利用RBD方法估計信道響應的實現步驟Fig.2 The steps for estimating CIR using RBD method

對圖1(b)中Anna Maersk的航跡上所有接收信號均做常規寬帶波束形成，所得結果隨距離的變化如圖3所示，其中分別用藍色和紅色加號標出了直達波θD和海底反射波θB的到達角。從圖3中可以看出清晰的直達波，并且隨著距離增加而逐漸向0?(水平方向)靠近。利用所有聲源-接收距離r上的直達波到達角θD的波束形成結果估計其對應的RBD-CIR，并提取其直達波在深度zj(第j個VLA陣元)上的到達時間TD(zj)，來計算其相對于參考陣元的到達時間差δT，其中參考陣元為第16個陣元，即δT(zj)=TD(zj)?TD(z16)。所得直達波相對到達時間差結果如圖4中黑色加號所示。為了便于比較，圖4中對每個距離上的時間差相對于上一個距離的結果整體向右平移了5 ms，從左至右的距離逐漸增加。可以觀察到，直達波的相對到達時間差δT隨距離的變化而變化，所以利用該提取量對距離進行估計是可行的。

圖3 接收信號常規寬帶波束形成圖Fig.3 Evolution of the beamformed receiving signals obtained from plane-wave beamforming

利用圖4中的直達波相對到達時間差δT對聲源-接收距離r進行估計。由于聲源深度zs也會對δT產生影響，所以這里分析δT對zs和r的敏感性。由于已知聲源為水面艦船，所以僅對zs=1～20 m的情況進行討論。利用Bellhop計算zs=5 m、r=3440 m(仿真圖1(b)中點P的情況)時的δT，以及zs和r分別在1～20 m和3000～4000 m的范圍內變化時的二者的差值在深度上的平均值(即公式(12))隨zs和r的變化分別如圖5(a)和圖5(b)所示，仿真時所使用的聲速剖面如圖1(d)所示。

圖4 不同距離上，由RBD-CIR提取的和利用EnKF 估計結果得到的直達波的相對到達時間差Fig.4 The direct ray-path arrival time differences extracted from RBD-CIRs and calculated by the EnKF estimation results,respectively,for different source-receiver distances

圖5 敏感性分析Fig.5 Sensitivity analysis

可以看出，相比于聲源-接收距離r，聲源深度zs的變化對δT的影響較小。所以在對r進行估計時，雖然也將對zs的估計包含在內，但由于δT對zs在1～20 m范圍內的變化不敏感，所以zs的估計值并不能作為最后的估計結果，而是為了對r的估計增加一定的寬限，所以本文并未對聲源深度的估計結果進行詳細研究。

將圖4中黑色加號所示的直達波相對到達時間差δT作為觀測值，利用M=20個集合成員的EnKF 濾波器對聲源深度zs和聲源-接收距離r進行追蹤，即xk=[zsr]T。所使用的初值、狀態方程的協方差和搜索區間如表1所示，并令狀態向量協方差的初值等于根據狀態向量的特性，將其協方差設為由于距離r的變化比較大，所以這里針對距離設定了較大的協方差參數。測量方程噪聲協方差取0.05 ms。由于zs對直達波相對到達時間差的影響不大，不妨設其初值為10 m，區間為1～20 m。對聲源-接收距離r的初值分別為1000 m、1600 m和3000 m的情況進行了討論。通常情況下，序貫方法不需要設定搜索區間，但為了使反演結果更符合海洋環境參數的一般設定，這里對狀態向量限定了上下限。但從圖6中聲源-接收距離的后驗PDF 在不同時間步長上的結果中看出，由于PDF 并未碰到搜索區間的上下限，所以上下限的設定對反演結果的影響很小。為了使每個距離的估計結果更加穩定，在每個距離點上均進行了10個時間步長的迭代，則對31個距離點估計的總時間步長為310。因此在每個距離點上需要計算200個點的直達波相對到達時間差。在10個時間步長以后，不同初值的測距估計結果非常接近。在10個時間步長之前，當r0=1600 m時(圖6(b))，由于其接近第一個測距點的測量值1640 m，其后驗PDF 快速收斂到測量值附近；而當r0=1000 m和3000 m時(圖6(a)和圖6(c))，二者后驗PDF 在10個步長后收斂到測量值附近。隨后，聲源-接收距離r的后驗PDF在每個距離點上(10個步長范圍內)均可以較快地過渡到收斂階段，并趨于穩定。

表1 序貫估計的初值及濾波器參數Table1 Initial settings of the sequential parameter estimation

圖6 聲源-接收距離初值不同時，聲源-接收距離的后驗PDF 隨時間步長的變化結果Fig.6 The estimated results of source-receiver distances as evolving marginal posterior PDFs with different initial values of the source-receiver distances

當聲源-接收距離初值r0為1600 m時，將聲源深度zs和距離r的后驗PDF(圖6(b))在每個距離點上最后一個步長的均值作為最終的估計結果，利用Bellhop[23]和圖1(d)中的聲速剖面模擬直達波的相對到達時間差δT，如圖4中紅色虛線所示。為了便于比較，采取了與RBD-CIR提取到的δT同樣的操作，即對每個距離上的時間差相對于上一個距離的結果整體向右平移了5 ms?？梢钥闯龆叻系煤芎?。圖7(a)中紅色叉號虛線表示r0=1600 m時，利用EnKF所得的最終測距結果，并將其與AIS系統測量的結果(藍色圓圈)以及幾何方法的計算結果(rG，灰色實線)進行比較。幾何方法是將信道視為理想波導(海水聲速為常數)，并將聲源深度視為0 m，通過直達波的到達角和接收陣所有陣元的深度，來估計聲源-接收距離，即rG(zj)=?zj/tan(θD),(j=1,2,···,N)。由于VLA 具有一定的長度，所以使用不同的陣元深度會得到不同結果的rG。圖7(b)中比較了當r0=1600 m時，EnKF測距結果(紅色叉號)和幾何方法結果rG(灰色實線)與AIS 測量值的誤差絕對值|δr|。EnKF的測距結果的誤差范圍為1～186 m，每個距離上的誤差小于6%。由于無法確定接收器深度，rG的誤差隨著所使用的接收器深度zj的變化而變化，誤差最大可達675.6 m。若使用所有深度的平均值(526.37 m)作為接收深度，即tan(θD)，此時的最大誤差530.1 m，相對誤差可達16%。

圖7 不同方法估計得到的聲源-接收距離及其與AIS 測量值的誤差Fig.7 The source-receiver distances obtained with different methods and the source-receiver distance differences between the estimations with different methods and the AIS measurements

比較不同距離上幾何方法估計的距離rG與AIS結果的誤差，可以發現隨著距離的增加，誤差也隨之增大。這是由于聲速剖面存在躍層，見圖1(d)。圖8是當聲源-接收距離r分別為1634 m(實線)和3440 m(虛線)時利用Bellhop模型仿真得到的在VLA 平均深度上(526.37 m)的直達波(D，由藍線表示)和海面反射波(S，由紅線表示)的路徑圖。聲源深度假設為10 m，所使用的聲速剖面如圖1(d)所示。從圖8可看出，當聲源-接收距離r=1634 m時，直達波和海面反射波經過躍層的長度較短，受其影響較小，直達波仍可視為聲源和接收器之間的直線，所以幾何方法估計的距離較為準確；但當聲源-接收距離r=3440 m時，直達波和海面反射波經過躍層的長度大大增加，根據Snell定律[22]，射線在此區域會發生較大彎曲，所以無法將直達波路徑視為聲源與接收器之間的直線，幾何方法估計結果的誤差大大增加。

圖8 聲源-接收距離不同時，利用Bellhop 仿真計算得到的直達波(D)和海面反射波(S)的路徑圖Fig.8 The direct(D)and the surface-bounced(S)ray-paths simulated by Bellhop with different source-receiver distances

3 與傳統匹配場方法比較

將本文方法的結果與傳統匹配場方法作比較。匹配場方法[24]利用代價函數在參數空間內進行搜索來估計聲源-接收距離r和聲源深度zs，這里代價函數表示為[24]

其中，F表示頻點數，本文在160～480 Hz中均勻選取了12個頻點；BMF是在頻率為fk時的Bartlett匹配濾波結果，并且有

其中，為第j個陣元的接收信號在fk時的復聲壓，*表示復轉置，為理論聲壓，在這里為利用Kraken聲場模型[25]根據圖9所示海洋環境模型所得聲壓。由海底回聲定位儀所測圖像(圖1(c))以及盲解卷積所得信道響應(圖2(c))可知，該區域海底情況較為復雜，且海底深度也在不斷變化。但為了簡化模型，這里將其假設為兩層分布的水平不變模型，上層為軟沉積層，下層為半無限大的硬基底。沉積層厚度為d，沉積層聲速cs和基底聲速cb、沉積層吸收系數αs和基底層吸收系數αb都假設為不隨深度變化。而沉積層密度ρs和基底密度ρb可以根據Hamilton 經驗公式(公式(15))由沉積層聲速cs和基底聲速cb計算得到(文獻[26]中也曾做過類似處理)

圖9 海洋環境模型Fig.9 The range-independent environment model

利用圖9中的海洋環境模型對r和zs進行估計，并且其搜索區間分別為900～4000 m和1～20 m，步長分別為10 m和1 m，因此在每個距離點上需要計算311個距離、20個聲源深度和12個頻點處的聲壓。由于該區域海底參數未知，所以根據已有的經驗對海底參數進行多種假設，其具體情況如表2所示。最后所估計得到的聲源-接收距離r如圖10(a)所示，其與AIS 測量值的誤差絕對值|δr|如圖10(b)所示。由圖可知，在海底參數取組合1和4時，r的匹配場估計結果在所有距離點上與測量值AIS 符合得較好，此時最大相對誤差分別為12%和9%，仍大于本文方法的誤差6%。若改變沉積層厚度以及其他海底參數(組合2和3)，此時r的估計結果在某些距離點上(例如第29、30個距離點)最大相對誤差均可達到70%；若排除這些點，r的最大相對誤差分別為12%和7%，仍大于本文方法的結果。由此可看出，海底參數的選取對匹配場方法的測距結果影響很大。

表2 匹配場方法所使用的海底參數組合Table2 Geoacoustic parameters for matched-field method

圖10 不同方法得到的聲源-接收距離估計結果及其與AIS 測量值的誤差Fig.10 The source-receiver distances estimated by different methods and the source-receiver distance differences between the estimations and the AIS measurements

4 結論

本文對2016年美國圣巴巴拉海峽實驗數據進行了處理，通過垂直陣在1.6～3.5 km 直達聲區范圍內接收到的水面艦船輻射噪聲數據，利用RBD方法來估計船和接收陣之間的信道響應，提取直達波在垂直陣不同陣元上的相對到達時間差，基于序貫方法對聲源-接收距離進行了估計，并將距離的估計結果與測量值、幾何方法和傳統匹配場方法的結果進行了比較。由于不需要對海底參數和陣列不變量進行估計，該方法在聲速剖面存在躍層和海底多層分布的復雜海洋環境下仍然適用，要優于傳統匹配場方法。本文方法距離估計值與測量值的相對誤差在6%以內，并且小于幾何方法的誤差。若根據已有經驗對海底參數進行假設，本文方法距離估計值與測量值的相對誤差仍小于傳統匹配場方法。由本文方法得到的估計值所仿真的直達波相對到達時間差與測量值吻合得很好。由于受到實驗條件的限制，直達聲區的范圍較小，下一步將對深海環境下該方法的有效性進一步進行驗證。

致謝由衷感謝佐治亞理工學院Karim G.Sabra教授提供的2016年美國圣巴巴拉海峽實驗數據，感謝參加海上實驗的全體工作人員，他們的辛勤勞動為本文提供了寶貴的實驗數據。