水聲探空動目標參數測量方法

石海杰, 李京華, 陳 剛

(1. 西北工業大學電子信息學院, 陜西 西安 710129;2. 中北大學信息與通信工程學院, 山西 太原 030051;3. 中國聯合網絡通信有限公司東營市分公司, 山東 東營 257000)

0 引 言

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

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

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

在相關領域現有研究基礎上,本文提出一種基于多普勒效應的單水聽器探測空中運動目標參數的方法。通過理論推導、數據仿真、實測驗證等途徑,設計出一種算法簡單,精度較高,能夠實現速度、水平正橫距離、通過時刻等運動參數測量的有效方法。

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

聲音由空氣傳入水中,會在空水界面產生反射和折射。由于二者聲阻抗的巨大差異,只有入射角很小的聲波才能傳入水中,水下聲信號好似由水面圓形小窗口引起,正如文獻[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。由于聲源與其對應的虛擬源具有相同的運動方式,因而具有相同的運動特征。

圖1 空中目標運動模型

假設τ時刻,A′點處的聲信號經過傳輸,在t時刻到達水聽器處,則有

(1)

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

(2)

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

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

(3)

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

(4)

由式(3)和式(4)聯合可得

(5)

式中,f0、cw和d是已知量;v、R0和tc是被測量。可采用最小二乘法進行參數估計,但這種方法計算復雜,特別是在多參數聯合估計時計算量比較大,本文將研究一種算法簡單、精度較高、適合于工程應用的方法來進行參數的測量。

2 空中動目標參數測量

2.1 參數測量模型

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

由式(5)可以推導得到

(6)

(7)

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

(8)

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

(9)

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

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

tc=t0-R1/cw

(10)

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

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

(11)

t時刻預測的t0為

(12)

t時刻預測得到的通過時刻tc為

(13)

2.2 參數測量方法

根據以上推導過程,可以將基于多普勒效應的運動目標參數估計方法描述如下:

步驟 1開始。一經發現目標,開始對目標連續瞬時測頻。測頻時間間隔為T,fn為本次測頻結果,fn-1為前次測頻結果;

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

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

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

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

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

步驟 7結束。

3 瞬時頻率估計方法

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

3.1 分段CSD方法理論基礎

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

采樣后的信號為

(14)

(15)

兩段信號的CSD可表示為

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

(16)

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

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

(17)

根據式(17),可得

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

(18)

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

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

(19)

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

根據式(19)和式(20),當給定兩段相鄰信號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與淋巴瘤患者臨床特征及生物學指標的相關性(張玲芳)(12):1143

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

表1 CSD方法測頻結果

圖2 CSD法測頻仿真結果

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

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

圖3 CSD測頻方法分析

4 仿真分析

4.1 數據仿真

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

假設空氣、海水和海底都是均勻液體,聲學參數為ρ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代表空氣參數、下標w代表海水參數、下標b代表海底沉積層參數。仿真參數如表2所示,其中,h為聲源高度,d是水聽器深度,lcpa是聲源的水平正橫距離,v為目標速度,f0是聲源頻率,H為海水深度。

表2 數據仿真參數

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

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

4.2 算法分析

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

圖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點的距離有關,距離越遠,fa越接近穩定值,估計誤差越小。因此當仿真的數據時間越長,速度估計誤差越小;而在實際工程中,如能在越遠處發現目標,測速誤差就會越小。

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

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

5 實驗驗證

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

圖6 實測信號

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

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

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

表3 測量結果

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

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

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

6 結 論

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

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