爆炸近場高精度P波角度提取方法

孟 銘，楊 明，李 劍，韓 焱

(中北大學信息探測與處理山西省重點實驗室，山西 太原 030051)

0 引言

基于偏振角度信息的DOA(Direct of Angle)定位方法是實現地下淺層震源近場定位的主要手段[1]，該方法采用三分量數據來計算P波的偏振方向，通過多個檢波器的P波偏振方向求取交匯點，從而實現震源定位。由于該方法無需求取初至時間，原則上僅需要兩組檢波器數據即可實現震源定位，因此是地下淺層空間定位領域研究的熱點。然而在爆炸點近場定位過程中，由于群波混疊嚴重、多徑干擾效應大，導致P波、S波偏振信息混疊，P波偏振角度 “提不準”，最終造成DOA定位模型“建不準”。

在直達P波角度提取算法的研究中，目前國內外學者主要采用MUSIC(Multiple Signal Classification)算法、ESPRIT(Estimating Signal Parameter via Rotational Invariance Techniques)算法[2-3]、協方差矩陣重構法和極化分析法等[4]。其中MUSIC算法、ESPRIT算法和協方差矩陣重構法均是基于均勻線性陣列構型的遠場角度估計方法，極化分析法最重要的工作就是時窗長度的選取，在信號強混疊的情況下，極化分析法難以實現時窗長度的準確選取[5-6]。

因此上述方法只適用于遠場大區域、大深度條件下P波偏振角度的提取，在近場信號強混疊情況下無法適用。針對這個問題，本文利用地下波場在Radon域的聚焦差異和P波的偏振特性，開展了基于高分辨率拋物Radon(High Resolution Parabolic Radon Transform，HRP-Radon)聯合ACM的直達P波偏振角度提取方法的研究

1 直達P波偏振角度提取原理

由地震波理論可知,同一類型的波引起空間質點的運動,其軌跡近似于橢球體，如圖1所示[7]。按照概率統計學方法，通過統計直達P波一段時間內主運動事件，即可得到該波主運動方向。

圖1 笛卡爾坐標系下的空間質點運動橢球示意圖Fig.1 Schematic diagram of spatial particle motion ellipsoid in Cartesian coordinate system

目前主要采用SCM(標準協方差極化濾波)提取直達P波角度信息。設直達P波對應三分量震動數據ui(t),(i=x,y,z)，t∈[t1～t2]。求取各分量數據的平均值μi(ξ)：

(1)

(2)

(3)

(4)

(5)

式中，θ0(t)為瞬時方位角，δ0(t)為瞬時傾角。基于SCM原理的偏振角度提取算法簡單，但是對時窗長度的選取較為敏感，并且在給定長度的時窗內求得的極化參數不具有時變特性。同時由于爆炸近場存在波形混疊嚴重、多徑干擾效應大等問題，導致無法有效和準確的提取時變信號的偏振角度信息。

2 爆炸近場直達P波偏振角度提取方法的改進

針對爆炸近場群波混疊條件下，SCM算法無法提取有效直達P波角度信息的難題，本文利用傳感器陣列的波場信息，通過HRP-Radon變換和遠程自然分離的P波振相特征在Radon域提取混疊信號中的P波信號。利用震動波的時變特性，自適應調整SCM算法中時窗長度，提取直達P波的瞬時極化特性，最終獲取精準的偏振角度信息。

2.1 基于HRP-Radon直達P波信號提取

本文提出的HRP-Radon變換結合P波波形振相特征的方法。HRP-Radon變換可根據橫縱波的特征參數(慢度、曲率等)差異將P波和S波映射為Radon域內不同位置的峰值點，實現縱、橫波的分離[9-11]。該方法首先提取P波的首振相、主振相和尾振相曲線，通過HRP-Radon變換得到這些振相曲線對應的Radon域投影，并將其作為Radon域投影區域P波提取邊界條件。然后參照信號特征邊界條件保留混疊波形中P波在Radon域有效區域，并根據P波有效投影區域實現P波反演重建。P波分離流程如圖2所示。

圖2 P波分離流程圖Fig.2 P-wave separation flow chart

首先對三分量混疊信號從時間域數據沿拋物線路徑進行積分變換到Radon域，并完成其逆變換(從Radon域變換到時間域)，具體公式如下：

(6)

(7)

對式(6)和式(7)兩端進行傅里葉變換，然后將公式轉變為矩陣形式：

m=LHd

(8)

d=Lm

(9)

采用高分辨率拋物線Radon變換，可以較大程度地降低由于空間采樣范圍有限導致的截斷效應對Radon域內信號分辨率的影響[12-13]。式(8)中的函數可優化為:

(10)

式(10)中的加權矩陣Qm與m的關系表達式如下：

(11)

共軛梯度算法占用內存小而且計算效率高，在式(10)中，使用共軛梯度算法完成矩陣(LHL+Qm(mk))的求逆過程。高分辨率拋物線Radon變換的數據重建誤差小，在變換過程中能量損耗特別低，變換后數據可以得到較完整的恢復，有效的實現P波的分離和反演重建。

在實現強混疊條件下直達P波有效分離的基礎上，還需對分離P波進行ACM極化分析，以提取P波偏振角度等極化參數，為構建高精度DOA定位模型提供精確的方向參量信息。

2.2 基于ACM的直達P波極化角度的提取

本文在SCM算法的基礎上，利用爆炸波時變特性瞬時極化，使時窗長度自適應調整到時窗內信號的最小周期，構建了ACM(自適應協方差極化濾波)偏振角度提取模型。

(12)

(13)

根據地震波理論，極化度T是評價質點在地下空間運動偏振特性的一個指標，如圖1所示。當極化度接近于1時，該質點在空間中呈線性偏振，由主特征值對應的特征向量得到的偏振角度信息可以最大程度反應該質點在地下空間的運動方向。當極化度接近于0時，該質點在空間中呈無序圓形運動，得到的偏振角度信息無法有效的衡量該質點在地下空間的運動方向。因此可以用極化度T來評價P波在地下空間中的偏振特性，即偏振角度信息的獲取精度。具體公式如下所示：

(14)

式(14)中，ρ表示主橢圓率極化參量，ρ1表示次橢圓率極化參量。

3 算法仿真

地下淺層爆炸近場的仿真模型如下：假設在地下介質為均勻水平層狀介質。建立XYZ三軸坐標系，地平面為XY平面，炸點位置為(0，30，-200)。以原點為中心，在X軸上布置21個傳感器，間距為31 m，P波傳播速度為5 500 m/s，S波傳播速度為3 800 m/s。采用雷克子波激勵方式產生相應信號，設信號采樣率為20 kHz，采樣時間為0.16 s。爆破震動波的初至時刻為0.01 s，頻率為160 Hz，衰減因子為e-i/160(i=0,1,…,n)，橫波在0.02 s產生，頻率為80 Hz，衰減因子為e-i/160(i=0,1,…,n)，噪聲源為高斯白噪聲，信噪比為35 dB,并抽取第5道數據為例進行算法仿真驗證。

圖3 仿真模型示意圖Fig.3 Schematic diagram of simulation model

通過矢量計算生成傳感器陣列的三軸數據，具體P波偏振角度信息提取方法如圖4所示。

圖4 混疊波形及Radon域投影Fig.4 Aliasing waveform and Radon domain projection

圖4(a)是模擬三分量傳感器接收到的地下爆炸產生的橫縱波混疊信號，圖4(b)為混疊信號經過HRP-Radon變換后在Radon域的投影圖。

圖5(a)是P波振相特征曲線，圖5(b)為P波振相位特征曲線在Radon域的投影區域，可作為P波Radon域有效區域提取的依據。

由圖4(b)和圖5(b)結合可以得出P波合理保留區域，如圖6(a)所示。圖6(b)為反變換重建后的P波。

圖5 P波振相特征曲線及Radon域投影Fig.5 P-wave vibration phase characteristic curve and Radon domain projection

圖7、圖8、圖9分別為混疊波形和分離P波在x，y，z軸的瞬時夾角對比圖。其中信號有效區域由信號初至波到時時刻開始到設定時間0.08 s結束，信號無效區域表示激勵產生的高斯白噪聲。實際數據處理時以信號初至波到時作為有效區域起始位置,以極化度瞬時變化時刻作為有效區域結束時間的劃分節點。通過圖10的極化度曲線對比圖可以看出分離P波偏振角度的極化度曲線[0.95，1]之間，并由表1可得，HRP-Radon變換聯合ACM算法能夠得到高精度直達P波的角度信息。

圖6 P波Radon域分離及時域重建Fig.6 P-wave Radon domain separation and time domain reconstruction

圖7 混疊波形與分離P波x軸夾角Fig.7 Inclusion angle between aliasing waveform and separated P-wave x-axis

圖8 混疊波形與分離P波y軸夾角Fig.8 Inclusion angle between aliasing waveform and separated P-wave y-axis

圖9 混疊波形與分離P波z軸夾角Fig.9 Inclusion angle between aliasing waveform and separated P-wave z-axis

圖10 混疊波形與分離P波極化度曲線Fig.10 Mixed Waveform and Separated P-Wave Polarity Curve

表1P波的偏振信息表
Tab.1 The polarization information table of P wave

理論值實際值絕對誤差提取精度vx/rad1.02081.0210.00020.9998vy/rad1.69761.6980.00040.9997vz/rad0.56780.56790.00010.9998

4 結論

本文提出了爆炸近場高精度P波角度提取方法。該方法首先利用P波、S波在Radon聚焦特性的不同，以HRP-Radon分離混疊波群中的P波信號。其次，利用爆炸近場震動信號的時變特性，自適應調整波束角度信息檢測的窗長，在SCM算法基礎上進行改進，構建基于ACM的偏振角度提取模型，以實現有效P波角度的提取。最后，利用極化度對波束角度信息的提取精度進行評價。數值仿真結果表明,上述方法能夠有效獲取爆炸近場混疊信號中P波角度信息，本文研究成果將為爆炸近場波場特征參數高精度提取提供有力的技術支撐，對地下空間近場震源定位研究具有一定的參考價值。但由于本算法利用遠場P波振相特征作為有效區域劃分邊界條件，因此在保證直達P波角度提取精度的前提下，如何優化近、遠場傳感器陣列布設是下一步研究的重點。