禚江浩 王 玲 許 可 萬建偉

(國防科技大學電子科學學院智能感知系, 湖南長沙 410073)

1 引言

被動聲納目標檢測分為窄帶檢測和寬帶檢測。窄帶檢測適用于具有穩(wěn)定線譜的水中目標,國內(nèi)外學者分別運用線譜的窄能量、頻率方差、空時穩(wěn)定性、幅相起伏性、功率譜熵等特征完成窄帶檢測[1,2]。但是窄帶檢測沒有考慮海洋傳播信道隨頻率的變化,而且待檢測的目標往往是非合作的,采用窄帶檢測可能導致穩(wěn)定性差的問題。寬帶檢測運用足夠?qū)挼膶拵盘栠M行目標檢測,往往可以得到更加穩(wěn)定的效果,對寬帶檢測的研究是近年來的一個熱門課題。

對于常規(guī)波束形成,能量檢測法(Conventional Energy Detection, CED)是常用的寬帶檢測方法,是理論上在非相干噪聲場中單一目標的最優(yōu)檢測器。M.Bono 提出了子帶峰值能量檢測(Subband Peak Energy Detection, SPED),用于解決寬帶能量檢測的目標軌跡模糊問題[3]。在SPED的基礎上,楊晨輝提出了波束域峰值能量檢測方法[4],性能優(yōu)于SPED。王聰提出了基于頻率著色的寬帶能量檢測方法[5],提高了方位時間歷程圖的顯示性能。

從匹配場波束形成的角度來看,寬帶檢測方法分為非相干處理方法和相干處理方法。非相干處理方法的處理單元為單個頻點的匹配處理,該方法的輸出為多個頻率點的匹配處理的均值。單個頻點的匹配處理考慮了單個頻點內(nèi)的空間相關(guān)性,但是均值的處理方法忽略了各個頻點之間的相干信息。目前的處理算法大多采用這種空間相干處理、頻率間非相干處理的方法[6]。為了降低旁瓣和提高定位精度,相干處理方法考慮了各頻點之間的相干性,但是這種處理方法的性能依賴于各頻點的歸一化系數(shù)的選擇,對噪聲比較敏感。

最早的寬帶相干處理方法是Clay實現(xiàn)的時域?qū)拵喔善ヅ鋱鎏幚韀7]。該方法將測量得到的脈沖響應與建模得到的脈沖響應相匹配,從而實現(xiàn)了對聲源的定位。Brienzo和Hodgkiss在爆炸聲源定位試驗中對該算法進行了驗證[8-9]。Clay提出的時域?qū)拵喔善ヅ鋱鎏幚硇枰缆曉吹念l譜,這對于被動聲納的應用來說是難以實現(xiàn)的。Hursky在研究高頻段的時域?qū)拵喔善ヅ鋱鎏幚頃r,利用雙水聽器的時域互相關(guān)實現(xiàn)了在聲源頻譜未知時的聲源定位[10]。Westwood實現(xiàn)了頻域的寬帶相干匹配場處理[11],其輸出結(jié)果為頻域互相關(guān)的相干累加。在單個頻點的接收數(shù)據(jù)的歸一化問題上,Michalopoulou提出了以第一個陣元為參考量進行歸一化的方法[12],這種歸一化方法適用于高信噪比的情況。孟華等對被動聲納寬帶相干處理技術(shù)進行了總結(jié),并且提出了一種改進的陣元歸一化方法,提高了寬帶相干處理器的檢測性能[13]。趙博將并行遺傳算法應用于匹配場聲源定位[14],提升了計算速度。周悅將匹配場處理建立在機器學習的框架下,基于信息理論準則對拷貝場和測量場進行距離度量[15],性能優(yōu)于傳統(tǒng)的匹配場估計方法。韋先聲將傳統(tǒng)的基于聲壓水聽器的匹配場處理擴展到基于矢量水聽器的匹配場處理[16],對比了不同匹配量下的相干處理器和非相干處理器的匹配性能。

匹配檢測處理算法的加權(quán)向量為聲源到水聽器的傳播損失聲壓。傳播損失聲壓可以通過建模得到。水下聲傳播模型按照建模的維度不同,可以分為時域模型和頻域模型[17]。直接對時域的波動方程的建模產(chǎn)生了時域模型,包括有限元法,有限差分法,邊界元法;此外,考慮到聲傳播的時間尺度遠小于海洋環(huán)境參數(shù)變化的時間尺度,因此可以將時域的波動方程轉(zhuǎn)化為頻域的亥姆霍茲方程,對亥姆霍茲方程的建模產(chǎn)生了頻域模型,包括射線法,簡正波法,拋物方程法。Jensen在基準問題的求解中,比較了基于簡正波模型的雙向耦合模式和單向耦合模式的建模精度,可以認為雙向耦合模式是精確建模,在后向散射可忽略的場景下,單向耦合模式與雙向耦合模式的精度相當[18]。從計算量上來看,單向耦合模式的計算量遠小于雙向耦合模式。隨著計算能力的提升,有限元法被用于基準問題的求解,這也進一步驗證了雙向耦合模式的準確度[19]。

本文基于Hursky提出的雙水聽器時域互相關(guān)方法[10],提出了一種基于多水聽器互相關(guān)的被動聲納寬帶目標檢測方法。與用于聲源定位的雙水聽器互相關(guān)方法不同的是,本文提出的方法考慮的是如何在低信噪比下利用多個水聽器的相干性實現(xiàn)寬帶目標檢測。本文的第二部分介紹了被動聲納寬帶檢測模型,第三部分介紹了多水聽器互相關(guān)方法,第四部分通過對42個水聽器的仿真實驗驗證了本文的方法可以在極低信噪比下實現(xiàn)寬帶目標檢測,檢測能力優(yōu)于雙水聽器。此外,寬帶濾波器系數(shù)的求解需要多次求解水聲模型,這會導致計算量非常大,本文提出了利用模型的互易性快速計算寬帶濾波器系數(shù)的加速方案,極大地縮短了計算時間。

2 問題描述

考慮M個水聽器對同一個寬帶聲源的接收聲壓信號。聲源的位置為r0,聲源發(fā)射信號為s(t)。用水聽器的深度和到聲源的距離表示水聽器的位置,M個水聽器的位置分別記為r1,r2…rM,接收聲壓信號分別為p1(t),p2(t)…pM(t)。被動聲納寬帶目標檢測問題可以描述為:在聲源位置r0和聲源發(fā)射信號s(t)未知的情況下,如何利用M個水聽器的接收信號p1(t),p2(t)…pM(t)判斷有無目標,即s(t)是否不為0。

通過傅里葉合成,我們可以進一步探討M個水聽器的接收信號。考慮聲波在海洋波導中的傳播,將海洋波導當做信道,考慮聲源-海洋波導-水聽器這一信道傳輸問題。聲源在第m個水聽器處激發(fā)的聲壓場為聲源頻譜中各個單頻信號激發(fā)的聲壓場的傅里葉合成,即

(1)

其中p(r0,rm,t)表示位置為r0的聲源在位置rm處激發(fā)的聲壓場,S(ω)為s(t)的頻譜,g(r0,rm,ω)表示當角頻率為ω時從r0到rm的信道頻率響應。通常g(r0,rm,ω)用傳播損失聲壓來表示,此時S(ω)為距離聲源位置1 m處測得的頻譜。

式(1)的傅里葉變換是無法直接進行數(shù)值計算的。我們可以考慮如下幾個因素,從而將式(1)轉(zhuǎn)化為可以用實數(shù)快速傅里葉變換(Real Fast Fourier Transform, RFFT)來計算的形式:

1)實際的水聽器都是有一定的接收帶寬的,這個接收帶寬可以用最大角頻率ωmax來表示,即接收聲壓信號的能量在角頻率ωmax以上是可以忽略的。

2)聲源發(fā)射信號是實信號,因此S(ω)是共軛對稱的。

3)g(r0,rm,ω)是通過水聲模型計算得到的,水聲模型推導的起點是亥姆霍茲方程,而亥姆霍茲方程關(guān)于角頻率是共軛對稱的,因此g(r0,rm,ω)關(guān)于角頻率也是共軛對稱的。

4)對時間軸和頻率軸進行離散化。將發(fā)射信號和接收信號的起始時刻都約定為0,如果要在時長為T的時間窗中計算接收聲壓信號,那么可以對時域和頻域做如下的離散化,

tj=jΔt,j=0,1…(N-1)

(2)

(3)

其中

(4)

Δt表示時域采樣間隔,Δω表示頻域的采樣間隔。對式(4)做簡單的移項,式(4)等價于Δf=1/T,這說明頻域的物理分辨率取決于信號長度,信號時長越長,在頻域可以達到更高的物理分辨率。

按照Jensen[12]的推導,用離散求和代替式(1)中的積分,同時考慮頻率軸的離散化導致的時間軸的周期延拓,可以得到以T為周期的接收聲壓信號為,

(5)

其中Re{·}表示取實部,并且

(6)

如果T足夠長,那么式(5)中的混疊項是可以忽略的,此時,

(7)

式(7)中的求和項可以由實數(shù)快速傅里葉變換同時計算tj的N個采樣點的值。

第m個水聽器的接收信號為聲源在rm處激發(fā)的聲場與環(huán)境噪聲的疊加,即,

pm(t)=p(r0,rm,t)+n(rm,t)

(8)

其中n(rm,t)表示在rm處的環(huán)境噪聲。

當各個水聽器的間距大于二分之一波長時,可以認為各個水聽器處的環(huán)境噪聲是不相關(guān)的。

3 多水聽器互相關(guān)方法

在本文提出的多水聽器互相關(guān)方法中,首先將M個水聽器的接收信號等分為兩組,分別對兩組中的每一個水聽器的接收信號通過寬帶濾波進行信道解調(diào),然后計算兩組水聽器累加后的解調(diào)信號的互相關(guān)系數(shù),最后通過互相關(guān)系數(shù)計算互相關(guān)指標作為檢測量,可以實現(xiàn)在接收端的信噪比的提升。

3.1 信道解調(diào)

信道解調(diào)是和信道傳輸相反的過程。假設聲源的位置為r,通過寬帶濾波可以計算得到在第m個水聽器處的接收信號為pm(t)的聲源發(fā)射信號,即

(9)

其中sD(r,rm,tj)表示第m個水聽器在聲源位置r處的解調(diào)信號,Pm(ωl)表示pm(t)的頻譜。

將式(9)中的Pm(ωl)表示為接收信號分量和環(huán)境噪聲分量之和,有,

(10)

其中P(r0,rm,ωl)是p(r0,rm,t)的頻譜,nD(r,rm,tj)表示環(huán)境噪聲n(rm,t)對于聲源位置r的解調(diào)信號。

根據(jù)式(10),當r=r0時,有,

sD(r0,rm,t)=s(t)+nD(r0,rm,t)

(11)

3.2 多水聽器互相關(guān)的計算

多水聽器互相關(guān)計算流程圖如圖1所示。

圖1 多水聽器互相關(guān)計算流程圖Fig.1 Multi-hydrophone correlation calculation flowchart

(12)

根據(jù)式(11),如果假定的聲源位置與實際的聲源位置一致,即r=r0,此時

(13)

(14)

3.3 用互易性原理加速計算

在式(14)中,我們需要對搜索網(wǎng)格中每一個聲源位置計算互相關(guān)系數(shù)R(r)。假設搜索網(wǎng)格中的聲源位置的數(shù)量是K,將需要搜索的聲源位置記為r(1),r(2)…r(K),那么對于每個頻點ωl,我們總共需要知道K×M個濾波器系數(shù),表示成集合的形式為,

Gl={g(r(k),rm,ωl)|k=1…K,m=1…M}

(15)

g(r(k),rm,ωl)是通過聲場模型計算得到的。聲場模型的計算通常分為兩部分:模型求解和聲場計算。其中模型求解是指根據(jù)海洋波導參數(shù)計算得到一組能夠表示聲場的參數(shù),聲場計算是指根據(jù)模型求解得到的參數(shù)計算聲場。在本文涉及到的雙向耦合模式和單向耦合模式這兩種聲場模型中,模型求解都是計算每個距離段上的前向散射系數(shù)和后向散射系數(shù),其中雙向耦合模式需要求解全局耦合方程,單向耦合模式需要求解推進式的耦合方程;兩種模型的聲場計算是將局部的簡正波解和散射系數(shù)進行加權(quán)求和,其計算量相比于模型求解是可以忽略的。例如,我們計算了聲源位置為r(k)的單向耦合模式的各個距離段的散射系數(shù),那么我們可以通過非常少的計算量得到g(r(k),r1,ωl),g(r(k),r2,ωl)…g(r(k),rM,ωl)總共M個濾波器系數(shù)。因此我們可以將式(15)中的Gl按照k的值劃分成K個子集,如下式,

(16)

(17)

顯然有,

(18)

每個子集的計算需要一次模型求解和M次聲場計算。用TGl表示計算Gl需要的計算時間,用Tmodel表示模型求解的計算時間,用Tfield表示聲場計算的計算時間,那么有,

TGl=KTmodel+KMTfield

(19)

當搜索網(wǎng)格很大,或者網(wǎng)格尺寸很小時,K的值會很大,此時計算Gl需要很長的時間,這給有限的計算資源帶來了很大的挑戰(zhàn),因此對Gl的計算加速是很有必要的。

利用聲場模型的互易性原理,可以對Gl的計算進行加速?；ヒ仔栽硎侵?交換聲源和水聽器的位置,信道傳輸函數(shù)不變。從物理上嚴格地講,交換前后的信道傳輸函數(shù)應該是完全一致的,但是由于建模過程中引入的假設可能會導致互易性的失效,比如絕熱近似假設了各個距離段之間沒有能量的耦合傳遞,這個假設違背了界面處的聲壓連續(xù)性和質(zhì)點振速連續(xù)性,因此不滿足互易性。此外,由于數(shù)值方法的引入,交換前后的信道傳輸函數(shù)也無法做到完全一致,只能達到非常高的相似度。在后續(xù)的仿真實驗中,我們對互易性原理進行了驗證。按照互易性原理,有

g(r(k),rm,ωl)=g(rm,r(k),ωl)

(20)

將式(20)代入(15),有

Gl={g(rm,r(k),ωl)|k=1…K,m=1…M}

(21)

按照類似的推導,我們將式(21)劃分為如下的M個子集,

(22)

(23)

并且有,

(24)

此時Gl的計算時間為

TGl=MTmodel+KMTfield

(25)

對比式(25)和式(19),我們通過互易性原理的運用,將模型求解次數(shù)從K次變成了M次。M表示水聽器的個數(shù),K表示搜索網(wǎng)格的網(wǎng)格點數(shù),在實際應用中通常有M?K。同時對于聲場模型而言,通常有Tfield?Tmodel,因此我們可以認為式(25)相比于式(19)在計算量上降低了K/M倍。

4 仿真實驗與分析

4.1 海洋波導參數(shù)

仿真的海洋波導是一個隨距離變化的淺海波導,如圖2所示。

圖2 淺海波導示意圖Fig.2 Schematic diagram of shallow sea waveguide

圖2所示的海洋波導包括一個海水層和一個沉積層。其中在2500 m到3000 m之間有一個海底斜坡。海水層的密度為1000 kg/m3,沉積層的密度為2000 kg/m3。海面的均方根粗糙度為0.5 m,海底的均方根粗糙度為0.5 m。計算域的深度為500 m,底部為理想剛性海底。海水層的介質(zhì)吸收系數(shù)為每波長0 dB,沉積層的介質(zhì)吸收系數(shù)為每波長1 dB。

4.2 互易性原理的驗證

本文通過仿真得到了聲源和水聽器的位置交換前后的接收信號,計算了兩者的互相關(guān)系數(shù)來驗證單向耦合模式的互易性。聲源和水聽器的深度都是50 m。水聽器到聲源的距離分別為1000 m,2000 m,3000 m。聲源的發(fā)射信號為線性調(diào)頻信號,其波形和頻譜如圖3。

圖3 聲源的波形和頻譜Fig.3 Waveform and frequency spectrum of sound source

圖4為水聽器和聲源位置交換前后的接收信號波形。

圖4 交換前后的接收信號的波形Fig.4 Waveform of the received signal before and after the exchange

表1是距離為1000 m,2000 m,3000 m時,交換前后的波形的互相關(guān)系數(shù)。從互相關(guān)系數(shù)上來看,交換前后的接收信號是非常相似的,這說明單向耦合模式滿足互易性。

表1 互相關(guān)系數(shù)

圖5是直接計算寬帶濾波器系數(shù)和用互易性原理計算寬帶濾波器系數(shù)的計算時間的對比。聲源位置的搜索點數(shù)是40000, 水聽器數(shù)量為42。從圖5中可以看出,在30 Hz到70 Hz的頻率范圍內(nèi),直接計算濾波器系數(shù)的計算時間大約為1000 s,運用互易性原理計算濾波器系數(shù)的計算時間大約為1 s,兩個計算時間的比值與我們在式(25)和式(19)得出的計算時間的比值是相符的。

圖5 直接計算濾波器系數(shù)與運用互易性原理 計算濾波器系數(shù)的計算時間Fig.5 Calculation time for direct calculation and calculation using the principle of reciprocity of filter coefficients

4.3 多水聽器互相關(guān)方法的目標檢測仿真實驗

本文通過仿真測試了水聽器數(shù)量為2和42兩種情況下的目標檢測結(jié)果。在兩種情況下,聲源和水聽器的深度都是50 m。水聽器數(shù)量為2時(雙水聽器),兩個水聽器到聲源的距離分別為2000 m和2500 m。水聽器數(shù)量為42時(多水聽器),前21個水聽器到聲源的距離分別為2000 m,2010 m…2200 m,后21個水聽器到聲源的距離分別為2500 m,2510 m…2700 m。水聽器數(shù)量為2時,本文提出的方法等價于Hursky提出的雙水聽器互相干方法。在距離維度上,搜索范圍是-2000 m到2000 m,每5 m設置一個搜索網(wǎng)格。在深度維度上,搜索的范圍是0 m到500 m,每10 m設置一個搜索網(wǎng)格。

聲源的波形和頻譜如圖3,為線性調(diào)頻信號。通過雙向耦合模式對海洋波導中的聲傳播進行建模,得到各個水聽器的接收信號。信道解調(diào)采用單向耦合模式,并且運用互易性原理對計算進行加速。在接收信號中加入互不相關(guān)的白噪聲用于模擬不同信噪比下的接收信號。在有目標和無目標兩種場景下,分別進行了500次仿真實驗。

圖6是雙水聽器的檢測量的分布,圖7是多水聽器的檢測量的分布。在圖6和圖7中,信噪比均指距離為2000 m處的第一個水聽器的接收信號的信噪比,其中離散點線表示有目標的檢測量,實線表示無目標的檢測量。從檢測量上看,多水聽器的檢測量具有很好的可區(qū)分性,在信噪比為-21 dB時仍有一定的可區(qū)分性,在信噪比高于-15 dB時是完全可區(qū)分的,此時選擇適當?shù)拈撝?可以達到0虛警率,100%檢測率。相比于多水聽器,雙水聽器的檢測量的可區(qū)分性較低,只有在信噪比為3 dB時才有一定的可區(qū)分性,在其他更低的信噪比的情況下,無法檢測目標。

圖6 雙水聽器的檢測量分布Fig.6 The distribution of dual hydrophones

圖7 多水聽器的檢測量分布Fig.7 The distribution of mutiple hydrophones

圖8是在不同的虛警率下的檢測概率隨信噪比的變化曲線。多水聽器的處理方法在信噪比為-15 dB時檢測概率為100%,在-21 dB時仍有一定的檢測概率。雙水聽器處理方法只有在信噪比為3 dB時才具有一定的檢測能力,此外檢測概率都很低。

圖8 檢測概率曲線Fig.8 Detection probability curve

圖9是雙水聽器的檢測率和虛警率隨著閾值的變化,圖10是多水聽器的檢測率和虛警率隨著閾值的變化。檢測率和虛警率兩條曲線離得越遠,檢測性能對于閾值的選擇更加魯棒。在信噪比為-21 dB時,多水聽器的兩條曲線很接近了,這時檢測器對閾值的選擇是很敏感的,而在信噪比高于-15 dB時,兩條曲線離得很遠。對于雙水聽器的檢測率和虛警率曲線,只有在信噪比為3 dB時離得較遠。

圖9 雙水聽器的檢測率和虛警率隨著閾值的變化Fig.9 Detection rate and false alarm rate of dual hydrophone varying with threshold

圖10 多水聽器的檢測率和虛警率隨著閾值的變化Fig.10 Detection rate and false alarm rate of multiple hydrophone varying with threshold

圖11是ROC曲線。從ROC曲線上看,隨著噪聲的增加,雙水聽器的ROC曲線變成了對角線,此時沒有檢測能力,而多水聽器的ROC曲線在信噪比高于-15 dB時都是直角,在信噪比為-21 dB時仍然有一定的檢測能力,在信噪比為-27 dB時其檢測能力仍略優(yōu)于雙水聽器的檢測能力。

圖11 ROC曲線Fig.11 ROC curve

從檢測量的分布,不同虛警率下的檢測概率,檢測率和虛警率隨閾值的變化,以及ROC曲線上看,多水聽器的檢測性能均優(yōu)于雙水聽器,增加水聽器的數(shù)量可以有效地抑制噪聲對檢測結(jié)果的影響,對多個水聽器的解調(diào)信號的求和操作實現(xiàn)了信號的相干疊加,噪聲的非相干疊加,實現(xiàn)了信噪比的提升。

5 結(jié)論

本文提出的多水聽器互相關(guān)方法,通過寬帶濾波計算每個水聽器的接收信號的解調(diào)信號,將多個解調(diào)信號相加從而提升了信噪比。此外,考慮到寬帶濾波器系數(shù)的計算量很大,本文提出了運用水聲模型的互易性原理的計算加速方案,極大地降低了濾波器系數(shù)的計算量。與雙水聽器相比,多水聽器互相關(guān)方法在低信噪比情況下的檢測性能更優(yōu)。仿真實驗驗證了本文方法的優(yōu)越性。