王 磊，陳昭男，張 君，周 唯，佟健飛

(1. 中國人民解放軍 91550部隊， 遼寧 大連 116023； 2. 西北工業(yè)大學 自動化學院， 陜西 西安 710129； 3. 中國科學院聲學研究所， 北京 100190)

亞音速彈道聲學測量以亞音速飛行彈體為目標，在亞音速飛行過程中，彈體與氣流之間會產(chǎn)生“氣動噪聲”。氣動噪聲是指空氣流動或者物體在空氣中運動引起空氣產(chǎn)生旋渦、沖擊或者壓力突變而導致空氣擾動所形成的噪聲，屬于寬帶非平穩(wěn)的噪聲源[1-2]。

由于亞音速彈體的信號顯著性特征不清晰，對于超音速彈體彈道聲學測量的相關研究較多[3-4]，關于亞音速彈道聲學測量分析的相關研究較少。對于低空飛行目標的小尺寸傳感器高精度定向方面，文獻[5]首次將壓差式矢量聲傳感器引入低空目標側向研究中，提出基于矢量傳感器陣的寬帶相干信號子空間最優(yōu)波束形成算法。

文獻[6]在亞音速末彈道聲學落點研究中，采用快速傅里葉變換與拉普拉斯小波分析技術對末端亞音速彈道噪聲進行了信號特征建模，提出了彈丸在亞音速飛行時，光滑彈體和氣流之間會產(chǎn)生氣動噪聲，并結合大量測試數(shù)據(jù)，得到了亞音速飛行彈丸氣動噪聲特點。但彈丸的氣動特性并非一致的，對于亞音速飛行的彈丸而言，彈體在飛行過程中氣動特性和彈體外形、體積和飛行速度等因素有關[7-8]，因此彈體和氣流作用產(chǎn)生的寬帶非平穩(wěn)氣動噪聲的聲源顯著特性不易獲取。

針對該類寬帶噪聲，經(jīng)典的譜分析技術不能顯著地刻畫其信號特征，合適的時頻分析手段，對于目標的信號特征分析會更有效[9]。本文采用時頻分析的方法提取氣動噪聲的聲學特性，并基于時頻分析模型對氣動噪聲的方位角進行估計。

1 亞音速氣動噪聲及聲學特性

圖1所示為一彈體發(fā)射后的飛行過程。

圖1 彈體亞音速飛行過程Fig.1 Subsonic fight process of projectile

彈體發(fā)射時，火藥在窄小的空間內燃燒產(chǎn)生高溫、高速的爆炸氣流將彈丸推出槍膛，在彈體出膛后，噴出膛口的氣流形成膛口激波，以聲速從槍口向四周傳播(近似為球面波)。隨后彈體作亞音速飛行，彈體和空氣氣流之間產(chǎn)生氣動噪聲，即高速氣流、不穩(wěn)定氣流以及彈體之間相互作用產(chǎn)生的噪聲。

(a) 原始波形時域圖(a) Time-domain diagram of original waveform

(b) 原始波形時頻譜圖(b) Time-frequency spectrum of original waveform圖2 實驗數(shù)據(jù)時域波形及時頻譜圖Fig.2 Time-domain waveform and time-frequency spectrum of experimental data

圖2為某次炮彈實驗采集的實際信號的時域圖和時頻譜圖。從時域圖中可以看出，0.1～0.2 s之間的沖擊波形為膛口激波，0.3～0.5 s之間的波形為彈體飛行過程中產(chǎn)生的氣動噪聲，其中膛口激波的強度大于氣動噪聲。從信號時頻譜圖中可以看出，膛口激波主要集中在100～600 Hz的頻帶范圍內，氣動噪聲則主要分布在1～14 kHz的頻帶范圍內。由于膛口激波和氣動噪聲處于不同頻段范圍內，對采集的信號使用不同的濾波器濾波，可實現(xiàn)兩種波形的分離。實驗數(shù)據(jù)濾波后的波形如圖3所示，當原始信號通過帶通濾波器后，氣動噪聲被削弱，從而保留了膛口激波；當原始信號經(jīng)過高通濾波器后，膛口激波被削弱，氣動噪聲被保留。

(a) 經(jīng)過帶通濾波器后的波形——膛口波(a) Waveform after bandpass filtering of muzzle wave

(b) 經(jīng)過高通濾波器后的波形——氣動噪聲(b) Waveform after high pass filtering of aerodynamic noise圖3 實驗數(shù)據(jù)濾波后的波形Fig.3 Waveform results after filtering of experimental data

2 基于小波變換的氣動噪聲波達方向估計

2.1 小波分析

氣動噪聲屬于寬帶非平穩(wěn)信號，其時域和頻域聲學特征不顯著，隨著時間變化，其統(tǒng)計特性也隨之發(fā)生變化。此時，傳統(tǒng)的傅里葉變換已經(jīng)不能夠滿足分析要求。而小波變換是一種信號的時間頻率分析方法，它具有多分辨特點，在時域和頻域同時具有良好的局部性質[10-11]。連續(xù)小波變換能夠同時對信號的低頻和高頻信息進行分析，因此更加適用于氣動噪聲信號的處理。

對于任意函數(shù)s(t)∈L2(R)的連續(xù)小波變換定義為：

WTs(a,b)=〈s(t),ψa,b(t)〉

(1)

在小波分析過程中，使用的小波函數(shù)不具有唯一性。同一個工程問題，使用不同的小波函數(shù)進行分析，結果不完全相同。因此，小波函數(shù)的選擇是小波分析應用到實際中的關鍵問題，目前小波函數(shù)的選擇常通過經(jīng)驗或不斷實驗來確定。主要選取的小波函數(shù)有Morlet小波、Mexh小波以及Gaussian小波，它們的時域、頻域數(shù)學形式見表1。

對圖2所示的炮彈實驗采集到的數(shù)據(jù)分別使用Morlet小波、Mexh小波及Gaussian小波進行分析，可得表2中不同小波函數(shù)下的小波時頻譜圖以及各尺度下分解出的信號，其中尺度和頻率成反比。

表1 不同類型小波函數(shù)

根據(jù)不同小波函數(shù)的時頻譜圖可以看出，以Morlet小波作為分析小波獲得的時頻譜圖與氣動噪聲特性對應，且基于Morlet小波下的時頻聚集性較好，能夠更好地提取出彈道波。根據(jù)不同尺度下分解的信號可以看出，Morlet小波在2尺度下的信號幅值較其他尺度下的小，且在高尺度下已經(jīng)濾除膛口波，而Mexh小波和Gaussian小波的效果并不如Morlet小波。除此之外，與圖3所示的實驗信號通過高通濾波器獲取的結果對比，Morlet小波依然分析較優(yōu)。根據(jù)以上結論，采用Morlet小波作為分析小波，對亞音速飛行氣動噪聲進行波達方向(Direction of Arrival, DOA)估計。

2.2 氣動噪聲DOA估計流程

設計四元十字陣，各陣元坐標分別為(d,0)、(0,d)、(-d,0)以及(0,-d)，考慮存在K個寬帶聲源分別以θi(i=1,…,K)的方位角入射到該十字陣，可得第m個陣元的接收信號為：

(2)

式中，τmk=(xmcosθk+ymsinθk)/c表示第k個聲源入射到第m個陣列所引起的時延，nm(t)為第m個陣元的量測噪聲。

對式(2)第m個傳聲器采集數(shù)據(jù)進行小波變換，借助小波變換的線性疊加性并根據(jù)Parseval定理，可將小波變換公式化簡為：

(3)

表2 不同小波函數(shù)下的小波時頻譜圖及分解信號對比

WT(ai,b)=A(ai,θ)Y(ai,b)+Yn(ai,b)

(4)

Y(ai,b)=[y1(ai,b),…,yk(ai,b)]T

因此可使用多信號分類(MUltiple SIgnal Classification， MUSIC)算法在小波域上進行聲源的來向估計。常規(guī)MUSIC算法主要利用信號和噪聲子空間的正交性進行參數(shù)估計，而對于本文所示的非平穩(wěn)信號，使用連續(xù)小波變換同時對信號時域和頻域進行分析，從而獲取多分辨時頻陣列信號模型，并引入MUSIC算法中以提高 DOA 估計性能。獲取式(4)所示陣列接收數(shù)學模型的協(xié)方差矩陣。

R=E{WT(ai,b)·WTH(ai,b)}

(5)

對協(xié)方差矩陣進行特征分解，獲得信號子空間UY(ai)與噪聲子空間UY(ai)，并借助信號子空間與噪聲子空間的正交性，根據(jù)MUSIC算法獲取聲源來向。將信號轉換到小波域，假設存在I個尺度，即可得到式(5)所示空間，遍歷整個角度范圍，搜索空間譜最大值所對應的角度即為入射聲源的方位角。

(6)

3 實驗驗證

2018年5月在靶場進行外場實驗，實驗采集火箭彈發(fā)射的聲信號，驗證本文算法的可行性。

實驗陣列為1 m孔徑平面?zhèn)髀暺麝嚵?BK4189)。具體實驗過程如圖4所示，火箭炮在距離立靶211 m處發(fā)射，四元傳聲器陣列位于彈道垂直距離30 m、與立靶的距離102 m處，其中令1號傳聲器對應方向為正x軸方向，2號傳聲對應方向為正y軸方向。

對實驗過程進行時間切片，截取彈體飛行到距離立靶14 m時刻的數(shù)據(jù)，估計方位角θ。實際上，該時刻彈體相對于陣列的角度為205°。使用常規(guī)寬帶MUSIC算法估計方位角θ，如圖5所示。常規(guī)寬帶MUSIC算法獲得的空間譜最大值所對應的方位角為223°，該結果與真實方位角相差較大，已經(jīng)不能估計出彈體相對于傳聲器的方位角。

圖4 實驗過程Fig.4 Experimental verification process

現(xiàn)分別基于Morlet小波、Mexh小波、Gaussian小波進行方位角估計，其估計結果如圖6所示?；贛exh小波的估計結果為225°，該結果已超過誤差允許范圍；基于Gaussian小波的空間譜圖結果，已經(jīng)出現(xiàn)了較多偽峰，完全估計不出方位角結果；基于Morlet小波的空間譜圖中可得截取時刻的角度為204°，該結果在誤差允許范圍內，從而驗證了Morlet小波作為分析小波時，可有效地估計彈體方位角。除此之外，使用基于時延估計的廣義互相關(Generialized Cross-Correlation， GCC)算法進行了方位角計算，用以參考，GCC算法獲取的角度結果為206°，與實際結果相差1°，估計結果較好。

由于Morlet小波作為分析小波時可有效地估計彈體方位角。因此選擇使用Morlet小波對整個彈體運動角度進行估計，獲取圖7所示的基于Morlet小波的彈體角度估計結果。除此之外，使用GCC算法進行角度估計，與基于Morlet小波算法的估計結果進行對比，驗證了本文所提算法的優(yōu)勢。

圖5 寬帶MUSIC算法空間譜估計結果Fig.5 Spatial spectrum estimation results of wide-band MUSIC algorithm

(a) 基于Morlet小波的空間譜估計結果(a) Spatial spectrum based on Morlet wavelet

(b) 基于Mexh小波的空間譜估計結果(b) Spatial spectrum based on Mexh wavelet

(c) 基于Gaussian小波的空間譜估計結果(c) Spatial spectrum based on Gaussian wavelet圖6 基于不同小波函數(shù)的MUSIC算法空間譜圖Fig.6 Spatial spectrum of MUSIC algorithm based on different wavelet functions

(a) 信號時域波形圖(a) Signal time-domain waveform

(b) 基于Morlet小波的估計結果(b) Estimation results based on Morlet wavelet

(c) GCC算法的估計結果(c) Estimated result of GCC algorithm圖7 基于Morlet小波的彈體角度估計結果Fig.7 DOA estimation results of projectile based on Morlet wavelet

由于整個實驗過程主要是對彈體運動過程產(chǎn)生的氣動噪聲進行估計獲取彈體運動角度軌跡，因此只估計傳聲器陣列采集到氣動噪聲的時間段的角度結果。圖7為信號時域波形圖，由圖7可得，氣動噪聲部分為0.3～0.5 s之間，理論上幅值最大時刻應該為傳聲器與彈道的垂直位置附近(相對于傳聲器225°附近)。GCC算法在0.3～0.4 s之間的估計結果較為平穩(wěn)地從180°運動到243°，之后的估計結果則遠離了真實軌跡。基于Morlet小波的角度估計結果可得：在采集到氣動噪聲的時間段內，所提算法的估計結果比較平穩(wěn)地從175°運動到275°，也就是傳聲器與彈道距離為30 m時，可有效估計出彈體與傳聲器的最遠距離為46.7 m。

4 結論

本文提出基于小波函數(shù)的MUSIC算法，將亞音速飛行的彈體產(chǎn)生的氣動噪聲轉換到小波域，并通過MUSIC算法實現(xiàn)彈體的方位角估計。該算法的處理過程中，需要對信號進行預處理，將其轉換到小波域，而在分析小波的選擇中，由實驗數(shù)據(jù)分析結果可得Morlet小波較Mexh小波、Gaussian小波在氣動噪聲處理更具優(yōu)勢。外場實驗驗證了所提算法的可行性和分析結果的正確性：使用不同小波下的MUSIC算法進行角度估計，可得Morlet小波下的估計結果更加準確，趨于真實結果。