殷 昊，龐新良，張震川

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

0 引 言

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

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

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

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

1 火箭發(fā)射信號時延估計

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

1.1 信號特征分析

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

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

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

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

1.2 火箭發(fā)射次聲信號短時能量突變

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

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

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

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

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

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

圖4為對火箭發(fā)射次聲信號的短時能量值進(jìn)行相鄰能量點(diǎn)之間計算差值后的能量突變幅度的曲線?？梢悦黠@地看出，次聲信號的能量變化存在多個極值，即火箭發(fā)射過程中，有多個階段產(chǎn)生了次聲信號。

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

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

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

1.3 火箭次聲信號時間差估計

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

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

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

計算t12+t23+t31=0，滿足條件。

2 近場次聲源的三站定位方法

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

此時聲源S處于近場范圍，需要按照球面波來計算。

2.1 基于時間差的三站定位原理

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

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

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

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

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

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

根據(jù)式（8）的陣元間距，再結(jié)合正弦定理即可得到三元陣的夾角以及相對于正北方向的夾角：∠ABC、∠ACB、∠BAC、LAB、LBC、LAC、α和β。假如次聲波最先到達(dá)3號傳感器即A點(diǎn)，A點(diǎn)和B點(diǎn)的時間差獲得來波方向?yàn)锽S1和BS4，A點(diǎn)和C點(diǎn)間的時間差獲得來波方向?yàn)镃S1和CS4。AB的來波方向和AC的來波方向相交得到4個點(diǎn)：S1、S2、S3和S4，這四個點(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。再將得到的其余三個可能發(fā)射點(diǎn)到三個陣元的距離長短與次聲到達(dá)時間進(jìn)行邏輯判斷，若某發(fā)射點(diǎn)到1號傳感器的距離大于到2號傳感器的距離，而t1＜t2，從邏輯上排除該發(fā)射點(diǎn)。從排除后剩下的點(diǎn)中計算發(fā)射點(diǎn)聲波傳輸?shù)紹點(diǎn)和C點(diǎn)的實(shí)測波程差與計算得到的波程差中選擇誤差最小的，即為實(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é)得出的時延量代入上述的定位方法進(jìn)行試驗(yàn)，根據(jù)最小波程差得到聲速 v=325.5 m·s-1，進(jìn)行定位計算距離A點(diǎn)為54.412 km，與正東方向的夾角為118.984 3°。與實(shí)際位置對比，定向誤差在2°以內(nèi)，距離誤差約為3.17%。

3 結(jié) 論

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