劉炎堃 郭永剛 李整林 李風華 張 波

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

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

0 引言

深海水下運動目標被動定位一直是水聲研究的重點內(nèi)容之一。匹配場算法已在聲源定位領域中獲得了廣泛應用，其基本原理是已知聲速和海底深度等參數(shù)，對聲場進行建模，通過比較測量場和拷貝場來對聲源進行定位[1]。匹配場算法的分辨率受海洋水文、海底底質(zhì)等環(huán)境參數(shù)影響[2]。當海洋的聲速剖面未知時且信號記錄時間足夠長的情況下，McCargar 等[3]、Kniffin 等[4]基于深度的信號分離技術(shù)，可以利用直達波和海面反射波的信號相關(guān)結(jié)構(gòu)來對聲源定位；該方法使用了改進的傅里葉變換來還原目標的位置。Lei 等[5]利用布放在深海的兩個水聽器之間的互相關(guān)函數(shù)，對所有可能的位置進行掃描，使互相關(guān)函數(shù)值達到最大的位置即被視為是聲源位置。匹配場技術(shù)和相關(guān)結(jié)構(gòu)定位技術(shù)均需要兩個或多個水聽器同時進行工作，且需對空間進行掃描，往往有很大的計算量。

海洋聲道的多徑特性同樣可以被用于定位聲源。孫梅等[6]研究了在射線模型下，水平振速與垂直振速的傳播損失與聲線到達接收點處的掠射角以及收發(fā)水平距離之間的關(guān)系，分析了深海直達波區(qū)域聲傳播特性。王夢圓等[7]在此基礎上，提出利用脈沖聲信號的直達波和海面反射波的到達時延，研究了估計水下聲源距離的方法。Gong 等[8]分析了使用拖曳水平陣進行被動目標定位的方法，并驗證了卡爾曼濾波器方法對于運動聲源的良好效果。Baggeroer等[9]與Duan等[10]研究了在深海可靠聲路徑中，位于臨界深度以下的水聽器可以在高信噪比的條件下接收到信號直達波和海面反射波，從而可以準確獲得聲傳播的多徑時延。Yang等[11?13]研究了在可靠聲路徑下，使用布設于深海的水聽器或水聽器陣列，利用直達波和海面反射波的到達時延和信噪比等信息對聲源進行定位。在使用寬帶信號時，僅僅使用單水聽器就能通過自相關(guān)函數(shù)獲得信號傳播的時延信息；利用擴展卡爾曼算法，使用時延信息進行位置估計可以大大減少計算量[13]。然而在該方法中，定位算法需要先定義聲源的初始狀態(tài)，再進一步對聲源進行定位算法的迭代。對聲源的初始狀態(tài)的定義不同將導致最后算法的定位效果的不同。

本文在多徑時延算法[13]的基礎上，提出了一種利用直達-海面反射波時延來對運動目標進行深度估計的算法，對聲源可能的運動路徑預先進行了選擇，從而避免了對初始狀態(tài)的定義，減少了需要的先驗信息。對于每條運動路徑，都可以看作是卡爾曼濾波的算法的觀測輸入，進而得到對運動的預測。通過比較預測位置的時延和理論計算出的時延，可以選擇出最優(yōu)的路徑組，從而完成對目標深度的估計。

1 路徑預選擇與估計

根據(jù)射線模型，由聲源位置到達接收位置的聲線將經(jīng)歷水體折射和邊界反射。圖1 為從聲源發(fā)出，到達間距為?d的兩個接收點的直達聲線與海面反射聲線的傳播示意圖。其中，S為聲源點，S′為聲源關(guān)于海面的鏡像點；R1、R2為接收點。利用聲線跟蹤技術(shù)可以確定到達目標參考點的一系列特征聲線。柱面坐標(r,z)系下的高斯射線方程[14]表示為

其中，r=r(s)及z=z(s)是射線的柱坐標，它們是弧長s的函數(shù)；c(r,z)為聲速。設聲線的曲率為p(s)，寬度為q(s)，單條高斯射線可以被表示為

其中，A為常數(shù)，n為與中心射線的垂直距離，ω為聲源信號的角頻率。τ(s)為傳播時間即相位延遲，滿足：

據(jù)此，可以從理論上得到射線的傳播時間，從而計算出直達波和海面反射波之間到達時間的時延。

圖1 聲源-接收器聲線示意圖Fig.1 Schematic diagram of source-receiver acoustic rays

圖2 為海深5500 m、位于深度4500 m 的水聽器對應的仿真時延圖，仿真中使用的是Munk 聲速剖面。從圖2 中可以看到，給定一個時延，有無數(shù)個可能的深度、距離與該時延值對應，這些深度、距離形成的位置組成了一條時延線；整個時延圖呈現(xiàn)出明顯的隨深度變化的趨勢。當聲源目標運動時，水聽器能接收到一系列的直達、海面反射波信號，從而可以計算得到一系列的時延線。每條時延線可以看作是所有點的一個子集，一個可能的路徑可以看作是在n條時延線上的點的連線；每條路徑都可以看作是時延線之間的一條連線，對運動路徑的預選擇即尋找可能表示目標運動的連線，這可以通過動態(tài)規(guī)劃的方法實現(xiàn)。

圖2 水聽器對應的直達-海面反射波時延圖Fig.2 Time-delay map of direct-to-surface reflection waves for hydrophones

1.1 路徑的預選擇

為了確定所有時延線之間的連線中哪些可以更好地表示目標的運動，假設某個由N條時延線組成的路徑為S，定義代價函數(shù)如下：

其中，Cost(S)為路徑S的代價函數(shù)，ak為包含一個點的距離Rk、深度Dk和速度vk的矩陣，即

其中，、分別為在ak狀態(tài)下的水平、垂直速度。函數(shù)F(ak?1,ak)為ak?1和ak兩個點之間的代價。為了限制運動的速度、方向變化，函數(shù)F(ak?1,ak)定義如下：

其中，F(xiàn)(ak?1,ak)中的第一項α|vk?1?vk|2是對速度矢量變化的限制，而第二項是對垂直方向位移的限制。參數(shù)α和β體現(xiàn)了對目標運動的大致推測：α設置的越大，趨于勻速的連線代價越?。沪略O置的越大，趨于水平的連線代價越小。這樣設置的目的是為了讓路徑趨于水平、勻速運動的路徑，以模擬真實的目標運動。

假設已經(jīng)獲取了M個時延，對應在時延圖上有著M個點集。設每個點集Am(m=1,2,···,M)上都有N個點，記為Am,km(km=1,2,···,N)。設以Am,km為終點的路徑為Sm,km，Sm,km上的所有的點可以表示為從m個點集上各取一個點后所形成的集合：

則代價函數(shù)的計算可以表示為

涉身的含義是什么？萊考夫涉身心智觀中的“‘涉身’是指人類‘依據(jù)生物性的能力以及身體與社會經(jīng)驗在外部環(huán)境中運作’，‘心智’則包括理性與概念，因此心智涉身性就是指人類理性和概念的建構(gòu)都是涉身的，它們通過人類身體、大腦以及與外部世界的互動式運作而形成?！盵3]30

公式(8)顯示了每一個以Am,km為終點的候補路徑的代價函數(shù)計算都僅與以前一個點集中的點為終點的候補路徑代價函數(shù)Cost(Sm?1,km?1)(km?1= 1,2,···,N)有關(guān)。至此，對于任意的km(km= 1,2,···,N)，都可以在O(N)的時間內(nèi)計算出Cost(Sm,km)。在深度上均勻選取q條候選路徑進行運動估計。q越大，路徑越密集，深度分辨率越高；q越小，路徑越稀疏，深度分辨率越低。在實際實驗中，q被設置為50，路徑終點之間的深度間隔約為3～5 m。在計算代價函數(shù)時記錄下該輪計算里的每一個點對應的上一輪計算中的點，即可直接回溯得到路徑。在給定km時，回溯的公式如下：

圖3 為通過回溯計算后得到的部分候選路徑，實際使用的候選路徑均勻分布在選定的深度范圍(從最淺的時延線深度到300 m)內(nèi)，為了圖片的清晰并未全部畫出。計算路徑代價的算法時間復雜度為O(MN)，路徑回溯的算法時間復雜度為O(qM)。由于qN，O(qM)O(MN)，因此，路徑的預選擇算法的時間復雜度為O(MN)，和整個空間中點的個數(shù)呈線性關(guān)系，而當每次迭代更新路徑時，需要的計算時間復雜度僅為O(N)。這種只與前一輪計算的結(jié)果相關(guān)的無后效性有利于減少大數(shù)據(jù)集和在線運算情況下的計算復雜度。

圖3 由路徑預選擇算法得到的部分候選路徑Fig.3 Candidate paths pre-selected by path choosing algorithm

1.2 基于卡爾曼濾波進行運動狀態(tài)估計

卡爾曼濾波器是一種線性系統(tǒng)中高效的自回歸濾波器[15]。在假設噪聲為高斯過程的情況下，卡爾曼濾波可以給出滿足最小均方誤差的結(jié)果。由于選擇的運動路徑中的點是離散的，且運動可以看作是一個線性的系統(tǒng)，卡爾曼濾波器可以很好地對運動進行估計。

設狀態(tài)向量為xk，xk的定義如下：

設狀態(tài)對應的白噪聲為ωk，則狀態(tài)方程為

其中，?t為接收到信號的時間間隔，ωk服從分布N(0,Q)，Q為狀態(tài)誤差的協(xié)方差為

其中，和為測量誤差在深度和距離方向的標準差，B的值為

以每一個預選擇路徑上的點ak作為目標狀態(tài)的觀測，設測量誤差為uk，則ak可以表示為

使用多條預選路徑組成的集合P 來進行聲源的深度估計。P 里的路徑有著所有預選路徑中最小的時延誤差，即對于任意的p ∈P和qP，有

圖4 卡爾曼濾波做運動估計示意圖Fig.4 Schematic diagram of Kalman filter results

圖5 算法流程圖Fig.5 Algorithm flow chart

隨著算法的進行，每一條預選路徑在每接收到一個新的時延時都會進行伸展，向外擴展出一個新的點。對新的點使用方程(10)進行計算，會得到新的路徑組。因此，每得到一個時延，都會對估計的深度進行一次迭代。完整的算法流程見圖5。

2 實驗介紹

實驗的數(shù)據(jù)來自2018年春季在南中國海某次海上聲學實驗。實驗采用單船結(jié)合接收陣垂直潛標的方式，接收陣為非均勻分布在120～3408 m深度范圍的24 陣元的水聽器陣列，如圖6 所示。實驗過程中，實驗船以大致3 n mile/h 的速度遠離接收陣(圖7)。拖曳換能器聲源的深度大致保持在150 m左右，聲源發(fā)射的信號中心頻率為300 Hz，帶寬為100 Hz。

圖6 實驗示意圖Fig.6 Schematic of the experiment

實驗過程中的距離水聽器2 km～35 km的海深分布如圖8(a)所示，可見海底較為平坦，平均海深約為3300 m，接收水聽器陣列所在海深為3500 m。計算使用的數(shù)據(jù)為船距離接收陣1.9 km～9.6 km 的信號數(shù)據(jù)。實驗期間測得的聲速剖面見圖8(b)。海底聲學參數(shù)被設置為聲速1565 m/s，密度1.6 g/cm3，吸收系數(shù)0.3 dB/λ[16]，用于本征聲線的計算。圖9 為基于聲速剖面計算得到的深度位于3020 m、3370 m 的水聽器的本征聲線，可見聲源到接收點的直達波和海面反射波的聲線路徑。在實驗過程中使用同水聽器固定在一起的壓力傳感器得到的位于2320 m、2370 m、3020 m、3370 m的水聽器深度較穩(wěn)定，可忽略接收深度變化的影響。在實驗中對接收到的信號使用脈沖壓縮方法進行截取。對4 個不同深度的水聽器同一時刻接收到的信號處理如圖10所示，可以看出信號到達的時間清晰的分為幾組，與利用射線模型仿真的結(jié)果一致，其中，4 條線上起伏的幅度代表了脈沖壓縮后的能量大小，直達波和一次海面反射的聲波位于圖10 中紅色方框標示部分。直達波與海面反射波的到達時延值可以從接收信號的自相關(guān)函數(shù)的峰值中獲得[13]。自相關(guān)函數(shù)的定義如下：

其中，T為選取的時間窗長，s(t)為接收信號。

圖7 航行中“實驗1 號”的速度變化Fig.7 Speed change of ship “Shiyan No.1”

圖8 實驗測得海深及聲速剖面Fig.8 Sea depth and sound profile

圖9 不同接收深度的本征聲線圖Fig.9 Acoustic rays in different receiver depth

圖10 4 個不同深度的潛標信號進行脈沖壓縮后的結(jié)果Fig.10 Pulse compression of four submersible signals at different depths

利用聲速剖面，預先計算得到2320 m、2370 m、3020 m、3370 m 深的水聽器對應的直達波與海面反射波的時延圖，參見圖11。獲得接收信號時延之后，即可對不同深度的水聽器在圖11上尋找對應的時延線。根據(jù)接收信號中得到的時延信息，利用公式(9)計算獲得路徑，并將路徑看作是卡爾曼濾波的一次觀測進行運動估計。將運動估計的結(jié)果與預先計算的時延圖中進行比對，可以得到最接近的路徑，從而迭代一次估計的聲源深度。在迭代的過程中，不同深度的水聽器接收到的信號在計算上是完全獨立的。由于少數(shù)時延值不一定能獲取到穩(wěn)定的最優(yōu)路徑組，在第一次迭代之前預先使用了10 個時延值進行了路徑選擇算法。

圖12 為實驗中使用的各個深度水聽器信號中得到的前40個時延值。從圖12中可以看出，隨著聲源逐漸遠離水聽器陣，各水聽器接收到的直達波和海面反射波之間的時延值逐漸下降，且越深的水聽器對應的時延值越高，與圖11一致。其中，由于噪聲干擾，位于2370 m的水聽器計算時延值出現(xiàn)了錯誤的峰值的情況，使時延值曲線產(chǎn)生了較大起伏。

圖11 不同接收深度下的直達波與海面反射波的時延圖Fig.11 Time delay map of direct waves and sea reflections under different receiver depth

圖12 使用互相關(guān)得到的實驗中各深度水聽器前40 個時延值Fig.12 The first 40 delay values of each depth hydrophone in the experiment using cross correlation

圖13 比較了用線性插值替換異常的時延值與異常時延值對應的時延線，可見異常值將導致目標可能的位置出現(xiàn)偏差，使路徑產(chǎn)生一定程度的畸變，導致公式(8)定義的代價函數(shù)值突增。圖14展示了算法對實驗數(shù)據(jù)進行深度估計的結(jié)果。圖14(a)為單獨使用某個水聽器的數(shù)據(jù)計算的深度迭代圖，其中藍色帶星號的曲線是由與拖曳聲源固定在一起的測深儀提供的深度數(shù)據(jù)；從圖14(a)可以看到，當使用深度位于2320 m、3020 m、3370 m處的水聽器接收到的信號進行聲源深度估計時，估計的相對誤差大致在10%以內(nèi)。使用深度為2370 m 的水聽器接收到信號進行深度估計時，由于時延數(shù)據(jù)的不準確，迭代初期的最優(yōu)路徑選擇并不穩(wěn)定，估計深度變化較大，導致在前20次深度迭代總深度估計的誤差大，最大值為第一次估計，相對誤差的絕對值達到70%；在第23 次迭代后，深度估計的誤差下降到了15%以內(nèi)。去除異常的時延值后，2370 m 深水聽器的深度估計曲線見圖14(b)。圖14(b)說明在時延值計算出現(xiàn)較大誤差時，將使估計的深度變得不準確；而后在時延值計算準確時，2370 m 深水聽器的深度估計結(jié)果誤差變小。因此，計算過程中可以在每次迭代前預先拋棄一定迭代次數(shù)以前的時延值來進行下一步計算，避免原異常數(shù)據(jù)對后續(xù)計算的影響。利用4 個水聽器的數(shù)據(jù)同時對目標深度進行估計，把每次估計迭代中代價函數(shù)值最低的前5 個深度都統(tǒng)計進結(jié)果中，則深度估計呈現(xiàn)的離散序列針狀圖見圖15。圖15 中每一列的深度寬度范圍為5 m。由于直達波與海面反射波的時延在時延圖上對應的是一條線，在數(shù)據(jù)量少的時候可能出現(xiàn)深度偏差比較大的估計。在圖15中，深度位于區(qū)間[0 m,100 m]和[200 m,300 m]的估計點數(shù)約占所有點數(shù)的25.9%。對整個估計的離散點序列進行高斯曲線擬合：

圖13 2370 m 水聽器的時延異常值產(chǎn)生的時延線偏移Fig.13 Time delay line offset due to the outliers value of 2370 m hydrophone

圖14 深度估計結(jié)果Fig.14 Depth estimation results

得到的置信水平在95%下的擬合參數(shù)見表1。從曲線的擬合結(jié)果來看，多次迭代中對聲源運動中的深度估計值大致分布在真實深度的附近，分布的均值為145.8 m。

圖15 各水聽器每次估計中代價值前五的深度統(tǒng)計結(jié)果Fig.15 Statistical results of hydrophones in top five depth estimation

表1 置信水平為95% 的深度估計曲線擬合參數(shù)Table 1 Curve fitting parameters for depth estimation with a confidence level of 95%

3 結(jié)論

本文提出了一種基于單水聽器進行水下運動目標深度估計方法。利用直達-海面反射波時延預先構(gòu)建出候選運動路徑，避免了對目標初始狀態(tài)的定義；使用候選路徑對目標運動的軌跡進行模擬，從而實現(xiàn)了對目標的深度估計。使用候選路徑的方法可以避免每次迭代在整個空間進行搜索，減少了計算復雜度。實驗結(jié)果驗證了該方法的有效性。由于需要獲取信號的到達時延值，該算法僅適用于深海直達聲區(qū)。若運動目標在深度或運動速度大小上變化較大，將導致路徑選擇算法的失效，后續(xù)研究將進一步分析路徑選擇對最終深度估計的影響。致謝 感謝參與2018年4月南中國海調(diào)查全體實驗人員，他們的辛勤勞動為本文提供了珍貴可靠的實驗數(shù)據(jù)。