提高台站地震资料信噪比的自适应极化滤波①

马见青，李庆春

（长安大学地质工程与测绘学院，陕西 西安 710054）

0 引言

地震波场在本质上是不同类型、不同偏振特性的振动相互干涉和叠加的结果［1］。极化滤波是在基于波的偏振特性基础上进行的一种空间滤波的信号处理方法，目前极化分析方法在电磁学、光学和地震勘探等领域已经取得了一些研究成果［2］。一方面通过对地震多分量记录的分析能够获得地震波的极化参数，可以实现各类噪声压制、特定波型的识别与分离，提高地震资料的信噪比，从而为后续的各种数据处理与解释服务［3－10］；此外还可以应用极化滤波分析横波分裂，进行各向异性研究［11－14］。对于充分利用地震资料中所包含的丰富的动力学和运动学信息，深化对地球本体或地下介质的认识具有重要的理论意义和实际应用价值［15］。

极化分析方法在地震信号处理中的应用研究始于上世纪60年代，Flinn［16］最早提出了极化分析方法的数学算法。以后 Montalbetti、Samson、Kanssewich、Vidale、Morozov、Zhang等［17－22］相继在时间域内对基于协方差矩阵的极化分析技术进行了大量研究。具体作法是对三分量地震数据sk（t）（k＝x，y，z）在一定长度的时窗内构造协方差矩阵M（ξ）；其特征值分别为λi（i＝1，2，3），且λ1≥λ2≥λ3；Vi为λi为 对应的特征向量（图1）；通过时窗内协方差矩阵的特征值和特征向量构造出一系列能够表征不同类型地震波极化特性的极化参数。这种方法原理简单，实现方便，但对时窗长度的选择依赖性较强，而且在给定长度的时窗内求得的极化参数不具有时变特性，由此设计的滤波器的滤波效果也不够理想。

图1 笛卡尔坐标系下的三分量极化椭球示意图Fig.1 schematic representation of an multi－component polarization ellipsoid in Cartesian coordinate system

鉴于上述方法的局限性，Diallo等［23］提出了基于自适应协方差矩阵的瞬时极化分析方法。该方法的时窗长度可以自适应于信号的瞬时频率，即协方差矩阵的计算时窗长度自适应调整到时窗内信号的最小周期，可以得到每一分量上各个采样点的瞬时极化参数，不需要进行插值（开始和结束处各1／2时窗除外）。这样我们就可以摆脱标准协方差矩阵极化分析中必须人为选择时窗长度的限制。最后，通过数值模拟实验及实际三分量台站资料处理进一步比较和解释了基于协方差矩阵的极化分析方法。

1 自适应协方差极化分析方法

对三分量地震数据sk（t），（k＝x，y，z）求其Hilbert变换谱（t），得到其解析信号

在时间t附近的局部信号可以利用解析信号ck（t）来近似表达：

利用多分量信号sk（t）的近似表达式（2）来构造协方差矩阵

协方差矩阵 M（t）与 M（ξ）不同，M（ξ）代表某一时窗内的平均极化分布，而M（t）代表的是地震信号中每一个时间点的瞬时极化分布，这正是自适应协方差矩阵方法（ACM）相比标准协方差矩阵方法（SCM）的优势所在。

该协方差矩阵M（t）的特征值分别为λ（i＝1，2，3），且λ1≥λ2≥λ3。Vi为λi对应的特征向量。因此可以得到极化椭球体的瞬时极化参数：

瞬时极化轴

可以证明，最小的特征向量V3（t）平行于这个平面向量 ，Vp（t），V3（t）的方向余弦θx（t），θy（t），θz（t）可以确定主椭圆面在三维空间的位置：

上述得到的极化参数都是时窗内的中心点位置的极化参数，通过将时窗沿时间轴逐个采样点向前移动，可以得到整个时间轴（除开始和结束处各1／2时窗）。根据三分量地震数据在空间某一时窗内的极化参数，我们就可以构造不同的极化滤波函数，保留具有某一组极化参数的波，从而达到波场分离或去噪的目的。

滤波器为［19］

2 自适应分析时窗的确定

Tkm表示t时刻M（t）的自适应时窗长度，定义为

或者可以由自适应协方差矩阵M（t）中每一个元素的瞬时频率来定义：

从方程（2）中我们可以看到，每一分量的信号都是由周期为2π／Ωk（t）的余弦函数近似得到。为了阐明自适应时窗长度Tkm（t）的原理，我们考虑Ωx（t）＝Ωy（t）＝Ωz（t）＝Ω（t）的情况，在这种情况下Tkm可以简化为Tkm（t）＝2πN／Ω（t），它是周期的整数倍，N为一个正整数。

通常情况下不同分量的瞬时频率是不同的，此时方程（8）和（9）所定义的Tkm（t）为一个平均周期。相比之下方程（9）更能提供最佳的自适应时间窗。

为了分析参数N的选取对本方法的影响，现设计一模型来阐述。图2（a）为一个三分量数据sx（t）、sy（t）和sz（t）（红色实线）以及利用方程（2）计算得到的相应近似表达式（黑色虚线）。在0～0.05s时间段内N＝1，在这个时间段内可以看到近似得到的信号和原始模拟信号的波形相当一致。同样在图2（b）内的矢端曲线图上也可以看到矢端曲线近似为一个完整的圆形。此外，这个矢端曲线图表明，当N＝1时极化分布不能刻画三维空间内比椭圆更复杂的结构。从图2（a）可以看出，随着N的增大信号的近似表示的精度在下降，但同时比椭圆更为复杂的曲线在矢端曲线图上出现（图2（b）～2（d））。因此在刻画平面内的极化属性时选取N＝1或2就足够了；当想刻画三维空间里的结构更为复杂的极化属性时，我们可以牺牲近似表达的精度，利用较大的N值。

图2 N值对三分量信号近似表示的影响Fig.2 Influence of Non the accuracy of the approximation for signals

3 数值模拟

为了进一步比较时间域内的标准协方差极化分析（SCM）和自适应瞬时极化分析（ACM），现人工合成一个三分量记录（图3），该记录由5个波段构成，A段表示平稳椭圆，B段表示旋转椭圆，C段表示线性极化的情况，D段表示三维空间的平稳椭圆，E段对应于旋转椭球的情况。前三段只在X－Y平面内分布，而后两段在整个三维空间均存在。图4为该三分量记录的空间极化轨迹图。

图3 三分量人工合成记录［23］Fig.3 Three components synthetic record［23］

图4 三分量人工合成记录的空间极化轨迹Fig.4 Hodograph of the three components synthetic record

图5 时间域瞬时极化参数Fig.5 Instantaneous polarization attributes in time domain

图5（a）为用标准协方差方法和自适应协方差方法得到的瞬时极化轴，我们可以看到，标准协方差方法与自适应协方差方法相比精度很低。这主要是由于标准协方差极化分析法的时窗长度很难准确地选择，也就是说时窗长度的选取对于标准协方差极化分析法的精度有很大的影响。图5（b）是极化率曲线，即Rmed（t）／Rmax）（t），从中可以很容易把线性极化波（rs（t）／R（t）≈0）从椭圆极化波和椭球极化波中区分出来。但标准协方差方法与自适应协方差方法相比，求得的极化率震荡很厉害。图5（c）是由方程（9）计算得到的自适应时间窗，可以看到时窗长度随频率的增大而变小，反之亦然。

值得注意的是，在信号周期和时窗长度接近的情况下，两种方法得到的极化轴曲线吻合程度很高，曲线形态差别较大的区域说明标准协方差方法的时窗长度选择不是最理想的。

4 实际三分量台站资料处理

图6是陕西省地震台网2002年记录的三分量实际地震数据，该台站位于纬度34.25°，经度108.95°。图7是分别通过SCM和ACM极化滤波方法计算得到的极化主轴、极化次轴（蓝线表示SCM，红线表示ACM），通过对比两种方法得到的极化轴曲线。可以看到，ACM对局部变化比较剧烈的信号更加敏感，而SCM实际上是ACM的平滑过程。而且，由于SCM的窗口长度固定，不适合用来刻画周期长度小于窗口长度的地震信号。在信号周期和窗口长度接近的地方，两条曲线的相似程度最高，在曲线形态相差较大的地方，表明该处采用SCM不是最佳的选择。

图6 陕西地震台记录的三分量地震数据Fig.6 The three components seismogram recorded by Shaanxi seismic net

图8是由ACM和SCM对三分量实际数据的滤波结果（红线表示SCM；蓝线表示 ACM），（a）、（b）、（c）分别对应于X、Y、Z分量，每一分量最上边是该分量的原始数据，接下来依次是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 窗口长度固定的不足，采用自适应于三分量信号瞬时频率的窗函数大大提高了极化滤波效果。

5 结语

本文实现了一种基于协方差矩阵近似表达的极化分析新方法——ACM极化分析方法。并与SCM极化滤波方法进行对比分析，可以得出以下结论：

（1）协方差矩阵在三分量地震极化滤波中发挥着极其重要的作用，可以通过协方差矩阵的特征值和特征向量构造各种能够描述不同类型地震波极化属性的极化参数，还可以将没有明显极化特性的噪音从地震信号中剔除。

（2）标准协方差极化分析方法的特点是计算结果相对比较稳定，对干扰不敏感。但时窗长度的选择对于该方法的效果有非常大的影响，并且由于时窗长度的影响无法确定地震记录的开始和结束部分的极化参数。该方法在目前的实际应用中有很大的局限性。

图7 由ACM和SCM计算的极化参数比较（蓝线表示SCM，红线表示ACM）Fig.7 Comparison of the polarization attributes obtained from the ACM and SCM （Blue lines：SCM；Red lines：ACM）

图8 由ACM和SCM对三分量实际数据的滤波结果比较Fig.8 Comparison of polarization filtering results of three components between ACM and SCM

（3）自适应协方差极化分析方法与标准协方差极化分析方法相比，优势在于时窗的长度自适应于三分量地震记录的瞬时频率，减小了在选择时窗长度时的人为影响，滤波精度大大提高。

（4）自适应协方差极化分析方法是在三分量地震记录的每一个时间采样点上求极化特征参数的，因此不需要进行插值处理。

（5）由于SCM和ACM极化分析方法都是是在时间域实现的，没有充分利用频率信息，对于在时间上有重叠的地震信号，仅仅在时间域或频率域做极化滤波处理，有时候并不能得到理想的分离及滤波效果。下一步研究将考虑把ACM极化分析方法推广到时频域中，利用时频分辨率较高的小波变换或S变换［24－25］，在时频域中求取信号的自适应瞬时极化参数，借助很高的时频分析能力，来提高各向异性介质中具有频散特性的面波极化分析的效果。

（References）

［1］加尔彼林EИ，著.何樵登，杨宝俊，译.地震勘探偏振法［M］.北京：石油工业出版社，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］.长春地质学院学报，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］黄中玉，高林，徐亦鸣，等.三分量数据的偏振分析及其应用［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－hua，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－hua，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］姚家骏，杨立明，冯建刚.常用时频分析方法在数字地震波特征量分析中的应用［J］.西北地震学报，2011，33（2）：105－110.YAO Jia－jun，YANG Li－ming，FENG 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－hong.Distribution Characteristics of Wavelet Relative Energy on Near－earthquake Wave［J］.China Earthquake Engineering Journal，2013，35（1）：166－170.（in Chinese）