三元次聲陣定位火箭發(fā)射的方法

殷 昊，龐新良，張震川

（國(guó)民核生化災(zāi)害防護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102205）

0 引 言

次聲源定位是一種典型的被動(dòng)式探測(cè)定位技術(shù)，根據(jù)利用的信號(hào)信息和參數(shù)的不同，可大致將定位技術(shù)分為四類：利用信號(hào)的到達(dá)時(shí)間（Time of Arrival,TOA）定位、利用信號(hào)的到達(dá)時(shí)間差（Time Difference of Arrival,TDOA）定位、利用信號(hào)的到達(dá)方位角（Azimuth of Arrival,AOA）定位、利用接收信號(hào)的強(qiáng)度（Received Signal Strength,RSS）定位。其中，基于到達(dá)時(shí)間差的定位方法定位精度較高，實(shí)用性較強(qiáng)，是目前應(yīng)用與研究最廣泛的一種聲源定位方法。

常用的是基于互相關(guān)的時(shí)延估計(jì)方法，其中最常用的方法主要是廣義互相關(guān)法（Generalized Cross Correlation,GCC）。該方法被Knapp等[1]于1976年提出，其基本原理是在求取信號(hào)互相關(guān)函數(shù)之前對(duì)功率譜進(jìn)行加權(quán)濾波，突出信號(hào)并抑制噪聲干擾部分，從而突出相關(guān)函數(shù)在時(shí)延處的峰值，取其峰值進(jìn)行時(shí)延估計(jì)。廣義互相關(guān)方法采用了眾多加權(quán)函數(shù)[2-6]，如由 Roth提出的Roth加權(quán)函數(shù)（The Roth Processor,ROTH）、平滑相關(guān)變換（The Smoothed Coherence Transform,SCOT）、相位變換（The Phase Transform,PHAT）、最大似然（Maximum Likelihood,ML）估計(jì)/HT（Hannan-Thomson）以及一些改善函數(shù)等，起到銳化峰值、抑制噪聲的作用。近幾年來，眾多學(xué)者致力于 GCC時(shí)延估計(jì)方法的改進(jìn)，小波變換[7]、時(shí)頻分析、希爾伯特變換[8]、神經(jīng)網(wǎng)絡(luò)[9]、經(jīng)驗(yàn)?zāi)B(tài)分解（Empirical Mode Decomposition,EMD）[10]等方法被引入到時(shí)延估計(jì)方法中，能有效地進(jìn)行非平穩(wěn)環(huán)境下的時(shí)延估計(jì)。為了更進(jìn)一步精確時(shí)延量，楊亦春等[11]針對(duì)震前異常次聲波一次GCC得到的時(shí)延量遠(yuǎn)遠(yuǎn)大于信號(hào)持續(xù)時(shí)間的情況，對(duì)信號(hào)采用二步法時(shí)延估計(jì)：第一步先根據(jù)能量分布粗略估計(jì)時(shí)延；第二步再用 GCC法計(jì)算精確時(shí)延，最終的時(shí)延量為兩次時(shí)延之和。考慮到一次互相關(guān)法受噪聲影響較大，文獻(xiàn)[12]提出了二次相關(guān)法，可在更低的信噪比環(huán)境下較準(zhǔn)確地估計(jì)時(shí)延。

因?yàn)榛鸺l(fā)射的時(shí)間與發(fā)射點(diǎn)的位置明確，對(duì)于檢測(cè)次聲時(shí)延估計(jì)、定位算法等具有重要意義，具有重要的研究?jī)r(jià)值[13]。火箭發(fā)射是一種運(yùn)動(dòng)次聲源，其產(chǎn)生的次聲波信號(hào)主要由火箭飛行時(shí)激發(fā)空氣產(chǎn)生，因此近距離采集到的次聲信號(hào)持續(xù)時(shí)間可長(zhǎng)達(dá) 3～4 min，甚至更久。文獻(xiàn)[14]考慮到次聲監(jiān)測(cè)網(wǎng)絡(luò)孔徑較大且信號(hào)相關(guān)性較小，獲取的時(shí)延估計(jì)值準(zhǔn)確度不高，因此采用最小方差法對(duì)陣列周邊的區(qū)域進(jìn)行掃描，進(jìn)而得出聲源位置的估計(jì)值，與實(shí)際位置相差數(shù)十千米。唐偉等[13]對(duì)朝鮮4.13“光明星3號(hào)”衛(wèi)星發(fā)射次聲信號(hào)進(jìn)行分析，分別采用單臺(tái)、兩個(gè)次聲臺(tái)站兩種方法對(duì)火箭爆炸區(qū)域進(jìn)行估算。

由于次聲波的遠(yuǎn)距離傳播特性，大多數(shù)學(xué)者考慮的是次聲源的遠(yuǎn)程定位。文獻(xiàn)[15]中考慮在小于100 km的范圍內(nèi)，次聲定位應(yīng)用的缺乏，開發(fā)了一種 srcLoc方法，該方法是一種近場(chǎng)、嚴(yán)格基于TDOA的聲學(xué)定位方法，其優(yōu)勢(shì)主要在于沒有限制性的大氣假設(shè)，試驗(yàn)結(jié)果表明該方法在96%的站點(diǎn)上比傳統(tǒng)的定位方法的定位精度提高了2倍以上。而火箭次聲信號(hào)的持續(xù)時(shí)間較長(zhǎng)，且聲源并非靜止的，難以識(shí)別判定不同發(fā)射階段所產(chǎn)生的次聲信號(hào)。因此，本文提出一種基于短時(shí)能量突變的時(shí)延估計(jì)方法對(duì)近場(chǎng)火箭發(fā)射次聲信號(hào)進(jìn)行定位。

1 火箭發(fā)射信號(hào)時(shí)延估計(jì)

基于 TDOA的聲源定位方法是利用傳感器陣列各個(gè)陣元接收聲信號(hào)因其傳輸過程中傳播路徑不同而引起的時(shí)間差，然后結(jié)合陣列的幾何關(guān)系和時(shí)延量估計(jì)出聲源的位置[16]，其包含時(shí)延估計(jì)和定位解算兩部分，該定位算法的關(guān)鍵在于獲取精確的時(shí)間差，利用時(shí)間差可以進(jìn)一步確定聲源的速度、方位和距離等參數(shù)。

1.1 信號(hào)特征分析

對(duì)某次三元陣采集到的火箭次聲信號(hào)先進(jìn)行去除大氣本底壓力處理后，再采用0.01～20 Hz的帶通濾波器對(duì)信號(hào)進(jìn)行去噪，得到的信號(hào)波形如圖1所示，其中圖1（a）、1（b）、1（c）分別表示1、2、3號(hào)傳感器采集到的次聲信號(hào)進(jìn)行濾波后的波形圖，橫坐標(biāo)表示的是次聲信號(hào)到達(dá)傳感器的時(shí)刻。從圖1可以看到，隨著時(shí)間的增加，次聲信號(hào)波形逐漸增強(qiáng)，有幾個(gè)起伏，然后慢慢消失，波形持續(xù)時(shí)間長(zhǎng)、強(qiáng)度高，這符合典型的液體燃料發(fā)動(dòng)機(jī)產(chǎn)生的次聲信號(hào)波形特征[17]，與實(shí)際情況相符。

次聲信號(hào)的頻譜圖如圖2所示，此次火箭發(fā)射次聲信號(hào)的頻率范圍是 4～12 Hz。次聲源到達(dá)三個(gè)傳感器的傳播路徑存在差異，在傳播過程中的衰減程度不一致，且火箭發(fā)射次聲源一直處于運(yùn)動(dòng)狀態(tài)，多普勒效應(yīng)會(huì)引起次聲信號(hào)的幅值和相位的偏移，每個(gè)傳感器接收到的次聲信號(hào)的頻帶范圍存在差異。

圖2 相應(yīng)的次聲信號(hào)的頻譜圖Fig.2 The corresponding spectrum diagrams of the infrasound signals

由于火箭次聲信號(hào)持續(xù)時(shí)間長(zhǎng)達(dá)3～4 min，且聲源位置相對(duì)于傳感器陣列位置一直在發(fā)生變化，繼續(xù)采用峰值檢測(cè)或者廣義互相關(guān)法求時(shí)間差的方法誤差較大，所以本文提出一種基于信號(hào)短時(shí)能量突變的時(shí)延估計(jì)算法，首先獲得接收信號(hào)的各個(gè)傳感器之間的時(shí)間差，再對(duì)其進(jìn)行位置解算。

1.2 火箭發(fā)射次聲信號(hào)短時(shí)能量突變

當(dāng)火箭發(fā)射產(chǎn)生的次聲信號(hào)到達(dá)各個(gè)傳感器時(shí)，傳感器輸出的次聲信號(hào)能量會(huì)發(fā)生突變。一般地，無次聲事件發(fā)生時(shí)，次聲傳感器接收到的信號(hào)能量幅度較低且波動(dòng)起伏不大。因此，當(dāng)傳感器接收到火箭發(fā)射產(chǎn)生的次聲信號(hào)時(shí)，傳感器輸出的信號(hào)能量幅值會(huì)發(fā)生突變且幅值較大，即火箭次聲信號(hào)到達(dá)傳感器的時(shí)間點(diǎn)就是信號(hào)能量發(fā)生突變的位置。實(shí)際上，傳感器接收到的次聲信號(hào)不僅包含火箭發(fā)射產(chǎn)生的次聲信號(hào)，還包含其他噪聲信號(hào)，如風(fēng)噪聲等，有用信號(hào)與噪聲信號(hào)的疊加導(dǎo)致不能準(zhǔn)確地確定火箭發(fā)射次聲信號(hào)能量突變的位置，即不能準(zhǔn)確確定次聲信號(hào)到達(dá)傳感器的時(shí)間點(diǎn)。次聲探測(cè)站的傳感器接收到的信號(hào)在一定程度上是經(jīng)過玫瑰狀降噪管過濾了一部分噪聲后的信號(hào)，因此雖然還有其他無用信號(hào)的存在，但仍然可以認(rèn)為火箭發(fā)射次聲信號(hào)到達(dá)各個(gè)傳感器的時(shí)間大約在信號(hào)能量發(fā)生突變的位置。信號(hào)能量發(fā)生突變位置的附近區(qū)域相對(duì)更容易找到，因此首先可以先找到包含信號(hào)到達(dá)傳感器時(shí)間點(diǎn)的該區(qū)域。火箭發(fā)射產(chǎn)生的次聲信號(hào)持續(xù)時(shí)間較長(zhǎng)，通過對(duì)次聲信號(hào)進(jìn)行加窗計(jì)算，計(jì)算窗內(nèi)的信號(hào)能量（即短時(shí)能量），先找出次聲信號(hào)能量發(fā)生突變的區(qū)域[18]。

設(shè)第i個(gè)時(shí)間窗內(nèi)次聲信號(hào)的短時(shí)能量計(jì)算公式[18]為

式中，s（·）表示火箭發(fā)射產(chǎn)生的次聲信號(hào)；w（·）為窗函數(shù)，此處窗函數(shù)選為矩形窗。本文采集到的火箭發(fā)射次聲信號(hào)數(shù)據(jù)較長(zhǎng)，因此選取窗長(zhǎng)為1 000。定義矩形窗函數(shù)的窗長(zhǎng)為N，其函數(shù)表達(dá)式為

根據(jù)式（2）計(jì)算得到次聲信號(hào)的短時(shí)能量序列{Ei,i=0,1,2,…}后，可以大致看出信號(hào)的能量變化情況。假設(shè)次聲信號(hào)首先在第I個(gè)時(shí)間窗內(nèi)發(fā)生能量突變，計(jì)算該時(shí)間窗內(nèi)的次聲信號(hào)的瞬時(shí)能量：

圖3為三個(gè)傳感器接收到的火箭發(fā)射產(chǎn)生的次聲信號(hào)根據(jù)上述過程計(jì)算的短時(shí)能量的結(jié)果。

圖3 相應(yīng)的次聲信號(hào)的短時(shí)能量曲線Fig.3 The corresponding short-time energy curves of the infrasound signals

圖4為對(duì)火箭發(fā)射次聲信號(hào)的短時(shí)能量值進(jìn)行相鄰能量點(diǎn)之間計(jì)算差值后的能量突變幅度的曲線。可以明顯地看出，次聲信號(hào)的能量變化存在多個(gè)極值，即火箭發(fā)射過程中，有多個(gè)階段產(chǎn)生了次聲信號(hào)。

圖4 短時(shí)能量突變幅度曲線Fig.4 The amplitude curves of short-time energy mutation

圖5為傳感器接收到的信號(hào)首次發(fā)生能量突變的時(shí)間窗內(nèi)的次聲信號(hào)能量變化幅值情況的曲線。即火箭發(fā)射次聲信號(hào)到達(dá)傳感器的時(shí)間點(diǎn)包含在該時(shí)間窗內(nèi)。

圖5 短時(shí)能量突變時(shí)間窗內(nèi)能量變化曲線Fig.5 Energy change curves within the time window of short-time energy mutation

1.3 火箭次聲信號(hào)時(shí)間差估計(jì)

設(shè)兩兩傳感器之間的時(shí)間差為tij（i,j=1,2,3），若tij＜0，說明次聲波信號(hào)先到達(dá)第i號(hào)傳感器，后到達(dá)第j號(hào)傳感器。理想情況下，閉環(huán)三元陣的時(shí)間差閉環(huán)和等于0，即傳感器之間的延時(shí)量滿足

而在實(shí)際情況中，由于噪聲和大氣傳播路徑的復(fù)雜性，閉環(huán)系統(tǒng)延時(shí)量之和只能趨于0[19-20]，即：

根據(jù)圖5可知，信號(hào)能量突變幅度在整個(gè)信號(hào)持續(xù)時(shí)間范圍內(nèi)出現(xiàn)多個(gè)極大值，根據(jù)各路信號(hào)能量突變幅度設(shè)置峰值檢測(cè)的條件獲取時(shí)間差。本文設(shè)置1號(hào)傳感器最小峰值高度為1.5，2號(hào)傳感器最小峰值高度為1.5，3號(hào)傳感器最小峰值高度為0.4，檢測(cè)峰值的時(shí)間間隔為 25個(gè)采樣點(diǎn)。可得到多個(gè)滿足條件的位置，采用第一個(gè)位置為信號(hào)到達(dá)時(shí)間點(diǎn)，進(jìn)而求解出本文中的次聲事件根據(jù)上述方法得到的時(shí)間差為

計(jì)算t12+t23+t31=0，滿足條件。

2 近場(chǎng)次聲源的三站定位方法

對(duì)采集到的次聲信號(hào)進(jìn)行頻譜分析，火箭次聲波頻率集中在4～12 Hz之間，相應(yīng)波長(zhǎng)λ約為27～81 m，由于聲源S與陣列之間的距離d滿足：

此時(shí)聲源S處于近場(chǎng)范圍，需要按照球面波來計(jì)算。

2.1 基于時(shí)間差的三站定位原理

基于時(shí)間差的次聲定位技術(shù)主要是利用到達(dá)任意兩個(gè)傳感器之間的時(shí)間差，并結(jié)合傳感器陣列之間的幾何關(guān)系[21]，進(jìn)而確定聲源位置。

三元陣中傳感器的位置關(guān)系如圖6所示。根據(jù)三元陣的經(jīng)緯度坐標(biāo)可計(jì)算出三個(gè)傳感器之間的距離：

圖6 三元陣位置示意圖Fig.6 The position diagram of tripartite array

建立如圖7所示的次聲三元陣，以其中一個(gè)陣元為原點(diǎn)，正北方向?yàn)閅軸正向、正東方向?yàn)閄軸正向建立平面直角坐標(biāo)系。本文中，A、B、C三點(diǎn)表示3號(hào)、1號(hào)、2號(hào)傳感器所在位置。

圖7 三元次聲陣定位原理圖Fig.7 Locating principle diagram of the ternate infrasound array

假設(shè)某一火箭發(fā)射事件發(fā)生后，三個(gè)傳感器接收到次聲信號(hào)的時(shí)間為：t1、t2和t3，若A點(diǎn)最先接收到次聲信號(hào)，以B點(diǎn)為圓心，聲速vt13為半徑畫圓，其中v為聲速，從A點(diǎn)向該圓畫兩條切線，連接B和兩個(gè)切點(diǎn)做法線，這兩個(gè)法線方向即為可能的來波方向；同理，以C點(diǎn)為圓心，vt23為半徑畫圓，從A點(diǎn)向該圓做兩條切線，連接C和兩個(gè)切點(diǎn)做法線，這兩個(gè)法線方向即為可能的另外兩個(gè)來波方向；兩兩來波方向相交于四個(gè)點(diǎn)S1、S2、S3和S4，這4個(gè)點(diǎn)即為可能的發(fā)射點(diǎn)。

根據(jù)式（8）的陣元間距，再結(jié)合正弦定理即可得到三元陣的夾角以及相對(duì)于正北方向的夾角：∠ABC、∠ACB、∠BAC、LAB、LBC、LAC、α和β。假如次聲波最先到達(dá)3號(hào)傳感器即A點(diǎn)，A點(diǎn)和B點(diǎn)的時(shí)間差獲得來波方向?yàn)锽S1和BS4，A點(diǎn)和C點(diǎn)間的時(shí)間差獲得來波方向?yàn)镃S1和CS4。AB的來波方向和AC的來波方向相交得到4個(gè)點(diǎn)：S1、S2、S3和S4，這四個(gè)點(diǎn)即為可能的發(fā)射點(diǎn)（v為聲速）：

由此可得S4的方位角為：φ=-1 80°+∠S4AB+α。此處φ表示目標(biāo)點(diǎn)與正北方向的夾角，設(shè)γ為目標(biāo)點(diǎn)與正東方向的夾角，得到最終的方位角為：γ=90°-φ。

同理，可得到其余可能的發(fā)射點(diǎn)的方位角和距離。從實(shí)際工程應(yīng)用方面考慮，可排除在陣內(nèi)的發(fā)射點(diǎn)S1。再將得到的其余三個(gè)可能發(fā)射點(diǎn)到三個(gè)陣元的距離長(zhǎng)短與次聲到達(dá)時(shí)間進(jìn)行邏輯判斷，若某發(fā)射點(diǎn)到1號(hào)傳感器的距離大于到2號(hào)傳感器的距離，而t1＜t2，從邏輯上排除該發(fā)射點(diǎn)。從排除后剩下的點(diǎn)中計(jì)算發(fā)射點(diǎn)聲波傳輸?shù)紹點(diǎn)和C點(diǎn)的實(shí)測(cè)波程差與計(jì)算得到的波程差中選擇誤差最小的，即為實(shí)際發(fā)射點(diǎn)。

2.2 定位算法實(shí)現(xiàn)框圖

采用編程語言進(jìn)行定位算法的實(shí)現(xiàn)，框圖如圖8所示。

圖8 三元次聲陣的定位算法流程圖Fig.8 Flowchart of locating algorithm of the ternate infrasound array

2.3 定位結(jié)果

采用第1節(jié)得出的時(shí)延量代入上述的定位方法進(jìn)行試驗(yàn)，根據(jù)最小波程差得到聲速 v=325.5 m·s-1，進(jìn)行定位計(jì)算距離A點(diǎn)為54.412 km，與正東方向的夾角為118.984 3°。與實(shí)際位置對(duì)比，定向誤差在2°以內(nèi)，距離誤差約為3.17%。

3 結(jié) 論

從火箭發(fā)射次聲信號(hào)波形看出，火箭發(fā)射次聲信號(hào)持續(xù)時(shí)間較其他次聲事件波形較長(zhǎng)，且由于火箭發(fā)射是一種運(yùn)動(dòng)次聲源，采用峰值檢測(cè)或廣義互相關(guān)等方法進(jìn)行時(shí)延估計(jì)誤差較大。因此本文采用了一種基于短時(shí)能量突變進(jìn)行時(shí)延估計(jì)，并采用近場(chǎng)的三站定位方法進(jìn)行定位解算，求解結(jié)果與實(shí)際情況方位角相差在2°以內(nèi)，距離誤差約為3.17%，定位精度較高。