CN117192609A - 一种基于Marchenko理论的克希霍夫偏移成像方法 - Google Patents
一种基于Marchenko理论的克希霍夫偏移成像方法 Download PDFInfo
- Publication number
- CN117192609A CN117192609A CN202310969035.5A CN202310969035A CN117192609A CN 117192609 A CN117192609 A CN 117192609A CN 202310969035 A CN202310969035 A CN 202310969035A CN 117192609 A CN117192609 A CN 117192609A
- Authority
- CN
- China
- Prior art keywords
- marchenko
- point
- seismic
- imaging
- wave field
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 238000003384 imaging method Methods 0.000 title claims abstract description 66
- 238000000034 method Methods 0.000 claims abstract description 55
- 238000004613 tight binding model Methods 0.000 claims abstract description 7
- 238000012546 transfer Methods 0.000 claims abstract description 3
- 238000013508 migration Methods 0.000 claims description 18
- 230000005012 migration Effects 0.000 claims description 18
- 230000001131 transforming effect Effects 0.000 claims description 7
- 238000011144 upstream manufacturing Methods 0.000 claims description 6
- 238000004088 simulation Methods 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 3
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 238000011161 development Methods 0.000 abstract description 2
- 238000005070 sampling Methods 0.000 description 6
- 238000001514 detection method Methods 0.000 description 2
- 239000011229 interlayer Substances 0.000 description 1
- 239000010410 layer Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012552 review Methods 0.000 description 1
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明属于地球物理勘探技术领域,尤其涉及一种基于Marchenko理论的克希霍夫偏移成像方法,包括:对初始速度模型采用快速行进法计算走时;根据走时估计出直达波场,结合地震记录用Marchenko方法迭代计算下行聚焦函数和上行聚焦函数,再通过下行聚焦函数和上行聚焦函数计算Marchenko格林函数;然后,将Marchenko格林函数代入克希霍夫方程计算反传波场;接着,对初始速度模型使用有限差分法计算震源波场;最后,对反传波场与震源波场应用成像条件成像。该方法能够提高克希霍夫偏移对多次波发育地区的适应性,利于复杂地质构造进行成像和解释,为油气勘探提供有力支持。
Description
技术领域
本发明属于地球物理勘探技术领域,尤其涉及一种基于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方程的时间梯度项表示为迎风差分格式来计算出网格点处的走时,采用窄带技术将计算区域内所有的地下成像点到地表的走时ti,j。
进一步地,所述步骤S2具体包括:
用地下成像点到地表的走时ti,j估计直达波场Td,并将直达波场Td的逆时波场作为初始的下行聚焦函数然后,使用Marchenko方法使用下行聚焦函数/>与地震记录R迭代求解下行聚焦函数f+和上行聚焦函数f-;最后根据下行聚焦函数f+、上行聚焦函数f-和地震记录R计算下行格林函数G+和上行格林函数G-,G+和G-直接相加得到Marchenko格林函数G。
进一步地,步骤S3具体包括:
用傅里叶变换将Marchenko格林函数G和地震记录R沿时间方向变换到频域得到频域格林函数Gω和频域地震记录Rω,再用克希霍夫方程将频域地震记录Rω反向传播到介质中,在成像点处生成反传波场
其中,*为复共轭运算,S为地表,为垂直与地表向内的单位向量。
进一步地,步骤S4、对初始速度模型使用有限差分法计算震源波场,包括:
使用有限差分法对初始速度模型做数值模拟得到地下点到地表震源的震源波场Psrc,用傅里叶变换将其变换到频域得到
进一步地,步骤S5、对反传波场与震源波场应用成像条件成像,包括:
成像条件是将反传波场和震源波场/>在频域相乘得到成像点xi处的反射系数I(xi):
计算出每个地下点的反射系数即可得到地震剖面,将全部地震记录偏移得到的地震剖面相叠加可以得到最终成像结果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方程为:
其中,为梯度算子,v(p)和t(p)分别为网格点处的速度和走时。在迎风差分法中,在网格中的一个离散的网格点p(i,j)处,Eikonal方程的时间梯度项可以用迎风差分格式为:
其中,ti,j表示点p(i,j)处的走时,vi,j表示点p(i,j)处的速度值,D表示点p(i,j)处走时t在x和z方向上的n阶向前或向后差分算子。然后,计算网格点(i,j)处的走时ti,j:
窄带技术将计算区域内所有的网格点分为完成点(走时计算已完成的节点)、窄带点(近似波前上的节点)、远离点(走时还未被计算的节点)。初始化之后,从窄带内选取走时最小的窄带点作为波前上的子震源点,并将其属性改为完成点;然后,对子震源邻近的8个网格点作判断,若为远离点则将其纳入窄带并通过迎风差分格式计算其走时,若为窄带点则采用迎风差分格式更新其走时,若为完成点则保持原有属性不变;最后,当窄带为空时计算结束,否则继续循环计算其余点的走时。
步骤S2、采用走时估计直达波场,再结合地震记录通过Marchenko方法计算Marchenko格林函数;
具体为:在步骤S1中计算出走时ti,j后,将走时作为直达波场的走时,再与震源子波做卷积得到该直达波场的振幅和相位,可以估计出直达波场Td。将直达波场Td的逆时波场作为初始的下行聚焦函数然后,使用Marchenko方法使用下行聚焦函数/>与地震记录R迭代求解下行聚焦函数f+和上行聚焦函数f-;最后根据下行聚焦函数f+、上行聚焦函数f-和地震记录R计算下行格林函数G+和上行格林函数G-,G+和G-直接相加得到Marchenko格林函数G。
步骤S3、用傅里叶变换将Marchenko格林函数G和地震记录R沿时间方向变换到频域得到频域格林函数Gω和频域地震记录Rω,再将Rω反向传播到介质中,在成像点处生成反传波
具体为:用傅里叶变换将G和R沿时间方向变换到频域得到频域格林函数Gω和频域地震记录Rω:
其中,为傅里叶变换,再用克希霍夫方程将频域地震记录Rω反向传播到介质中,在成像点处生成反传波场/>
其中,*为复共轭运算,S为地表,为垂直与地表向内的单位向量。
步骤S4、对初始速度模型使用有限差分法计算震源波场。
使用步骤S1中的方法对初始速度模型做数值模拟得到地下点到地表震源的震源波场Psrc,用傅里叶变换将其变换到频域得
步骤S5、对反传波场与震源波场应用成像条件成像得到成像点的反射系数,计算出每个地下点的反射系数得到地震剖面,将全部地震记录偏移得到的地震剖面相叠加得到最终成像结果。
成像条件是将反传波场和震源波场/>相乘得到点xi处的反射系数I(xi):
计算出每个地下点的反射系数即可得到地震剖面,将全部地震记录偏移得到的地震剖面相叠加可以得到最终成像结果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所示。对感兴趣的局部区域进行成像与全局成像相比可以成比例地缩短成像时间。
以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。
Claims (6)
1.一种基于Marchenko理论的克希霍夫偏移成像方法,其特征在于,该方法包括:
步骤S1,对初始速度模型采用快速行进法计算走时;
步骤S2,采用走时估计直达波场,再结合地震记录通过Marchenko方法计算Marchenko格林函数;
步骤S3,将Marchenko格林函数代入克希霍夫方程计算反传波场;
步骤S4,对初始速度模型使用有限差分法计算震源波场;
步骤S5,对反传波场与震源波场应用成像条件得到成像点的反射系数,计算出每个地下点的反射系数得到地震剖面,将全部地震记录偏移得到的地震剖面相叠加得到最终成像结果。
2.按照权利要求1所述的基于Marchenko理论的克希霍夫偏移成像方法,其特征在于,所述步骤S1具体包括:
在离散的网格点处,将Eikonal方程的时间梯度项表示为迎风差分格式来计算出网格点处的走时,并采用窄带技术将计算区域内所有的网格点分为完成点、窄带点、远离点。初始化之后,从窄带内选取走时最小的窄带点作为波前上的子震源点,并将其属性改为完成点;然后,对子震源邻近的8个网格点作判断,若为远离点则将其纳入窄带,并通过迎风差分格式计算走时,若为窄带点则采用迎风差分格式更新其走时,若为完成点则保持原有属性不变,当窄带为空时计算结束,否则继续循环计算其余点的走时,最终得到走时ti,j。
3.按照权利要求1所述的基于Marchenko理论的克希霍夫偏移成像方法,其特征在于,所述步骤S2具体包括:
用地下成像点到地表的走时ti,j估计直达波场Td,并将直达波场Td的逆时波场作为初始的下行聚焦函数然后,使用Marchenko方法使用下行聚焦函数/>与地震记录R迭代求解下行聚焦函数f+和上行聚焦函数f-;最后根据下行聚焦函数f+、上行聚焦函数f-和地震记录R计算下行格林函数G+和上行格林函数G-,G+和G-直接相加得到Marchenko格林函数G。
4.按照权利要求1所述的基于Marchenko理论的克希霍夫偏移成像方法,其特征在于,步骤S3具体包括:
用傅里叶变换将Marchenko格林函数G和地震记录R沿时间方向变换到频域得到频域格林函数Gω和频域地震记录Rω,再用克希霍夫方程将频域地震记录Rω反向传播到介质中得到反传波场
其中,*为复共轭运算,S为地表,为垂直与地表向内的单位向量。
5.按照权利要求1所述的基于Marchenko理论的克希霍夫偏移成像方法,其特征在于,步骤S4、对初始速度模型使用有限差分法计算震源波场,包括:
使用有限差分法对初始速度模型做数值模拟得到地下点到地表震源的震源波场Psrc,用傅里叶变换将其变换到频域得到
6.按照权利要求1所述的基于Marchenko理论的克希霍夫偏移成像方法,其特征在于,步骤S5、对反传波场与震源波场应用成像条件成像,包括:
成像条件是将反传波场和震源波场/>在频域相乘得到点xi处的反射系数I(xi):
计算出每个地下点的反射系数即可得到地震剖面,将全部地震记录偏移得到的地震剖面相叠加可以得到最终成像结果I。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310969035.5A CN117192609A (zh) | 2023-08-03 | 2023-08-03 | 一种基于Marchenko理论的克希霍夫偏移成像方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202310969035.5A CN117192609A (zh) | 2023-08-03 | 2023-08-03 | 一种基于Marchenko理论的克希霍夫偏移成像方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN117192609A true CN117192609A (zh) | 2023-12-08 |
Family
ID=88989532
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202310969035.5A Pending CN117192609A (zh) | 2023-08-03 | 2023-08-03 | 一种基于Marchenko理论的克希霍夫偏移成像方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN117192609A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN118033746A (zh) * | 2024-04-07 | 2024-05-14 | 吉林大学 | 一种面向目标的海洋拖缆地震数据Marchenko成像方法 |
-
2023
- 2023-08-03 CN CN202310969035.5A patent/CN117192609A/zh active Pending
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN118033746A (zh) * | 2024-04-07 | 2024-05-14 | 吉林大学 | 一种面向目标的海洋拖缆地震数据Marchenko成像方法 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US8437998B2 (en) | Hybrid method for full waveform inversion using simultaneous and sequential source method | |
US9625593B2 (en) | Seismic data processing | |
Verschuur et al. | Removal of internal multiples with the common-focus-point (CFP) approach: Part 2—Application strategies and data examples | |
US20120275267A1 (en) | Seismic Data Processing | |
CN111158049B (zh) | 一种基于散射积分法的地震逆时偏移成像方法 | |
Hindriks et al. | Reconstruction of 3-D seismic signals irregularly sampled along two spatial coordinates | |
GB2547942A (en) | Method for deghosting and redatuming operator estimation | |
CN111856577B (zh) | 一种降低逆时偏移地表炮检距道集计算量的方法 | |
CN111257939B (zh) | 一种时移地震虚拟震源双向波场重构方法和系统 | |
Song et al. | Efficient wavefield inversion with outer iterations and total variation constraint | |
CN117192609A (zh) | 一种基于Marchenko理论的克希霍夫偏移成像方法 | |
CA2808083C (en) | Hybrid method for full waveform inversion using simultaneous and sequential source method | |
Ovcharenko et al. | Extrapolating low-frequency prestack land data with deep learning | |
Yang et al. | Mini-batch optimized full waveform inversion with geological constrained gradient filtering | |
EA038811B1 (ru) | Способ и система для формирования геофизических данных | |
CN113687417B (zh) | 一种三维叠前地震数据层间多次波预测和压制方法 | |
CN110967734B (zh) | 基于快速傅立叶变换的虚源重构方法及系统 | |
Jia et al. | Subbasalt Marchenko imaging with offshore Brazil field data | |
US8010293B1 (en) | Localized seismic imaging using diplets | |
CN113866821B (zh) | 一种基于照明方向约束的被动源干涉偏移成像方法和系统 | |
Biondi | Prestack exploding-reflectors modeling for migration velocity analysis | |
Staring et al. | R-EPSI and Marchenko equation-based workflow for multiple suppression in the case of a shallow water layer and a complex overburden: A 2D case study in the Arabian Gulf | |
Chen et al. | Poststack internal multiples attenuation based on virtual events | |
Kim et al. | Least-squares reverse time migration using analytic-signal-based wavefield decomposition | |
CN118033746B (zh) | 一种面向目标的海洋拖缆地震数据Marchenko成像方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination |