槽波地震技术对于探测小尺度的地质构造比地面地震勘探更具优势，因此它作为矿井物探的一种重要手段而被广泛应用。槽波地震勘探是利用煤层作为低速波导，用以探查煤层不连续性的一种地球物理方法。它具有能量强、精度高、传播距离远、波形特征易于识别等特点[1-5]，这些特点对于探测小断层、陷落柱、采空区及煤层厚度变化带等，具有地面地震勘探无可比拟的优势。地面地震勘探的反射界面通常平行或近似平行于地面，而槽波勘探对象与测线之间夹角可以是0°～90°任意一种情况，加上井下环境较为恶劣，有众多的干扰信号，且能布置探测的空间与通道极其有限。因此，具有频散特性的槽波信号经常湮没于干扰信号中，而槽波的频散特性又使其难以应用常规去噪手段，这就使得去噪成为槽波数据处理的一个重要环节。

槽波资料频带较宽，主频较高，对其影响严重的随机噪声在频谱上与槽波信号相互交叠，同时随机噪声不具有特定的方向信息。因此，无论是频率域还是空间域都难以有效地将两者分离。然而，有效信号和干扰波在偏振方向和偏振度上存在不同，基于此，MONTALBETTI和KANASEWICH[6]将极化分析的方法引入地震数据处理中，通过时间域的极化滤波去除远震资料中的噪声。这种方法假定时间域的某个信号在整个频段上具有相同的极化特征，但这与实际情况存在不符;SAMSON和OLSON[7]设计了频率域的自适应极化滤波，它对噪声较弱的情况有着明显效果，却难以对强噪声背景下的信号进行有效提取。为了改进强噪声背景下的极化滤波效果，SAMSON在频率域使用谱矩阵对信号进行了极化特征分析;NGUYEN[9]和ROBERTO[10]等采用极化分析方法用于多分量地震数据的滤波处理;PARK在频域采用多窗谱分析的方法，结合谱密度矩阵得到一个自适应的极化滤波器[11-12]。由此，在前人极化分析的基础上，综合考虑槽波信号与背景噪声的特点，对频率域的极化滤波函数进行改进。同时采用时间域的滑动时窗，由谱密度矩阵计算时窗内不同频率信号的偏振系数，提出了一种槽波资料频率域的极化滤波方法。

1 频率域极化滤波原理

偏振分析能够定量描述质点的振动轨迹，极化滤波是根据质点的不同偏振特性进行分离的一种滤波方法。有效波和干扰波具有不同的偏振特性，有效波通常具有特定的偏振方向，如Rayleigh波为椭圆偏振，Love波为线性偏振，且偏振度较高;而随机噪声没有固定的偏振方向，并且偏振度较低(图1)，据此可以设计偏振滤波器保留或去除具有某种偏振特性的地震波。

根据SAMSON和OLSON，频率域的极化滤波可表达为

(1)

式中，g为极化滤波因子，是用来调节随机噪声压制程度的正整数，取值为1≤g≤6，值越大压制程度越高;z(ω)为x(t)傅里叶变换后的频域数据;P(ω)为偏振系数，0≤P≤1，它反映了极化程度，P为0时表明质点偏振状态是随机状态，P为1时表明高度极化。通过频率域的极化滤波，能够对具有较小偏振系数的随机噪声进行压制，而保留具有较大偏振系数的有效信号，并通过极化滤波因子取不同的值，调整压制程度。P的表达式为

(2)

其中，S为谱密度矩阵;Tr代表它的迹。对于一个三分量数据来说;S为一个3×3的矩阵。S(ω)可由式(3)构造

(3)

其中

(4)

(5)

其中，为x(t)多窗谱分析得到的结果;为k阶窗函数，k=0,1,…,K-1,W为分辨率带宽，N为窗口内数据的数目。NW作为一个参数控制多窗谱分析的分辨率，它决定了多窗谱分析中窗的个数。

2 多窗谱分析

计算偏振系数首先需要求取谱密度矩阵，常规傅里叶窗函数求取谱密度矩阵时会产生奇异值，增加了谱密度矩阵估计的畸变。因此，在求取谱密度矩阵时引入多窗谱分析来减少频谱泄漏，克服常规傅里叶变换中高频成分因频谱泄漏导致的偏倚，并消除频率成分的剧烈变化。多窗谱分析是一种低方差、高分辨的谱分析方法，适用于短序列、高噪声背景下的准周期信号分析，其优势在信噪比较低的序列中尤为突出[13]。多窗谱法的原理是构造一组正交序列窗口，使频率泄漏最小，称为Slepian序列[14-16]，它在傅里叶变换时作用于数据，对其加权平均，获得一个低方差的平滑谱估计。Slepian窗是从最优化的观点出发，寻求一组能够使感兴趣的波段能量最大化的窗。Slepian窗的构建由下式Toeplitz矩阵的特征向量得到

(6)

其中，矩阵D的特征值为0<λN-1< …<λ1<λ0<1，对应的特征向量即为Slepian序列。Toeplitz矩阵的前2NW-1个特征值接近于1，这些特征值对应的向量是一组能够有效减少频率泄漏的窗。因此，多窗谱法中通常选取前2NW-1阶的Slepian序列作为窗函数。

3 极化滤波函数的改进

由谱密度矩阵结合多窗谱分析得到的偏振系数，在有效信号频率处具有接近于1的值，噪声频率处接近于0。式(1)所表示的极化滤波函数形式决定了偏振系数从0到1为逐渐过渡的，过渡过程难免引起有效信号的泄露及噪声压制不足，而希望的是有效信号具有尽量大的系数，同时噪声频段得到充分压制。为此，在式(1)中引入一个由正弦函数代替的权函数，有

(7)

式(7)由P(ω)计算的正弦函数作为新的权函数代替偏振系数，这样在偏振系数较小的噪声频率处(图2实线接近0的位置处)对其有了更高的压制，有效信号偏振系数(图2实线接近1的位置处)则被提升，同时减小了噪声与有效信号之间的过渡，使2者之间的压制与提升作用更明显。

另一方面，在有噪声的情况下，有效信号频段难免混入噪声，从而影响到偏振系数。为此，还需消除噪声对有效信号偏振系数的影响。谱密度矩阵S为一个非负的Hermitian矩阵，它可转换为

(8)

其中，λj为谱密度矩阵S的特征值，uj为对应的特征向量。当仅有一个非零特征值(λ1≠0，λ2=λ3=0)时

(9)

式(9)代表了完全极化的状态，此时由式(2)计算得到的偏振系数P=1，但在含有强背景噪声的情况下，完全极化状态并不存在，有效信号的偏振系数难以达到1。为了最大限度地发挥权函数(式(7))的作用，需要去除有效信号中噪声的影响。

对于时间信号x(t)的傅里叶变换z，在含有噪声时可表示为有效信号与噪声的结合

z=av+e

(10)

其中，a为与频率相关的复系数;v为有效信号;e为噪声频谱。此时，谱密度矩阵可表示为

S=E{a2}vv*+E{ee*}

(11)

式(11)右端第1项即代表式(9)所表示的完全极化状态，第2项则表示噪声谱密度矩阵，由N(ω)表示

N(ω)=E{e(ω)e*(ω)}

(12)

对于有效信号中N(ω)的估计较为困难，但N(ω)的理想模型是可以求取的。对噪声白化的过程可用于清除噪声对有效信号的影响[12]，通过用N-1/2(ω)z(ω)代替z(ω)来重新构建S(ω)，得到

A(ω)=N-1/2(ω)S(ω)N-1/2(ω)

(13)

通过式(11)可知，A满足:仅有噪声时P(A)=0(S=ee*)，同时P(A)=1(S={a2}vv*)，这与P(S)的理论极值一致。因此，通过求取噪声的谱密度矩阵N(ω)，并用N-1/2(ω)z(ω)代替z(ω)即可消除有效信号中噪声的影响，而噪声谱密度矩阵又可由初至前的背景噪声求取。最终，对消除了噪声影响的有效信号应用式(7)表示的极化滤波权函数，得到提高了信噪比的滤波输出。

4 滤波效果分析

4.1 模型试算

设计煤层厚度为5 m的均匀介质3层模型，时间采样间隔dt=0.25 ms，纵波震源雷克子波主频120 Hz。观测系统为超前探测的方式，接收线沿x轴及y轴展布，炮线分布于x轴，炮间距分别为10 m，道间距为2 m(图3)。通过交错网格高阶有限差分技术，采用弹性波波动方程数值模拟进行正演计算，得到三分量的槽波地震记录(图4(a))。合成记录上，瑞雷型槽波主要在x和z分量上，勒夫型槽波分布于x和y分量。

加入随机噪声的记录如图4(b)所示，在噪声的影响下，反射槽波能量变得不明显，尤其远道处，反射槽波甚至淹没于随机噪声中无法识别。采用前述频率域的极化滤波方法对其进行去噪处理，并与时间域的极化滤波进行对比，结果如图4(c),(d)所示。可知，频率域极化滤波后的记录整体信噪比得以明显提升，直达槽波和反射槽波能量被最大限度的保留。相比时间域的极化滤波，频率域的极化滤波在槽波弱信号的保留方面有较好的效果，图4(c),(d)中x及z分量反射槽波淹没于噪声、能量较弱的部分，在频率域极化滤波后得以恢复，而时间域的极化滤波仅能提高能量较强的槽波信号信噪比。同时，去噪前后的差值记录(图4(e),(f))也显示出，时间域的极化滤波去除的噪声包含较多的有效信号，而频率域的极化滤波几乎不含有效信号，保持了较高的保真度。

这是由于随机噪声的偏振系数较小，在频率域极化滤波过程中，根据初至前的背景噪声求取噪声谱密度矩阵，重新构建记录的谱密度矩阵，在一定程度上消除了噪声对有效信号的影响;再经过改进后的极化滤波函数进行频率域的滤波，偏振方向为随机且偏振系数较小的噪声被削弱，而在固定方向上偏振的槽波具有较大的偏振系数，由此在极化滤波过程中被加强。因此，远道处的弱反射槽波能量在频率域极化滤波后能够被有效识别。

4.2 实际资料处理

采用三分量检波器在某矿得到的槽波资料如图5(a)所示，记录中的背景噪声较强，槽波能量仅在近道处可以识别，而远道处以及整体的槽波频散特征都被强背景噪声覆盖。采用频率域的极化滤波对资料进行处理，根据数据特征，极化滤波的参数设置为ρ(t，g，N)∈[50，2，4]，其中的t为滤波时窗长度(ms)。图5(b)为极化滤波后的槽波记录，与滤波前的记录相比，背景噪声得到有效压制，槽波能量分布情况也更加容易分辨，记录的信噪比明显提升。

5 结 论

(1)相比时间域的极化滤波，频率域求取的偏振系数具有时变特性，得到的滤波结果也会更加自然。利用这样的优势，通过改进极化滤波权函数，提升了对噪声信号的压制效果。

(2)多窗谱分析采用的Slepian窗，应用于高噪声背景下的信号分析，使得频率泄漏较小，能够对有效信号有一个更准确的谱估计结果。

(3)在该方法中极化滤波因子g的取值决定了对噪声的压制程度，g越大，对噪声的压制越明显;g越小，对输入数据的改变作用越弱。因此，在实际应用中要根据输入数据选择合适的g值，提高噪声压制效果。

参考文献(References):

[1] 胡国泽，滕吉文，皮娇龙，等.井下槽波地震勘探-预防煤矿灾害的一种地球物理方法[J].地球物理学进展,2013,28(1):439-451.

HU Guoze,TENG Jiwen,PI Jiaolong,et al.In-seam seismic exploration techniques-A geophysical method predicting coal-mine disaster[J].Progress in Geophys,2013,28(1):439-451.

[2] 陆斌，程建远，胡继武，等.采煤机震源有效信号提取及初步应用[J].煤炭学报,2013,38(12):2202-2207.

LU Bin,CHENG Jianyuan,HU Jiwu,et al.Shearer source signal extraction and preliminary application[J].Journal of China Coal Society,2013,38(12):2202-2207.

[3] 姬广忠，程建远，朱培民，等.煤矿井下槽波三维数值模拟及频散分析[J].地球物理学报,2012,55(2):645-654.

JI Guangzhong,CHENG Jianyuan,ZHU Peimin,et al.3-D numerical simulation and dispersion analysis of in-seam wave in underground coal mine[J].Chinese Journal Geophysics,2012,55(2):645-654.

[4] 姬广忠，程建远，朱培民.煤层Love型槽波数值模拟及其频散特征分析[J].煤炭科学技术,2011,39(6):106-109.

JI Guangzhong,CHENG Jianyuan,ZHU Peimin.Numerical simulation of seam love type channel-wave and analysis on dispersion features[J].Coal Science and Technology,2011,39(6):106-109.

[5] 杨思通，程久龙.煤巷小构造Rayleigh型槽波超前探测数值模拟[J].地球物理学报,2012,55(2):655-622.

YANG Sitong,CHENG Jiulong.The method of small structure prediction ahead with Rayleigh in-seam wave in coal roadway and seismic wave field numerical simulation[J].Chinese Journal Geophysics,2012,55(2):655-622.

[6] MONTALBETTI J F,KANASEWICH E R.Enhancement of teleseismic body phases with a polarization filter[J].Geophys,1970(21):119-129.

[7] SAMSON J C,OLSON J V.Data-adaptive polarization filters for multichannel geophysical data[J].Geophysics,1981(46):1432-1431.

[8] SAMSON J C.The reduction of sample-bias in polarization estimators for multichannel geophysical data with anisotropic noise[J].Geophys,1983(75):289-308.

[9] NGUYEN D T,BROWN R D.Polarization filter for multi-component seismic data[J].Exploration Geophysics,2001:93-101.

[10] ROVERTO D F,GEMMA M.Polarization filter with singular value decomposition[J].Geophysics,2001,66(3):932-938.

[11] PARK J,VERNON F L,LINDBERG C R.Frequency dependent polarization analysis of high-frequency seismograms[J].Geophys,1987(92):12664-12674.

[12] DU Z,FOULGER G R,MAO W.Noise reduction for broad-band,three-component seismograms using data-adaptive polarization filters[J].Geophysical Journal International,2000,141(3):820-828.

[13] 王云专，王珊，董相杰，等.多窗谱分析在Q值估算中的应用[J].地球物理学进展,2009,24(6):2156-2162.

WANG Yunzhuang,WANG Shan,DONG Xiangjie,et al.The application of the multiple tapers spectral analysis in the Q estimation[J].Progress in Geophys,2009,24(6):2156-2162.

[14] THOMSON D J.Spectrum estimation and harmonic analysis[J].Proc.IEEE,2005,70(9):1055-1096.

[15] RIEDEL K S,SIDORENKO A.Minimum bias multiple taper spectral estimation[J].IEEE Trans.Signal Processing,1995,43(1):188-195.

[16] HU Y,LOIZOU P C.Speech enhancement based on wavelet thresholding the multitaper spectrum[J].IEEE Trans.Speech and Audio Processing,2004,12(1):59-67.