淺層震源定位中高精度時間差測量方法

李 劍，賀 銘，韓 焱，雍順寧

(1.中北大學信息探測與處理山西省重點實驗室，山西 太原 030051 ；2.中國兵器工業試驗測試研究院，陜西 華陰 714200)

0 引言

地下淺層分布式震源定位技術是在地下震源近場監測區域內，將大量傳感器節點埋設在地下不同深度，通過自組織和多跳方式構成的無線網絡，協作地感知、監測、采集、處理以及傳輸震源產生的震動信號，并通過分析各個節點收集到的震動信息的特征，以實現震源的定位[1-3]。與地震、煤礦深層開采、石油勘探等大區域、大當量、大深度震動監測相比較，地下淺層分布式震源定位具有以下特點：1)分布的范圍相對較小，屬于小區域監測[4]；2)地下震動的深度較淺，一般不超過100 m，屬于地下淺層震動；3)傳輸介質的密度更加不均勻，往往包括浮土、巖石、沙石、密實土壤等多種類型[5]；4)對震動位置的定位精度要求更高，往往在1 m以內[6]。地下震源定位技術屬于近場震源定位技術，即在震源的近區進行定位。該技術可以有效地解決工程爆破定向拆除、地下毀傷、兵器試驗場炮彈炸點定位等各種震源定位問題，是目前地下空間定位問題中研究的一個熱點。

地下淺層分布式震源定位目前主要借鑒地震定位方法，利用地震檢波器拾取震動信號，通過提取震動波到達各節點的時間差，構建TDOA(Time Difference of Arrival)定位方程實現定位[7-9]。其中，時間差的提取精度直接關系到震源定位精度[10]，在空中爆炸、水中爆炸炸點定位中，由于沖擊波信號強度大、上升速度快、頻率響應高，一般通過對節點間獲取的第一個脈沖信號進行相關算法即可得到有效地時差信息。然而在地下爆炸定位中，與上述的沖擊波相比，爆炸近場產生的震動波具有峰值上升速度慢、波形復雜、多波形混疊現象嚴重等特點，最初的脈沖峰值點并不一定是全波最強的峰值點[11]，因此在地下炸點定位過程中，采用傳統的時差測量算法，即利用節點間獲取地震動信號全波形的最高峰值點或最初峰值點獲取時差信息時誤差較大[12-13]。本文針對此問題，提出了淺層震源定位中高精度時間差測量方法。

1 高精度時間差信息測量原理

地下爆炸產生的震動波是一種寬頻帶、多譜能的信號。由于土介質的物理特性近似帶通濾波器，隨著震動波在土壤介質的傳播，高頻成分逐漸衰減，某些低頻成分處于土壤通帶內，可以實現遠距離傳播[14]，這部分頻率成分稱為爆破震動信號的優勢頻率。

由傅里葉定理可知，一個單分量的時間域信號可以看作是一系列不同頻率正弦信號的疊加，對于爆破震動信號的優勢頻率的某一單頻信號而言，傳播過兩個節點的相位差信息能夠表征該頻率成分下的時差信息，雖然利用初至波到達時間可以得到節點間的時差信息，但相位差的分辨率是對應信號周期的360倍，因此結合到達時間信息，通過相位差能夠極大地提高時差分辨率。時間差的提取流程分為三個部分，如圖1所示，包括提取節點間的相位差信息、校準相位差信息，提取時間差信息。

圖1 時間差提取的流程圖Fig.1 Time difference extraction flow chart

1.1 基于譜分解技術的相位差信息提取方法

相位差的提取步驟如圖2所示，分為基于ST的時頻分析、譜分解和求取相位差。

圖2 相位差的提取流程圖Fig.2 Phase difference extraction flow chart

第一步 通過時頻域分析得到節點群共有的P波優勢頻帶。

通過ST變換將各節點獲取的時域信號轉換到時頻域中，利用震動波中P波、S波以及面波等波形成分的特點，識別出P波的優勢頻帶，同時得到節點群共有的P波優勢頻帶，取最高優勢頻率作為下一步譜分解的對象。

第二步 對P波最高優勢頻率進行譜分解。

假設P波最高優勢頻率為f0，根據傅里葉變換，其頻域表達式為式(2)：

x(t)=sinw0t

(1)

F(sinw0t)=-jπδ(w+w0)+jπδ(w-w0)

(2)

由式(2)可知，頻率為w0離散時間正弦信號的頻譜為點w0和-w0的兩處沖激脈沖，為了抽取最高優勢頻率w0的譜信號及時域波形，結合式(2)得出了式(3)：

(3)

第三步 采用相關法提取一個周期內的相位差信息。

設P波最高優勢頻率f0，經式(3)譜分解后，對應的信號為x1(t)，x2(t)：

x1(t)=A1sin(wt+φ1)

(4)

x2(t)=A2sin[w(t+τ)+φ2]

(5)

(6)

由式(6)可以得到兩節點在優勢頻率f0處的相位差信息，但相位差的范圍僅限制在[0,π]之間，存在相位模糊的問題，同時相位差信息不能夠表示兩信號間相位的超前和滯后關系，因此采用初至時間和系統相位信息等先驗知識對相位差信息進行修正。

1.2 基于初至波到達時間和系統相位信息的相位差校準方法

相位差的校準主要分為兩個步驟，包括基于初至波到達時間的相位校準方法和基于系統固有相位特征的相位校準方法，如圖3所示。

圖3 相位差校準示意圖Fig.3 Phase difference calibration diagram

第一步 基于初至波到達時間的相位校準方法。

初至波到達時間，能夠表征P波到達各個節點的先后順序即相位關系，通過到達時間差能夠解決相位模糊的問題，從而修正相位差信息。

設節點i的初至波到達時間為ti，節點j的初至波到達時間為tj，通過1.1節的步驟得了對應兩節點的相位差為φij，φij∈[0，π]。設φij′為修正后的相位差，tφij為修正后的相位差對應的時間差信息。按照波動理論可得式(7)、式(8)，其中N為該頻率成分周期重復的次數，是主要修正的變量；節點i、j對應的時間差信息為式(9)。

φij′=φij+2πN

(7)

(8)

Δtij=ti-tj

(9)

按照幾何關系和波場傳播理論，修正后的相位差對應的時差信息應滿足式(10)：

(10)

結合式(7)—式(10)可得N的約束條件：

(11)

同時根據兩節點的初至波到達時間，可得出兩節點在P波最高優勢頻率的相位關系：

第二步 基于系統相位特征的二次相位校準方法在無線傳感器網絡中，各節點的相頻特性并非完全一致，因此在求取節點間的相位差時，存在不可避免的系統誤差。為了提高相位差的測量精度，利用傳感器的相頻特性進行校準。

上一節中，利用兩節點的到達時間對相位差φij進行修正得到了φij′。設φi為第i個節點在優勢頻率f0處的固有相位值，φj為第j個節點在優勢頻率f0處的固有相位值，φij″為二次修正后的相位差信息，則：

(12)

通過基于初至波到達時間的相位校準和基于系統相位特征的相位二次校準方法，解決了相位模糊的問題，同時修正了系統的相位誤差。

1.3 利用修正后相位信息計算時間差

通過式(12)，得到時間差：

tφij=φij″/2πf0

(13)

式(13)中，tφij為修正后的時間差，f0為所述的優勢頻率，φij″為的修正后的相位差。

2 算法仿真及分析

為了進一步驗證本方法的可信性，通過人工模擬的方式產生兩組爆破震動波，如圖4所示。設信號采樣率為20 kHz，采樣時間為0.6 s。第一組爆破震動波的初至時刻為0.1 s，頻率為150 Hz，衰減因子為e-i/160(i=0,1，…，n)，縱波波束模型與x軸夾角為π/6，與y軸夾角為π/3，與z軸夾角為π/3。橫波在0.12 s產生，頻率為80 Hz，衰減因子為e-i/160(i=0,1，…，n)，橫波波束模型與x軸夾角為π/3，與y軸夾角為π/3，與z軸夾角為π/6；第二組波群初至波到達時間為0.115 s，縱波頻率為150 Hz，衰減因子為e-i/150(i=0,1，…，n)，縱波波束模型與x軸夾角為-π/6，與y軸夾角為π/3，與z軸夾角為π/6。橫波在0.135 s開始線性疊加，頻率為80 Hz，衰減因子為e-i/160(i=0,1，…，n)，橫波波束模型與x軸夾角為-π/6，與y軸夾角為π/6，與z軸夾角為π/3，兩波群的縱波初始相位差為14.13 rad，噪聲源為高斯白噪聲，信噪比為35 dB。

圖4 模擬爆破震動信號Fig.4 Simulated blasting vibration signal

采用STA /LTA算法提取兩個節點的初至波到達時間，圖5為基于改進STA/LTA的識別因子圖。

圖5 識別因子圖Fig.5 Recognition factor graph

按照最大峰值原則，提取峰值因子及對應的初至波到達時間，得到了表1，由表1可知，節點1、2對應的時間差為0.015 10 s。

表1 兩節點的初至波達到時間表

通過時頻分析，抽取150 Hz信號進行譜分析，如圖6所示。

圖6 P波優勢頻率下的時域信號Fig.6 Time domain signal at P wave dominant frequency

按照式(7)得到了節點間的相位差φij=1.574 36 rad，結合初至波到達時間信息，通過式(8)得到了波形重復周期N=2，修正后的相位差信息φij′=14.134 36 rad，對應的時間差為0.015 004 5 s，具體參數如表2所示。

表2 基于相位差信息的相關測量參數

結合表1和表2，通過相位差測量的方法將節點1，2的時差測量誤差由0.6%提高到0.3%，將時間分辨率由5×10-5提高到5×10-7。采用相位測量的方法雖然增加了程序運行時間，但時間差的測量精度卻提升了2倍，有助于提高后續基于TDOA的地下炸點的定位精度。

3 試驗驗證

在中國兵器試驗測試研究院華陰基地進行了地下震源定位試驗，將3 kg TNT埋設在地下3 m處起爆，采用中北大學信息探測與處理山西省重點實驗室研發的分布式地下震源定位系統進行定位，采樣率為20 kHz,以傳感器節點6與傳感器節點10獲取的震動信號為例，驗證提取時間差的流程。

由圖7可知，節點6處的三分量震動信號初至波呈脈沖式，波群持續時間短，上升速度快；節點10的三分量震動信號初至波上升速度緩慢，波群持續時間長。分別對上述兩個傳感器的X軸信號進行ST變換方式，在時頻域進行信號分析。

由圖8可知，針對第一發信號而言，直達P波的時頻域特征比較明顯，其中節點6的直達P波的有效頻帶在180～410 Hz，能量最強，相比之下，節點10的直達P波有效頻帶開始縮減，能量開始衰減，有效頻帶下降到了180～280 Hz，瞬時能量下移到S波區。通過圖8獲取節點6，10的瞬時優勢頻率，結合STA/LTA方法，提取信號的初至波到達時間[15]，如圖9所示。

圖7 兩節點時域波形圖Fig.7 Two-node time domain waveform

圖8 節點群的時頻圖Fig.8 Time-frequency diagram of a node group

圖9 初至波識別因子圖Fig.9 First arrival wave identification factor graph

通過分析識別因子圖，得到了初至波到達時間如表3所示。

提取初至波到達時間后，結合時頻域分析的結果，抽取兩節點在280 Hz對應的時域波形，如圖10所示。

表3 初至波到達時間相關參數表

圖10 兩節點在優勢頻率下的時域波形圖Fig.10 First arrival wave identification factor graph

利用譜分解技術和相關算法，得到了兩節點的相位差，同時利用初至波到時消除相位模糊，得到表4。

表4 基于相位差信息的時差測量表

同時采用直接互相關法、Roth脈沖沖擊響應法、平滑相關變換(SCOT)法、相位變換法(PHAT)、HB加權法和最大似然法(Hannan-Thomson)等6種廣義互相關法進行時間差測量[16-18]，結果如圖11所示，其中橫坐標代表采樣點數，縱坐標代表歸一化后的頻譜峰值。

圖11 基于廣義相關法的時間差測量結果圖Fig.11 Based on the generalized correlation method of time difference measurement result map

對圖11中的各種廣義相關譜峰值進行搜索，并將峰值對應的采樣點數轉換為時間差信息，如表5所示。

表5 基于廣義相關法的時差測量結果

由表5可知，采用廣義互相關中的不同方法，得出的時差信息變化較大，主要是因為爆炸近場產生的震動波具有峰值上升速度慢、波形復雜、多波形混疊現象嚴重等特點，最初的脈沖峰值點并不一定是全波最強的峰值點，因此，檢測出的峰值存在誤判現象，雖然平滑相關變換(SCOT)法、相位變換法(PHAT)、HB加權法和最大似然法(Hannan-Thomson)等檢測到的峰值較為明顯，但時差信息偏差較大。

結合表4和表5可知，雖然基于初至波到時的時差提取方法和基于廣義互相關的時差提取方法，時差測量分辨率均為ADC采樣周期，由于初至波到時提取方法利用的是震動波最初起跳點，不存在時差信息的誤判問題，其性能要優于基于廣義互相關的時差提取方法。

而本文提出的相位測量方法是在初至波到時提取方法上的改進，其將時間差的分辨率提高了2個數量級，有效地提高了時間差測量的精度，為后續實現高精度的微震定位提供了可靠的數據基礎。

4 結論

本文提出了淺層震源定位中高精度時間差測量方法，該方法利用爆炸群波的寬譜特征信息，提高時差的測量精度。數值仿真和試驗驗證表明，與基于初至波到時的時差測量方法相比，基于群波相位差信息的時差測量方法得到的時間差分辨率、精度明顯提高，在地下空間定位領域具有一定的工程應用價值。

由于在時差測量過程中，確定有效的頻率成分是實現高精度時差測量的關鍵，但在震源定位時，使用的傳感器數量較多，頻率成分的提取速度較慢，因此，如何快速準確地提取震動場的優勢頻率是實現下一步實時震源定位的重要研究方向。