高飛， 孫磊,2

(1.海軍研究院 海洋環(huán)境研究所, 天津 300061； 2.哈爾濱工程大學(xué) 水聲工程學(xué)院, 黑龍江 哈爾濱 150001)

0 引言

水下移動目標(biāo)定位一直以來都是水聲學(xué)研究的熱點問題，廣泛應(yīng)用于海洋漁業(yè)、水下搜救、水聲制導(dǎo)等領(lǐng)域，極具軍事和民用價值。傳統(tǒng)的水聲定位通常使用多個水聽器陣元組合成一定陣型，例如垂線陣、水平陣、L形陣、球面陣等，基于陣列信號采用Bayes匹配場算法處理[1-3]，可達到良好的目標(biāo)定位效果。

單水聽器聲源定位主要利用寬帶聲信號的頻散特征及多途聲波的聲到達時差，可有效節(jié)約觀測成本，并提高水下探測效率。對于水下寬帶聲信號，基于聲波能量傳輸與群速度的相關(guān)性，可利用簡正波理論模型計算得到寬帶信號的模內(nèi)離散、模間離散信息，并依據(jù)各頻率分量時域上的差異構(gòu)建代價函數(shù)，進行聲源定位及聲傳播環(huán)境參數(shù)的反演[4-6]。

多途波是指聲源單次輻射信號沿不同傳播路徑到達同一接收點的聲波，由于傳播路徑總長度的差異，通常可從時域?qū)ζ溥M行區(qū)分。Aubauer等[7]通過聲線幾何關(guān)系推導(dǎo)，在不考慮聲速垂向變化的前提下，利用直達波與海底反射波的聲到達時差對淺海魚群進行定位；Nosal等[8]則通過坐底水聽器數(shù)據(jù)得到直達波、海面反射波聲到達時差，融合實測聲速剖面反演淺海魚群位置；Gebbie等[9]和Duan等[10]則采用連續(xù)聲源信號，基于直達波、海面反射波聲到達時差反演移動聲源位置，并驗證了定位結(jié)果的準(zhǔn)確性。

事實上，海洋環(huán)境的不確定性、信號接收的偶然性及數(shù)據(jù)處理的方法差異都會影響多途聲到達時差的計算結(jié)果，直接關(guān)系到聲源定位的精度。本文采用擴展卡爾曼濾波(EKF)算法進行移動聲源定位，以首達波與次達波聲到達時差為觀測值進行連續(xù)濾波，可有效克服觀測值與實際值間的偏差對定位結(jié)果的影響，得到穩(wěn)定、高精度的聲源定位參數(shù)。

1 定位算法理論

構(gòu)建水平不變深海多途聲傳播環(huán)境如圖1所示。移動聲源(深度zs，水平距離rs)和水聽器(深度zr，水平距離0)位于淺層水域，其深度遠(yuǎn)小于水深H. 圖1中GPS、SOFAR分別為差分全球定位系統(tǒng)、深海聲道。

分別設(shè)直達波、海面反射波、海底反射波、海面海底反射波傳播路徑長度分別為mD、mSR、mBR、mSBR，傳播時間分別為τD、τSR、τBR、τSBR，則有

(1)

(2)

rs通常取104m以上量級，對于深海淺層收發(fā)分置的水聲設(shè)備，zs、zr深度值約為102m量級，故可得mD≈mSR≈rs. 當(dāng)目標(biāo)信號脈沖長度較大時(約500 ms)，直達波與海面反射波重合，難以從時域?qū)ζ溥M行分離，故下文將直達波和海面反射波統(tǒng)稱為首達波。

(3)

(4)

值得一提的是，在實際海洋環(huán)境中，受聲速垂向變化的影響，聲波并非沿直線傳播，上文聲線路徑長度是一種近似推導(dǎo)過程。

1.1 多途波識別

受海洋環(huán)境噪聲影響，通常難以直接識別原始水聽器時間序列數(shù)據(jù)(見圖2(a))中的目標(biāo)信號。對于頻率已知的拖曳式移動聲源，利用功率譜計算或帶通濾波處理原始信號，可分別從頻域和時域識別目標(biāo)信號。

選取水聽器與移動聲源相對距離為10.5 km的一段水聽器原始數(shù)據(jù)(見圖2(a))，移動聲源間隔10 s發(fā)射一次中心頻率為500 Hz、脈沖長度為500 ms的正弦波激勵信號，分段計算其功率譜，其中頻率為400～600 Hz的功率譜計算結(jié)果如圖2(b)所示。從中不難發(fā)現(xiàn)，聲波多路徑到達特征明顯，圖2(b)中方框、圓圈、虛線方框分別標(biāo)注不同路徑到達波。

為確認(rèn)多途到達波的具體傳播路徑及相對能量大小，通過Butterworth帶通濾波處理時域信號(見圖2(c))。對比電壓幅度可知，到達聲波C、A、B能量依次增加。采用Bellhop射線模型(模型代碼下載地址：http:∥oalib.saic.com)模擬試驗數(shù)據(jù)海洋環(huán)境中的聲線分布特征(見圖2(d))，前3個到達波束分別為首達波、次達波、第2次海底反射波，其相對到達時間與聲波A、B、C吻合較好，故可判定聲波A、B、C分別為首達波、次達波、第2次海底反射波，且次達波能量最大，首達波次之，第2次海底反射波最小，故下文選擇直達波與第1次海底反射波的到達時差進行聲源定位。詳盡多途分類及各到達波的能量相對大小可參照文獻[11]。

1.2 多途波時延算法

相關(guān)性函數(shù)用于描述隨機過程中不同時段或空間內(nèi)的相互關(guān)系，對于連續(xù)水聲信號時間序列，時域相關(guān)性函數(shù)[12-13]可表示為

(5)

式中：相關(guān)系數(shù)Γ為時延Δτ的函數(shù)；τU、τL分別為分段數(shù)據(jù)的時間上限、下限(s)；pm、pe為匹配、實測聲壓時域波形數(shù)據(jù)。通過pm=sin (2πfτ)計算匹配聲壓，其中f=500 Hz為頻率，τ=τL∶1/fs∶τU，fs=20 000 Hz為采樣頻率。

本文取匹配聲壓長度為0.5 s(與移動聲源激勵信號長度相同)，即τU=τL+0.5，沿原始采樣時域數(shù)據(jù)逐個采樣點平移，并獲取各組最大相關(guān)系數(shù)。以計算得到的與首達波、次達波的最大相關(guān)系數(shù)分別確定其相對到達時間τF、τS，進而計算聲到達時差ΔτF-S(見圖3)。

1.3 EKF算法

EKF算法通過建立非線性系統(tǒng)模型，主要為描述過程變化動態(tài)的狀態(tài)方程和系統(tǒng)參量與測量值間關(guān)系的測量方程，以實現(xiàn)其最小方差條件下的遞推，從而達到多維平穩(wěn)、非平穩(wěn)隨機過程的最優(yōu)估計。對于速度為v(m/s)的移動聲源，其狀態(tài)方程可表述為

xk|k-1=Φk-1|k-1+qk-1，

(6)

式中：狀態(tài)矢量x=[rs,vr,zs,vz]T，vr、vz分別為水平、垂直方向的聲源移動速度；xk|k-1為k-1時刻得到的k時刻狀態(tài)預(yù)測結(jié)果；k-1|k-1為k-1時刻的狀態(tài)最優(yōu)估計值；qk-1～N(0，αk-1)為過程噪聲，αk-1為過程噪聲均方差；狀態(tài)轉(zhuǎn)移矩陣Φ可表示為

(7)

Δt為迭代時間步長(s)。在水平距離不變的深海環(huán)境中，k時刻的首達波與次達波聲到達時差ΔτF-S可表示為聲速剖面c(z)、水聽器位置以及H、rs、zs的函數(shù)，假設(shè)c(z)、H、水聽器位置已知，則測量方程可表示為

ΔτF-S,k=h(rs,k|k-1,zs,k|k-1)+pk,

(8)

式中：pk～N(0，βk)為測量噪聲，βk為測量噪聲均方差；函數(shù)h(·)為聲線傳播模型，無具體表達形式。于是測量方程可用Taylor級數(shù)展開，并取1階近似：

ΔτF-S,k≈h(s,k-1|k-1,s,k-1|k-1)+s,k-1|k-1)+s,k-1|k-1),

(9)

式中：上標(biāo)“^”表示已得到的參數(shù)估計值。利用Bellhop模型計算得到的首達波與次達波聲到達時差如圖4所示。

圖4中水聽器深度位于30 m，聲速剖面如圖5所示，海表聲速混合層厚度約45 m，45～200 m水層為強躍層。從中不難發(fā)現(xiàn)ΔτF-S隨聲源深度、水平距離呈負(fù)梯度變化特征，即(9)式中的?h/?zs、?h/?rs均小于0，說明隨著移動聲源深度、距離的增加，ΔτF-S逐漸減小。

依據(jù)(9)式，測量方程的雅克比矩陣為

(10)

(10)式中的?h/?rs、?h/?zs雖無具體表達式，但可依據(jù)Bellhop射線模型計算各空間點處的聲到達時差(見圖4)，進而得到聲到達時差隨聲源距離、深度變化的梯度值。

按照非線性離散系統(tǒng)下的EKF進行狀態(tài)更新迭代：

Pk|k-1=Φk-1|k-1ΦT+αk-1，

(11)

(12)

k|k=xk|k-1+Kk(ΔτF-S,k-Hkxk|k-1)，

(13)

k|k=Pk|k-1-KkHkPk|k-1，

(14)

式中：Pk|k-1為基于k-1時刻預(yù)測的k時刻協(xié)方差矩陣；k|k為k時刻的協(xié)方差最優(yōu)估計值矩陣；Kk為k時刻的EKF增益矩陣。(11)式、(12)式、(13)式、(14)式分別為預(yù)測協(xié)方差方程、EKF增益方程、狀態(tài)更新方程、協(xié)方差更新方程。EKF聲源定位的算法具體步驟如下：

1)k=0，狀態(tài)矢量0|0、協(xié)方差0|0初始化；

2)基于(6)式，利用k-1時刻狀態(tài)k-1|k-1預(yù)測k時刻的狀態(tài)xk|k-1；

3)基于(11)式，計算預(yù)測協(xié)方差矩陣Pk|k-1；

4)依據(jù)預(yù)測狀態(tài)xk|k-1和Bellhop射線模型計算聲到達時差梯度Hk，進而基于(12)式，計算EKF增益Kk；

5)基于(13)式，輸入聲到達時差ΔτF-S，得到更新狀態(tài)k|k；

7)單次迭代計算結(jié)束，返回第2步進行k+1時刻的EKF.

2 聲源定位與結(jié)果驗證

試驗海域水深H≈4 150 m，移動聲源間隔10 s發(fā)射中心頻率f=1 000 Hz的正弦波，聲速剖面如圖5所示，利用位于30 m水深的水聽器數(shù)據(jù)計算首達波與次達波聲到達時差(見圖6)。

利用1.2節(jié)中的時域相關(guān)性算法計算得到首達波與次達波聲到達時差，所用數(shù)據(jù)時長為30 min，ΔτF-S的測量值分布于1.4～1.8 s之間，且隨時間呈減小的趨勢。根據(jù)圖4，可初步推測移動聲源的水平距離逐漸增加，或者所處深度逐步增大。

本文中深海水聲調(diào)查區(qū)域，海表聲速混合層約45 m，受海洋聲速剖面影響，混合層以下直達波迅速減弱消失；海表至約20 m水層受散射衰減影響，直達波幅度較弱，故選擇30 m水層數(shù)據(jù)。

采用壓力傳感器采集記錄移動聲源實時深度，移動聲源與水聽器間的水平距離通過GPS信號計算得到，水平移動速度通過測量船航速共享可得。由于無移動聲源垂向速度實測數(shù)據(jù)，這里不對vz進行分析。水聲試驗過程中，連接移動聲源的拖纜放纜長度保持不變，聲源深度約49.1 m，且起伏較小。移動聲源以3.40 m/s沿直線近似勻速離開水聽器垂線陣。

分析可知，EKF得到的聲源位置與實測數(shù)據(jù)總體吻合較好。受海洋環(huán)境及試驗過程中的不確定性影響，試驗獲取的首達波與次達波聲到達時差隨距離和深度呈現(xiàn)嚴(yán)格單調(diào)的變化趨勢(見圖6)，利用(5)式計算的ΔτF-S具有一定誤差，這也是導(dǎo)致濾波結(jié)果起伏的重要因素。

計算得到EKF迭代穩(wěn)定后的聲源深度(見圖7(b))和速度(見圖7(c))的平均值分別為48.8 m、3.42 m/s，實測數(shù)據(jù)的均值分別為49.1 m、3.38 m/s，二者差別較小，驗證了本文聲源定位結(jié)果的準(zhǔn)確性。

3 結(jié)論

1)深海淺層水域是海洋魚群、水下航行器等活動頻繁的水層，也是聲納等水下探測設(shè)備經(jīng)常布放的深度。本文利用深海水聲調(diào)查數(shù)據(jù)，并依據(jù)其建立了聲源定位模型，采用時域相關(guān)性分析計算多途聲到達時差，以聲源深度、水平距離、速度為變量參數(shù)構(gòu)建了EKF狀態(tài)方程((6)式)，并結(jié)合Bellhop射線模型得到EKF測量方程((8)式)，考慮聲速剖面，進而進行了移動聲源定位。

2)通過將EKF定位結(jié)果與壓力傳感器、GPS及母船航速對比，結(jié)果表明迭代穩(wěn)定后的定位結(jié)果與實測數(shù)據(jù)吻合較好，可驗證本文算法的可靠性。

3)準(zhǔn)確的識別多途信號中的直達波與第1次海底反射波并計算其聲到達時差，是影響EKF算法精度的重要因素，需要考慮深海聲速剖面、波導(dǎo)截止頻率對淺表層直達波的影響。