提高臺站地震資料信噪比的自適應極化濾波①

馬見青，李慶春

（長安大學地質(zhì)工程與測繪學院，陜西 西安 710054）

0 引言

地震波場在本質(zhì)上是不同類型、不同偏振特性的振動相互干涉和疊加的結果［1］。極化濾波是在基于波的偏振特性基礎上進行的一種空間濾波的信號處理方法，目前極化分析方法在電磁學、光學和地震勘探等領域已經(jīng)取得了一些研究成果［2］。一方面通過對地震多分量記錄的分析能夠獲得地震波的極化參數(shù)，可以實現(xiàn)各類噪聲壓制、特定波型的識別與分離，提高地震資料的信噪比，從而為后續(xù)的各種數(shù)據(jù)處理與解釋服務［3－10］；此外還可以應用極化濾波分析橫波分裂，進行各向異性研究［11－14］。對于充分利用地震資料中所包含的豐富的動力學和運動學信息，深化對地球本體或地下介質(zhì)的認識具有重要的理論意義和實際應用價值［15］。

極化分析方法在地震信號處理中的應用研究始于上世紀60年代，F(xiàn)linn［16］最早提出了極化分析方法的數(shù)學算法。以后 Montalbetti、Samson、Kanssewich、Vidale、Morozov、Zhang等［17－22］相繼在時間域內(nèi)對基于協(xié)方差矩陣的極化分析技術進行了大量研究。具體作法是對三分量地震數(shù)據(jù)sk（t）（k＝x，y，z）在一定長度的時窗內(nèi)構造協(xié)方差矩陣M（ξ）；其特征值分別為λi（i＝1，2，3），且λ1≥λ2≥λ3；Vi為λi為 對應的特征向量（圖1）；通過時窗內(nèi)協(xié)方差矩陣的特征值和特征向量構造出一系列能夠表征不同類型地震波極化特性的極化參數(shù)。這種方法原理簡單，實現(xiàn)方便，但對時窗長度的選擇依賴性較強，而且在給定長度的時窗內(nèi)求得的極化參數(shù)不具有時變特性，由此設計的濾波器的濾波效果也不夠理想。

圖1 笛卡爾坐標系下的三分量極化橢球示意圖Fig.1 schematic representation of an multi－component polarization ellipsoid in Cartesian coordinate system

鑒于上述方法的局限性，Diallo等［23］提出了基于自適應協(xié)方差矩陣的瞬時極化分析方法。該方法的時窗長度可以自適應于信號的瞬時頻率，即協(xié)方差矩陣的計算時窗長度自適應調(diào)整到時窗內(nèi)信號的最小周期，可以得到每一分量上各個采樣點的瞬時極化參數(shù)，不需要進行插值（開始和結束處各1／2時窗除外）。這樣我們就可以擺脫標準協(xié)方差矩陣極化分析中必須人為選擇時窗長度的限制。最后，通過數(shù)值模擬實驗及實際三分量臺站資料處理進一步比較和解釋了基于協(xié)方差矩陣的極化分析方法。

1 自適應協(xié)方差極化分析方法

對三分量地震數(shù)據(jù)sk（t），（k＝x，y，z）求其Hilbert變換譜（t），得到其解析信號

在時間t附近的局部信號可以利用解析信號ck（t）來近似表達：

利用多分量信號sk（t）的近似表達式（2）來構造協(xié)方差矩陣

協(xié)方差矩陣 M（t）與 M（ξ）不同，M（ξ）代表某一時窗內(nèi)的平均極化分布，而M（t）代表的是地震信號中每一個時間點的瞬時極化分布，這正是自適應協(xié)方差矩陣方法（ACM）相比標準協(xié)方差矩陣方法（SCM）的優(yōu)勢所在。

該協(xié)方差矩陣M（t）的特征值分別為λ（i＝1，2，3），且λ1≥λ2≥λ3。Vi為λi對應的特征向量。因此可以得到極化橢球體的瞬時極化參數(shù)：

瞬時極化軸

可以證明，最小的特征向量V3（t）平行于這個平面向量 ，Vp（t），V3（t）的方向余弦θx（t），θy（t），θz（t）可以確定主橢圓面在三維空間的位置：

上述得到的極化參數(shù)都是時窗內(nèi)的中心點位置的極化參數(shù)，通過將時窗沿時間軸逐個采樣點向前移動，可以得到整個時間軸（除開始和結束處各1／2時窗）。根據(jù)三分量地震數(shù)據(jù)在空間某一時窗內(nèi)的極化參數(shù)，我們就可以構造不同的極化濾波函數(shù)，保留具有某一組極化參數(shù)的波，從而達到波場分離或去噪的目的。

濾波器為［19］

2 自適應分析時窗的確定

Tkm表示t時刻M（t）的自適應時窗長度，定義為

或者可以由自適應協(xié)方差矩陣M（t）中每一個元素的瞬時頻率來定義：

從方程（2）中我們可以看到，每一分量的信號都是由周期為2π／Ωk（t）的余弦函數(shù)近似得到。為了闡明自適應時窗長度Tkm（t）的原理，我們考慮Ωx（t）＝Ωy（t）＝Ωz（t）＝Ω（t）的情況，在這種情況下Tkm可以簡化為Tkm（t）＝2πN／Ω（t），它是周期的整數(shù)倍，N為一個正整數(shù)。

通常情況下不同分量的瞬時頻率是不同的，此時方程（8）和（9）所定義的Tkm（t）為一個平均周期。相比之下方程（9）更能提供最佳的自適應時間窗。

為了分析參數(shù)N的選取對本方法的影響，現(xiàn)設計一模型來闡述。圖2（a）為一個三分量數(shù)據(jù)sx（t）、sy（t）和sz（t）（紅色實線）以及利用方程（2）計算得到的相應近似表達式（黑色虛線）。在0～0.05s時間段內(nèi)N＝1，在這個時間段內(nèi)可以看到近似得到的信號和原始模擬信號的波形相當一致。同樣在圖2（b）內(nèi)的矢端曲線圖上也可以看到矢端曲線近似為一個完整的圓形。此外，這個矢端曲線圖表明，當N＝1時極化分布不能刻畫三維空間內(nèi)比橢圓更復雜的結構。從圖2（a）可以看出，隨著N的增大信號的近似表示的精度在下降，但同時比橢圓更為復雜的曲線在矢端曲線圖上出現(xiàn)（圖2（b）～2（d））。因此在刻畫平面內(nèi)的極化屬性時選取N＝1或2就足夠了；當想刻畫三維空間里的結構更為復雜的極化屬性時，我們可以犧牲近似表達的精度，利用較大的N值。

圖2 N值對三分量信號近似表示的影響Fig.2 Influence of Non the accuracy of the approximation for signals

3 數(shù)值模擬

為了進一步比較時間域內(nèi)的標準協(xié)方差極化分析（SCM）和自適應瞬時極化分析（ACM），現(xiàn)人工合成一個三分量記錄（圖3），該記錄由5個波段構成，A段表示平穩(wěn)橢圓，B段表示旋轉橢圓，C段表示線性極化的情況，D段表示三維空間的平穩(wěn)橢圓，E段對應于旋轉橢球的情況。前三段只在X－Y平面內(nèi)分布，而后兩段在整個三維空間均存在。圖4為該三分量記錄的空間極化軌跡圖。

圖3 三分量人工合成記錄［23］Fig.3 Three components synthetic record［23］

圖4 三分量人工合成記錄的空間極化軌跡Fig.4 Hodograph of the three components synthetic record

圖5 時間域瞬時極化參數(shù)Fig.5 Instantaneous polarization attributes in time domain

圖5（a）為用標準協(xié)方差方法和自適應協(xié)方差方法得到的瞬時極化軸，我們可以看到，標準協(xié)方差方法與自適應協(xié)方差方法相比精度很低。這主要是由于標準協(xié)方差極化分析法的時窗長度很難準確地選擇，也就是說時窗長度的選取對于標準協(xié)方差極化分析法的精度有很大的影響。圖5（b）是極化率曲線，即Rmed（t）／Rmax）（t），從中可以很容易把線性極化波（rs（t）／R（t）≈0）從橢圓極化波和橢球極化波中區(qū)分出來。但標準協(xié)方差方法與自適應協(xié)方差方法相比，求得的極化率震蕩很厲害。圖5（c）是由方程（9）計算得到的自適應時間窗，可以看到時窗長度隨頻率的增大而變小，反之亦然。

值得注意的是，在信號周期和時窗長度接近的情況下，兩種方法得到的極化軸曲線吻合程度很高，曲線形態(tài)差別較大的區(qū)域說明標準協(xié)方差方法的時窗長度選擇不是最理想的。

4 實際三分量臺站資料處理

圖6是陜西省地震臺網(wǎng)2002年記錄的三分量實際地震數(shù)據(jù)，該臺站位于緯度34.25°，經(jīng)度108.95°。圖7是分別通過SCM和ACM極化濾波方法計算得到的極化主軸、極化次軸（藍線表示SCM，紅線表示ACM），通過對比兩種方法得到的極化軸曲線。可以看到，ACM對局部變化比較劇烈的信號更加敏感，而SCM實際上是ACM的平滑過程。而且，由于SCM的窗口長度固定，不適合用來刻畫周期長度小于窗口長度的地震信號。在信號周期和窗口長度接近的地方，兩條曲線的相似程度最高，在曲線形態(tài)相差較大的地方，表明該處采用SCM不是最佳的選擇。

圖6 陜西地震臺記錄的三分量地震數(shù)據(jù)Fig.6 The three components seismogram recorded by Shaanxi seismic net

圖8是由ACM和SCM對三分量實際數(shù)據(jù)的濾波結果（紅線表示SCM；藍線表示 ACM），（a）、（b）、（c）分別對應于X、Y、Z分量，每一分量最上邊是該分量的原始數(shù)據(jù)，接下來依次是SCM和ACM的濾波結果。對于SCM，窗口長度分別采用2、5、10以及50個采用點的長度，對于ACM方法，N分別給不同的值（N＝2、5、10、50）。通過改變窗口長度以及N值大小來比較、分析兩種方法的實際濾波效果。

首先來分析SCM的濾波效果。我們注意到，當窗口長度很小時（N＝2），在有效信號的附近布滿了高頻干擾；隨著窗口長度的增加，高頻干擾逐漸降低；但窗口長度選擇過大（N＝50），又會將有效信號削弱。因此，對于SCM方法有一個非常大的局限性——如何選擇合適的窗口長度，對信號進行極化濾波，重構有效信號。

對于ACM極化濾波方法，從圖中可以看到該方法大大減小了有效信號在濾波前后的波形差異。同樣，對比原始信號和濾波后的信號，ACM在有效信號附近幾乎沒有高頻干擾存在。隨著N值逐漸增大，主要影響極化分析的敏感程度，但不會明顯影響濾波器中前兩個特征值的比值，這也解釋了N取不同值時濾波結果沒有明顯差異的原因。

總之，在SCM基礎上提出的ACM 極化濾波方法克服了SCM 窗口長度固定的不足，采用自適應于三分量信號瞬時頻率的窗函數(shù)大大提高了極化濾波效果。

5 結語

本文實現(xiàn)了一種基于協(xié)方差矩陣近似表達的極化分析新方法——ACM極化分析方法。并與SCM極化濾波方法進行對比分析，可以得出以下結論：

（1）協(xié)方差矩陣在三分量地震極化濾波中發(fā)揮著極其重要的作用，可以通過協(xié)方差矩陣的特征值和特征向量構造各種能夠描述不同類型地震波極化屬性的極化參數(shù)，還可以將沒有明顯極化特性的噪音從地震信號中剔除。

（2）標準協(xié)方差極化分析方法的特點是計算結果相對比較穩(wěn)定，對干擾不敏感。但時窗長度的選擇對于該方法的效果有非常大的影響，并且由于時窗長度的影響無法確定地震記錄的開始和結束部分的極化參數(shù)。該方法在目前的實際應用中有很大的局限性。

圖7 由ACM和SCM計算的極化參數(shù)比較（藍線表示SCM，紅線表示ACM）Fig.7 Comparison of the polarization attributes obtained from the ACM and SCM （Blue lines：SCM；Red lines：ACM）

圖8 由ACM和SCM對三分量實際數(shù)據(jù)的濾波結果比較Fig.8 Comparison of polarization filtering results of three components between ACM and SCM

（3）自適應協(xié)方差極化分析方法與標準協(xié)方差極化分析方法相比，優(yōu)勢在于時窗的長度自適應于三分量地震記錄的瞬時頻率，減小了在選擇時窗長度時的人為影響，濾波精度大大提高。

（4）自適應協(xié)方差極化分析方法是在三分量地震記錄的每一個時間采樣點上求極化特征參數(shù)的，因此不需要進行插值處理。

（5）由于SCM和ACM極化分析方法都是是在時間域實現(xiàn)的，沒有充分利用頻率信息，對于在時間上有重疊的地震信號，僅僅在時間域或頻率域做極化濾波處理，有時候并不能得到理想的分離及濾波效果。下一步研究將考慮把ACM極化分析方法推廣到時頻域中，利用時頻分辨率較高的小波變換或S變換［24－25］，在時頻域中求取信號的自適應瞬時極化參數(shù)，借助很高的時頻分析能力，來提高各向異性介質(zhì)中具有頻散特性的面波極化分析的效果。

（References）

［1］加爾彼林EИ，著.何樵登，楊寶俊，譯.地震勘探偏振法［M］.北京：石油工業(yè)出版社，1989.ГальперинЕИ，Edit.HE Qiao－deng，YANG Bao－jun，Trans－late.Polarization Method of Seismic Exploration［M］.Beijing：Petroleum Industry Press，1989.（in Chinese）

［2］Nguyen D T，Brown R J，Lawton D C.Polarization filter for Multi－conponent Seismic Data［J］.Exploration Geophysics，2001：93－101.

［3］Jurkevice A.Polarization Analysis of Three－component Array Data［J］.Bull Seis Soc Am，1988，78：1725－1743.

［4］Jackson G M，MasonI M.Principal Components Transforms of Triaxial Recordings by Singular Value Decomposition［J］.Geophyiscs，1991，56：528－533.

［5］Reading A M.Polarization Filtering for Automatic Picking of Seismic data and Improved Converted Phase Detection［J］.Geophysical Journal Lnternational，2001，147：227－234.

［6］Schimmel M，Gallart J.The use of Instantaneous Polarization Attributes for Seismic Signal Detection and Image Enhancement［J］.Geophysical Journal Lnternational，2003，155：653－668.

［7］李英康，崔作舟.分離縱波和橫波的偏振旋轉法［J］.地球物理學報，1994，37（增刊Ⅱ）：372－382.LI Ying－kang，CUI Zuo－zhou.P and S－waves Separated by Polarization Revolving Method［J］.Chinese J Geophys，1994，37（Supp.2）：372－382.（in Chinese）

［8］葛勇，韓立國，韓文明，等.極化分析研究及其在波場分離中的應用［J］.長春地質(zhì)學院學報，1996，26（1）：83－88.GE yong，HAN Li－guo，HAN Wen－ming，et al.Study and Application of Polarization Analysis in Separation of P and S Waves［J］.Journal of Changchun University of Earth Sciences，1996，26（1）：83－88.（in Chinese）

［9］黃中玉，高林，徐亦鳴，等.三分量數(shù)據(jù)的偏振分析及其應用［J］.石油物探，1996，35（2）：9－16.HUANG Zhong－yu，GAO Lin，Xu Yi－ming，et al.The Analysis of Polarization in Three－component Seismic Data and its Application［J］.Geophysical Prospecting for Petroleum，1996，35（2）：9－16.（in Chinese）

［10］劉建華，劉福田，胥頤.三分量地震資料的偏振分析［J］.地球物理學進展，2006，21（1）：6－10.LIU Jian－h(huán)ua，LIU Fu－tian，XU yi.The Analysis of Polarization in Three－component Seismic Data［J］.Progree in Geophysics，2006，21（1）：6－10.（in Chinese）

［11］Meissner R，Hegazy M A .The ratio of the PP－to the SS－reflection Coefficients a Possible Future Method to Estimate oiland Gas Reservoirs［J］.Geophysical Prospecting，1981，29：533－540.（in Chinese）

［12］Helbig K，Mesdag C S.The Potential of Shear Wave Observations［J］.Geophysical Prospecting，1982，30：413－431.

［13］Li X L，Crampin S.Complex Component Analysis of Shear Wave Splitting：Theory［J］.Goephysical Journal International，1991，107：597－604.

［14］黃忠賢，陳虹，王貴華，等.面波偏振與中國大陸巖石層橫向不均勻性［J］.地球物理學報，1994，37（4）：456－468.HUANG Zhong－xian，CHEN Hong，WANG Gui－h(huán)ua，et al.Surface Wave Polarization and Lateral Heterogeneities of the Lithosphere in China［J］.Chinese J Geophys，1994，37（4）：456－468.（in Chinese）

［15］張中杰.多分量地震資料的各向異性處理與解釋方法［M］.哈爾濱：黑龍江教育出版社，2002.ZHANG Zhong－jie.Methods of Anisotropic Processing and Lnterpretation for Multi－Component Seismic Data［M］.Harbin：Heilongjiang Education Press，2002.（in Chinese）

［16］Flinn E A.Signal Analysis Using Rectilinearity and Direction of Particle Motion［J］.Proceedings of the IEEE，1965，53（12）：1874－1876.

［17］Montalbetti J F，Kanasewich E R.Enhancement of Teleseismic Body Phase with a Polarization Filter［J］.Geophy J R Astr Soc，1970，12：19－129.

［18］Samson J C.Matrix and Stokes Vector Representations of Detectors for Polarized Waveforms：Theory，with Some Applications to Teleseismic Waves［J］.Geophys，1977，51：583－603.

［19］Kanasewich E R.Time Sequence Analysis in Geophysics［D］.Edmonton：University of Alberta Press，1981.

［20］Viddale J E.Complex Polarization Analysis of Particle Motion［J］.Bull Seis Soc Am，1986，76：1393－1405.

［21］Morozov I B，Smithson S B.Lnstantaneous Polarization Attributes and Direction Filtering［J］.Geophysics，1996，15：872－881.

［22］Zhang J，Walter W R，et al.Time－domain Pure－state Polarization Analysis of Surface Waves Traversing California［J］.Pure Appl Geophys，2003，160（8）：1447－1478.

［23］Diallo M S，Kulesh M，Holschneider M，et al.Lnstantaneous Polarization Attributes Based on Adaptive Covariance Method［J］.Geophysics，2006，71（5）：99－109.

［24］姚家駿，楊立明，馮建剛.常用時頻分析方法在數(shù)字地震波特征量分析中的應用［J］.西北地震學報，2011，33（2）：105－110.YAO Jia－jun，YANG Li－ming，F(xiàn)ENG Jian－gang.Application of Common Time－frequency Analysis Methods in Analyzing Characteristic Quantity of Digital Seismic Wave［J］.Northwestern Seismological Journal，2011，33（2）：105－110.（in Chinese）

［25］許康生，李英，李秋紅.近地震波的小波相對能量分布特征分析［J］.地震工程學報，2013，35（1）：166－170.XU Kang－sheng，LI Ying，LI Qiu－h(huán)ong.Distribution Characteristics of Wavelet Relative Energy on Near－earthquake Wave［J］.China Earthquake Engineering Journal，2013，35（1）：166－170.（in Chinese）