1 引言

现今，随着地震勘探的深入，一些沙漠、丘陵、盆地、湖泊等工区勘探条件越显复杂.地震勘探过程中得到的地震记录混有各种随机干扰^{[1, 2]}.若随机干扰强度到达一定程度，必将影响后期地震资料的动校正速度的分析、静校正量的拾取等，最终影响数据成像效果，从而导致地震勘探工作者对地下地质结构及地质体的解释变得困难.

到目前为止，人们设计了多种不同的随机噪声的消除方法，例如：f-x域反褶积滤波方法、多项式拟合、径向预测滤波、反Q滤波^[3~5]等等.这些方法在不同程度上解决了一些随机噪声的压制问题，但是仍然存在难以滤除强随机噪声的问题.最为经典的维纳滤波方法，由于其适用于任何平稳信号，并且能在一定程度上滤除含噪信号中的随机噪声而得到较为广泛的应用^{[6, 7]}.将维纳滤波应用在地震勘探中，在一定条件下可以滤除随机噪声，但当采集到的地震信号信噪比低到一定程度，维纳滤波处理过程中矩阵方程求解将会出现矩阵病态问题，使得求解滤波结果不精确，以至信号幅值产生衰减.

维纳滤波求解过程中出现的矩阵病态问题可以通过数据降维解决，主分量分析^{[8, 9]}（Principal Component Analysis，PCA）与Mercer核函数^{[10, 11]}结合是一种较为理想的方法.主分量分析是一种线性降维方法，它把高维数据投影到低维子空间，是一种从高维数据空间提取数据特征的极好工具，目前应用于图像识别^[12]、化学^[13]、雷达目标识别^[14]等各个领域.为了使主分量分析能提取出非线性的数据特征，在PCA方法的基础上加入Mercer核函数，称为核主分量分析（KPCA）^[15~18].它的主要原理是将输入样本数据空间通过非线性映射到新的特征空间，使得数据在高维特征空间中线性可分，利用核函数代替内积计算，从而极大降低了计算的复杂度.

综合以上，本文提出一种基于核函数主分量的维纳滤波压制地震资料中强随机噪声方法（K-WPC），利用线性核函数对要处理的数据进行核变换，使其映射到高维数据空间，再利用主分量分析进行降维.核函数与主分量分析的优势在于：核函数在一定程度上减弱了滤除随机噪声时带来的幅值衰减问题，同时使得低信噪比条件下维纳滤波求解过程中出现的矩阵病态问题得到改善.此滤波方法对正演地震信号和实际地震信号做消噪处理，都得到令人满意的效果.

2 K-WPC方法的基本原理 2.1 核函数主分量维纳滤波

经典维纳滤波中的维纳霍夫方程为^[8]

(1)

其中Rxx、Rxd分别表示数据序列x的自相关函数E{xx^T}、数据序列x和理想输出数据序列d的互相关函数E{dx}，其中E表示数学期望.求得滤波因子序列w_n为

(2)

为解决求解维纳滤波因子时R_xx^-1易出现病态的问题，本文采取核主分量维纳滤波方法.K-WPC方法首先通过一个非线性核函数ϕ，把观测数据映射到更高维的特征空间中，然后在高维数据空间里采用标准的PCA方法降维并求解维纳滤波^[18].这里PCA采用点积形式表示，利用Mercer核函数，新特征空间的点积可以直接在输入数据空间中计算，而不需要明确计算出非线性核函数ϕ.

一个数据样本集x_i（i=1，2，…，N，x_i∈Rⁿ，N为样本数）在非线性映射ϕ的作用下映射到高维空间，即得到ϕ(x_i).为了确保ϕ(x_i)的均值为0，即保证采取各项减去其均值得到新的ϕ(x_i).为了避免在高维特征空间求解维纳滤波计算量大和病态问题，对ϕ(x_i)进行PCA降维处理. ϕ(x_i)的协方差矩阵C为

(3)

C的特征值λ及其特征向量ν满足λv=Cv，其中λ₁≥λ₂≥…≥λn.

由于ν为ϕ（x₁），ϕ（x₂），…，ϕ（x_N）的线性组合，所以一定存在一组相应的系数矩阵{α_l}，使得

(4)

将式（3）、（4）代入到λν=Cν中，有

(5)

可得

(6)

现定义两矩阵：

（1）n×n矩阵K称为核矩阵，它的第ij位元素称为内积核，记作K(x_i，xj)，i，j=1，2，…，N.根据再生核空间理论，向量ϕ(x_i)，ϕ(xj)在特征空间中的内积＜ϕ(x_i)，ϕ(xj)＞可以用内积核K(x_i，xj)表示，即K(x_i，x_j)=＜ϕ(x_i)，ϕ(x_j)＞，而不用精确地计算映射到特征空间的函数ϕ（·），其中〈·〉表示内积.

（2）n×1矩阵α，第l个元素记作α_l.对式（6）进行化简整理得：

其中K²表示矩阵自身相乘，两端均有K，所以上式等同于：

(7)

由式（7）可知，α为核矩阵K的特征向量，与C的特征值λ对应.假设特征值是降序排列，即λ₁≥λ₂≥…≥λ_N，λ_Z为最小的非零值.为了确保协方差矩阵C的特征向量ν是归一化的，即ν_m^Tν_m=1，m=1，2，…，Z.根据式（4）和式（7）可以得到ν_m^Tν_m=1的等价归一化条件，α_m^Tα_m=1/λ_m，m=1，2，…，Z.计算ϕ(x_i)向协方差矩阵C的Z个非零特征向量张成的特征空间V=[ν₁，ν₂，…，ν_Z]上的投影P_α，i，

(8)

其中P_α，i为ϕ(x_i)的主分量，向量α_m，l为核矩阵K中第m个特征值对应的特征向量αm的第l个系数，m=1，2，…，p，核矢量K_i=[K（x_i，x₁），K（x_i，x₂），…，K（x_i，x_N）]T，矩阵A=[α₁，α₂，…，α_Z].

由上述过程可知，KPCA方法只需要在原空间通过核函数计算核矩阵，而不需要明确地计算出映射ϕ.通过定义核矩阵，主分量分析也可以当做是求解该核矩阵特征值和特征向量的问题.同时，在选取Z个特征值的时候，选取较大的特征值.因为原特征值越大，对重构时的贡献也越大，所以只需要选取原特征值中前大部分占绝对优势的特征值，使得维数Z小于原先矩阵维数，从而达到了降维的目的，以使矩阵的病态问题得到缓减.

利用线性核函数将地震数据映射到高维特征空间后，使用PCA方法得到所需要的主分量，利用所得到的主分量P_α，i，代替原始输入的地震信号x_i.利用提取出的主分量P_α=[P_α，1，P_α，2，…，P_α，Z]和理想参考信号d=[d1，d2，…，dN]^T，得到

K-WPC算法的滤波因子求解表达式为

(9)

综上所述，得到K-WPC算法的滤波表达式：

(10)

2.2 消减地震勘探随机噪声的K-WPC

采用K-WPC方法来解决地震勘探资料中的随机噪声消减问题.地震信号x（n）包含有效信号s（n）和随机噪声ω（n），即

(11)

为了消减地震信号中的随机噪声，针对低信噪比地震数据，基于噪声抵消的思想实现K-WPC对地震信号逐道进行处理.其实现过程如下：

首先，对第K道地震记录x_K（n），构建特征信号向量：

(12)

M是特征信号向量长度，相应地有x_kn=s_kn+ω_kn.

第二，选择适合于地震信号的齐次核函数K(x，y)=(〈x，y〉)^q，q∈IN^*，以噪声作为维纳滤波的待处理信号，计算得到核矩阵K，每一个元素为K(ω_i，ω_K)，i，j=1，2，…，N，并根据公式(8)计算第K道记录的主分量P_{α, k}.

第三，以观测数据x_kn作为参考数据d，采用公式（9）求K-WPC算法的滤波因子w_kn=（P_α^TP_α）^-1P_α^Tx_kn，利用主分量代替原待处理信号得到该道记录的维纳滤波y_kn=w_kn^TPα, k.

第四，从观测数据减去维纳滤波器的输出得到残差数据e_n应为待检测有效信号的估计ŝ_kn=x_kn-y_kn.

采用K-WPC消减地震勘探随机噪声，由于K-WPC经过核变换，在高维特征空间求解维纳滤波，提高了维纳滤波对输入空间非线性系统求解的能力，有利于更好地恢复地震信号；同时，通过结合核主分量分析降低数据维数，有效地解决了维纳滤波求逆过程中的病态问题.利用噪声抵消的思想实现K-WPC，有助于更好地防止信号在时域和频域的畸变.

3 实验仿真 3.1 正演信号噪声处理实验

本文在正演信号仿真实验时采用同源噪声为背景噪声，进行信号的滤波处理，同源噪声如图 1所示.

利用人工合成的共炮点记录进行试验.地震波是Ricker子波，偏移距为零，不考虑距离衰减.其中t₀为双程垂直走时，取值为300ms，主频为30Hz，D为道间距，取值为50m，K为道序号，取值为40道，v为该反射同相轴相应的叠加波速，取值为1900m/s，所满足的双曲时距方程为：

(13)

将核主成分维纳滤波算法（K-WPC）应用于上述模拟地震记录.采样间隔T_s=1/500s，每一道数据均取801个采样点，雷克子波的主振幅度为2.5，滤波器阶数取12，线性核函数矩阵维数取12.

单道试验数据处理结果如图 2所示.滤波前SNR=-16.4622dB，维纳滤波后SNR=-11.2764dB，K-WPC滤波后SNR=10.655dB.多道正演信号处理结果见图 3.

根据单道数据和多道正演信号的处理结果（如图 2，3所示）可知，维纳滤波方法对于强噪声环境下滤除随机噪声能力有限，输出结果中仍含有很强的随机噪声，输出数据同相轴不明显（如图 3c所示），信噪比只是稍稍提高，无法达到去噪处理的要求.对于K-WPC滤波结果，恢复出来的共炮点记录在宏观上可以很好地表现出原始记录中同相轴的位置，在其输出的信号中，随机干扰滤除较彻底，同时同相轴幅值保真性好（如图 3d所示）.

对第二道记录滤波前后含同相轴信号处的对比分析见图 4.混有强随机噪声的合成记录见图 4a，图 4a中标有星号的曲线为Ricker子波，信噪比较低，有效信号在波峰和波谷处因噪声干扰都存在较严重的畸变；对于传统维纳滤波结果，如图 4b所示，此方法只能滤除部分随机噪声，并且在有效信号处产生一定的畸变，且有效信号的峰值点和谷值点误差较大；K-WPC滤波的结果如图 4c所示，有效信号的峰值点和谷值点误差较小，几乎不存在毛刺现象，滤波后的记录接近于原始合成记录.可见，核主分量维纳滤波与传统的维纳滤波方法相比可以更好地消减随机噪声.在信号保幅和减少信号畸变方面得到较大改善.

3.2 实际地震数据处理结果

为了更好地说明K-WPC方法的适用性，将此方法应用于实际地震资料.

选取某炮区共炮点记录，如图 5a所示，该记录存在较强的随机噪声，同相轴不清晰，甚至被随机噪声截断.

分别采用传统维纳滤波和核主分量维纳滤波对图 5a共炮点记录进行处理，滤波效果分别如图 5b，5c所示.由于噪声干扰，传统维纳滤波方法使得数据矩阵求解不精确，随机噪声仍有残留.核主分量维纳滤波处理后，随机噪声消减效果比传统维纳滤波的效果令人满意.图 5d为原始信号与K-WPC处理之后信号的差值，差值图像中主要为被K-WPC滤除的随机噪声.可见，与一般维纳滤波相比，核主分量维纳滤波由于存在一系列核变换、主分量提取、降维处理使得计算得到的滤波因子较一般维纳方法的精准，最终表现为在消减随机噪声方面效果更为明显.从差值图像可以看出，K-WPC在消减随机噪声的同时，具有较好的保幅性.

4 结论

将K-WPC算法应用于地震资料处理、理论研究与仿真实验，证明了该方法对地震资料中的随机噪声有很好的压制作用.将核函数与PCA算法相结合进行矩阵变换和降维，降低了在求解过程中产生病态矩阵的几率，使得矩阵求解精确，滤波性能明显提高.正演模型处理前后信噪比可以从-16. 4622dB左右提高到10.6550dB左右.同时，处理后信号同相轴连续性较好，有效信号的幅值衰减和畸变方面都较一般维纳方法改善很多.