大氣次聲源定位算法及誤差分析

尚楚翔，滕鵬曉，呂 君，楊 軍，程 巍

(1.中國科學院聲學研究所，北京100190；2.中國科學院大學，北京101408)

0 引 言

自然事件和人類活動會伴有次聲信號產生。大氣中的次聲源主要包括3種[1]：人為制造的脈沖信號；隕石、火山噴發、地震等自然災害；對流層中的大氣擾動。由于大氣層中存在溫度和風向的分層結構，次聲因其頻率低、衰減小的特點，可以在大氣中形成次聲波導，傳播上千千米而沒有明顯的能量損失[2]，因此逐漸成為隕石、火山、核爆炸次聲源定位的有效技術和實用方法。

文獻[2-3]中根據已有的地震探測技術，提出逐步多通道相關方法(Progressive Multi-channel Correlation Method, PMCC)。通過次聲信號檢測方法和聲達時間差(Time Difference of Arrival, TDOA)算法，可以計算出不同頻段下次聲波到達陣列的方位角和俯仰角。1996年，全面禁止核試驗條約組織(Comprehensive Nuclear Test Ban Treaty Organization, CTBTO)成立，構建了國際監測系統(International Monitoring System, IMS)，將次聲作為核爆炸監測定位和當量估計的手段之一。IMS在全球設有60個次聲監測站，保證任何位置發生當量大于1 kt以上的核爆炸至少能被2個監測站探測到，并使用PMCC方法對其定位。每個監測站采用4元次聲陣列，陣列形狀為中心一點，其余三點組成三角形[4]。1996年，謝金來[5-6]建立了中國科學院聲學研究所核爆炸次聲波監測系統，包括北京香山次聲站、昆明次聲站和杭州次聲站。監測系統采用4元方陣，實現了四點定向、兩陣定位的次聲監測網絡。Kennett等[7]通過增加陣元數目解決了IMS中4元中心三角陣列對高頻信號出現混疊現象的問題，此后監測站逐步升級為8元或9元陣列。次聲陣列監測系統不僅可以在最優置信度、實時性和靈敏度三方面顯著提高核爆炸監測能力，還可以用于全球環境、人為災害和自然災害的監測[6]。

目前，CTBTO成員國研究人員廣泛采用兩個或多個陣列的地理坐標及其測得的多個方位角，通過在地球坐標系上進行測繪，定位次聲源的經緯度。由于在陣列布置時，各陣元無法保證高度一致，忽略陣元間的海拔高度差，方位角的計算結果會存在一定的誤差。本文基于最小二乘法，通過分析PMCC方法的主要流程，確定其次聲源定位的主要誤差來源，推導了在考慮和忽略陣元間高度差時計算次聲波入射角的方法，針對特定陣型仿真模擬實現各方向入射次聲波的計算誤差可視化，研究最大海拔高度差和入射角計算誤差之間的關系，并討論了入射角的計算誤差對后續定位的影響。

1 理論模型

PMCC方法主要包含次聲信號檢測方法和定位算法。其中定位算法采用TDOA算法，包括時間延遲估計(Time Delay Estimation, TDE)和基于時延的定位方法兩部分。次聲波在大氣遠距離傳播過程中，不同頻率的聲波路徑存在差異，到達陣列的入射角度不同[8]，因此在定位前會先將各陣元的接收信號通過多個帶通濾波器進行分頻處理，以實現對不同頻率次聲波的定位計算。

1.1 時間延遲估計

TDE方法大致可包括最小均方(Least Mean Square, LMS)自適應濾波器法、自適應特征值分解法、聲學傳遞函數比法等[9]。PMCC方法采用廣義互相關(Generalized Cross- Correlation, GCC)[10]算法估計次聲信號到達陣元的時間延遲。其原理相當于將接收信號進行前置濾波，適當消除噪聲及干擾的影響后再進行相關運算，取其峰值進行時延估計，以提高時延估計的精度，具體流程如圖1所示。圖1中FFT表示快速傅里葉變換，IFFT表示快速傅里葉逆變換。

圖1 廣義互相關估計時延流程圖Fig.1 Flow chart of generalized cross-correlation estimation of time delay

對于兩個陣元接收到的信號x1(n)和x2(n)，簡化模型后可以表示為

其中：s(n)為次聲源信號，αi(i=1,2)表示次聲信號到達陣元i的衰減系數，vi表示陣元i處的噪聲信號。假設噪聲信號與聲源信號互不相關，噪聲信號之間互不相關，互功率譜表示為

廣義互相關定義為

其中：ψ(ω)為加權函數，相當于前置濾波作用。文獻[10]給出了多種通用的加權函數，PMCC方法主要使用HB(Hassab-Boucher)加權函數，其表達式為

在實際應用中，使用信號頻譜相乘估計功率譜，相關函數的最大峰會被弱化，同時噪聲的相關性也會對結果產生影響。因此布置次聲陣列通常選擇背景噪聲較小的位置，并且較大的陣元間距可以減小噪聲相關性。

1.2 信號一致性閾值檢測

目前CTBTO主要采用8元或9元陣列構建IMS監測站，由大尺寸的5元陣列嵌套小尺寸的3元或4元陣列組成。根據次聲波長，陣元間距通常設定在1～3 km。較大的陣元間距可以減少空間噪聲的相關性，提高廣義互相關加權的濾波性能。同時可以增加聲波到達陣元的時延，減小讀取時延數據的測量誤差對后續定位的影響。但陣元間距太大會使信號本身的互相關性下降，并且當陣元間距超過信號高頻部分的半波長時，時延計算會出現空間混疊現象?？梢钥醋鲀蓚€陣元接收了不同周期內的高頻信號，互相關函數會出現多個峰值，導致時延估計錯誤。

次聲陣列Rn的任意3元子陣列(i,j,k)接收次聲源信號的時延滿足閉合時間關系：

由于背景噪聲的相位具有隨機性，并且讀取時延數據存在測量誤差，閉合時間關系無法嚴格滿足。修正閉合時間關系更改為[11]

當陣列的修正閉合時間關系在閾值cthreshold范圍內，說明陣元位置有效，否則說明某陣元位置干擾過大導致時延計算出現了較大誤差。

為檢測次聲源信號并且避免高頻信號出現混疊，檢測過程會從小尺寸的3元陣列出發，逐步擴大陣型以提高計算精度。這種方式體現PMCC方法的逐級性，如圖2所示。

圖2 PMCC方法的逐級檢測Fig.2 Progressive detection by PMCC method

在監測過程中，需先排除連續零信號的傳聲器。當初始三元陣的接收信號時延滿足修正閉合時間關系時，說明存在次聲源，生成次聲源監測子網絡。逐個添加其他陣元，計算子網絡的修正閉合時間關系并進行閾值判斷，如果滿足閾值范圍則保留該陣元。遍歷全部陣元后利用監測網絡進行后續定位。

1.3 基于時延的定位方法

在獲取不同頻率次聲信號到達各陣元的時間延遲后，可通過幾何定位法或目標函數搜索法計算入射角?？紤]到次聲遠距離傳播的路徑遠大于陣列尺寸，聲波到達陣列可視為平面波，入射角采用最小二乘法利用多通道時延進行計算。

圖3 任意陣列示意圖Fig.3 Schematic diagram of arbitrary array

圖3為任意陣列示意圖。設入射信號方向為r=-[sinθcosφ,sinθsinφ,cosθ]，其中方位角為φ(與x軸夾角)，俯仰角為θ(與z軸夾角)；陣元m的坐標為pm=(xm,ym,zm)；c表示當地聲速；τij表示陣元i相對于陣元j的延遲，正值表示信號先到陣元j、后到陣元i，τji可以表示為[12]

在實際測量中，陣元的地理坐標、采集設備和環境因素都有可能造成式(7)不嚴格相等，用誤差表示為

n個陣元的誤差線性方程組可寫為

其中：e為誤差向量，A為位置矩陣，b為聲程差向量。在最小二乘意義下，估計誤差的模的平方和

式(10)取最小值，可得到：

視速度V是表示次聲信號的重要參數，由本地聲速和俯仰角決定：

為了進一步濾除干擾信號，在時頻域對計算結果進行聚類。

接收信號經過連續時間窗和分頻處理后，通過時延定位可以得到由獨立單元組成的時頻圖[13]，每個獨立單元成為聚類的檢測點，如圖4(a)所示。每個檢測點Pi由時間ti、頻率fi、視速度Vi和方位角φi表示，應用權重歐幾里得距離表示相鄰時間窗或頻帶之間檢測點的距離[11]：

其中：σ表示對應參數的權值。設定閾值dmax，當時頻圖中兩個相鄰檢測點的距離在閾值范圍內，這兩個檢測點同屬于一類信號。當類中檢測點數量大于建立族所需要的最小數目時，建立一個PMCC族。圖4(a)為聚類前的時頻圖，灰色檢測點距離相鄰的黑色檢測點較遠，在進行聚類后形成只包含黑色檢測點的PMCC族，如圖4(b)。建立PMCC族有效地排除了外界擾動信號，提高定位算法的抗干擾能力。

圖4 聚類前后檢測信號的變化Fig.4 Changes in detection signals before and after clustering

2 誤差分析

在實際應用中，陣元的地理位置使用經緯度表示，不考慮其海拔高度，陣元m的坐標取為pm=(xm,ym)。此時式(9)中A變為n×2維矩陣：

此時誤差向量表示為

同樣利用最小二乘法得到入射角度：

顯然，式(16)與式(11)的結果不完全相等，這是因為忽略陣元海拔高度差會影響入射角的計算結果。

仿真采用的4元次聲陣列示意圖如圖5所示，陣元間距為1～3 km，最大海拔高度差為0.04 km。由于誤差只與陣元位置坐標和信號入射角度有關，利用式(11)計算各方向次聲信號到達各陣元的真實時延，同時是次聲陣列采集到的時延。利用式(16)求解在該時延下的入射角，進而對比各方向的誤差結果?？梢灶A測，信號入射角度越接近垂直入射，陣元海拔高度差所產生的時延在總時延中所占的比例越大，忽略高度差所產生的誤差也就越大。

圖5 仿真采用的中心一點其余三點組成三角形的4元次聲陣列Fig.5 The triangular quaternary infrasound array composed of one element at the center plus other three elements

圖6為取方位角0°～360°，步長1.8°的次聲波入射到圖5陣列，不考慮陣元高度差時產生的方位角誤差。在圖6(a)、6(b)中，當入射信號俯仰角接近0°，即信號接近垂直入射時，方位角誤差極大，甚至相差180°，導致定位反向。隨著俯仰角增加，方位角誤差逐漸減小。接近水平入射的信號方位角誤差范圍可保證在2°以內。考慮到大氣次聲遠距離傳播到達陣列的俯仰角基本大于30°，將俯仰角取值范圍改為30°到90°。圖6(c)中，俯仰角為30°，方位角為90°或270°附近時存在誤差極大值，與當前陣型構造有關。各方向方位角誤差在3°以內。3°的方位角誤差在次聲傳播水平距離500 km時，對應定位偏移誤差約為25 km，次聲傳播水平距離1 000 km對應50 km的定位偏移誤差。

圖6 不考慮陣元高度差時次聲波入射到圖5陣列產生的方位角誤差Fig.6 The azimuth errors when the infrasound wave is incident to the array shown in Fig.5 without considering element height difference

圖7為取方位角0°～360°，步長1.8°的次聲波入射到圖5陣列，不考慮陣元高度差時產生的俯仰角誤差。由圖7可見，各方向俯仰角誤差全部在4°以內。隨著俯仰角的增加，俯仰角誤差逐漸減小。在入射信號方位角為344°、俯仰角為86°附近存在誤差極大值，與當前陣型構造有關，不同陣型的極大值位置和大小不同。4°的俯仰角誤差在俯仰角為57°、聲速為344 m·s-1時，會產生14 m·s-1左右的視速度誤差。

圖7 不考慮陣元高度差時次聲波入射到圖5陣列產生的俯仰角誤差Fig.7 The pitch angle errors when the infrasound wave is incident to the array shown in Fig.5 without considering element height difference

顯然，入射角計算誤差只取決于陣元海拔的相對高度差，而非絕對高度。仍然采用圖5中的4元次聲陣列，保持陣元水平位置不變，陣元的海拔高度服從0到設定最大海拔高度差值的均勻隨機分布，在經過多次仿真試驗后，可以得到對于同一入射角度的次聲波，最大海拔高度差與入射角計算誤差的關系，圖8為方位角30°時平均俯仰角誤差和平均方位角誤差。隨著陣元最大海拔高度差的增加，入射角計算誤差呈線性增大。在最大海拔高度差為0.2 km的情況下，方位角誤差可達到5°左右。對比圖8(a)與圖8(b)可以發現，當俯仰角較小時，方位角誤差的增速高于俯仰角誤差的增速。隨著俯仰角的增加，俯仰角誤差幾乎沒有影響，與圖7相對應；但方位角誤差會急劇減小，對應圖6。

圖8 方位角為30°時平均俯仰角誤差和平均方位角誤差Fig.8 The average pitch angle error and the average azimuth angle error with the azimuth angle of 30°

在忽略陣元海拔高度差的情況下，入射信號俯仰角為0°時，方位角誤差最大，俯仰角誤差存在極大值。隨著俯仰角的增加，誤差逐漸減小。次聲陣列的陣型會對特定的入射方向產生誤差極大值。隨著陣元間最大海拔高度差的增加，入射角計算誤差呈線性增大趨勢。對于仿真的4元中心三角陣型，在限定入射方向與最大海拔高度差后，計算誤差均可保持在5°以內，但由于次聲傳播距離較遠，5°的入射角誤差水平偏移距離相當于水平傳播距離的9%，視速度誤差可達到18 m·s-1左右。

3 結 論

本文依據PMCC方法，分析了基于HB加權的廣義互相關時間延遲估計算法、信號一致性閾值檢測方法和基于最小二乘法的平面波時延定位方法，確定了忽略陣元間的海拔高度差是PMCC方法產生誤差的主要原因。推導了在忽略陣元間高度差的情況下次聲波入射角的計算方法。通過假設存在高度差的4元中心三角陣型，對各方向入射聲波因忽略陣元間海拔高度差而產生的誤差進行模擬與分析。同時在假設陣型水平位置不變的情況下，討論了最大海拔高度差與計算誤差的關系。結果表明，在次聲遠距離傳播的條件下，忽略陣元間的海拔高度差對后續定位誤差有較大影響，PMCC定位次聲源時需要考慮陣元海拔高度差的影響。