2016, Vol. 33| 基于行波分离和角度域衰减的地震波叠前逆时成像条件 |  [PDF全文] |
由于早期计算机硬件条件不足和计算成本的大幅增加，限制了地震波逆时偏移方法的大规模应用，但该技术是基于双程地震波动方程，具有诸多优势，因此，越来越受到地球物理工作者的广泛重视.而后随着计算机技术的快速发展，特别是基于CPU/GPU高性能并行计算平台的出现和快速发展，较大程度地改善了逆时偏移技术的工业化应用现状，使得大规模逆时偏移快速计算成为现实，并在墨西哥湾、北海等地区的盐下和岩丘侧翼成像方面取得了突破[1-2]，与此同时，该技术也逐步向更为复杂的多波多分量弹性波逆时成像等前沿领域方向发展[3].
针对国内陆地地震资料，在其逆时偏移技术的工业化应用过程中，仍然存在如下几个重要的问题亟待完善[4-8]：1)计算效率；2)偏移噪声；3)针对逆时偏移的保幅预处理；4)复杂构造区速度精细建模等等.本文主要围绕逆时偏移噪声问题开展研究.当前工业化软件中最常用的就是互相关逆时成像条件，其缺点是存在较强能量的低频背景噪声，目前国内外学者针对低频噪声的成因问题已形成一些认识，例如该噪声是有炮检点同方向波场沿着传播路径相关成像而引入等等[9-12].目前压制逆时偏移噪声主要有如下2种思路：1)在成像过程中进行优化，即在炮检点波场延拓成像过程施加一定的约束.例如Fletcher[13]在地震波动方程中引入衰减项压制背向散射波；Xie等[14]利用单程波算子分离不同方向的地震波场；Yan等[15]采用倾斜叠加方法得到不同方向波场，并在局部角度域区分方向进行成像；Yoon[16]通过采用波印廷矢量计算波场传播方向，并选择传播方向相反的炮检点波场进行成像；Baysal等[17]采用无反射波动方程，但其只能去除垂直入射方向形成的噪声；Loewenthal等[18]采用大于波长长度的窗函数对速度模型慢度进行平滑.2)在成像后进行去噪.例如高通滤波方法，计算效率高，但其缺点是存在吉布斯效应，并且由于其是逐道滤波，横向不保幅；导数滤波方法可以去除近零频率成分，但相位发生改变，且高频噪声增强；拉普拉斯去噪方法[19-20]可以压制入射角接近90度时的噪声，但其对其他角度的入射波压制效果较差，且存在不保幅和对噪声敏感、相位发生改变的问题；扩散滤波方法[21]可以实现多尺度、相对保幅、相位保持的低频噪声压制，缺点是与上述方法相比其计算量相对较大，同时参数确定困难等.
综上可知，在逆时成像后的数据上再进行去噪，存在诸多问题.本文的研究重点转为在逆时成像过程中实现逆时噪声的有效压制，并提高成像质量.以二维理论模型为例，采用波印廷矢量实现方向行波(左行波、右行波、上行波和下行波)的分解和地层界面反射角度的计算，首次开展炮点和检波点的上行波和下行波、左行波和右行波两两作互相关计算，优选成像结果并进行组合叠加，同时剔除形成逆时偏移噪声的波场分量，再将反射角度的余弦函数作为炮点和检波点延拓波场的权函数，进一步实现角度域逆时噪声的衰减，最终实现陡倾角复杂构造模型的准确成像，计算结果验证了逆时成像条件的准确有效性.
1 基本理论以二维地震波动方程为例，其偏微分方程[23]
| $ \frac{{{\partial ^2}P}}{{\partial {t^2}}} = v_0^2\left( {\frac{{{\partial ^2}P}}{{\partial {x^2}}} + \frac{{{\partial ^2}P}}{{\partial {z^2}}}} \right) + f\left( {{x_s},{z_{\rm{s}}}} \right), $ | (1) |
式中，P为质点振动位移波场，v0为地震波传播速度，t为传播时间，xs和zs分别为震源在模型x和z空间方向的坐标位置，f为震源函数.在具体的数值离散公式推导时，我们将式(1)变换为一阶双曲型波动方程形式，便于应用时间2阶、空间10阶差分精度的交错网格有限差分法[22]，且根据空间网格大小和数值频散关系[23]来优选最大频率，从而保证计算结果具有较高的数值精度；同时采用正余弦型吸收衰减因子[24]和内外侧镶边思想[25](外侧15个网格点，内侧5个网格点)构建的完全匹配层(PML)吸收边界条件可以有效消除或衰减人为截断边界处的反射干扰波场，从而保证有效模拟区域计算结果的信噪比.本文重点在于优化逆时成像条件实现噪声的压制，因此，波动方程相关的数值离散公式推导在此不再赘述.
Yoon等[16]采用坡印廷矢量来确定波前面的法线方向，即地震波场的传播方向，
常规的相关法逆时成像条件[22]为
| $ {\mathop{\rm Im}\nolimits} \;age1 = \int\limits_t {{P_{\rm{s}}}{P_{\rm{r}}}} , $ | (2) |
式(2)中，Ps和Pr分别代表炮点波场和检波点波场.当引入地层界面反射角度的余弦函数作为炮点和检波点延拓波场的权函数时，式(1)可以进一步修改为[7]
| $ {\mathop{\rm Im}\nolimits} \;age2 = \int\limits_{t,\beta } {\cos \beta \;{P_{\rm{s}}}{P_{\rm{r}}}} . $ | (3) |
上述处理方式是有效的，其原因是低频噪声主要分布于大角度区域[27-31](90°附近)，采用余弦函数形式的权函数可以有效衰减该角度区域的低频偏移噪声.接着，将基于波印廷矢量的方向行波波场分离算子引入到式(3)中，得到如下2种方向行波的相关法逆时成像条件：
上行波和下行波两两作互相关的逆时成像条件
| $ {\mathop{\rm Im}\nolimits} \;age2,1 = \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,up}}}}{P_{{\rm{r,up}}}}} + \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,up}}}}{P_{{\rm{r,down}}}} + \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,down}}}}{P_{{\rm{r,up}}}}} } +\\ \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,down}}}}{P_{{\rm{r,down}}}}} ; $ | (4) |
左行波和右行波两两作互相关的逆时成像条件
| $ {\mathop{\rm Im}\nolimits} \;age2,2 = \int\limits_{t,\beta } {\cos \left( \beta \right){P_{{\rm{s,left}}}}{P_{{\rm{r,left}}}}} + \int\limits_{t,\beta } {\cos \left( \beta \right){P_{{\rm{s,left}}}}{P_{{\rm{r,right}}}}} + \int\limits_{t,\beta } {\cos \left( \beta \right){P_{{\rm{s,right}}}}{P_{{\rm{r,left}}}}} +\\ \int\limits_{t,\beta } {\cos \left( \beta \right){P_{{\rm{s,right}}}}{P_{{\rm{r,right}}}}} . $ | (5) |
式(4)和式(5)中，满足如下关系式：Im age2, 1=Im age2, 2.根据第2小节中的理论模型试验分析表明，不同传播方向炮点和检波点延拓波场的互相关成像是形成低频逆时偏移噪声的主要原因，因此，这部分成像波场分量应该从改进的逆时成像条件(式(3))中去除.而相同传播方向炮点和检波点上行波场和下行波场、左行波场和右行波场两两作互相关成像并组合叠加，可以在有效压制偏移噪声的同时保持对陡倾角的刻画能力，于是，笔者重构式(4)和式(5)中的逆时成像条件并进行组合叠加，得到如下新的逆时成像条件，通过理论模型测试验证了式(6)的准确有效性：
| $ {\mathop{\rm Im}\nolimits} \;age3 = \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,up}}}}{P_{{\rm{r,up}}}}} + \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,down}}}}{P_{{\rm{r,down}}}}} + \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,left}}}}{P_{{\rm{r,left}}}}} +\\ \int\limits_{t,\beta } {\cos \beta {P_{{\rm{s,right}}}}{P_{{\rm{r,right}}}}} . $ | (6) |
为了验证所提优化的逆时成像条件的应用效果，设计了含有高陡倾角的复杂界面速度模型(图 1).模型总大小为2.5 km×1 km，横纵向空间步长均为5 m.模型顶层①和底层⑧厚度分别为150 m和250 m，对应的地震波速度分别为1 800 m·s-1和3 500 m·s-1.模型中间部分包含一组不同倾角的倾斜界面，且倾斜界面自右向左的倾角依次为300、450、600和900，模型中的序号及其对应的介质速度分别为②：2 000 m·s-1、③：2 300 m·s-1、④：2 600 m·s-1、⑤：3 000 m·s-1、⑥：2 700 m·s-1、⑦：2 500 m·s-1，模型中密度为常数，其大小为1.2×103 kg·m-3.采用最大频率为60 Hz的Ricker子波，在地表位置距离模型最左侧布置震源，炮间距为25 m，共激发100炮.正演记录长度为1.5 s，采用的时间步长为0.2 ms，满足计算所需的稳定性条件.同时在整个地表布置检波器，间距为5 m，一共500个检波器.对采集到的正演炮集记录采用第1小节所提的逆时成像条件进行叠前逆时偏移处理研究.
 |
| 图 1 复杂介质速度模型 Fig.1 Velocity model of complex medium |
图 2(a)为采用常规相关法逆时成像条件(式(2))的逆时偏移叠加剖面，该剖面实现了较好的陡倾角成像效果，但剖面中存在较强能量的低频干扰和偏移噪声(黑色箭头所示)，掩盖了地层细节，同时降低了剖面的信噪比.图 2(b)为在图 2(a)的基础上，采用波印廷矢量计算的反射角度的余弦函数作为权函数(式(3))实现角度域逆时噪声衰减后的逆时偏移叠加剖面.对比图 2(a)和图 2(b)可知，后者的低频逆时噪声得到了有效压制，特别是介质②的底界面和介质⑧的顶界面均得到了清晰的成像，信噪比得到了有效提高，断面刻画也更加清晰.但剖面中仍存在一些低频干扰能量(黑色箭头所示).鉴于上述两种逆时成像条件仍无法完全有效地压制逆时偏移噪声，笔者将基于波印廷矢量的行波分离算子引入到上述两种逆时成像条件.
 |
| 图 2 采用两种逆时成像条件的全波场叠前逆时偏移剖面 Fig.2 Full wave field pre-stack reverse-time migration with different imaging conditions |
图 3为炮点和检波点波场中相同传播方向的行波波场(上下行波场)两两作互相关成像后的叠前逆时偏移叠加剖面.其中图 3(b)是在图 3(a)的逆时成像条件基础上实现角度域衰减后的逆时偏移叠加剖面.分析图 3可知，此时逆时偏移结果的信噪比较高，地层刻画清晰，与图 1中速度模型的地层界面位置对应较好，偏移噪声较弱，地层水平界面刻画准确，其中图 3(b)进一步提高了介质②的底界面和介质⑧的顶界面的成像精度.分析还可知，图 3中的逆时偏移结果对直立界面的刻画较差，其能量较弱.
 |
| 图 3 采用两种逆时成像条件的上下行波分离叠前逆时偏移剖面 Fig.3 Up/down one-way wave pre-stack reverse-time migration with different imaging conditions |
综上可知，同传播方向的炮检点上下行波波场两两作互相关成像，可以提高对水平界面的成像能力，但对陡倾角界面的成像贡献较小.对比图 2和图 3还可知，不同传播方向的炮检点上下行波波场作互相关成像是形成偏移噪声的主要原因.
于是，笔者根据方向行波波场两两作互相关成像结果，筛选出有利于提高对复杂构造的成像能力及信噪比的成像波场分量，实现左右行波和上下行波互相关成像对陡倾角和水平界面成像能力的互补，并去除形成偏移噪声的波场分量，得到了新的逆时成像条件(式(6))及其逆时偏移结果(图 4).应用新的逆时成像条件后，采用方向行波分离的相关法逆时成像条件对应的剖面(图 4(a))，其低频噪声已得到了有效压制(对比图 2(a)和图 4(a))，陡界面刻画清晰；而在图 4(a)的基础上，采用反射角度的余弦函数作为权函数进行角度域逆时噪声衰减后的剖面(图 4(b))，偏移噪声得到了进一步压制(对比图 2(b)和图 4(b))，信噪比得到有效提高；对比图 4(a)和图 4(b)还可知，图 4(b)在保持陡倾角和水平界面成像能力的基础上，还能够进一步提高对介质②的底界面和介质⑧的顶界面的成像精度.由此可见，笔者提出新的逆时成像条件可以有效提高对复杂构造的成像精度.
 |
| 图 4 采用两种优化的逆时成像条件后的叠前逆时偏移剖面 Fig.4 Pre-stack reverse-time migration with different optimization imaging conditions |
以松辽盆地某区块深层火山岩三维地震资料为例，将经保幅预处理后的炮集数据和深度域速度模型作为逆时偏移的输入，采用文中方法与国内某专业公司的逆时偏移成像结果进行对比，最终的深度域处理结果均比例到时间域(图 5).分析图 5(a)可知，火山岩边界、内幕、断点等刻画均不清晰(椭圆位置)，且偏移划弧噪声较为严重，信噪比较低；而经文中方法逆时偏移后，火山岩机构的刻画更加清晰，同时信噪比得到有效提高，这就验证了文中方法的准确有效性.
 |
| 图 5 实际地震资料逆时偏移处理效果对比 Fig.5 Pre-stack reverse-time migration processing on practical seismic data |
提出一种新的逆时成像条件，采用波印廷矢量来实现上行波、下行波、左行波和右行波方向行波波场的分离，并实现地层界面反射角度的计算，同时对炮检点方向行波波场两两作互相关计算并成像，优选出提高成像质量的波场分量，去除形成逆时偏移噪声的分量，同时将反射角度的余弦函数作为权函数实现角度域逆时噪声的衰减，进一步提高复杂构造逆时成像质量，得到如下结论：
1)无论是上下行波还是左右行波，不同传播方向炮检点方向行波波场作互相关成像是形成偏移噪声的主要原因，需要对其进行有效压制和去除；
2)同传播方向的炮检点上下行波波场作互相关成像可以提高对水平界面和倾斜界面的刻画能力，但对近似直立界面的成像效果较差；而同传播方向的炮检点左右行波波场作互相关成像可以提高对近似直立界面和倾斜界面的刻画能力，但对水平界面的成像效果较差.于是结合同传播方向的炮检点上下行波和左右行波波场作互相关成像可以在提高成像结果信噪比的前提下，同时提高对水平界面、倾斜界面和直立界面的成像能力.
3)将反射角度的余弦函数作为权函数引入到最终的逆时成像条件中，可以进一步压制低频偏移噪声，并提高逆时成像精度.
4)模型试验和实际地震资料验证了文中方法的准确有效性，这对于实际地震资料复杂构造逆时成像具有重要的指导意义.
| [1] | YOON K, MARFURT K, STARR E. Challenges in reverse-time migration[C]//Expanded Abstracts of 74th Annual Internat SEG Mtg, 2004:1057-1060. |
| [2] | 陈可洋.高阶弹性波波动方程正演模拟及逆时偏移成像研究[D].大庆:大庆石油学院, 2009. http://www.oalib.com/references/14739568 |
| [3] | CHEN K, CHEN S, LI L, et al. Numerical experiments on the elastic wave united prestack reverse-time migration[J]. Geophysical Prospecting for Petroleum , 2014, 53 (1) :8–16. |
| [4] | 康玮, 程玖兵. 叠前逆时偏移假象去除方法[J]. 地球物理学进展 , 2012, 27 (3) :1163–1172. |
| [5] | 许璐, 孟小红, 刘国峰. 逆时偏移去噪方法研究进展[J]. 地球物理学进展 , 2012, 27 (4) :1548–1556. |
| [6] | 刘红伟, 刘洪, 邹振, 等. 地震叠前逆时偏移中的去噪与存储[J]. 地球物理学报 , 2010, 53 (9) :2171–2180. |
| [7] | 丁亮, 刘洋. 逆时偏移成像技术研究进展[J]. 地球物理学进展 , 2011, 26 (3) :1085–1100. |
| [8] | 王娟, 李振春, 陶丽. 逆时偏移成像条件研究[J]. 地球物理学进展 , 2012, 27 (3) :1173–1182. |
| [9] | 徐兴荣, 王西文, 王宇超, 等. 基于波场分离理论的逆时偏移成像条件研究及应用[J]. 地球物理学进展 , 2012, 27 (5) :2084–2090. |
| [10] | LIU F, ZHANG G, MORTON S, et al. An effective imaging condition for reverse-time migration using wave-field decomposition[J]. Geophysics , 2011, 76 (1) :S29–S39. doi:10.1190/1.3533914 |
| [11] | CHEN K, WU Q, FAN X, et al. Study on the impulse response of seismic wave prestack reverse-time migration and its application[J]. Geophysical Prospecting for Petroleum , 2013, 52 (2) :163–170. |
| [12] | COSTA J, SILVA NETO F, et al. Obliquity-correction imaging condition for reverse time migration[J]. Geophysics , 2009, 74 (3) :S57–S66. doi:10.1190/1.3110589 |
| [13] | FLETCHER R F, FOWLER P, KITCHENSIDE P. Suppressing artifacts in prestack reverse-time migration[C]//75th Annual International Meeting, SEG Expanded Abstracts, 2006:2049-2053. |
| [14] | XIE X, WU R. A depth migration method based on the full-wave reverse-time calculation and local one-way propagation[C]// 76th Annual International Meeting, SEG Expanded Abstracts, 2006:2333-2337. |
| [15] | YAN R, XIE X. A new angle-domain imaging condition for prestack reverse-time migration[C]//79th Annual International Meeting, SEG Expanded Abstracts, 2009:2784-2788. |
| [16] | YOON K, MARFURT K. Reverse-time migration using the Poynting vector[J]. Exploration Geophysics , 2006, 37 (1) :102–107. doi:10.1071/EG06102 |
| [17] | BAYSAL E, KOSLOFF D, SHERWOOD J. Reverse-time migration[J]. Geophysics , 1983, 48 (11) :1514–1524. doi:10.1190/1.1441434 |
| [18] | LOEWENTHAL D, MUFTI R. Reverse time migration in spatial frequency domain[J]. Geophysics , 1983, 48 (5) :627–635. doi:10.1190/1.1441493 |
| [19] | 陈可洋. 基于拉普拉斯算子的叠前逆时噪声压制方法[J]. 岩性油气藏 , 2011, 23 (5) :87–95. |
| [20] | 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 Conference Technical Program Expanded Abstracts, Beijing:CPS & SEG, 2009:5. |
| [21] | CHEN K, FAN X, WU Q, et al. Multi-scale decomposition and reconstruction method based on diffusion filtering and its preliminary application[J]. Chinese J Comput Phys , 2013, 30 (6) :855–861. |
| [22] | 陈可洋. 基于高阶有限差分的波动方程叠前逆时偏移方法[J]. 石油物探 , 2009, 48 (5) :475–478. |
| [23] | 陈可洋. 地震波数值模拟中差分近似的各向异性分析[J]. 石油物探 , 2010, 49 (1) :19–22. |
| [24] | 陈可洋. 完全匹配层吸收边界条件研究[J]. 石油物探 , 2010, 49 (5) :472–477. |
| [25] | 陈可洋. 边界吸收中镶边法的评价[J]. 中国科学院研究生院学报 , 2010, 27 (2) :170–175. |
| [26] | 王保利, 高静怀, 陈文超, 等. 逆时偏移中用Poynting矢量高效地提取角道集[J]. 地球物理学报 , 2013, 56 (1) :262–268. |
| [27] | 陈可洋, 吴清岭, 范兴才, 等. 地震波逆时偏移中不同域共成像点道集偏移噪声分析[J]. 岩性油气藏 , 2014, 26 (2) :118–124. |
| [28] | CHEN K, CHEN S, Li L, et al. Seismic data repair technology based on diffusion filtering method[J]. Chinese J Comput Phys , 2014, 31 (4) :465–470. |
| [29] | CHEN K, WU Q, FAN X, et al. Dual elastic wave wavefield separating simulation method and related theory derviation[J]. Chinese J Comput Phys , 2013, 30 (6) :843–853. |
| [30] | CHEN K. First-order velocity-stress elastic wave field separation numerical simulating scheme of Biot two-phase isotropic medium[J]. Chinese J Comput Phys , 2011, 28 (3) :404–412. |
| [31] | CHEN K, FAN X, WU Q, et al. Seismic wave pre-stack interpolation processing for improving the precision of reverse-time migration and its application[J]. Geophysical Prospecting for Petroleum , 2013, 52 (4) :409–416. |
引用本文 |
