一种相对保幅的低频逆时噪声压制方法及其应用

陈可洋

（中国石油大庆油田有限责任公司勘探开发研究院，黑龙江 大庆 163712）

在复杂构造高精度地震波成像方法中，逆时偏移是理论较为成熟且最精确的成像方法之一。该方法采用双程地震波波动方程，对方程的近似最少，适合于任意陡倾角、速度在纵横向变化较为剧烈的情况，可以解决单程波偏移方法受地层倾角限制的缺陷和克希霍夫偏移方法的多值走时和计算盲区等问题，同时还可实现多次波、回转反射波等通常认为是干扰波类型波场的准确成像（其他方法均将多次波等波场当作一次反射波进行处理，这必然引入较大误差），且能够准确地处理焦散面的多值走时和相位变化问题[1-5]。同时，随着计算机技术的快速发展，特别是基于CPU/GPU高性能集群并行计算技术和大容量磁盘的快速存储技术的出现，较大程度地改善了逆时偏移技术工业化应用的现状[6]。鉴于上述诸多优点，地震波叠前逆时成像技术受到地球物理学界的高度重视和广泛应用。

众所周知，逆时偏移方法区别于其他现行地震成像方法的显著特点：经其处理后的结果上存在较强能量的低频背景噪声。若不对该噪声进行有效压制，则将掩盖有效的、真实的地层细节。地球物理界对此开展了多项研究工作，并取得了丰硕的研究成果，主要可归纳为以下两个方面：一方面对逆时成像条件进行改进。例如：Yoon等[7]提出了在零延迟互相关条件中加入波印廷矢量来消除低频成像噪声。Liu等[8]把全波场分解成上下行单程波波场，并运用分解出的单程波波场结合相关逆时成像条件达到消除低频成像噪声的目的。陈可洋等[9]首次分析了低频噪声成因是由以震源和检波点及它们在速度分界面的投影为椭圆焦点的三条椭圆曲线引入，上下行波场分离逆时成像条件是消除这种噪声的有效方法，并指出速度平滑可以削弱低频背景噪声，但它以损失成像精度为代价。康玮等[10]比较了拉普拉斯算子、波印廷矢量法和上下左右行波分解法消除低频背景噪声的效果。Antoine等[11]采用最小二乘预测误差滤波器消除成像噪声。另一方面是对逆时成像结果进行后续噪声压制处理。例如：Zhang等[12]指出拉普拉斯算子滤波相当于成像波场角度域衰减。刘红伟等[13]指出在叠前对数据进行相位和振幅校正，在成像后运用拉普拉斯算子滤波法消除成像噪声。陈可洋[14]研究了低频噪声是由上行反射波参与相关逆时成像引起，在叠加域、成像点域、共炮点域进行了拉普拉斯算子去噪分析，同时对比了带通滤波方法。陈康等[15]采用四阶拉普拉斯算子有效地去除了低频成像噪声，并使滤波前、后子波的振幅和相位保持相对不变。以上低频噪声压制处理均在深度域实现，低频背景噪声均得到了有效压制。在前人研究基础上，开展了一种全新的低频噪声压制方法探索研究，它是基于热力学传导方程的扩散滤波方法，通过参数试验和实际资料处理，以达到相对保幅低频噪声压制的目的，对复杂构造地震资料高精度逆时成像具有重要的指导意义。

1 基本原理

在数字图像处理领域，Perona和Malik[16]首次提出了用于热力学传导的偏微分方程，即PM方程：

式中：t为扩散时间；div为散度算子；∇是梯度算子；是扩散函数，为一个有界非负的递减函数；U为t时刻的扩散滤波结果，其中U0为t=0时刻的原始数据，即扩散滤波迭代计算的初始条件。

为了将式（1）应用于低频逆时噪声压制处理，对其进行简化。令扩散函数为某一常数g0，此时有，然后对其进行有限差分离散化，得到如下数值计算公式：

式（2）为三维扩散滤波方程，其中，i和Δx分别为x方向的离散网格节点号和网格步长；j和Δy分别为y方向的离散网格节点号和网格步长；k和Δz分别为z方向的离散网格节点号和网格步长；n和Δt分别为时间t方向的离散网格节点和离散步长；N为扩散滤波迭代次数。

对式（2）还可作进一步简化，此时，令空间步长Δx=Δy=Δz，方程系数，忽略空间网格和时间网格步长及扩散滤波系数的影响，得：

式（2）中仅有一个参数，即扩散滤波迭代次数N。现定义经扩散滤波N次迭代后的信号变为UN，原始信号为U0，两者的差为Uerr=U0-UN。经扩散滤波后的信号UN是在原始信号U0基础上进行平滑处理，主要反映原始信号的低频分量，即低频背景噪声，而残差信号Uerr主要反映原始信号的高频分量，即掩盖在低频背景噪声之下的有效地层细节信息。所提方法的相对保幅性将在后续的实际处理试验中进行详细探讨。

2 应用实例

2.1 低频噪声压制效果分析

以FZ地区三维地震资料叠前逆时偏移处理剖面为例，采用3种低频噪声压制方法和2组文中扩散滤波参数（迭代次数N）进行三维数据体去噪试验和应用效果对比分析。图1（a）为经照明补偿后的叠前逆时偏移叠加剖面，其深度域频谱如图2（a）所示。图1（b1）和图1（b2）分别为图1（a）经邻点加权平滑滤波处理后的结果及分离出的低频噪声，对应的深度域频谱分别为图2（b1）和图2（b2）。图1（c1）和图1（c2）分别为图1（a）经文中扩散滤波方法处理后（迭代次数为30次）的结果及分离出的低频噪声，对应的深度域频谱分别为图2（c1）和图2（c2）。图1（d1）和图1（d2）分别为图1（a）经文中扩散滤波方法处理后（迭代次数为200次）的结果及分离出的低频噪声，对应的深度域频谱分别为图2（d1）和图2（d2）。图1（e1）和图1（e2）分别为图1（a）经高通滤波处理后（斜坡为0.1～2 Hz）的结果及分离出的低频噪声，对应的深度域频谱分别为图2（e1）和图2（e2）。

分析图1（a）可知，逆时偏移剖面上存在较强能量的低频背景噪声，掩盖了有效的地层细节信息，且主要地震能量集中于频谱的低频端（图2（a）中粉红色箭头所示）。从图1（b1）和图1（c1）分析可知，这两种方法的低频噪声压制效果相当，地层细节特征得到有效刻画，低频能量也得到有效压制（对比图2（b）和图2（c）），但在去噪后剖面里均存在一定能量的有效信号的残留（图1（b2）和图1（c2）），其中，在经邻点加权平滑滤波处理后的剖面里存在更多的高频信息（图1粉红色椭圆）。对比图1（c1）和图1（d1）可知，形成这两种去噪效果差异的原因仅与文中扩散滤波方法的迭代次数大小有关，迭代次数越大，低频噪声压制程度越轻，这一点也可从图1（c2）和图1（d2）的残差剖面对比分析得到。其中，图1（d2）主要为压制掉低频背景逆时噪声，且无明显的有效信号残留。从频谱分析可知，图2（d1）与图2（c1）具有相似的主频，但前者频带更宽，这一点也可从对应扩散滤波剖面的频谱分析得到，此时图2（d2）中的频谱能量要比图2（c2）更集中于低频端，因此，保留的高频信息更多。分析图1（d1）和图1（e1）可知，文中扩散滤波方法与高通滤波方法具有相似的噪声压制效果。对比频谱可知，这两种去噪方法的主频相当，其中图2（d1）的频谱能量分布要比图2（e1）更广，频带也更宽（图2鲜绿色椭圆）。分析图1（d2）和图1（e2）可知，由于高通滤波方法采用的是逐道处理，该方法处理剖面存在横向能量不均过度的问题（图1蓝色和红色椭圆），这正是带通滤波方法模糊陡倾角成像波场的主要原因，而文中所提的扩散滤波方法则是全局去噪方法。因此，不存在上述问题，且去噪后的地震剖面横向能量过渡自然，有效信号的保真性和保幅性更好。

图1 实际地震资料叠前逆时偏移低频噪声压制前后及其残差剖面Fig.1 Practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression sections and their residual sections

图2 实际地震资料叠前逆时偏移低频噪声压制前后及其残差的深度域归一化频谱Fig.2 Depth normalized amplitude spectrum of practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression and their residual section spectrum

2.2 相对保幅性分析

为了进一步验证文中所提方法的相对保幅性，从第2.1小节3种方法的三维处理数据体满覆盖位置提取出一道逆时偏移记录（图3）和1 km深度处的水平切片，同时从水平切片中提取出一条满覆盖位置的主测线幅值（图4）进行相对保幅性分析。

图3 实际地震资料叠前逆时偏移低频噪声压制前后深度方向波形曲线及其残差对比Fig.3 Wave in the depth direction of practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression and their residual curves comparison

图3（a）和图4（a）分别为原始逆时偏移深度域信号波形和水平切片振幅曲线（com，黑色）和经不同去噪方法处理后的有效信号，其中，smo（粉红色）代表邻点加权平滑滤波处理结果，dif（蓝色）代表扩散滤波方法处理后（迭代次数为30次）的结果，dif1（红色）代表扩散滤波方法处理后（迭代次数为200次）的结果，bp（深黄色）代表高通滤波处理后（斜坡为0.1～2 Hz）的结果。图3（b）和图4（b）分别为原始逆时偏移深度域信号波形和与水平切片振幅曲线不同去噪方法处理后的低频噪声，其中com-smo（粉红色）代表邻点加权平滑滤波法去除的低频噪声，com-dif（蓝色）代表扩散滤波方法（迭代次数为30次）去除的低频噪声，com-dif1（红色）代表扩散滤波方法（迭代次数为200次）去除的低频噪声，com-bp（深黄色）代表高通滤波法（斜坡为0.1～2 Hz）去除的低频噪声。

图4 实际地震资料叠前逆时偏移低频噪声压制前后水平切片中某inline方向曲线及其残差对比Fig.4 Wave along an inline direction in horizontal slice of practical seismic data pre-stack reverse-time migration before and after low frequency noise suppression and their residual curves comparison

分析图3（a）可知：由于原始信号振幅较大，原始逆时偏移结果的波形曲线与经不同去噪方法处理后的结果差异较小，但从去除的低频噪声曲线中可以明显地看出不同去噪方法的压制效果（图3（b）），其中文中扩散滤波方法（迭代次数为30次）去噪效果与邻点加权平滑滤波法效果相当，前者更加平滑，与文中扩散滤波方法（迭代次数为200次）效果相比，迭代次数越大，地层细节特征将得到越多的保留，低频波形曲线过度自然光滑，但如果该参数取得过大，则达不到低频噪声压制的目的，此时的去噪结果将保留更多的低频能量。因此，迭代次数的选择需根据试验和处理要求来实现优化选择。而高通滤波法（斜坡为0.1～2 Hz）去除结果存在一定的波形抖动（图3（b）天蓝色椭圆），这与频域斜坡的选取存在一定的影响，若斜坡选取过陡，则吉普斯效应就会出现在处理结果中。虽然高通滤波方法在横向相对保幅性方面存在缺陷，但仍是目前低频噪声压制中的一种高效的方法，因其采用逐道处理方式，需求内存较小，不需要多次迭代计算。

分析图4（a）可知，原始水平切片的振幅曲线与不同去噪方法处理后的结果差异较明显，去噪后的横向能量均比去噪前减弱，其中扩散滤波方法（迭代次数为30次）和邻点加权平滑滤波法的低频噪声压制效果较强，扩散滤波方法（迭代次数为200次）和高通滤波法（斜坡为0.1～2 Hz）较弱，但彼此间的地层细节差异不明显。从去除的低频噪声水平切片振幅曲线中可以明显地看出不同去噪方法的压制效果（图4（b）），其中文中扩散滤波方法的（迭代次数为30次）去噪效果与邻点加权平滑滤波法效果相当，但与文中扩散滤波方法（迭代次数为200次）效果相比，后者的地层细节特征将得到更大程度的保留，其低频波形曲线过度更光滑且自然，而在高通滤波法（斜坡为0.1～2 Hz）的去除结果里存在一定的波形抖动（图4（b）绿色椭圆），这与其逐道处理方式有关，导致了处理结果在横向能量不保幅的问题。

3 结束语

1）引入一种基于热力学传导的扩散滤波方法用于有效削弱或消除低频逆时噪声，在离散方程系数一定情况下，其简化计算公式仅受迭代次数一个参数影响，不同迭代次数实现不同低频噪声压制程度的去噪效果。因此，根据参数试验和处理要求可实现迭代次数的优化选择，同时该方法简单易实现，且可以推广应用于二维或一维数据的低频噪声压制中。

2）分析三维地震资料逆时偏移结果的低频噪声压制效果表明：三种方法均可以有效压制低频噪声。文中所提的扩散滤波法具有最佳的相对保幅性，去噪结果的频带也更宽，主频段谱能量更均匀。邻点加权平滑滤波法的低频噪声压制程度较强，其去噪效果与较小扩散滤波迭代次数情况等效。

3）高通滤波法具有逐道快速低频噪声压制的特点，仍是一种有效的低频噪声压制方法，其低频噪声压制程度可控（由滤波参数控制），但在三维数据体的横向和纵向上的相对保幅性较差，易受吉普斯效应的影响，同时对陡倾角成像波场存在一定的模糊性。

[1]陆基孟.地震勘探原理（上、下册）[M].北京：中国石油大学出版社，2004.

[2]凌云.地震数据采集、处理解释一体化实践与探索[M].北京：石油工业出版社，2007.

[3]渥·伊尔马滋.地震资料分析——地震资料处理、反演和解释（上、下册）[M].刘怀山，王克斌，童思友，等译.北京：石油工业出版社，2006.

[4]陈可洋.高阶弹性波波动方程正演模拟及逆时偏移成像研究[D].大庆：大庆石油学院,2009.

[5]陈可洋.地震波逆时偏移方法研究综述[J].勘探地球物理进展，2010，33（3）：153-159.

[6]刘红伟，李博，刘洪，等.地震叠前逆时偏移高阶有限差分算法及GPU实现[J].地球物理学报，2010，63（7）：1725-1733.

[7]Yoon K,Marfurt K J.Reverse-time migration using the poynting vector[J].Exploration Geophysics,2006,37(1):102-107.

[8]Liu F Q,Zhang G,Morton S A,et al.Reverse-time migration using one way wavefield imaging condition[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:2170-2174.

[9]陈可洋，吴清岭，范兴才，等.地震波叠前逆时偏移脉冲响应研究与应用[J].石油物探，2013，52（2）：163-170.

[10]康玮，程玖兵.叠前逆时偏移假象去除方法[J].地球物理学进展，2012，27（3）：1163-1172.

[11]Antoine G,Bruno K,Biondo B.Least square attenuation of reverse time migration artifacts[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006:2348-2352.

[12]Zhang Y,James S.Practical issues of reverse time migration:true amplitude gathers,noise removal and harmonic source encoding[C].CPS/SEG Beijing 2009 International Geophysical Cofference&Exposition,2009:5.

[13]刘红伟，刘洪，邹振，等.地震叠前逆时偏移中的去噪与存储[J].地球物理学报，2010，63（9）：2171-2180.

[14]陈可洋.基于拉普拉斯算子的叠前逆时噪声压制方法[J].岩性油气藏，2011，23（5）：87-95.

[15]陈康，吴国忱.逆时偏移拉普拉斯算子滤波改进算法[J].石油地球物理勘探，2012，47（2）：249-255.

[16]Perona P,Malik J.Scale-space and edge detection using anisotropic diffusion[J].IEEE Tran Pattern Analysis and Machine Intelligence,1990,12(7):629-639.