CN102043165A - 基于基追踪算法的面波分离与压制方法 - Google Patents

Abstract

本发明涉及一种基于基追踪算法的面波分离与压制方法；该方法首先根据地震剖面的实际频率范围及地震信号的结构特征自适应地建立过完备原子库，利用基追踪算法对地震剖面的每道地震信号进行稀疏分解，并从中筛选出能近似重建地震信号的重要原子，然后，根据筛选出的重要原子自适应地设置能量阈值和主频阈值范围，以从每道地震信号分解结果中找出组成面波的原子，用地震信号减去面波原子重建的面波成分得到压制面波后的地震信号；该方法具有自适应匹配信号的特点，能够同时体现地震信号时域和频域的局部特性，与一维傅里叶变换滤波压制面波的方法相比，能有效地压制面波，同时对有效信号造成的损伤较小。

Description

基于基追踪算法的面波分离与压制方法

技术领域

本发明涉及一种用于压制地震资料面波干扰的面波基追踪算法分离与压制方法。

背景技术

面波是陆上地震勘探遇到最多的相干干扰波。它的特点是：能量强；频率低，一般为几赫兹～30赫兹；速度低，为横波传播速度的0.9倍，视速度一般为100～1000m/s，其中以200～500m/s的视速度最为常见。面波的时距曲线为直线。由于介质的频散效应和吸收衰减作用，面波的能量会随着深度和传播距离的增加而迅速衰减，同时，其频率由高向低变化，形成“扫帚状”，即发生频散。在地震资料处理中通常利用面波与反射波的频率差异和速度差异来去除面波。过去采用的模拟滤波器进行滤波存在严重的缺点，如结构复杂、改变特性较困难、存在不需要的相移等。后来广泛采用了数字滤波方法来消除它，如一维滤波(康冶，于承业，贾卧等.f-x域去噪方法研究[J].石油地球物理勘探，2003，38(2)：136-138)、F-K滤波(闫立志，景新义，李刚.F-K滤波在噪音减去法中的作用[J].海洋地质动态，2006，22(10)：28-32)、τ-p变换(曾有良.Radon变换波场分离技术研究.中国石油大学硕士研究生学位论文，2007.4)和小波变换(罗国安，杜世通.小波变换及信号重建在压制面波中的应用[J].石油地球物理勘探，1996，31(3)：337-349)等。这些方法虽然见到了一定的处理效果，但仍然存在两方面的不足：①经滤波处理后，受面波干扰影响的地震信号与未受影响的地震信号能量相差悬殊。若在叠前不进行必要的处理，则对速度参数的提取、对振幅随炮检距变化关系的处理以及叠加等均会带来较大的误差，而最终影响地震剖面的质量。②在消除面波干扰的同时也消掉了与面波干扰同频带的有效信号。这对于分析中、深层反射信号和提高地震记录的分辨率都是极为不利的。

S.Mallat和Z.Zhang于1993年首次提出了信号在过完备库上分解的思想，并引入了匹配追踪算法(Matching Pursuit，MP)(Mallat S，Zhang Z.Matching pursuits withtime-frequency dictionaries[J].IEEE Trans Signal Process.1993，41(12)：3397-3415)，1998年，D.Donoho等人又另辟蹊径，提出了基追踪算法(Basis Pursuit，BP)(S.Chen，D.Donoho and M.Saunders，“Atomic decomposition by basis pursuit”，SIAM Journal onScientific Computing，1998，20(1)：33-61)，此方法采用表示系数的范数作为信号稀疏性的度量，通过最小化l1范数将信号稀疏表示问题定义为一类有约束的极值问题，进而转化为线性规划问题进行求解。此方法提出后，随之就出现了相关的应用，如：信号去噪、特征识别和医学图像处理等，尽管已经有学者将基追踪算法应用到了地震资料去噪方面，但大多是去除随机噪声，对于将其应用到面波的识别及提取就少之甚少。根据地震信号的特性将此方法应用到地震资料面波压制中，为获得高分辨率的地震资料提供了新的研究方向。

发明内容

本发明的目的是提供一种基于基追踪算法的面波分离与压制方法，该方法能够有效的压制面波，并且在压制面波的同时对有效信号的损伤较小。

本发明所述的基于基追踪算法的面波分离与压制方法，包括如下步骤：

(1)Ricker原子过完备原子库的形成：

按时间采样间隔ΔT得到大小为M×N的地震剖面F，其中M为每道地震信号的采样点数，N为地震剖面的道数。对地震剖面的每道地震信号{F_i，i＝1，2，3，K，N}做离散傅里叶变换，利用傅里叶变换频谱图并根据ΔT确定地震剖面F的模拟频率范围为0～f_max。选用Ricker原子构造稀疏分解用的过完备原子库D，一个Ricker原子的表达式如下：

$g_r\left(t\right)={(1-2{\mathrm{\π}}^2{f}^2{(t-u)}^2)e}^{-{\mathrm{\π}}^2{f}^2{\left(\text{t-u}\right)}^2}$

其中，t＝0，ΔT，2ΔT，...，(M-1)ΔT，主频参数f为区间(0，f_max]上的整数，位移参数u为区间[0，M-1]上的整数。n＝M×f_max，由n个原子g_γ构成过完备原子库：

D＝{g_γ|γ＝(f，u)，f＝l，2，L，f_max，u＝0，1，L，M-1}

(2)利用基追踪算法对地震剖面F的每道地震信号进行稀疏分解

对地震剖面F的每道地震信号{F_i，i＝1，2，3，K，N}利用基追踪算法进行稀疏分解，并记录每道地震信号Fi分解后得到的原子所对应的分解系数和参数：

$\{(P_{\mathrm{\γij}},{\mathrm{\γ}}_{\mathrm{ij}})=(P_{\mathrm{\γij}},f_{\mathrm{ij}},u_{\mathrm{ij}})|i=\mathrm{1,2,3},K,N;j=\mathrm{1,2,3},K,n\}$

其中，j为过完备原子库中原子的序号，是原子所对应的分解系数；

(3)重要原子筛选

对第(2)步的分解结果进行重要原子筛选，筛选过程如下：设置系数阈值初始值

找出并记录下每道地震信号Fi分解结果中分解系数满足

的原子所对应的分解系数和参数：

$\{(P_{\mathrm{\γij}}^{\′},{\mathrm{\γ}}_{\mathrm{ij}}^{\′})=(P_{\mathrm{\γij}}^{\′},f_{\mathrm{ij}}^{\′},u_{\mathrm{ij}}^{\′})|\left|P_{\mathrm{\γij}}^{\′}\right|>P_r^{\′},j=\mathrm{1,2},\mathrm{\Λ},n_i^{\′}\}$

其中n′_i为从第i道分解结果中筛选出来的原子个数。计算

F′i为重建后

第i道地震信号，设置重建误差

如果mis＜0.01％，

结束筛选过程，分解系数满足

的原子作为第i道地震信号F_i筛选出的重要原子；否则

重复以上过程。

(4)设置主频阈值f_λ和能量阈值p_λ

从第(3)步筛选出来的重要原子中求出绝对值最大的分解系数对应原子的主频值f_m，根据面波的频谱分布特性设置主频阈值范围f_λ＝(2±1)f_m。再求出第三步筛选出的分解系数

绝对值的平均值

设置能量阈值p_λ＝p_a；

(5)面波的分离与压制

根据(4)步中设定的主频阈值f_λ和能量阈值p_λ，对(3)步中筛选出来的每一道地震信号F_i的重要原子进行面波对应原子的筛选，取原子

为组成第i道地震信号中面波的原子，其中对应的分解系数和主频参数满足

由原子

重建建第i道面波成份

用原始剖面的第i道地震信号减去该道的面波成份F″_i得到压制面波后的第i道地震信号F′″_i＝F_i-F″_i。

与现有技术相比，本发明的有益效果是：

一、本发明将基追踪算法应用到面波的分离与压制中，地震记录主要是由许多子波叠加而成，而零相位的Ricker子波具有很好的分辨率，因此选用零相位的Ricker子波来构造原子库。由于基追踪算法能将信号分解为原子库中一组原子的最佳线性组合，而面波具有低频高能量的特性，因此将地震信号利用基追踪算法进行稀疏分解时，分解得到的每个原子都有其实际意义，原子的主频参数较小且对应的系数大的原子为组成面波的原子，原子的主频参数较大且对应的系数小的原子为组成有效信号的原子。

二、本发明在设置区分面波与有效波的阈值时，考虑面波低频高能量的特性，根据基追踪算法分解结果中绝对值最大的分解系数对应原子的主频值来设置主频阈值fλ，而能量阈值pλ的设置，则通过求筛选出来的重要原子对应的分解系数绝对值的平均值来确定，以达到既能有效的压制面波，又能减小对有效信号的损失的效果。

三、本发明在构造原子库时，根据实际地震剖面的频率范围来确定原子库中原子的主频参数范围，既符合实际情况，又避免了不必要的计算量。

四、本发明是基于稀疏分解的原理，与传统的一维傅里叶滤波的方法相比，它具有自适应匹配信号的特点，能根据信号的结构特征同时体现信号在时域和频域的局部特性。稀疏分解与传统的傅里叶变换有着根本的区别，傅里叶变换使用的基函数是频谱单一、时域无限的复函数e^-jwt，与被分析信号作内积，而稀疏分解使用的基函数是一个具有时频局域性的时频窗函数，其时宽和频宽均可调，而且用基函数组成的原子库是过完备的(这一点有别于小波)。利用传统的一维傅里叶滤波的方法会损失有效信号的低频成份，而利用本方法进行面波的分离与压制时，则会减小有效信号的损失。

总之，本发明方法根据地震信号的结构特征建立自适应原子库，同时考虑面波的频率及能量两方面特性，根据基追踪算法分解结果设置阈值，能有效地压制面波，同时对有效信号造成的损伤较小。

下面结合附图和具体实施方式对本发明作进一步详细说明。

附图说明

图1(a)一个主频为20Hz的Ricker子波的时域波形示意图.

图1(b)一个主频为20Hz的Ricker子波的频谱示意图.

图2合成二维地震记录。

图3本发明对合成地震记录进行面波分离与压制结果

其中：图3(a)为压制面波后的记录；图3(b)为分离出的面波。

图4一维傅里叶滤波对合成记录进行面波分离与压制结果

其中：图4(a)为压制面波后的记录；图4(b)为分离出的面波。

图5实际地震剖面及子剖面频谱图

其中：图5(a)为实际二维地震剖面(3501×354)，图5(b)为从实际二维地震剖面中截取的地震子剖面，大小为256×80，5(c)为地震子剖面的傅里叶变换频谱图。

图6本发明不同阈值下对实际地震剖面进行面波分离与压制结果

其中：图6(a)，图6(b)，图6(c)为能量阈值均取p_m，主频阈值分别取2f_m，2.5f_m，3f_m时压制面波后的地震剖面；图6(d)为主频阈值取2.5f_m时分离出的面波。

图7一维傅里叶滤波对实际地震剖面进行面波分离与压制结果

其中：图7(a)为压制面波后的剖面；图7(b)为分离出的面波

具体实施方式

基于基追踪算法的面波分离与压制方法，包括如下步骤：

(1)Ricker原子过完备原子库的形成：

本实施例合成了二维地震记录以及从实际的地震剖面(3501×354)中截取的地震子剖面(256×80，256代表采样点数、80代表道数)，合成的二维地震记录及地震子剖面时间采样间隔均为ΔT＝2ms。对于合成的二维地震记录，可以依据先验值来设置原子主频参数，对于截取的地震子剖面，先对地震子剖面的每道地震信号做傅里叶变换，将结果显示在一个频谱图中，图6(a)显示了实际地震剖面80道地震信号的傅里叶变换频谱图，从图中确定地震剖面的模拟频率范围0～50Hz，根据此频率范围来设置原子库中原子的主频参数。

确定频率范围后，后续工作就是形成Ricker原子过完备原子库D。Ricker子波的基本形式表示如下：

$g_r\left(t\right)={(1-2{\mathrm{\π}}^2{f}^2{(t-u)}^2)e}^{-{\mathrm{\π}}^2{f}^2{\left(\text{t-u}\right)}^2}$

其中，t＝0，ΔT，2ΔT，...，(M-1)ΔT，根据上述确定的频率范围及截取的地震子剖面的长度，取主频参数f为区间(0，50]上的整数，取位移参数u为区间[0，255]上的整数。n＝256×50，由n个原子g_γ构成过完备原子库：

D＝{g_γ|γ＝(f，u)，f＝1，2，Λ，50；u＝0，1，Λ，255}

(2)利用基追踪算法对地震剖面F的每道地震信号进行稀疏分解

对地震剖面F的每道地震信号{F_i，i＝1，2，3，K，N}利用基追踪算法进行稀疏分解，并记录每道地震信号F_i分解后得到的原子所对应的分解系数和参数：

$\{(P_{\mathrm{\γij}},{\mathrm{\γ}}_{\mathrm{ij}})=(P_{\mathrm{\γij}},f_{\mathrm{ij}},u_{\mathrm{ij}})|i=\mathrm{1,2,3},K,N;j=\mathrm{1,2,3},K,n\}$

其中，j为过完备原子库中原子的序号；

(3)重要原子筛选

对第(2)步的分解结果进行重要原子筛选，筛选过程如下：设置系数阈值初始值

找出并记录下每道地震信号Fi分解结果中分解系数满足的原子所对应的分解系数和参数：

$\{(P_{\mathrm{\γij}}^{\′},{\mathrm{\γ}}_{\mathrm{ij}}^{\′})=(P_{\mathrm{\γij}}^{\′},f_{\mathrm{ij}}^{\′},u_{\mathrm{ij}}^{\′})|\left|P_{\mathrm{\γij}}^{\′}\right|>P_r^{\′},j=\mathrm{1,2},\mathrm{\Λ},n_i^{\′}\}$

其中n′_i为从第i道分解结果中筛选出来的原子个数。计算

F′_i为重建后第i道地震信号，设置重建误差

如果mis＜0.01％，结束筛选过程，分解系数满足

的原子作为第i道地震信号F_i筛选出的重要原子；否则

重复以上过程。

(4)设置主频阈值f_λ和能量阈值p_λ

从第(3)步筛选出来的重要原子中求出绝对值最大的分解系数对应原子的主频值f_m，根据面波的频谱分布特性设置主频阈值范围f_λ＝(2±1)f_m。再求出第三步筛选出的分解系数

绝对值的平均值

设置能量阈值p_λ＝p_a；

(5)面波的分离与压制

根据(4)步中设定的主频阈值f_λ和能量阈值p_λ，对(3)步中筛选出来的每一道地震信号F_i的重要原子进行面波对应原子的筛选，取原子

为组成第i道地震信号中面波的原子，其中

对应的分解系数和主频参数满足

由原子重建建第i道面波成份

用原始剖面的第i道地震信号减去该道的面波成份F″_i得到压制面波后的第i道地震信号F′″_i＝F_i-F″_i。

仿真实验

本发明的效果可以通过以下仿真实验结果进行验证：

为比较本发明与基于一维傅里叶滤波的方法分离和压制面波的效果，我们对合成的二维地震记录以及实际的地震剖面分别做了实验。

实验中采用合成的二维地震记录(256×80)和实际的地震剖面(256×80)；过完备原子库采用Ricker子波构建；实验中与传统的一维傅里叶滤波做对比，压制面波后剖面的质量根据面波的去除效果以及信号的损失程度来衡量。

图1为一个主频为20Hz的Ricker子波的时域波形和频谱示意图，从图1(a)中可以看出，Ricker子波由一个波峰(或称作主瓣)和两个波谷(或称做旁瓣)组成，其两个波谷之间时间间隔ΔT的倒数称为视主频或主频。Ricker子波的波形取决于它的主频，即它的振幅谱的峰值频率；图1(b)为其频谱图，从图中看出其峰值频谱为20Hz。Ricker子波是零相位的，它的波形是对称的，在相同的带宽条件下，与其他相位的子波相比，Ricker子波主瓣窄，旁瓣幅度小，能量集中在较窄的时间范围内，具有较高的分辨率。

图2为合成二维地震记录，第一层和第二层介质中选用反射波的主频分别为50Hz和35Hz，三组面波同相轴从上到下主频分别为12Hz、10Hz、8Hz。

表1为对合成二维地震记录压制面波时选取不同的主频阈值和能量阈值下压制面波后信号的信噪比(dB)。

表1对合成二维地震记录压制面波时选取不同的主频阈值和

能量阈值下压制面波后信号的信噪比(dB)。

从表中可以看出当主频阈值取2.5f_m，能量阈值取p_a时，去除面波后信号的信噪比最高，即去除的效果最好。

图3为本发明对合成地震记录进行面波的分离与压制结果，图3(a)和3(b)为主频阈值取2.5f_m，能量阈值取p_a时压制面波后的记录和分离出的面波，从视觉效果来看，面波基本上被提取干净，且对有效信号的损失较小。

图4为一维傅里叶滤波对合成记录进行面波分离与压制结果，图4(a)为压制面波后的记录；图4(b)为分离出的面波，比较图3(a)和图4(a)中压制面波后的记录，图4(a)中第40道左右的有效信号受到了损失，再比较图3(b)和图4(b)，很明显4(b)中虽然提取出了面波，但同时也能提取出了有效信号，而图3(b)中却没有对有效信号造成任何损失。从对合成记录的处理效果可以看出，基于基追踪算法的面波分离与压制方法在处理效果上明显优于一维傅里叶滤波。

表2为本方法与传统的FFT方法对合成二维地震记录压制面波后的信噪比(这里认为去面波前除面波外的信号为有效信号)。

表2本方法与传统的FFT方法对合成二维地震记录

压制面波后的信噪比

去面波前SNR	基追踪算法去面波后SNR	FFT去面波后SNR
			-3.6370dB	16.7133dB	7.9543dB

从表2可以看出，本发明方法去除面波后的信噪比要比传统的FFT方法高了8.7316dB，这说明了本发明方法与传统的FFT方法相比，更能有效的压制面波，减小有效信号的损失。

图5为实际地震剖面及子剖面频谱图，图5(a)为实际二维地震剖面(3501×354)，图5(b)为从实际二维地震剖面中截取的地震子剖面，大小为256×80，5(c)为地震子剖面的傅里叶变换频谱图，从图5(b)中可以看出，截取的地震子剖面的前十五道受面波的干扰很小，其它道均受到了低频面波的干扰。图5(c)中可以看出地震子剖面的模拟频率范围为0～50Hz。

图6为本发明不同阈值下对实际地震剖面进行面波分离与压制结果，图6(a)，6(b)，6(c)为能量阈值为p_a，主频阈值分别为2f_m，2.5f_m，3f_m时压制面波后的地震剖面，这几幅图中被面波淹没的有效信号同相轴均显示出来了；表3为不同主频阈值下未受面波干扰部分(前十五道)损失信号的均方差。从表3中可以看出，当主频阈值取2f_m时损失的信号最少，但比较图6(a)，6(b)，6(c)，图6(a)中压制面波后的地震剖面质量最差，图6(b)和6(c)中压制面波后的地震剖面质量差不多，但图6(b)对应损失的信号较6(c)要少，故选主频阈值为2.5f_m时能达到最佳效果；6(d)为主频阈值取2.5f_m时分离出的面波。

表3不同主频阈值下信号损失的均方差

fλ	2fm	2.5fm	3fm
				MSE	4.2669e+005	2.1939e+007	3.5839e+007

图7为一维傅里叶滤波对实际地震剖面进行面波分离与压制结果，图7(a)为滤波后的剖面；图7(b)为分离出的面波，比较图6(b)和图7(a)，在第38道至46道之间，图7(a)中明显有效信号损失严重。原始剖面图5(b)中的前十五道地震信号几乎没有受面波干扰，比较图6(d)和图7(b)中分离出的面波的前十五道，图6(d)中提取了少量的信号，而7(b)中却提取了大量的信号，这进一步说明了用一维傅里叶滤波处理时会严重的损失有效信号。

由仿真结果可以看出：基于基追踪算法的面波分离与压制方法，在处理效果方面优于传统的一维傅里叶滤波方法，本发明方法具有自适应匹配信号的特点，在压制面波的同时能够减少对有效信号的损失。

Claims (1)

1.一种基于基追踪算法的面波分离与压制方法，其特征在于：包括如下步骤：

(1)Ricker原子过完备原子库的形成：

按时间采样间隔ΔT得到大小为M×N的地震剖面F，其中M为每道地震信号的采样点数，N为地震剖面的道数。对地震剖面的每道地震信号{F_i，i＝1，2，3，K，N}做离散傅里叶变换，利用傅里叶变换频谱图并根据ΔT确定地震剖面F的模拟频率范围为0～f_max。选用Ricker原子构造稀疏分解用的过完备原子库D，一个Ricker原子的表达式如下：

$g_r\left(t\right)={(1-2{\mathrm{\π}}^2{f}^2{(t-u)}^2)e}^{-{\mathrm{\π}}^2{f}^2{\left(\text{t-u}\right)}^2}$

其中，t＝0，ΔT，2ΔT，...，(M-1)ΔT，主频参数f为区间(0，f_max]上的整数，位移参数u为区间[0，M-1]上的整数。n＝M×f_max，由n个原子g_γ构成过完备原子库：

D＝{g_γ|γ＝(f，u)，f＝1，2，L，f_max，u＝0，1，L，M-1}

(2)利用基追踪算法对地震剖面F的每道地震信号进行稀疏分解

对地震剖面F的每道地震信号{F_i，i＝1，2，3，K，N}利用基追踪算法进行稀疏分解，并记录每道地震信号Fi分解后得到的原子所对应的分解系数和参数：

$\{(P_{\mathrm{\γij}},{\mathrm{\γ}}_{\mathrm{ij}})=(P_{\mathrm{\γij}},f_{\mathrm{ij}},u_{\mathrm{ij}})|i=\mathrm{1,2,3},K,N;j=\mathrm{1,2,3},K,n\}$

其中，j为过完备原子库中原子的序号；

(3)重要原子筛选

对第(2)步的分解结果进行重要原子筛选，筛选过程如下：设置系数阈值初始值

找出并记录下每道地震信号Fi分解结果中分解系数满足

的原子所对应的分解系数和参数：

$\{(P_{\mathrm{\γij}}^{\′},{\mathrm{\γ}}_{\mathrm{ij}}^{\′})=(P_{\mathrm{\γij}}^{\′},f_{\mathrm{ij}}^{\′},u_{\mathrm{ij}}^{\′})|\left|P_{\mathrm{\γij}}^{\′}\right|>P_r^{\′},j=\mathrm{1,2},\mathrm{\Λ},n_i^{\′}\}$

其中n′_i为从第i道分解结果中筛选出来的原子个数。计算

F′_i为重建后第i道地震信号，设置重建误差

如果mis＜0.01％，结束筛选过程，分解系数满足

的原子作为第i道地震信号Fi筛选出的重要原子；否则

重复以上过程。

(4)设置主频阈值f_λ和能量阈值p_λ

从第(3)步筛选出来的重要原子中求出绝对值最大的分解系数对应原子的主频值f_m，根据面波的频谱分布特性设置主频阈值范围f_λ＝(2±1)f_m。再求出第三步筛选出的分解系数

绝对值的平均值

设置能量阈值p_λ＝p_a；

(5)面波的分离与压制

根据(4)步中设定的主频阈值f_λ和能量阈值p_λ，对(3)步中筛选出来的每一道地震信号F_i的重要原子进行面波对应原子的筛选，取原子

为组成第i道地震信号中面波的原子，其中

对应的分解系数和主频参数满足

由原子

重建建第i道面波成份

用原始剖面的第i道地震信号减去该道的面波成份F″_i得到压制面波后的第i道地震信号F′″_i＝F_i-F″_i。

Images