采用單水聽器的勻速圓弧運動直升機三維參數(shù)估計算法

張華霞 孫偉濤 王惠剛*② 榮少巍

①(西北工業(yè)大學(xué)航海學(xué)院 西安 710072)

②(西北工業(yè)大學(xué)深圳研究院 深圳 518057)

1 引言

跨介質(zhì)場景中運動聲源的參數(shù)估計是目標定位領(lǐng)域的重要問題[1]。直升機線頻輻射噪聲透射入水，折射波具有余弦平方方向性，水下聲壓場隨水平距離增大呈指數(shù)衰減趨勢，信噪比也隨之迅速降低，難以獲得穩(wěn)定的長時間多普勒頻移信號[2,3]。同時，水聽器接收的多普勒頻率與動聲源位置密切相關(guān)，做非直線運動的直升機會進一步增加反演聲源運動參數(shù)的復(fù)雜度[4]。因此，如何利用局部多普勒頻移信息反演做圓弧運動聲源飛行參數(shù)具有重要意義。

目前現(xiàn)有理論模型大多研究估計直升機沿直線航行的飛行參數(shù)。例如在2維域方面，即空中目標過頂直線飛行模型，已有大量文獻分析了空氣-水兩種介質(zhì)的聲傳播模型，并推導(dǎo)出該場景聲源運動參數(shù)的方法[5–7]。在3維域方面，一種方法是使用單矢量水聽器和恰當數(shù)學(xué)模型來估計目標方位和其他飛行參數(shù)[4]；另一種方法是采用寬孔徑水聽器陣列獲取空中目標的方位信息，對空氣中飛行的直升機進行參數(shù)估計和定位[8,9]。文獻[10]中已對勻速直線飛行的參數(shù)估計方法和結(jié)果進行總結(jié)。然而，相比直升機直線航行，沿圓弧航行時直升機與聲源之間距離函數(shù)更為復(fù)雜，導(dǎo)致以上算法在估計圓弧航行的飛行參數(shù)時會失效。

針對勻速圓弧運動的3維場景，本文提出單水聽器估計直升機勻速圓弧運動參數(shù)的算法，可實現(xiàn)固有頻率、飛行高度、速度和圓弧半徑的同步解算，從而解決了單水聽器對勻速圓弧運動聲源定位信息以及狀態(tài)參數(shù)的估計問題。首先，在兩種等聲速介質(zhì)的3維空間模型中，本文建立了水聽器接收動聲源聲波的幾何模型。然后根據(jù)聲源與水聽器相對運動的幾何規(guī)律，推導(dǎo)出多普勒頻移與聲源頻率、速度、高度之間的函數(shù)關(guān)系。最后，應(yīng)用聲源和水聽器的幾何關(guān)系、多普勒頻移曲線，可以估計出聲源的運動參數(shù)。本文通過仿真勻速圓弧運動聲源的多普勒信號，驗證了本方法估計直升機的固有頻率、速度、高度等航行參數(shù)性能。

2 勻速圓弧飛行的多普勒頻移模型

在空氣-水介質(zhì)環(huán)境中，如圖1所示，直升機聲源在水聽器節(jié)點上空以O(shè)為圓心、r為半徑做勻速圓弧運動。靜止的水聽器節(jié)點記為點H，位于水下深度d處。恒定頻率為f0的直升機輻射聲源記為點S，它以恒定亞聲速v(vcacw)恒定高度h飛過點H上方附近。投射到靜止節(jié)點H所在的水平面上動聲源記為S′，該投影動點S′形成的軌跡記為O′。動點S′和 節(jié)點H之間距離記為ρ(t)， 即點聲源S和節(jié)點H直線距離r(t)的 水平投影，ρ(t)隨 動點S′的位置改變而變化。入射角為θI(t)的聲源發(fā)出的聲線在點T處發(fā)生折射進入水下，到達水聽器節(jié)點H。

圖1 空中沿曲線運動的點聲源與靜止水聽器節(jié)點的3維圖

根據(jù)文獻[10]中對勻速直線運動直升機飛行參數(shù)估計方法，水聽器所接收到的多普勒頻移曲線fd(t) 由 參數(shù){f0,v,β(t),α(t)}決定。對于固有頻率為f0的 動點聲源S，水聽器在t時刻所接收到的多普勒頻率為

其中，偏向角α(t)是 聲源S在圓弧處切線與直線S′H的夾角，空氣中的視線角β(t) 與 入射角θI(t)互為余角，即β(t)+θI(t)=π/2。

式(1)給出多普勒頻移fd(t)與聲源運動參數(shù)之間的關(guān)系。由于難以直接求解偏向角α(t)與空氣視線角β(t) 隨 時間t的變化關(guān)系，因此需要將該角度函數(shù)轉(zhuǎn)化為聲源S與水聽器節(jié)點H之間位置關(guān)系，以方便求解多普勒頻移。

首先，根據(jù)圖2所示的俯視圖中動聲源S與水聽器節(jié)點H的位置關(guān)系，可以求解偏向角α(t)。在時間t時刻，聲源S與節(jié)點H的 坐標分別為(xS(t),yS(t)),(xH,yH)， 因此偏向角可由聲源S與水聽器節(jié)點H的坐標位置求解

圖2 空中動點聲源與靜止水聽器節(jié)點的俯視圖

然后，根據(jù)圖3所示的聲線傳播路徑幾何關(guān)系，可以求解空氣中視線角β(t) 。 動點聲源S和靜止節(jié)點H之間的側(cè)視圖，展示了入射角為θI(t)和折射角為θT(t)的聲線傳播路徑關(guān)系。

圖3 聲波傳播路徑的側(cè)視圖

在 以∠STH為頂角△STH中，由 余 弦定理可 知入射角θI(t)和 折射θT(t)角滿足關(guān)系為

其中，聲波在空氣介質(zhì)與水介質(zhì)的傳播距離分別為ra(t)和rw(t)，可寫為

點聲源S和水聽器節(jié)點H的直線距離r(t)為

入射角θI(t)和 折射角θT(t)之間滿足的斯涅爾(Snell)定理，且與空氣視線角β(t)的互補關(guān)系和水下視線角的互補關(guān)系可表示為

將式(6)代入式(3)后，整理可得水下視線角β(t)與點聲源高度h、聲源速度v、聲傳播路徑投影ρ(t) 、水聽器深度d之間的隱函數(shù)方程

其中，水下視線角β(t)與 折射角θT(t)互為余角，即βˉ(t)+θT(t)=π/2。

最后，將式(2)和式(7)代入式(1)中，化簡后可以得到方程組

由于聲源S在水平方向沿曲線航跡運動，靜止水聽器H觀測到聲波，水聽器H觀測聲源S的偏向角α(t) 與 空 氣 視 線 角β(t) 會 隨 聲 源S與 水 聽 器H相對位置改變而變化。在圖2所示坐標系中，設(shè)定聲源S與水聽器H之間距離變小的階段，速度v為正；設(shè)定聲源S與水聽器H之間距離變大的階段，速度v為負。

由式(2)可知，聲源S在運動過程中，偏向角α(t)變 化范圍為[ 0,π]rad 。其中，聲源S朝向水聽器H運動時，偏向角α(t)∈[0,π/2]rad ，遠離水聽器H運動時，偏向角α(t)∈[π/2,π]rad 。特別地，若聲源S相距水聽器H無限遠，偏向角為α(t)=0 rad 或α(t)=π rad ；速度v的瞬時方向垂直于聲源S與水聽器H連線時，偏向角α(t)為 π /2 rad 。由式(7)可知，空氣視線角β(t)變化范圍為[0,arctanh/ρmin(t)]rad 。其中，若聲源S相距水聽器H無限遠，空氣視線角為β(t)=0 rad，聲源S經(jīng)至水聽器H最近點時，空氣視線角β(t)為arctanh/ρmin(t)rad 。聲源S朝向水聽器H運動過程中，觀測到的多普勒頻率大于聲源S的聲波頻率，即fd(t)>f0(t)， 聲源S遠離水聽器H運動過程中，觀測到的多普勒頻率小于聲源S的聲波頻率，即fd(t)f0(t)， 速度v的瞬時方向垂直于聲源S與水聽器H連線時，觀測到的多普勒頻率等于聲源S的聲波頻率，即fd(t)=f0(t)。

3 估計飛行參數(shù)

若已知空氣介質(zhì)聲速ca、 水介質(zhì)聲速cw、水聽器深度d，本節(jié)使用單水聽器所接收的水聲數(shù)據(jù)來估計多普勒頻移曲線fd(t)，并根據(jù)式(8)、聲源與水聽器幾何關(guān)系可以估計直升機的固有頻率f0、速度v、飛行高度h。根據(jù)幾何常識，確定一段圓弧需要選取3個不共線點處的瞬態(tài)頻率來估計運動參數(shù)。

圖4展示了圓弧上選取3個點的幾何關(guān)系俯視圖，聲源S在時刻tA,tM和tC依次經(jīng)過不共線的3點A,M,C， 其中點M為弧線lAC的 中點，即lAM=lMC。由于聲源S在點A發(fā)出的聲波傳播至水聽器H具有一定時間差 ?tAB，水聽器H接收到聲波信號時，聲源S已運動至點B處。若時間差?tAC大于聲波傳至水聽器時間，則點B位于弧lAC上，如圖4(a)所示；若時間差 ?tAC小于聲波傳至水聽器時間，則點B位于弧lAC外，如圖4(b)所示。

在圖4所示坐標系中，聲源經(jīng)過點{A,M,C}是等間隔的，令點{A,M,C}的坐標分別為

圖4 圓弧上選取3個測量點的幾何關(guān)系圖

這里

為求解聲源飛行參數(shù)信息，本節(jié)需確定輔助點B的信息，其為聲源由點A處航行?tAB時間后到達的位置，然后根據(jù)點{A,B,C,M}的幾何關(guān)系來求解參數(shù)。

首先，令點B坐標為

這里

其中，s gn是符號函數(shù)。

聲源位于點A時，聲線由聲源S傳至水聽器節(jié)點H的傳播時間為

根據(jù)弧長公式，聲源S一定時間內(nèi)經(jīng)過的弧長與其弧度有關(guān)系為

其中，聲源到達點B的時刻為

然后，根據(jù)點{A,B,C,M}的位置關(guān)系，建立幾何關(guān)系。以圓心點O為頂點、選取圓弧上的不同的2點連線PQ(即{P,Q|P ?=Q}?{A,B,C,M})為底邊所構(gòu)成的三角形應(yīng)用余弦定理，有

對于∠POQ，有如式(17)的關(guān)系式

以節(jié)點投影H′為頂點、圓弧上的任2 點{P,Q|P ?=Q}?{A,B,C,M}為底邊所構(gòu)成的三角形應(yīng)用余弦定理，有

∠PH′Q

對于 ，有如式(19)的關(guān)系式

對于ρP，有如式(20)的關(guān)系式

以節(jié)點投影O為頂點、H′P為底邊所構(gòu)成的三角形應(yīng)用余弦定理，有

聲源位于P ∈{A,B,C,M}時，以折射點TP為頂點、聲源S與水聽器H連線為底邊所構(gòu)成的三角形，其幾何關(guān)系如圖3所示，應(yīng)用余弦定理可以得到

根據(jù)弧長公式，聲源S一定時間內(nèi)經(jīng)過的弧長與其弧度有如式(23)的關(guān)系

P ∈{A,B,C,M}

聲源位于 時，其瞬時多普勒頻率為

4 實驗與分析

4.1 多普勒頻移特征

3維空間內(nèi)，假定聲波在空氣和水中的速度為1500 m/s和340 m/s。如圖5所示，水聽器靜止于(0,0,?100)m處 ，螺旋槳基頻頻率f0=68 Hz直升機分別在高度150 m和100 m處的水平面以v=123 m/s沿逆時針圓弧航行，聲源在t時刻3維坐標為

圖5 仿真選取的5組測量點與靜止水聽器幾何示意圖

考慮到時間間隔、飛行朝向以及飛行高度，仿真中分別取了5組觀測點，以說明選取聲源坐標對飛行參數(shù)估計精度的影響，其中第1,2組高度為h1=150 m、時刻間隔為0.5 s，第3,4組高度為150 m、時刻間隔為0.2 s；第5組高度為100 m、時刻間隔為0.5 s。其中，第1,3,5組聲源朝向水聽器運動，第2,4組聲源遠離水聽器運動。

由式(25)可知，任意時刻聲源距離水聽器直線的投影ρ(t)為

將聲源運動參數(shù){f0,v,d,h,ρ(t)}代入式(8)后，可以得到聲源做勻速圓周圓弧的多普勒頻率曲線fd(t)，如圖6所示，展示了2個不同高度下(50 m, 150 m)的多普勒偏移曲線。聲源做勻速圓周運動時，多普勒頻移曲線也會呈現(xiàn)周期性。

4.2 估計多普勒頻移

聲源沿圓周做圓弧運動時，根據(jù)單水聽器采集的多普勒頻移信號及其深度信息，本節(jié)使用所提出的水下探空算法解算空中勻速圓周航行直升機的全部飛行參數(shù)。

由圖6所示的多普勒頻移曲線fd(t)可以得到時域信號，即

設(shè)置采樣范圍t ∈[1,6]s ，采樣頻率fs=500 Hz 。本文使用加性Alpha穩(wěn)定噪聲模擬包含大量沖擊的水下噪聲場。根據(jù)式(27)得到的時域信號加入α=1.6的Alpha穩(wěn)定噪聲后，得到如圖7所示的時域圖，其廣義信噪比為0 dB。

圖7 水聽器接收的時域信號

估計直升機飛行參數(shù)前，需獲取實測信號的時頻估計曲線。由文獻[10]中的3種估計多普勒頻移曲線方法可以看出，以短時傅里葉變換算法為基礎(chǔ)的傳統(tǒng)時頻分析方法受到時間帶寬積制約，時間與頻率分辨率不會同時達到最佳，這會導(dǎo)致估計的頻率曲線產(chǎn)生一定理論誤差。因此本文采用了基于多重Log-sum累積相位功率(Accumulated Phase-difference Power on Multiple Logarithm Sum, APPMLS)的時頻估計算法[11,12]。該方法在時延估計、壓縮感知等方面都有廣泛應(yīng)用[13–15]。

圖8顯示了基于多重Log-sum累積相位功率針對圖6中在兩種不同高度下的理想頻率變化曲線所估計的多普勒頻移曲線，信噪比為0 dB，并與理想的曲線對比。將該音頻文件每次取125個樣本點作為1幀數(shù)據(jù)，相鄰兩幀重疊10個樣本點。圖8中多普勒頻移曲線具有較高的頻率分辨率，并且變化趨勢較為平滑，未受到噪聲影響而發(fā)生明顯波動。此外，該時頻曲線已經(jīng)有效抑制噪聲，可以直接據(jù)此做參數(shù)估計，無需做進一步濾波處理。

圖6 勻速圓弧運動聲源的多普勒頻移曲線

4.3 估計3維參數(shù)

根據(jù)第3節(jié)所述，根據(jù)式(9)、式(10)設(shè)定3點坐標，并利用式(11)—式(15)設(shè)定輔助點B的位置，并聯(lián)立式(16)、式(18)、式(21)和式(22)所示的余弦定理，然后將選定直升機運行軌跡中等間隔3點{A,M,C}和輔助點B的時刻和多普勒頻率代入式(24)。聯(lián)立該非線性方程組，即可以求解出直升機航行參數(shù)。

在圖8中，選取3個等間隔時刻及其多普勒頻移信息，即 (tA,f(tA)),(tM,f(tM))和 (tC,f(tC))。表1列出了5組數(shù)據(jù)信息以及相應(yīng)的直升機飛行參數(shù)估計結(jié)果。從表1可以看出，由于前兩組中相鄰時刻間隔比后兩組長，時頻變化特征相對較明顯，因此第1,2,5組的估計結(jié)果較為準確。圖5(b)顯示了直升機運動軌跡與靜止水聽器的位置關(guān)系，由于水聽器在非圓心位置，第1,3組的直升機位置比第2, 4組距離水聽器較近，因此第1,2組的估計結(jié)果偏差相對較??；而第5組直升機高度較低，多普勒頻移特征明顯，估計結(jié)果偏差相對第1組較小。整體而言，可以看出5組估計結(jié)果基本一致且誤差均在合理范圍之內(nèi)，因此可以說明基于單水聽器的勻速圓弧運動直升機3維參數(shù)估計算法獲取飛行參數(shù)信息的有效性。

表1 參數(shù)估計結(jié)果

圖8 時頻曲線的估計結(jié)果

5 結(jié)論

針對在水下估計反潛機的運動參數(shù)問題，本文從兩層等聲速介質(zhì)的聲傳播幾何模型出發(fā)，利用直升機的螺旋槳諧波頻率來計算水下聲學(xué)的多普勒頻移曲線。從3維空間幾何關(guān)系上推導(dǎo)了直升機在空中做圓弧運動時的航行參數(shù)，包括固有頻率、速度、高度、圓弧半徑等信息，與多普勒之間的非線性函數(shù)關(guān)系。利用4個時刻動聲源單水聽器之間的幾何關(guān)系，聲線由聲源傳至水聽器時間以及多普勒頻移曲線，可以計算勻速圓弧運動聲源的飛行參數(shù)。分析仿真數(shù)據(jù)的結(jié)果表明，本文提出的算法可以利用單水聽器準確地估計勻速圓弧運動直升機的飛行參數(shù)。