石海杰, 李京華, 陳 剛

(1. 西北工業(yè)大學(xué)電子信息學(xué)院, 陜西 西安 710129;2. 中北大學(xué)信息與通信工程學(xué)院, 山西 太原 030051;3. 中國聯(lián)合網(wǎng)絡(luò)通信有限公司東營市分公司, 山東 東營 257000)

0 引 言

水下探測空中聲源在軍民領(lǐng)域都得到廣泛關(guān)注,特別在軍事領(lǐng)域更具有重要意義[1-5]。潛艇是現(xiàn)代海軍重要裝備,其在近岸防衛(wèi)、突破封鎖、偵查掩護以及戰(zhàn)略威懾等方面都有重要作用。潛艇在國防領(lǐng)域的突出作用必然伴隨各種反潛技術(shù)的出現(xiàn)和發(fā)展,航空反潛是現(xiàn)代反潛技術(shù)中的重要方式。由于其具有機動靈活、通訊方便、協(xié)作迅捷的特點,航空反潛一直是空潛對抗中優(yōu)勢一方,為了提高潛艇的生存能力,水下對抗反潛飛機成為值得研究的課題。

航空反潛多采用螺旋槳飛機或直升機,這類飛機噪聲中含有螺旋槳轉(zhuǎn)動產(chǎn)生的線譜信號。線譜信號具有頻率較低、功率集中、穩(wěn)定性強的特點,可以傳輸?shù)捷^遠距離,是聲探測設(shè)備檢測和識別目標的主要信息載體。當(dāng)反潛機以一定高度勻速通過水面上方時,水下聲探測設(shè)備接收的信號將會產(chǎn)生多普勒頻移,其中包含目標速度、距離、通過時刻等參數(shù)的信息,這為目標運動參數(shù)的測量提供了可能。

利用聲信號多普勒頻移特征,在單一介質(zhì)中進行參數(shù)測量的方法已有較多研究[6-11],其中,Gomez -Tejedor[6]利用多普勒效應(yīng)分析了4種直線運動,分析結(jié)果與實際運動狀態(tài)有較好的一致性;Lindgren[7]利用聲傳感網(wǎng)絡(luò)獲取空氣中目標的多普勒信號,采用改進的高斯牛頓迭代法解決最小二乘最優(yōu)解問題,取得較好的定位效果;Timlelt[8]利用地面?zhèn)髀暺鹘邮者^頂直升機多普勒聲信號,估計聲源頻率和飛行速度等參數(shù),在10 dB信噪比條件下取得了較高的估計精度;Statman[9]推導(dǎo)得到運動目標多普勒瞬時頻率的解析表達,并采用最小二乘法進行參數(shù)估計,分析得出瞬時頻率變化率與正橫距離有較強相關(guān)性的結(jié)論;Xu[10]利用瞬時頻率估計方法實現(xiàn)了水中目標的定位;Feng[11]利用時頻分析方法測量水聲多普勒信號的瞬時頻率,該方法精度較高,但計算量較大。近年來,傳感網(wǎng)絡(luò)、數(shù)據(jù)融合、人工智能等技術(shù)在水聲探測方法中的應(yīng)用也是層出不窮,但多適用于復(fù)雜環(huán)境,實用性還有待驗證[12-15]。以上方法均沒有考慮聲音跨界傳輸?shù)膯栴},其結(jié)論亦不適用于水下探空的情況。Ferguson[16]和Lo[17]采用寬孔徑水聽器陣列,對空氣中過頂飛行的直升機進行參數(shù)估計和定位。Buckingham[18]實測直升機過頂飛行數(shù)據(jù),分別獲取了空氣中、海水中和海底沉積層中的聲音信號,實測數(shù)據(jù)表明過頂信號具有時變特性,驗證了在水中可以有效接收到空氣中聲源產(chǎn)生的多普勒頻移信號,這滿足了水下探空的前提條件。同時,水聲探測設(shè)備的不斷進步,也為水下探空提供了必要條件[19-20]。

在相關(guān)領(lǐng)域現(xiàn)有研究基礎(chǔ)上,本文提出一種基于多普勒效應(yīng)的單水聽器探測空中運動目標參數(shù)的方法。通過理論推導(dǎo)、數(shù)據(jù)仿真、實測驗證等途徑,設(shè)計出一種算法簡單,精度較高,能夠?qū)崿F(xiàn)速度、水平正橫距離、通過時刻等運動參數(shù)測量的有效方法。

1 空中目標水聲信號瞬時頻率模型

聲音由空氣傳入水中,會在空水界面產(chǎn)生反射和折射。由于二者聲阻抗的巨大差異,只有入射角很小的聲波才能傳入水中,水下聲信號好似由水面圓形小窗口引起,正如文獻[21-24]中所述的虛擬聲源。

空中點聲源等高勻速運動模型如圖1所示,圖中描繪了空水界面之上頻率為f0的聲源,以速度v沿平行于空水界面的直線AB勻速運動,A′B′為AB在空水界面上的投影,也就是虛擬聲源的運動軌跡,P為水聽器所在位置,A″B″為AB在水聽器所在平面上的投影,S為運動聲源距離水聽器的最近距離點(closest point of approach, CPA),S′為S在水面上的投影,水聽器深度為d,聲源高度為h,起始時刻聲源與S點距離為l,虛源與S′點距離同樣為l。由于聲源與其對應(yīng)的虛擬源具有相同的運動方式,因而具有相同的運動特征。

圖1 空中目標運動模型

假設(shè)τ時刻,A′點處的聲信號經(jīng)過傳輸,在t時刻到達水聽器處,則有

(1)

式中,cw是水中聲速;R是τ時刻虛擬聲源與水聽器的距離,則有

(2)

式中,R0是目標的水平正橫距離;d是水聽器深度;R1是虛擬聲源的正橫距離;tc是目標通過時刻;v是目標速度。

將式(2)代入式(1)中,可解得

(3)

t時刻水聽器接收的信號相位用z(t)表示,則τ時刻聲源發(fā)射的信號相位2πf0τ再加上一個相位與z(t)相等。則接收端的瞬時頻率f(t)[19]可表示為

(4)

由式(3)和式(4)聯(lián)合可得

(5)

式中,f0、cw和d是已知量;v、R0和tc是被測量?？刹捎米钚《朔ㄟM行參數(shù)估計,但這種方法計算復(fù)雜,特別是在多參數(shù)聯(lián)合估計時計算量比較大,本文將研究一種算法簡單、精度較高、適合于工程應(yīng)用的方法來進行參數(shù)的測量。

2 空中動目標參數(shù)測量

2.1 參數(shù)測量模型

本節(jié)將研究在已知聲源基準頻率f0的條件下,利用水聽器所接收的信號,根據(jù)運動目標瞬時頻率模型,測量目標運動速度v、通過時刻tc和水平正橫距離R0(正橫距離在水平方向上的投影)的方法。

由式(5)可以推導(dǎo)得到

(6)

(7)

式中,fa是聲源從較遠處接近水聽器運動時信號的瞬時頻率;fr為聲源遠離水聽器運動時信號的瞬時頻率。當(dāng)目標距離CPA點較遠時,可以認為fa和fr為固定值,由式(6)可得

(8)

由于fa可以在目標到達CPA點之前測量得到,因此參數(shù)v可以用來估計正橫距離和通過時刻。將式(5)對時間求導(dǎo),并假設(shè)v?cw(直升機反潛作業(yè)時,其速度一般不會很高,而水中聲速為1 500 m/s,這種假設(shè)是合理的),此時可得到等效源正橫距離的解析表達式。

(9)

再根據(jù)水聽器深度d可以計算目標水平正橫離R0。

由式(5)可推知,目標在CPA點的頻率fd(t)|t=tc=f0,假設(shè)t0時刻是水聽器接收到信號頻率f0的時刻,則得到目標通過時刻計算公式:

tc=t0-R1/cw

(10)

目標勻速通過CPA點,其多普勒頻率變化是連續(xù)的,也就是在多普勒頻率-時間曲線上,取CPA點附近的小段曲線可以近似認為是斜率為常數(shù)的直線。這就為在目標到達CPA點前一小段時間預(yù)測正橫距離和通過時刻提供了可能。

在目標到達CPA點(tc時刻)前的某時刻t進行預(yù)測,得到正橫距離為

(11)

t時刻預(yù)測的t0為

(12)

t時刻預(yù)測得到的通過時刻tc為

(13)

2.2 參數(shù)測量方法

根據(jù)以上推導(dǎo)過程,可以將基于多普勒效應(yīng)的運動目標參數(shù)估計方法描述如下:

步驟 1開始。一經(jīng)發(fā)現(xiàn)目標,開始對目標連續(xù)瞬時測頻。測頻時間間隔為T,fn為本次測頻結(jié)果,fn-1為前次測頻結(jié)果;

步驟 2估計fa。當(dāng)滿足一定條件時,fa=fn,判斷準則為|fn-1-fn|/T

步驟 3估計速度v。根據(jù)式(8)估計目標運動速度v;

步驟 4進行運動參數(shù)估計。啟動條件為|fn-f0|/T

步驟 5估計正橫距離。根據(jù)式(11)計算正橫距離R1,并根據(jù)水聽器深度,計算水平正橫距離R0;

步驟 6估計通過時刻。根據(jù)式(13)計算通過時刻tc;

步驟 7結(jié)束。

3 瞬時頻率估計方法

由以上的推導(dǎo)可知,本文所提出的方法依賴于對瞬時頻率的估計,因此本節(jié)討論適合在此應(yīng)用的頻率估計方法?；诳焖俑道锶~變換(fast Fourier transform, IFFT)類方法有諸多優(yōu)勢,但其頻率分辨率受到采樣時長的限制,在時變信號測量中不再適用。而以維格-威利分布(wigner-ville distribution, WVD)為代表的時頻分析類方法由于算法復(fù)雜,計算量較大,無法滿足實時性要求。針對以上問題,本文采用分段互譜密度(cross spectral density, CSD)的方法進行瞬時頻率的估計。

3.1 分段CSD方法理論基礎(chǔ)

對兩段時間長度為T的相鄰信號x1(t)、x2(t)采樣,采樣頻率為fs,采樣點數(shù)為2N。假設(shè)信號頻率為f1+f2,其中f1=k0fs/N=k0/T(k0為整數(shù)),是通過傅里葉變換能夠分辨的頻率;f2

采樣后的信號為

(14)

(15)

兩段信號的CSD可表示為

Y(k)=X2(k)·[X1(k)]*

(16)

式中,[·]表示點乘;*表示共軛;X1(k)是x1(n)的傅里葉變換;X2(k)是x2(n)的傅里葉變換。根據(jù)維納-辛欽定理,變換相關(guān)運算與傅里葉變換的順序,可得

Y(k)=DFT{x2(n)·[x1(n)]*}=DFT[e-j2π(f1+f2)T]=DFT[e-j2π(k0+f2T)]=δ(2πf2T)

(17)

根據(jù)式(17),可得

f2=arg max[Y(k)]/2πT

(18)

f1為傅里葉變換可分辨的頻率,因此

f1=ser max[X1(k)]fs/N

(19)

式中,ser max[]表示取序列最大值對應(yīng)的序號。

根據(jù)式(19)和式(20),當(dāng)給定兩段相鄰信號x1(n)和x2(n)時,其測頻表達式為

f=ser max[X1(k)]fs/N+arg max[Y(k)]/2πT

(20)

3.2 分段CSD算法仿真

為了驗證CSD算法的有效性,選取仿真頻率為f1+f2=100.2 Hz的信號,采樣頻率fs=2 000 Hz,單次測頻采樣時長T=1 s,即傅里葉變換頻率分辨率為1 Hz。其中f1=100 Hz是可以通過FFT分辨并測量的頻率,f2=0.2 Hz是FFT無法分辨的頻率,必須采用CSD方法測量。

18F-FDG PET/CT SUVmax與淋巴瘤患者臨床特征及生物學(xué)指標的相關(guān)性(張玲芳)(12):1143

分別在無噪聲和信噪比SNR=0 dB的條件下進行頻率估計仿真實驗,結(jié)果如表1和圖2所示。表1所示為30次、50次、100次蒙特卡羅實驗的估計均值和均方根誤差。結(jié)合表1中f1估計結(jié)果和圖2(a)與圖2(c)可以看出,無論有無噪聲,FFT均能測出精確到1 Hz的頻率,說明FFT方法能夠較好地估計信號的頻率,但是其測頻精度受到采樣時長的限制,無法分辨出f2的存在。

表1 CSD方法測頻結(jié)果

圖2 CSD法測頻仿真結(jié)果

結(jié)合表1中f2估計結(jié)果和圖2(b)與圖2(d)可以看出,CSD方法能夠有效估計出f2=0.2 Hz的頻率。同時也能看出CSD方法測頻精度是受到噪聲影響的,由圖2(d)可以看出,當(dāng)信噪比為0 dB時,測頻結(jié)果在真實頻率0.2 Hz附近有一定起伏。表1中f2估計結(jié)果表明,隨著蒙特卡羅次數(shù)的增加,其估計均值逐漸接近于真實值,均方根誤差(root mean square error, RMSE)基本維持在0.01 Hz量級,說明在0 dB信噪比條件下,CSD測頻方法具有無偏一致性,測頻精度可以達到0.01 Hz量級,是一種有效的測頻方法。

分段CSD測頻方法除了受到噪聲影響外,還與單次測頻采樣時長有關(guān),當(dāng)采樣時間過短時,噪聲對于估計結(jié)果的影響將會嚴重。為了分析噪聲和采樣時間對于估計結(jié)果的影響,仿真了信噪比SNR=0 dB的條件時,不同采樣時間對估計結(jié)果分布情況的影響。仿真結(jié)果如圖3(a)所示?？梢钥闯?采樣時間越短,估計結(jié)果偏離真實值越大,估計結(jié)果越分散。同時仿真了在不同信噪比的條件下,估計結(jié)果的均方根誤差隨著采樣時間變化的情況,仿真結(jié)果如圖3(b)所示,可以看出,均方根誤差隨著采樣時間增大而減小,隨著信噪比降低而增大。從圖3(b)中還可以看出,在采樣時間T=1 s時,信噪比大于0 dB(包含0 dB)的信號,估計誤差不大,明顯好于信噪比小于0 dB的信號。

圖3 CSD測頻方法分析

4 仿真分析

4.1 數(shù)據(jù)仿真

由于實測數(shù)據(jù)受到海域、海況、季節(jié)、時間等多種因素影響,其信噪比等參數(shù)不穩(wěn)定,難于準確測量,不適合定量分析。因而,采用仿真方法可以進行算法性能分析?？焖賵瞿Ｐ鸵圆〝?shù)積分方法為理論基礎(chǔ),是聲波動方程的全波解,是一種適用于本文聲場計算的數(shù)值仿真方法[25-26]。

假設(shè)空氣、海水和海底都是均勻液體,聲學(xué)參數(shù)為ρa=1.293 kg/m3,ca=340 m/s;ρw=1 023 kg/m3,cw=1 480 m/s;ρb=1 620 kg/m3,cb=1 806 m/s。其中ρ代表密度,c代表聲速,下標a代表空氣參數(shù)、下標w代表海水參數(shù)、下標b代表海底沉積層參數(shù)。仿真參數(shù)如表2所示,其中,h為聲源高度,d是水聽器深度,lcpa是聲源的水平正橫距離,v為目標速度,f0是聲源頻率,H為海水深度。

表2 數(shù)據(jù)仿真參數(shù)

仿真參數(shù)和前述環(huán)境聲學(xué)參數(shù)如表2所示,利用Zhang[15]所述的基于快速場模型的運動目標水聲信號產(chǎn)生方法,仿真空氣中運動目標激發(fā)的水聲信號,仿真過程中加入適當(dāng)?shù)脑肼?最終的水聽器接收信噪比為0 dB。仿真結(jié)果如圖4所示,其中4(a)為時域信號,可以看出聲源在中間時刻(30 s附近)有明顯的過頂飛行過程;圖4(b)是該仿真信號的時頻圖,可以看出在聲源過頂飛行過程中產(chǎn)生了明顯的多普勒頻移。

圖4 波數(shù)積分方法仿真水中接收信號

4.2 算法分析

在第2.1節(jié)中,對運動目標參數(shù)測量方法做理論分析時,做了小段多普勒頻率-時間曲線近似直線的假設(shè),這種假設(shè)成立的條件受到正橫距離和目標速度關(guān)系的影響。在單次測頻采樣時長確定的條件下,理論上來說,有效通過時長(R0/v)越長,相當(dāng)于測頻點越密集,小段頻率-時間曲線更接近于直線,參數(shù)估計結(jié)果的精度就越高。但參數(shù)估計結(jié)果還依賴于相鄰兩個測點的頻率差值,測點相鄰過近時,單個頻點的測頻誤差對估計結(jié)果影響增大。本小節(jié)將分析信噪比SNR=0 dB、采樣時間T=1 s的條件下,估計誤差與各參數(shù)間的關(guān)系。

圖5所示為信噪比SNR=0 dB、采樣時長T=1 s的條件下,速度分別為50 m/s、75 m/s和100 m/s時,速度v、通過時刻tc、水平正橫距離R0的估計誤差隨著R0/v的變化情況。其中圖5(a)是速度相對誤差,圖5(b)是通過時刻絕對誤差,圖5(c)是水平正橫距離相對誤差。

圖5 估計誤差分析

從圖5(a)中可以看出測速誤差受到有效通過時長(R0/v)的影響不大,這與實際情況相符。由式(8)可知,測速誤差和目標與CPA點的距離有關(guān),距離越遠,fa越接近穩(wěn)定值,估計誤差越小。因此當(dāng)仿真的數(shù)據(jù)時間越長,速度估計誤差越小;而在實際工程中,如能在越遠處發(fā)現(xiàn)目標,測速誤差就會越小。

從圖5(b)中可以看出通過時刻估計誤差受R0/v影響,并且當(dāng)R0/v的取值在5～15范圍時,估計誤差較小,超出這個范圍,估計誤差變大。正是由于估計誤差同時受到直線近似誤差、單點測頻誤差和頻率估計間隔共同影響造成的。

從圖5(c)中可看出水平正橫距離估計誤差同樣受到有效通過時長(R0/v)的影響,并且當(dāng)R0/v的取值范圍為10～15時,估計誤差有較小值。超出這個范圍時,估計誤差變大,特別是當(dāng)R0/v的取值小于2時,誤差顯著變大。這是由于目標速度過快,有效通過時間過短,信號的多普勒頻率變化在很短的時間內(nèi)發(fā)生,而測頻方法對瞬時頻率分辨有限所致。

5 實驗驗證

某次實測實驗兩組數(shù)據(jù)的時域圖和時頻圖如圖6所示。

圖6 實測信號

實驗采用的是某型號螺旋槳飛機,飛機飛行高度為100 m,以150 km/h的速度等高勻速飛行,通過水聽器所在水域的上方,水平正橫距離為200 m,水聽器位于水下深度為15 m,兩組實測數(shù)據(jù)命名為實測數(shù)據(jù)1和實測數(shù)據(jù)2。同時,按照第4.1節(jié)的仿真方法,假設(shè)飛機飛行速度v=41.7 m/s,水平正橫距離為200 m,其他參數(shù)與第4.1節(jié)所述相同,仿真一組待測數(shù)據(jù),命名為仿真數(shù)據(jù)1。

實測實驗數(shù)據(jù)表明(見圖6),沒有飛機過頂飛行時的環(huán)境噪聲信號幅度(見圖6(a)0～20 s的時域信號)明顯小于有飛機過頂飛行的信號幅度(見圖6(a)20～40 s的時域信號),說明0 dB信噪比的假設(shè)是適合本文應(yīng)用背景的。因此,本文在0 dB信噪比條件下對算法進行分析和驗證,采樣時間T取1 s。

按照本文的參數(shù)估計方法,對以上3組數(shù)據(jù)進行目標參數(shù)估計,單次估計采樣時長T=1 s,估計結(jié)果如表3所示。

表3 測量結(jié)果

其中，v、tc和R0分別代表目標速度、通過時刻和水平正橫距離的真實值,M代表估計值,E代表誤差。由于實測數(shù)據(jù)的通過時刻沒有參照,因此未做誤差分析。對估計結(jié)果進行對比分析可得出如下結(jié)論。

仿真數(shù)據(jù)1的有效通過時長(R0/v)為4.79。仿真結(jié)果中目標速度、通過時刻和水平正橫距離的測量誤差分別為0.4%、0.48 s和3.6%;圖5分析結(jié)果中,R0/v=4.79處3個參數(shù)的誤差分別約為0.25%、0.3 s和1.0%。兩者對比可以發(fā)現(xiàn),估計結(jié)果與分析結(jié)果是相符合的,水平距離估計誤差略大,這是由單次估計的隨機性造成的。

兩組實測數(shù)據(jù)的有效通過時長(R0/v)為4.79(受到實驗電纜長度、實驗飛機安全速度條件限制)。目標速度和水平正橫距離兩組數(shù)據(jù)的估計誤差均值分別為8.2%和18.1%。相同條件下的仿真數(shù)據(jù)兩參數(shù)估計誤差分別為0.4%和3.6%。兩者對比可以發(fā)現(xiàn),實測數(shù)據(jù)估計誤差大于仿真數(shù)據(jù)估計誤差,這是由于實際海洋環(huán)境受到海域、海況和人類活動噪聲等多種因素的影響,在沒有做環(huán)境具體分析和校正補償?shù)那闆r下,其誤差必然要大于受控的仿真情況。因此,海域、海況、人類活動等因素對于估計精度的影響是該方法應(yīng)用于工程時必須要考慮的問題,也是作者在此后研究中將要關(guān)注的問題之一。

6 結(jié) 論

本文建立了一種利用水下聲信號探測空中運動目標的模型,提出了行之有效的運動目標參數(shù)估計方法。通過仿真分析,得出參數(shù)估計精度受到有效通過時長(R0/v)影響的結(jié)論。利用運動目標瞬時頻率模型,把正橫距離點處的多普勒-時間曲線等效成小段直線,把目標參數(shù)估計問題轉(zhuǎn)化為瞬時頻率估計問題,使問題得以采用簡單、有效的方法加以解決;采用分段CSD方法,提高瞬時頻率估計精度,在滿足實時性要求的條件下,提高了運動目標參數(shù)的估計精度。

本文方法在仿真的條件下取得了較高的估計精度,并通過實測實驗對方法進行了可行性的驗證。實際海洋環(huán)境存在各種噪聲,如何提高在復(fù)雜海洋環(huán)境下本文所述方法的估計精度,是需要進一步深入研究的工作。