采用單水聽器的勻速圓弧運(yùn)動(dòng)直升機(jī)三維參數(shù)估計(jì)算法

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

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

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

1 引言

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

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

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

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

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

圖1 空中沿曲線運(yùn)動(dòng)的點(diǎn)聲源與靜止水聽器節(jié)點(diǎn)的3維圖

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

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

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

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

圖2 空中動(dòng)點(diǎn)聲源與靜止水聽器節(jié)點(diǎn)的俯視圖

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

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

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

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

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

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

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

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

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

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

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

3 估計(jì)飛行參數(shù)

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

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

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

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

這里

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

首先，令點(diǎn)B坐標(biāo)為

這里

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

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

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

其中，聲源到達(dá)點(diǎn)B的時(shí)刻為

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

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

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

∠PH′Q

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

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

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

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

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

P ∈{A,B,C,M}

聲源位于 時(shí)，其瞬時(shí)多普勒頻率為

4 實(shí)驗(yàn)與分析

4.1 多普勒頻移特征

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

圖5 仿真選取的5組測(cè)量點(diǎn)與靜止水聽器幾何示意圖

考慮到時(shí)間間隔、飛行朝向以及飛行高度，仿真中分別取了5組觀測(cè)點(diǎn)，以說(shuō)明選取聲源坐標(biāo)對(duì)飛行參數(shù)估計(jì)精度的影響，其中第1,2組高度為h1=150 m、時(shí)刻間隔為0.5 s，第3,4組高度為150 m、時(shí)刻間隔為0.2 s；第5組高度為100 m、時(shí)刻間隔為0.5 s。其中，第1,3,5組聲源朝向水聽器運(yùn)動(dòng)，第2,4組聲源遠(yuǎn)離水聽器運(yùn)動(dòng)。

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

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

4.2 估計(jì)多普勒頻移

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

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

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

圖7 水聽器接收的時(shí)域信號(hào)

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

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

圖6 勻速圓弧運(yùn)動(dòng)聲源的多普勒頻移曲線

4.3 估計(jì)3維參數(shù)

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

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

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

圖8 時(shí)頻曲線的估計(jì)結(jié)果

5 結(jié)論

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