劉炎堃 郭永剛 李整林 李風(fēng)華 張 波

(1 中國科學(xué)院聲學(xué)研究所 聲場聲信息國家重點(diǎn)實(shí)驗(yàn)室 北京 100190)

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

0 引言

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

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

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

1 路徑預(yù)選擇與估計(jì)

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

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

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

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

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

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

圖2 水聽器對應(yīng)的直達(dá)-海面反射波時(shí)延圖Fig.2 Time-delay map of direct-to-surface reflection waves for hydrophones

1.1 路徑的預(yù)選擇

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中，和為測量誤差在深度和距離方向的標(biāo)準(zhǔn)差，B的值為

以每一個(gè)預(yù)選擇路徑上的點(diǎn)ak作為目標(biāo)狀態(tài)的觀測，設(shè)測量誤差為uk，則ak可以表示為

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

圖4 卡爾曼濾波做運(yùn)動估計(jì)示意圖Fig.4 Schematic diagram of Kalman filter results

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

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

2 實(shí)驗(yàn)介紹

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

圖6 實(shí)驗(yàn)示意圖Fig.6 Schematic of the experiment

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

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

圖7 航行中“實(shí)驗(yàn)1 號”的速度變化Fig.7 Speed change of ship “Shiyan No.1”

圖8 實(shí)驗(yàn)測得海深及聲速剖面Fig.8 Sea depth and sound profile

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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