基于淺海聲強(qiáng)干涉的雙水聽(tīng)器運(yùn)動(dòng)聲源速度估計(jì)

姚 遠(yuǎn) ,孫 超 ,謝 磊 ,劉雄厚

(1.西北工業(yè)大學(xué) 航海學(xué)院,陜西 西安,710072;2.陜西省水下信息技術(shù)重點(diǎn)實(shí)驗(yàn)室,陜西 西安,710072)

0 引言

運(yùn)動(dòng)聲源速度估計(jì)是水聲學(xué)的重要研究課題之一,在目標(biāo)定位、跟蹤以及辨識(shí)等方面有著重要的研究?jī)r(jià)值與意義。以往多采用大孔徑基陣[1-2]進(jìn)行運(yùn)動(dòng)聲源速度估計(jì),但存在布放難、造價(jià)高以及陣型失配等問(wèn)題。單水聽(tīng)器具有布放簡(jiǎn)便、成本低等優(yōu)點(diǎn),非常適合在小平臺(tái)探測(cè)模塊上使用。因此,基于單水聽(tīng)器的運(yùn)動(dòng)聲源速度估計(jì)是一個(gè)重要的研究方向。

現(xiàn)有的單水聽(tīng)器測(cè)速方法大體分為多普勒類(lèi)、聲壓互譜類(lèi)和干涉特性類(lèi)等。基于多普勒特性的測(cè)速方法主要利用線譜噪聲產(chǎn)生的多普勒頻移測(cè)速,首先提取時(shí)頻分布中線譜噪聲的瞬時(shí)頻率,再利用瞬時(shí)頻率與聲源運(yùn)動(dòng)速度間的函數(shù)關(guān)系估計(jì)出聲源速度[3-5]。基于聲壓互譜類(lèi)的測(cè)速方法主要利用運(yùn)動(dòng)聲源相對(duì)水聽(tīng)器徑向運(yùn)動(dòng)信息實(shí)現(xiàn)聲源速度估計(jì)[6-7]。第一類(lèi)方法在聲源高速運(yùn)動(dòng)且線譜頻率較高時(shí)測(cè)速結(jié)果較為準(zhǔn)確,而水下聲源在巡航工況時(shí)往往處于低航速運(yùn)動(dòng)狀態(tài);第二類(lèi)方法需要對(duì)聲源徑向運(yùn)動(dòng)速度隨時(shí)間的變化趨勢(shì)進(jìn)行長(zhǎng)時(shí)間擬合測(cè)速。

基于干涉特性的測(cè)速方法主要利用聲源激發(fā)各階簡(jiǎn)正波之間的干涉效應(yīng)實(shí)現(xiàn)聲源測(cè)速。研究人員利用運(yùn)動(dòng)聲源輻射噪聲低頻分析記錄(low frequency analysis and recording,LOFAR)譜中表現(xiàn)出的明暗相間的連續(xù)譜條紋狀干涉圖樣,開(kāi)展了一系列聲源測(cè)速研究[8-9],不過(guò)理論模型都是基于對(duì)聲場(chǎng)貢獻(xiàn)最大的前4 條聲線簡(jiǎn)化而來(lái)的。與單水聽(tīng)器相比,雙水聽(tīng)器同樣具有布放簡(jiǎn)便、成本低等優(yōu)點(diǎn),并且增加了水聽(tīng)器之間相對(duì)距離這一先驗(yàn)信息。有學(xué)者對(duì)時(shí)間頻率域干涉條紋圖進(jìn)行Hough 變換,提取出聲源運(yùn)動(dòng)參數(shù),通過(guò)構(gòu)建雙水聽(tīng)器的距離與聲源運(yùn)動(dòng)參數(shù)之間的空間關(guān)系,實(shí)現(xiàn)了聲源測(cè)速[10-12]。然而,在一些場(chǎng)景下,運(yùn)動(dòng)聲源輻射噪聲LOFAR 譜中寬帶連續(xù)譜干涉條紋非常微弱,上述測(cè)速方法難以應(yīng)用。而由于聲源輻射噪聲中線譜成分比連續(xù)譜成分平均高10～20 dB[13],相較連續(xù)譜而言,能量高的線譜更易觀測(cè)。

已有研究主要是將輻射噪聲線譜聲強(qiáng)干涉起伏看作是對(duì)距離頻率域干涉譜在特定頻率上的調(diào)制,即它的起伏對(duì)應(yīng)干涉譜條紋的亮暗,再通過(guò)構(gòu)建兩頻率線譜聲強(qiáng)干涉起伏相似度代價(jià)函數(shù),實(shí)現(xiàn)聲源距離和波導(dǎo)不變量參數(shù)估計(jì)[14-16]。目前尚未見(jiàn)到有關(guān)利用線譜聲強(qiáng)干涉起伏實(shí)現(xiàn)聲源運(yùn)動(dòng)速度估計(jì)方面的報(bào)道。文中基于雙水聽(tīng)器利用線譜聲強(qiáng)干涉起伏分別獲得相對(duì)1 號(hào)水聽(tīng)器和2 號(hào)水聽(tīng)器的聲源運(yùn)動(dòng)參數(shù),再結(jié)合聲源與雙水聽(tīng)器之間的距離位置關(guān)系,解算出聲源運(yùn)動(dòng)速度。其中正確的聲源運(yùn)動(dòng)參數(shù)是時(shí)間維聲強(qiáng)經(jīng)過(guò)插值變換后的不同頻率線譜聲強(qiáng)干涉起伏滿(mǎn)足正比關(guān)系的必要條件。因此,通過(guò)在一定區(qū)間范圍內(nèi)進(jìn)行參數(shù)搜索,構(gòu)建不同頻率線譜聲強(qiáng)干涉起伏間的相關(guān)系數(shù)代價(jià)函數(shù),獲得聲源最近通過(guò)時(shí)間和距速比參數(shù),結(jié)合聲源與雙水聽(tīng)器距離位置關(guān)系,可以實(shí)現(xiàn)聲源速度估計(jì)。文中方法不受聲源運(yùn)動(dòng)速度約束,適合在小平臺(tái)探測(cè)模塊上使用,對(duì)背景噪聲下低速運(yùn)動(dòng)聲源有較好的測(cè)速效果。

1 聲源相對(duì)水聽(tīng)器運(yùn)動(dòng)模型

勻速直線運(yùn)動(dòng)聲源相對(duì)水聽(tīng)器的運(yùn)動(dòng)示意圖如圖1 所示。聲源運(yùn)動(dòng)到最近通過(guò)位置(the closest position of approach,CPA)處的時(shí)刻稱(chēng)為最近通過(guò)時(shí)間,記為tCPA,CPA 處的聲源與水聽(tīng)器之間的距離稱(chēng)為最近通過(guò)距離,記為rCPA。

圖1 聲源相對(duì)水聽(tīng)器的運(yùn)動(dòng)示意圖Fig.1 Illustration of relative position between hydrophone and moving source

在t時(shí)刻運(yùn)動(dòng)聲源與水聽(tīng)器間的位置關(guān)系滿(mǎn)足

式中,v0表示聲源運(yùn)動(dòng)速度。

當(dāng)聲源依次通過(guò)雙水聽(tīng)器,相對(duì)1 號(hào)水聽(tīng)器和2 號(hào)水聽(tīng)器的最近通過(guò)時(shí)間分別記為tCPA1和tCPA2,相對(duì)1 號(hào)水聽(tīng)器和2 號(hào)水聽(tīng)器的最近通過(guò)距離分別記為rCPA1和rCPA2。聲源運(yùn)動(dòng)軌跡與雙水聽(tīng)器之間的位置關(guān)系如圖2 所示。

圖2 聲源運(yùn)動(dòng)軌跡與雙水聽(tīng)器之間的位置關(guān)系示意圖Fig.2 Illustration of relative position between dual hydrophones and moving source trajectory

雙水聽(tīng)器的水平間距可以表示為

進(jìn)一步化簡(jiǎn),有

式中,b1=rCPA1/v0,b2=rCPA2/v0,分別表示聲源分別相對(duì)1 號(hào)和2 號(hào)水聽(tīng)器的最近通過(guò)距離與速度之比,簡(jiǎn)稱(chēng)為距速比?？梢?jiàn),只要估計(jì)出聲源相對(duì)雙水聽(tīng)器的聲源運(yùn)動(dòng)參數(shù)tCPA1、tCPA2以及b1和b2,便能計(jì)算出聲源運(yùn)動(dòng)速度v0。下文將給出利用線譜聲強(qiáng)干涉起伏估計(jì)聲源運(yùn)動(dòng)參數(shù)具體過(guò)程。

2 淺海波導(dǎo)聲源速度估計(jì)

2.1 接收信號(hào)模型

以自身旋轉(zhuǎn)部件的機(jī)械噪聲和螺旋槳推進(jìn)噪聲為主的噪聲源可建模為K個(gè)低頻線譜分量的疊加[16]

式中:k為線譜分量的號(hào)數(shù);Ak和 ωk分別為第k個(gè)線譜的幅度和角頻率;?k為線譜的相位。

由于聲源特性以及傳播過(guò)程引起的幅度起伏,一般聲源線譜幅度具有起伏性。文中主要關(guān)注聲源與接收器傳播距離變化引起的幅度起伏,暫不考慮聲源起伏,即假設(shè)聲源在線譜頻率處輻射能量穩(wěn)定。

根據(jù)簡(jiǎn)正波理論,波導(dǎo)中位于r0=(0,zs)的遠(yuǎn)場(chǎng)聲源到接收器(r,zr)的信道傳遞函數(shù)G(r,ω)可表示為[17]

式中:krm為第m階簡(jiǎn)正波的水平波數(shù);Ψm為第m階簡(jiǎn)正波的模態(tài)深度函數(shù);Bm為聲壓的幅度項(xiàng)。

水聽(tīng)器接收信號(hào)頻譜表示為

式中,S(ω)為聲源s(t)的頻譜。

2.2 線譜聲強(qiáng)干涉起伏

定義第k個(gè)線譜聲強(qiáng)為[17]

式中,Δkmn為第m階與第n階簡(jiǎn)正波的水平波數(shù)差,即Δkmn=krm-krn。

式(7)省略了隨著距離和頻率慢變的非相干項(xiàng),保留了簡(jiǎn)正波兩兩干涉疊加的相干項(xiàng),該項(xiàng)決定了線譜聲強(qiáng)干涉起伏[18]。通常在多號(hào)簡(jiǎn)正波頻散特性非常一致的淺海波導(dǎo)中,波導(dǎo)不變量與相互干涉的簡(jiǎn)正波號(hào)數(shù)m,n無(wú)關(guān),根據(jù)文獻(xiàn)[17],寬帶連續(xù)譜輻射噪聲的接收信號(hào)聲強(qiáng)在距離-頻率平面上的恒定強(qiáng)度條紋軌跡滿(mǎn)足

式中:r0和 ω0分別表示參考距離和參考頻率;β表示波導(dǎo)不變量。軌跡上任意兩頻率 ω1和ω2處聲強(qiáng)I1(r1)和I2(r2)滿(mǎn)足如下關(guān)系[19]

式中:A1和A2分別為頻率 ω1和ω2處聲源線譜頻譜幅度;r2=r1(ω2/ω1)1/β,r1和r2為軌跡上兩頻率點(diǎn)聲強(qiáng)對(duì)應(yīng)的距離。

由圖1 可知,由于運(yùn)動(dòng)聲源靠近CPA 的過(guò)程中,聲源與水聽(tīng)器之間的距離r相對(duì)時(shí)間t是非線性變化的。水聽(tīng)器接收信號(hào)聲強(qiáng)干涉起伏Ik在時(shí)間維度t是均勻采樣的,而對(duì)應(yīng)的距離維度是非均勻的,其不同頻率之間的線譜聲強(qiáng)干涉起伏I1和I2滿(mǎn)足式(9)。對(duì)于隨時(shí)間變化的線譜聲強(qiáng)干涉起伏,可以通過(guò)對(duì)式(1)兩端同時(shí)除以v0得到,即

依據(jù)式(10)對(duì)t時(shí)間域進(jìn)行重采樣,并對(duì)聲強(qiáng)Ik(t)進(jìn)行插值變換,得到在t′時(shí)間域均勻采樣的聲強(qiáng)(t′),從而滿(mǎn)足類(lèi)似式(9)的關(guān)系。根據(jù)式(9)和式(10),重采樣后時(shí)間域的不同頻率線譜聲強(qiáng)干涉起伏之間滿(mǎn)足[20]

2.3 聲源運(yùn)動(dòng)速度估計(jì)

綜上,利用雙水聽(tīng)器獲得的聲源運(yùn)動(dòng)參數(shù)估計(jì)結(jié)果及聲源與雙水聽(tīng)器之間的距離位置關(guān)系,從而實(shí)現(xiàn)聲源速度估計(jì),基本處理流程如圖3 所示。

圖3 聲源速度估計(jì)的基本處理流程Fig.3 Basic flowchart of acoustic source velocity estimation

3 數(shù)值實(shí)驗(yàn)與分析

文中通過(guò)數(shù)值仿真驗(yàn)證基于雙水聽(tīng)器利用多線譜聲強(qiáng)干涉起伏實(shí)現(xiàn)聲源速度估計(jì)方法的有效性。仿真使用波導(dǎo)環(huán)境為淺海等聲速波導(dǎo)環(huán)境,水深為40 m,水體聲速為1 500 m/s,半空間海底聲速為1 750 m/s,密度為1 750 kg/m3,吸收系數(shù)為0.2 dB/λ,λ 表示波長(zhǎng)。波導(dǎo)環(huán)境參數(shù)如圖4 所示。數(shù)值仿真利用Kraken 聲場(chǎng)仿真軟件。

圖4 波導(dǎo)環(huán)境參數(shù)Fig.4 Waveguide environment parameters

3.1 數(shù)值仿真

雙水聽(tīng)器水平間距為1 000 m。聲源以速度2 m/s 作勻速直線運(yùn)動(dòng),運(yùn)動(dòng)時(shí)間區(qū)間為0～1 800 s,采樣間隔為1 s。聲源在第1 號(hào)水聽(tīng)器CPA 處的時(shí)刻為950 s,該時(shí)刻與第1 號(hào)水聽(tīng)器的最近通過(guò)距離為800 m;聲源在第2 號(hào)水聽(tīng)器CPA 處的時(shí)刻為1 350 s,該時(shí)刻與第2 號(hào)水聽(tīng)器的最近通過(guò)距離為1 400 m。聲源線譜頻率分別為300 Hz 和350 Hz。聲源深度為10 m,水聽(tīng)器深度為35 m。波導(dǎo)以反射類(lèi)簡(jiǎn)正波為主,波導(dǎo)不變量視作1。

第1 號(hào)水聽(tīng)器接收信號(hào)的信噪比(signal-tonoise ratio,SNR)設(shè)為0 dB,由于運(yùn)動(dòng)聲源輻射噪聲LOFAR 譜中線譜成分相較連續(xù)譜成分能量要高,這里將300 Hz 和350 Hz 線譜聲強(qiáng)提高15 dB作為已知線譜聲強(qiáng)起伏,水聽(tīng)器接收信號(hào)的LOFAR譜如圖5 所示。SNR 定義為帶寬內(nèi)信號(hào)功率與高斯白噪聲功率的比值。

圖5 第1 號(hào)水聽(tīng)器接收信號(hào)的LOFAR 譜Fig.5 LOFAR spectrum of the receiving signal with hydrophone No.1

給出時(shí)間區(qū)間[100,950]的頻率300 Hz 和350 Hz 線譜聲強(qiáng)干涉起伏如圖6 所示。其中,左側(cè)圖中兩頻率線譜聲強(qiáng)的干涉起伏,通過(guò)式(12)對(duì)時(shí)間進(jìn)行重采樣且對(duì)線譜聲強(qiáng)進(jìn)行插值變換,變換后聲強(qiáng)隨時(shí)間的變化干涉起伏如圖6 右側(cè)圖所示,圖中黑色實(shí)線表示兩頻率線譜聲強(qiáng)峰值對(duì)應(yīng)位置。可以看出,變換后的兩線譜聲強(qiáng)干涉起伏滿(mǎn)足式(11)對(duì)應(yīng)關(guān)系,即時(shí)間變量t′在時(shí)間區(qū)間對(duì)應(yīng)頻率 ω1線譜聲強(qiáng)干涉起伏與經(jīng)過(guò)尺度伸縮的時(shí)間區(qū)間內(nèi) ω2線譜聲強(qiáng)干涉起伏滿(mǎn)足正比關(guān)系,此時(shí)和兩者之間的相關(guān)系數(shù)達(dá)到最大。

圖6 第1 號(hào)水聽(tīng)器100～950 s 時(shí)間區(qū)間內(nèi)300 Hz 和350 Hz線譜聲強(qiáng)干涉起伏Fig.6 Interference fluctuation of 300 Hz and 350 Hz acoustic intensity with hydrophone No.1 from 100 s to 950 s

若距速比和最近通過(guò)時(shí)間搜索值bs和tCPAs偏離真實(shí)值,則相關(guān)系數(shù)減小。在不同bs和tCPAs假設(shè)值條件下,分別對(duì)300 Hz 和350 Hz 線譜聲強(qiáng)在時(shí)間區(qū)間[350,900]重采樣,得到插值變換后的線譜聲強(qiáng)干涉起伏,兩者間的相關(guān)系數(shù)代價(jià)函數(shù)如圖7所示。由圖可看出,代價(jià)函數(shù)最大值對(duì)應(yīng)最近通過(guò)時(shí)間和距速比估計(jì)結(jié)果分別為961 和408。

圖7 根據(jù)第1 號(hào)水聽(tīng)器頻率300 Hz 和350 Hz 聲強(qiáng)干涉起伏計(jì)算得到代價(jià)函數(shù)Fig.7 Cost functions of acoustic intensity interference fluctuation of hydrophone No.1 at 300 Hz and 350 Hz

時(shí)間區(qū)間從[200,750]到[350,900],各區(qū)間間隔為10 s。給出不同時(shí)間區(qū)間的最近通過(guò)時(shí)間和距速比估計(jì)結(jié)果,如圖8 所示。紅色方框?yàn)檎鎸?shí)值,藍(lán)色星號(hào)表示各參數(shù)估計(jì)結(jié)果,可見(jiàn),其穩(wěn)定分布在真實(shí)值附近。通過(guò)K-means 聚類(lèi)算法求取各區(qū)間參數(shù)估計(jì)結(jié)果的質(zhì)心[22],如紅色星號(hào)所示,對(duì)應(yīng)的最近通過(guò)時(shí)間=956 和距速比=407 作為第1 號(hào)水聽(tīng)器的聲源運(yùn)動(dòng)參數(shù)估計(jì)結(jié)果。

圖8 根據(jù)第1 號(hào)水聽(tīng)器頻率300 Hz 和350 Hz 聲強(qiáng)干涉起伏計(jì)算得到不同時(shí)間區(qū)間的參數(shù)估計(jì)結(jié)果Fig.8 Motion parameter estimation results in different time intervals which obtained by acoustic intensity interference fluctuation of hydrophone No.1 at 300 Hzand 350 Hz

同理,時(shí)間區(qū)間從[500,1 150]到[650,1 300],各區(qū)間間隔為10 s,第2 號(hào)水聽(tīng)器不同時(shí)間區(qū)間的最近通過(guò)時(shí)間和距速比參數(shù)估計(jì)結(jié)果如圖9 所示。由圖可知,估計(jì)結(jié)果穩(wěn)定分布在真實(shí)值附近。質(zhì)心值對(duì)應(yīng)的最近通過(guò)時(shí)間=1 345 和距速比=708 作為第2 號(hào)水聽(tīng)器的聲源運(yùn)動(dòng)參數(shù)估計(jì)結(jié)果。

圖9 根據(jù)第2 號(hào)水聽(tīng)器頻率300 Hz 和350 Hz 聲強(qiáng)干涉起伏計(jì)算得到不同時(shí)間區(qū)間的參數(shù)估計(jì)結(jié)果Fig.9 Motion parameter estimation results in different time intervals which obtained by acoustic intensity interference fluctuation of hydrophone No.2 at 300 Hz and 350 Hz

根據(jù)式(14),利用雙水聽(tīng)器的聲源運(yùn)動(dòng)參數(shù)估計(jì)結(jié)果計(jì)算出速度為2.03 m/s。圖10 是在不同聲源速度條件下,利用所提方法實(shí)現(xiàn)的聲源速度估計(jì)。可以看出,所提方法能有效實(shí)現(xiàn)不同速度聲源的速度估計(jì)。存在誤差的原因一方面是噪聲的存在對(duì)估計(jì)結(jié)果產(chǎn)生一定影響;另一方面,文中根據(jù)反射類(lèi)波導(dǎo)條件,在常規(guī)意義上將波導(dǎo)不變量看作經(jīng)驗(yàn)值1,實(shí)際上波導(dǎo)不變量并不嚴(yán)格等于1,這導(dǎo)致估計(jì)得到的距速比和最近通過(guò)時(shí)間存在偏差,影響聲源速度估計(jì)結(jié)果。

圖10 不同聲源速度估計(jì)結(jié)果Fig.10 Estimation results of velocity with different acoustic sources

關(guān)于雙水聽(tīng)器的布放間距,這里分別以聲源平行和垂直于雙水聽(tīng)器的運(yùn)動(dòng)場(chǎng)景進(jìn)行說(shuō)明。通常情況下,聲源速度約在1～5 m/s 之間。對(duì)于聲源平行于雙水聽(tīng)器運(yùn)動(dòng)情況:若雙水聽(tīng)器布放距離為100 m,那么最近通過(guò)時(shí)間差為20～100 s。由圖8 和圖9 可看出,最近通過(guò)時(shí)間估計(jì)值存在10 s 量級(jí)誤差,此時(shí)由式(14)計(jì)算得到聲源速度估計(jì)的誤差最高將達(dá)到50%;而若雙水聽(tīng)器布放距離為1 000 m,那么最近通過(guò)時(shí)間差為200～1 000 s,此時(shí)聲源速度估計(jì)的誤差最高為5%。

同理,對(duì)于聲源垂直于雙水聽(tīng)器運(yùn)動(dòng)情況:根據(jù)式(14)可以看出,雙水聽(tīng)器的最近通過(guò)時(shí)間相消,此時(shí)速度估計(jì)結(jié)果由距速比決定。由圖8 和圖9可以看出,距速比估計(jì)結(jié)果存在大概10 s 量級(jí)誤差。隨著雙水聽(tīng)器布放距離增大,距速比之差增大,距速比誤差量級(jí)一定,聲源速度估計(jì)誤差將降低。因此在滿(mǎn)足波導(dǎo)環(huán)境平坦條件前提下,盡量保證雙水聽(tīng)器間距遠(yuǎn)距離布放。該結(jié)論對(duì)于聲源非平行和非垂直于雙水聽(tīng)器的運(yùn)動(dòng)場(chǎng)景同樣適用。

4 結(jié)束語(yǔ)

針對(duì)低SNR 時(shí)難以利用LOFAR 譜中寬帶連續(xù)譜條紋狀干涉圖樣估計(jì)聲源運(yùn)動(dòng)速度的問(wèn)題,提出一種基于雙水聽(tīng)器利用線譜聲強(qiáng)干涉起伏估計(jì)聲源速度方法。利用聲源與雙水聽(tīng)器的距離位置關(guān)系,通過(guò)不同頻率線譜聲強(qiáng)干涉起伏獲得最近通過(guò)時(shí)間和距速比參數(shù),實(shí)現(xiàn)淺海反射類(lèi)簡(jiǎn)正波主導(dǎo)的波導(dǎo)環(huán)境下的低速聲源運(yùn)動(dòng)速度估計(jì)。仿真實(shí)驗(yàn)結(jié)果表明,所提方法對(duì)背景噪聲下低速運(yùn)動(dòng)聲源有較好的測(cè)速效果,不需要采用大孔徑基陣,非常適合在小平臺(tái)探測(cè)模塊上使用。

未來(lái)研究中可以考慮在夏季典型負(fù)躍層聲速剖面時(shí),聲場(chǎng)可能由反射類(lèi)簡(jiǎn)正波主導(dǎo)變?yōu)橛烧凵漕?lèi)簡(jiǎn)正波主導(dǎo),此時(shí)所提方法應(yīng)用性將受到限制。如何在折射類(lèi)簡(jiǎn)正波主導(dǎo)情況下實(shí)現(xiàn)聲源速度等參數(shù)估計(jì)將是下一步要開(kāi)展的工作。