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

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

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

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

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

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

1 亞音速氣動(dòng)噪聲及聲學(xué)特性

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

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

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

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

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

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

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

(b) 經(jīng)過高通濾波器后的波形——?dú)鈩?dòng)噪聲(b) Waveform after high pass filtering of aerodynamic noise圖3 實(shí)驗(yàn)數(shù)據(jù)濾波后的波形Fig.3 Waveform results after filtering of experimental data

2 基于小波變換的氣動(dòng)噪聲波達(dá)方向估計(jì)

2.1 小波分析

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

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

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

(1)

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

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

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

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

2.2 氣動(dòng)噪聲DOA估計(jì)流程

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

(2)

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

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

(3)

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

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

(5)

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

(6)

3 實(shí)驗(yàn)驗(yàn)證

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

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

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

圖4 實(shí)驗(yàn)過程Fig.4 Experimental verification process

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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