一种基于Marchenko理论的克希霍夫偏移成像方法

技术领域

本发明属于地球物理勘探技术领域，尤其涉及一种基于Marchenko理论的克希霍夫偏移成像方法。

背景技术

克希霍夫偏移是地球物理学中广泛使用的地下结构成像方法，相关的文献包括：Fehler M C,Huang L.Modern imaging using seismic reflection data[J].Annualreview of earth and planetary sciences,2002,30(1):259-284.；Keho T H,Beydoun WB.Paraxial ray Kirchhoff migration[J].Geophysics,1988,53(12):1540-1546.；Thorbecke J,Slob E,Brackenhoff J,et al.Implementation of the Marchenko method[J].Geophysics,2017,82(6):WB29-WB45.；Wapenaar K,Thorbecke J,Van Der Neut J,etal.Marchenko imaging[J].Geophysics,2014,79(3):WA39-WA57，现有方法常使用射线理论来计算格林函数实现快速成像(Keho TH and Beydoun W B,1988)。与基于波动方程的方法相比，克希霍夫偏移计算速度快、易于解释成像剖面、可以实现单个共炮点道集的偏移等。但是，该方法通常是一个基于射线追踪的高频近似方法，不能完全解释复杂结构中的波动现象，并且忽略了地震波在地下的多次散射，不能适应易发育层间多次波的地下结构成像(Fehler M C and Huang L,2002)。采用克希霍夫偏移直接对包含有多次波的地震数据进行成像时，层间多次波会在地震剖面中产生假象。Marchenko成像是一种数据驱动的成像方法，能估计出包含有准确的多次波信息的格林函数，从而避免在成像剖面中出现多次波引起的假象(Wapenaar K.etal.,2014；Thorbecke J.et al.,2017)。

发明内容

本发明所要解决的技术问题在于提供一种基于Marchenko理论的克希霍夫偏移成像方法，来解决克希霍夫偏移对层间多次波适应性差的问题。

本发明是这样实现的，

一种基于Marchenko理论的克希霍夫偏移成像方法，该方法包括：

步骤S1，对初始速度模型采用快速行进法计算走时。

步骤S2，采用走时估计直达波场，再结合地震记录通过Marchenko方法计算Marchenko格林函数；

步骤S3，将Marchenko格林函数代入克希霍夫方程计算反传波场；

步骤S4，对初始速度模型使用有限差分法计算震源波场；

步骤S5，对反传波场与震源波场应用成像条件成像得到成像点的反射系数，计算出每个地下点的反射系数得到地震剖面，将全部地震记录偏移得到的地震剖面相叠加得到最终成像结果。

进一步地，所述步骤S1具体包括：

在网格中的一个离散的网格点处，Eikonal方程的时间梯度项表示为迎风差分格式来计算出网格点处的走时，采用窄带技术将计算区域内所有的地下成像点到地表的走时t

进一步地，所述步骤S2具体包括：

用地下成像点到地表的走时t

进一步地，步骤S3具体包括：

用傅里叶变换将Marchenko格林函数G和地震记录R沿时间方向变换到频域得到频域格林函数G

其中，*为复共轭运算，S为地表，

进一步地，步骤S4、对初始速度模型使用有限差分法计算震源波场，包括：

使用有限差分法对初始速度模型做数值模拟得到地下点到地表震源的震源波场P

进一步地，步骤S5、对反传波场与震源波场应用成像条件成像，包括：

成像条件是将反传波场

计算出每个地下点的反射系数即可得到地震剖面，将全部地震记录偏移得到的地震剖面相叠加可以得到最终成像结果I。

本发明与现有技术相比，有益效果在于：

本发明将Marchenko理论引入克希霍夫偏移中，为克希霍夫偏移领域带来了重大突破，能使克希霍夫偏移适应多次波发育区域的地下结构成像。通过该方法，能够更好地理解地下结构，利于复杂地质构造的成像和解释。

附图说明

图1本发明实施例提供的方法流程图；

图2本发明实施例提供的真实速度模型，其中，三角形表示炮点，三角形表示检波点。

图3本发明实施例提供的观测系统记录的地震记录，其中，斜下箭头指示一次波，斜上箭头指示多次波，第1炮、121炮和241炮的地震记录分别3(a)、3(b)、3(c)；

图4本发明实施例提供的初始速度模型；

图5本发明实施例提供的快速行进法计算的走时，网格点(181,1)、(181,121)和(181,241)的走时分别为5(a)、5(b)和5(c)；

图6本发明实施例提供的Marchenko格林函数，网格点(181,1)、(181,121)和(181,241)的Marchenko格林函数为6(a)、6(b)和6(c)；

图7本发明实施例提供的克希霍夫方法计算的反传波场，网格点(181,1)、(181,121)和(181,241)处的反传波场分别7(a)、7(b)和7(c)；

图8本发明实施例提供的部分地震记录偏移得到的全局地震剖面，网格点(181,1)、(181,121)和(181,241)到地表检波点的震源波场分别8(a)、8(b)和8(c)；

图9本发明实施例提供的叠加后的全局地震剖面；

图10本发明实施例提供的部分地震记录偏移得到的局部地震剖面；

图11本发明实施例提供的叠加后的局部地震剖面。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

一种克希霍夫偏移成像模拟方法，如图1所示。在对真实速度模型正演得到地震记录后，该方法包括如下步骤：

步骤S1，对初始速度模型采用快速行进法计算走时；

步骤S2，采用走时估计直达波场，再结合地震记录通过Marchenko方法计算Marchenko格林函数；

步骤S3，将Marchenko格林函数代入克希霍夫方程计算反传波场；

步骤S4，对初始速度模型使用有限差分法计算震源波场；

步骤S5，对反传波场与震源波场应用成像条件成像得到成像点处的反射系数，计算出每个地下点的反射系数得到地震剖面，将全部地震记录偏移得到的地震剖面相叠加得到最终成像结果。

其中，步骤S1、对初始速度模型采用快速行进法计算走时；

具体为：快速行进法计算的Eikonal方程为：

其中，

其中，t

窄带技术将计算区域内所有的网格点分为完成点(走时计算已完成的节点)、窄带点(近似波前上的节点)、远离点(走时还未被计算的节点)。初始化之后，从窄带内选取走时最小的窄带点作为波前上的子震源点，并将其属性改为完成点；然后，对子震源邻近的8个网格点作判断，若为远离点则将其纳入窄带并通过迎风差分格式计算其走时，若为窄带点则采用迎风差分格式更新其走时，若为完成点则保持原有属性不变；最后，当窄带为空时计算结束，否则继续循环计算其余点的走时。

步骤S2、采用走时估计直达波场，再结合地震记录通过Marchenko方法计算Marchenko格林函数；

具体为：在步骤S1中计算出走时t

步骤S3、用傅里叶变换将Marchenko格林函数G和地震记录R沿时间方向变换到频域得到频域格林函数G

具体为：用傅里叶变换将G和R沿时间方向变换到频域得到频域格林函数G

其中，

其中，*为复共轭运算，S为地表，

步骤S4、对初始速度模型使用有限差分法计算震源波场。

使用步骤S1中的方法对初始速度模型做数值模拟得到地下点到地表震源的震源波场P

步骤S5、对反传波场与震源波场应用成像条件成像得到成像点的反射系数，计算出每个地下点的反射系数得到地震剖面，将全部地震记录偏移得到的地震剖面相叠加得到最终成像结果。

成像条件是将反传波场

计算出每个地下点的反射系数即可得到地震剖面，将全部地震记录偏移得到的地震剖面相叠加可以得到最终成像结果I。

应用实施例：

在本实施例中，真实速度模型如图2所示。模型网格点数为201*241，网格间距为8m*8m，模型大小为1600m*1920m。数值模拟时，震源子波的主频为20Hz，观测系统为均匀布设在地表的241个炮点和检波点，炮间距和道间距均为8m、采样时间间隔为0.002s、采样时间3s、采样点数为1501。其中，第1炮、121炮和241炮的地震记录分别如图3(a)、3(b)、3(c)所示。

步骤S1、对初始速度模型采用快速行进法计算走时。

本实施例中，采用如图4所示的初始速度模型作为初始速度模型。初始速度模型的网格点数为201*241，网格间距为8m*8m，模型大小为1600m*1920m。

针对初始速度模型(图4)采用快速行进法计算每个网格点到地表的走时。其中，网格点(181，1)、(181，121)和(181，241)的走时分别如图5(a)、5(b)和5(c)所示。

步骤S2、采用走时估计直达波场，再结合地震记录通过Marchenko方法计算Marchenko格林函数；

根据地震记录和步骤S1中计算的走时，用Marchenko方法计算每个网格点处的Marchenko格林函数。其中，网格点(181,1)、(181,121)和(181,241)的Marchenko格林函数如图6(a)、6(b)和6(c)所示，可以看出，Marchenko方法能在Marchenko格林函数中将一次波和多次波刻画出来。

步骤S3、将Marchenko格林函数代入克希霍夫方程计算反传波场。

在步骤S2中计算出模型中每个网格点的Marchenko格林函数后，将地震记录和Marchenko格林函数代入公式(6)计算出每一炮到地下每个点的反传波场。其中，第121炮(图3b)在网格点(181,1)、(181,121)和(181,241)处的反传波场分别如图7(a)、7(b)和7(c)所示。

步骤S4、对初始速度模型使用有限差分法计算震源波场。

对初始速度模型(图4)采用有限差分法计算震源波场时，震源子波的主频为20Hz，观测系统的炮点布设在网格点处，241个检波点以8m道间距均匀地布设在地表，采样时间间隔为0.002s、采样时间为3s、采样点数为1501。使炮点遍历所有网格点，计算每个网格点到地表检波点的震源波场。其中，网格点(181,1)、(181,121)和(181,241)到地表检波点的震源波场分别如图8(a)、8(b)和8(c)所示。

步骤S5、对反传波场与震源波场应用成像条件成像得到成像点的反射系数，计算出每个地下点的反射系数得到地震剖面，将全部地震记录偏移得到的地震剖面相叠加得到最终成像结果。

在每个网格点处，对步骤S4计算的反传波场和步骤S5计算的震源波场应用成像条件，可以将一个炮集偏移为1600m*1920m大小的全局地震剖面，部分炮集偏移的全局地震剖面如图9所示。将241个炮集的全局地震剖面直接叠加就可以得到最终的全局地震剖面，如图10所示。

在步骤S2到步骤S5中，如果将一个炮集偏移只针对深度为480m～1440m、距离为240m～960m内的局部区域进行成像，可以得到该区域的局部地震剖面，部分炮集偏移的局部地震剖面如图11所示。将241个炮集的局部地震剖面进行叠加，可以得到最终的局部地震剖面，如图12所示。对感兴趣的局部区域进行成像与全局成像相比可以成比例地缩短成像时间。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。