孟 銘，楊 明，李 劍，韓 焱

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

0 引言

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

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

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

1 直達(dá)P波偏振角度提取原理

由地震波理論可知,同一類型的波引起空間質(zhì)點(diǎn)的運(yùn)動,其軌跡近似于橢球體，如圖1所示[7]。按照概率統(tǒng)計學(xué)方法，通過統(tǒng)計直達(dá)P波一段時間內(nèi)主運(yùn)動事件，即可得到該波主運(yùn)動方向。

圖1 笛卡爾坐標(biāo)系下的空間質(zhì)點(diǎn)運(yùn)動橢球示意圖Fig.1 Schematic diagram of spatial particle motion ellipsoid in Cartesian coordinate system

目前主要采用SCM(標(biāo)準(zhǔn)協(xié)方差極化濾波)提取直達(dá)P波角度信息。設(shè)直達(dá)P波對應(yīng)三分量震動數(shù)據(jù)ui(t),(i=x,y,z)，t∈[t1～t2]。求取各分量數(shù)據(jù)的平均值μi(ξ)：

(1)

(2)

(3)

(4)

(5)

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

2 爆炸近場直達(dá)P波偏振角度提取方法的改進(jìn)

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

2.1 基于HRP-Radon直達(dá)P波信號提取

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

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

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

(6)

(7)

對式(6)和式(7)兩端進(jìn)行傅里葉變換，然后將公式轉(zhuǎn)變?yōu)榫仃囆问剑?/p>

m=LHd

(8)

d=Lm

(9)

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

(10)

式(10)中的加權(quán)矩陣Qm與m的關(guān)系表達(dá)式如下：

(11)

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

在實現(xiàn)強(qiáng)混疊條件下直達(dá)P波有效分離的基礎(chǔ)上，還需對分離P波進(jìn)行ACM極化分析，以提取P波偏振角度等極化參數(shù)，為構(gòu)建高精度DOA定位模型提供精確的方向參量信息。

2.2 基于ACM的直達(dá)P波極化角度的提取

本文在SCM算法的基礎(chǔ)上，利用爆炸波時變特性瞬時極化，使時窗長度自適應(yīng)調(diào)整到時窗內(nèi)信號的最小周期，構(gòu)建了ACM(自適應(yīng)協(xié)方差極化濾波)偏振角度提取模型。

(12)

(13)

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

(14)

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

3 算法仿真

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

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

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

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

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

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

由圖4(b)和圖5(b)結(jié)合可以得出P波合理保留區(qū)域，如圖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軸的瞬時夾角對比圖。其中信號有效區(qū)域由信號初至波到時時刻開始到設(shè)定時間0.08 s結(jié)束，信號無效區(qū)域表示激勵產(chǎn)生的高斯白噪聲。實際數(shù)據(jù)處理時以信號初至波到時作為有效區(qū)域起始位置,以極化度瞬時變化時刻作為有效區(qū)域結(jié)束時間的劃分節(jié)點(diǎn)。通過圖10的極化度曲線對比圖可以看出分離P波偏振角度的極化度曲線[0.95，1]之間，并由表1可得，HRP-Radon變換聯(lián)合ACM算法能夠得到高精度直達(dá)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 結(jié)論

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