基于弱反射方程和特征值分解的叠前逆时偏移噪声压制
张大洲，熊章强
(中南大学 地球科学与信息物理学院，有色金属成矿预测教育部重点实验室，湖南 长沙，410083)
摘要：采用双程波动方程的逆时偏移可以使成像不受构造倾角限制且精度较高，但双程波动方程的采用使地震波在分界面处产生反射，在偏移过程中产生低频干扰，尤其在浅层更为严重，从而使得成像剖面中出现噪声。通过研究弱反射方程和特征值分解法，采用基于特征值分解的一阶压力-速度弱反射方程压制叠前逆时偏移噪声。通过理论模型分析表明：该方法不但保持了偏移成像不受构造倾角限制且精度较高的优点，而且有效地压制了界面反射波，从而消除逆时偏移过程中产生的低频干扰，使成像效果有显著的改善。
关键词：叠前逆时偏移；弱反射波方程；特征值分解；噪声压制
中图分类号：P631 文献标志码：A 文章编号：1672-7207(2013)08-3363-07
Suppress noise in pre-stack reverse time migration base on weak reflection equation and eigenvalue decomposition
ZHANG Dazhou, XIONG Zhangqiang
(Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education,
School of Geosciences and Info-Physics, Central South University, Changsha 410083, China)
Abstract: Pre-stack reverse-time migration using two-way wave equation can bring high image accuracy and is not restricted by structural obliquity. But owing to using two-way wave equation, seismic waves in the medium have the reflection at the interface, engendering low-frequency disturbance during migration, especially in shallow lays, and producing some noise in migration section. We propose the method, that is, the first-order weak reflection equation of pressure-speed is combined with eigenvalue decomposition, through the research of both the weak reflection equation and eigenvalue decomposition methods. The theoretical model analyses show that, this approach can not only get high accuracy imaging without being restricted by structural obliquity, but also effectively suppress the interface reflected wave in the process of reverse-time migration, consequently eliminating low-frequency disturbance generated by the interface reflected wave in the imaging. So it has a very marked improvement to the effect of imaging.
Key words: pre-stack reverse-time migration; weak reflected wave equation; eigenvalue decomposition; suppress noise
叠前逆时偏移作为一种重要的地震偏移方法被广泛的应用于地震勘探数据处理。与其他偏移方法相比，该方法是在时间轴上实现波场外推，可以看作是沿时间反向的正演过程。传统的偏移方法是利用单程波沿深度方向延拓，而叠前逆时偏移则是基于双程波沿时间轴外推，由于这种方法允许地震波在全方位传播，因而不受构造倾角的限制。逆时偏移最早是由Whitemore[1]提出的，此后Baysal等[2]和 Loewenthal等[3]分别从不同角度对逆时偏移的原理和应用进行详细的论述，并在叠后逆时偏移应用中取得了较好的效果。Chang 和McMechan[4]以及Sun 和 McMechan[5]分别对声波和弹性叠前逆时偏移进行了研究，Dong和McMechan[6]研究了叠前逆时偏移方法,并应用到实际资料处理。由于叠前逆时偏移采用双程波动方程，从而保证了偏移有较高精度且不受倾角限制，这也是叠前逆时偏移的优势所在。但在波场外推过程中正是由于采用双程波动方程，因而在界面处会产生反射。由叠前逆时偏移原理可知，在波场外推过程中仅需要透射波，如果界面反射波存在，则在偏移过程中就会产生低频干扰，从而使得成像剖面中出现噪声。这种干扰在浅层以及波阻抗对比强烈处尤为严重。近年来针对这一问题有许多不同的处理方法，如Baysal等[7]，王山山等[8]和吴清岭等[9]利用无反射波动方程来实现压制反射波的目的，其实这里称为无反射波动方程是不严谨的，因为采用这种波动方程并不能完全消除界面反射，仅仅是对界面反射波能量有一定的削弱，因此称为弱反射波方程更为准确；Guitton等[10]对偏移结果做最小平方滤波来削弱这种噪声，但并不能完全消除；Mulder等[11]通过结合高通滤波与迭代偏移等处理手段使偏移效果得到改善，但低频干扰与有效信号的频带相互叠合，要想在不损失有效信号的前提下做滤波处理来压制噪声是非常困难的，而且与迭代偏移的结合又会增加处理成本；Fletcher等[12]通过在产生噪声的分界面处增加一衰减层以达到去除噪声的目的，但这一处理需要人工拾取该分界面；Xie等[13]通过用一个单程波算子把全波场分离成不同方向传播的波来消除噪声，但由于单程波算子基于射线理论构造，因而不能适应复杂地质情况；Yoon等[14]和杜启振[15]通过在成像过程中引入Poynting矢量来压制界面反射波，这种方法是通过计算Poynting矢量来判断波传播方向，然后利用方向滤波来消除反射波，如果模型简单时效果较好，但是当模型复杂时反射波和透射波相互重叠，因此不能将两者很好的分离。本文作者采用一种基于特征值分解的方法来压制反射波，该方法在波场传播过程中完全遵守全波动方程，仅在成像过程中利用特征值分解来提取压制反射波后的波场值，在此基础之上将该方法与一阶压力-速度弱反射波动方程结合得到一种新的基于特征值分解的弱反射方程，这样就使得压制效果得到更大的提高。通过理论模型分析，该方法可以较大地压制叠前逆时偏移过程中波场外推时界面反射波，从而减小偏移过程中低频干扰对于逆时偏移结果的影响。从理论偏移模型分析可知，应用此方法后叠前逆时偏移剖面中的噪声被压制，成像效果得到增强，从而证明了该方法的有效性。
1 方法原理
1.1 一阶压力-速度弱反射波动方程推导
二维非均匀介质声波波方程为：
(1)
其中：u=(ux，uz)为位移场，Vp为介质的纵波速度，ρ为密度。
由于波阻抗Zp=ρVp，如果要在界面上不产生反射，则波阻抗Zp应为常数，当Zp为常数时，式(1)变为：
(2)
将方程(2)两边除ρ并展开为一阶压力-速度方程如下：
(3)
式中：p为压力场；vx和vz为速度场。式(3)就是一阶压力-速度弱反射方程。
1.2 基于特征值分解的反射波压制方法
二维情况下，非均匀各向同性介质中声波一阶压力-速度方程可表示为：
(4)
式中：弹性模量
；p(x，z)为压力场；vx和vz为质点在x和z方向上的速度分量。
式(4)可表示为一阶双曲型标量方程，即：
(5)
其中：
；系数矩阵A和B分别为：
，
(6)
由于式(5)为双曲型方程，因此一定存在一个非奇异矩阵LA，使得下式成立
(7)
式中：N为矩阵A的特征值；LA为特征值对应的特征向量；
为特征向量的逆。
对于一维声波方程
可变换为
(8)
对式(8)进一步变换得到：
(9)
将式(9)写成分量形式可得：
(10)
式(10)描述的是以特征速度
传播的波，至于传播方向则由的正负号来决定，其振幅为
[16]。由于波传播的方向可以由特征值的符号决定，这给压制反射波提供了条件，只要根据反射波的方向使
，这样就削弱了反射波。
当需要削弱z方向的反射波时，仅对B系数矩阵进行变换。矩阵B的特征值依次由大到小为
，
，
，由特征值计算得到的特征向量及其逆为：
(11)
这样式(5)可以改写为：
(12)
设N=N1+N2，其中 
，
则
，
就地面地震而言，反射波是上行波，若要吸收上行波，则波动方程为
，将上述两式展开并仅取计算压力的一项后如下：
(13)
式中：p1(x，z)为成像所用的压力场；vx，vz和p(x，z)为前一个时间间隔的波场。这样对于整个波场外推的波动方程为：
(14)
通过前一个时间间隔的波场，利用式(13)就可以提取到用于成像的压制了反射波的波场值。
1.3 弱反射方程结合特征值分解的反射波压制方法
将式(3)表示一阶双曲型标量方程可得：
(15)
其中
，吸收矩阵A1和B1分别为：
(16)
通过计算特征值及其向量并整理后得到分量公式如下：
(17)
利用式(17)就可以得到削弱反射波后的波场，式中前3个为弱反射方程，第4个为利用前3个方程计算出的弱反射波场进一步计算得到用于成像的波场值，通过2种方法的结合使得反射波得到较大的压制。
2 反射波压制效果分析
为了测试反射波的压制效果，本文特设计一个两层介质模型：模型尺寸长和高分别为400 m和300 m，分界面在120 m处，上层波速度为3 000 m/s，下层波速为4 000 m/s。空间步长为2 m，时间步长为0.2 ms，震源为100 Hz的雷克子波，震源位置在(200，0)处。图1(a)所示为未加任何压制方法的波场快照(t=65 ms)，从图中可以看出当波传播到分界面时发生透射和反射，相对来说透射波能量比反射波强；图1(b)所示为利用弱反射方程进行波场外推时的波场快照，从图1可以看出：当波传播到分界面时仍然发生透射和反射，但是与图1(a)相比反射波能量有所减弱；图1(c)所示为利用特征值分解方法来压制反射波后的波场快照，从图中可以看出同样反射波被压制，且吸收效果比弱反射方程好；图1(d)所示为利用弱反射方程结合特征值分解对反射波压制后的波场快照，从图1还可以看出：利用这种方法后反射波几乎被完全消除，在波场中仅有透射波，这样就有非常利于叠前逆时偏移成像。
对于以上3种方法压制反射波的效果通过波场快照做了初步分析，下面从波场快照中提取一道数据进行进一步分析。图2所示为X=200 m处4种不同情况下的波形对比，从图2(a)可以看出：4种情况下透射波的能量是相同的，但是反射波能量得到不同程度的压制。图2(b)所示为X=200 m处的反射波能量对比，从图2(b)可以看出未采用任何反射波压制时最大振幅为0.080 23；采用弱反射方程时反射波最大振幅为0.016 14，80.39%的反射波被压制；采用特征值分解吸收反射波后反射波最大振幅为0.006 664 5，91.93%的反射波被压制；当采用弱反射方程结合特征值分解吸收反射波后反射波最大振幅为0.001 677，97.96%的反射波被压制。从以上的数值分析来看，利用弱反射方程结合特征值分解可以很好的压制界面反射波。
图1 不同方法压制反射波波场快照对比(t=65 ms)
Fig. 1 Snapshot of different methods to suppress reflected wave (t=65 ms)
图2 X=120 m处不同反射波压制方法单道波形对比
Fig. 2 Comparison of single-channel waveform which use different reflection wave suppression methods at X=120 m
3 理论模型偏移分析
3.1 两层模型
为了说明反射波在逆时偏移中的影响以及压制反射波后成像效果，特设计一两层模型，模型长和高都为2 000 m，空间采用间隔为10 m，时间采样间隔0.5 ms，震源采用50 Hz的雷克子波。图3所示为两层模型逆时偏移结果对比，图3(a)所示为没有压制反射波时逆时偏移剖面，从图3(a)可以看出尽管界面被很好地成像，但在界面上层有很强的噪声(图中显示白色部分)。图3(b)所示为采用弱反射方程结合特征值分解来压制反射波后的偏移剖面，从图3(b)可以看出成像剖面中不但界面被很好地成像，噪声也被很好的消除了。从以上模型分析可以得出：反射波在叠前逆时偏移过程中会产生低频干扰使得偏移成像剖面中出现噪声，利用弱反射方程结合特征值分解的方法可以很好地压制偏移噪声，从而得到较好的成像剖面。
图3 两层模型叠前逆时偏移成像剖面对比
Fig. 3 Comparison of pre-stack reverse-time migration section in two-layer model
3.2 多层断层模型
多层断层模型如图4(a)所示，模型长和高都为1 500 m，空间采样间隔为5 m，时间采样间隔为0.5 ms，震源为30 Hz的雷克子波。图4(b)所示为没有压制反射波时逆时偏移剖面，从图4(b)可以看出尽管各界面被成像，但在模型上半部分存在噪声，尤其在浅部噪声更大。图4(c)所示为采用弱反射方程结合特征值分解来压制反射波后的偏移剖面，从图4(c)可以看出无论是浅部还是深度都没有噪声干扰，地层各界面非常清晰，尤其是大倾角的断层面也被很好的成像(模型下部断面由于采集时反射波不能被地面检波器接收，故未能成像)，这也充分说明本文采用的反射波压制方法没有改变地震波全方位传播的特性，因此能够对陡倾角界面成像。与模型相比，偏移剖面的构造位置准确，取得了较好的成像效果。
图4 多层断层模型叠前逆时偏移
Fig. 4 Multi-layer fault model pre-stack reverse migration
4 结论
(1) 在叠前逆时偏移过程中界面反射波会产生低频干扰，从而使得成像剖面出现偏移噪声，尤其在浅部更为严重。
(2) 在波场外推过程中利用所推导的弱反射方程结合特征值分解方法可以有效地压制界面的反射波。
(3) 在叠前逆时偏移过程中引入弱反射方程结合特征值分解的方法可以很好的压制界面反射波，从而使得成像剖面不受低频干扰，减少了偏移噪声，大大改善了逆时偏移成像的效果。
参考文献：
[1] Whitmore N D. Iterative depth migration by backward time propagation[C]//53th SEG Annual Meeting. Las Vegas, 1983： 382-385.
[2] Baysal E, Kosloff D D. Sherwood J W C. Reverse time migration[J]. Geophysics, 1983, 48(11): 1514-1524.
[3] Loewenthal D, Mufti I R. Reverse time migration in spatial frequency domain[J]. Geophysics, 1983, 48(5): 627-635.
[4] Chang W F, McMechan G A. Reverse migration of offset vertical seismic profiling data using the excitation time imaging condition[J]. Geophysics, 1986, 51(1): 67-84.
[5] Sun R, McMechan G A. Pre-stack reverse-time migration for elastic with application to synthetic offset vertical seismic profiles[J]. Proceedings of the IEEE, 1986, 74(3): 457-465.
[6] DONG Zhengxin, McMechan A. 3-D prestack migration in anisotropic media[J]. Geophysics, 1993, 58(1): 79-90.
[7] Baysal E, Kosloff D D, Sherwood J. A two-way nonreflecting wave equation[J]. Geophysics, 1984, 49(2): 132-141.
[8] 王山山, 聂勋碧, 邓玉琼. 混合法双程无反射波动方程偏移[J]. 石油地球物理勘探, 1993, 28(5): 537-544.
WANG Shanshan, NIE Xunbi, DENG Yuqiong. Hybrid migration using two-way nonrefleeting wave equation[J]. Oil Geophysical Prospecting, 1993, 28(5): 537-544.
[9] 吴清岭, 李伟, 陈可为, 等. 无反射波动方程正演模型及其应用[J]. 石油地球物理勘探, 2000, 35(2): 147-154.
WU Qingling, LI Wei, CHEN Kewei, et al. Reflection-free wave equation modeling and its application[J]. Oil Geophysical Prospecting, 2000, 35(2): 147-253.
[10] Guitton A, Kaelin B, Biondi B. Least-squares attenuation of reverse-time-migration artifacts[J]. Geophysics, 2007, 72(1): S19-S23.
[11] Mulder W A, Plessix R E. A comparison between one-way and two-way wave-equation migration[J]. Geophysics, 2006, 69(6): 1491-1504.
[12] Fletcher R F, Fowler P, Kitchenside P, et al. Suppressing artifacts in prestack reverse time migration[C]//75th SEG Annual Meeting. 2005: 2049-2052.
[13] XIE Xiaobi, WU Rushan. A depth migration method based on the full-wave reverse-time calculation and local one-way propagation[C]//76th SEG Annual Meeting. New Orleans, 2006: 2333-2337.
[14] Yoon K, Marfurt K J. Revesrse-time migration using the pointing vector[J]. Exploration Geophysics, 2006, 59(1): 102-107.
[15] 杜启振, 秦童. 横向各向同性介质弹性波多分量叠前逆时偏移[J]. 地球物理学报, 2009, 52(3): 501-807.
DU Qizheng, QIN Tong. Multicomponent prestack reverse-time migration of elastic waves in transverse isotropic medium[J]. Chines J Geophys, 2009, 52(3): 501-807.
[16] 董良国. 弹性波数值模拟中的吸收边界条件[J]. 石油地球物理勘探, 1999, 34(1): 45-56.
DONG Liangguo. Absorptive boundary condition in elastic-wave numerical modeling[J]. Oil Geophysical Prospecting, 1999, 34(1): 45-56.
(编辑 陈爱华)
收稿日期：2012-09-16；修回日期：2012-12-20
基金项目：国家自然科学基金资助项目(41004047，41274123)；中南大学自由探索项目
通信作者：张大洲(1979-)，男，甘肃榆中人，博士，讲师，从事面波勘探及地震偏移成像等方面研究；电话：0731-88830451；E-mail：dazhou2005@163.com
