CN114325829A - 一种基于双差思想的全波形反演方法 - Google Patents
一种基于双差思想的全波形反演方法 Download PDFInfo
- Publication number
- CN114325829A CN114325829A CN202111574138.9A CN202111574138A CN114325829A CN 114325829 A CN114325829 A CN 114325829A CN 202111574138 A CN202111574138 A CN 202111574138A CN 114325829 A CN114325829 A CN 114325829A
- Authority
- CN
- China
- Prior art keywords
- model
- waveform inversion
- double
- difference
- seismic
- 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.)
- Granted
Links
Images
Landscapes
- Geophysics And Detection Of Objects (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明涉及一种基于双差思想的全波形反演方法,该方法包括以下步骤:1)对原始地震数据进行预处理;2)建立初始速度模型:3)利用双差波形反演方法进行速度迭代更新,获得最终的速度模型。与现有技术相比,本发明能够在光滑模型下减弱地震子波误差对波形反演的影响,获得较高精度与分辨率的速度模型。
Description
技术领域
本发明涉及勘探地震学中的速度建模领域,尤其是涉及一种基于双差思想的全波形反演方法。
背景技术
近地表的地下速度模型对于地震数据处理(静校正、偏移成像等)具有关键影响,然而由于地表、地下介质复杂,地下速度建模是一个极具挑战的问题,地球物理学发展了一系列基于反演的地下速度模型构建方法,包含了走时层析与波形反演方法,传统走时层析方法对地震波进行高频近似,利用走时信息,可高效的提供地下速度模型,但是分辨率低,在复杂速度环境下应用效果较差,随着勘探技术的提高,逐渐从利用走时信息走向利用波形信息来刻画地下速度结构,上世纪80年代,Tarantola等率先提出全波形反演理论框架(Full Waveform Inversion,简称FWI)(Tarantola等,1984),根据全波形反演的定义(PannG S等,1986),利用合成数据匹配观测数据指导更新速度模型,建立L2范数下观测地震记录与模拟地震记录误差的目标函数,利用优化算法降低目标函数来更新初始模型,因其利用波形信息进行反演,地震运动学与动力学保真,具有较高的精度与分辨率,成为近年来研究的热点,也广泛应用于各种复杂构造的速度建模。
然而全波形反演也存在一定局限性,一方面,它依赖于初始模型的建立(GauthierO等,1986;Mora P,1987),另一方面如果引入系统误差,如地震子波选取不当,也会影响全波形反演的结果。
Lee,Kim等提出了不依赖子波的全波形反演方法,然而计算量大,需计算Jacobian矩阵(Lee K H等,2003)。Choi等在频率域提出利用褶积法和反褶积法,先对数据与参考道振幅相除再建立目标函数来消除子波影响(Choi Y等,2005)。随后也成功应用于时间域全波形反演(Choi Y等,2011)。敖瑞德等提出基于包络的不依赖子波方法,通过包络对数的目标函数建立全波形反演(敖瑞德等,2015)。除了直接消除子波影响,另一种思路在于尝试估计地震子波进行反演。刘立彬等介绍了一种针对时间域波形反演的子波估计方法,对一给定的试探地震子波,通过频率域对比观测数据来反演预测子波,可较好的吻合真实子波,在全波形反演中取得了不错的效果(刘立彬,2020)。但是与Pratt等的实验相似,需要观测数据良好的信噪比并建立在较准确的初始模型上。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于双差思想的全波形反演方法。
本发明的目的可以通过以下技术方案来实现:
一种基于双差思想的全波形反演方法,该方法包括以下步骤:
1)对原始地震数据进行预处理;
2)建立初始速度模型:
3)利用双差波形反演方法进行速度迭代更新,获得最终的速度模型。
所述的步骤1)中,预处理包括去噪、道均衡和去多次波。
所述的步骤2)中,根据先验信息,采用层析成像建立初始速度模型。
所述的步骤3)具体包括以下步骤:
31)基于二维常密度声波方程,进行数值正演模拟计算,获得模拟地震记录dcal;
32)根据真实地震记录dobs与模拟地震记录dcal的残差,计算每一炮的核函数,即伴随震源;
33)基于二维常密度声波方程,利用伴随震源进行波场反传,得到反传波场;
34)利用伴随状态法,计算基于双差波形反演目标函数的梯度;
35)更新速度模型;
36)开始下一次迭代,重复步骤31-35),直到获得最终的速度反演结果。
所述的步骤31)中,二维常密度声波方程的表达式为:
其中,u(x,z,t)表示声波波场,v(x,z)为速度场,表示质点在(x,z)位置的纵波速度,s(t)表示震源函数,(xs,zs)表示震源所处位置,x表示横向位置,z表示深度位置,t表示时间,δ(·)为0/1函数,当自变量取值0时,其值为1,当自变量取值不为0时,则其值为0。
所述的步骤34)中,目标函数E的表达式为:
所述的目标函数E的梯度具体为:
所述的步骤35)中,采用LBFGS法更新速度模型。
所述的步骤35)具体包括以下步骤:
351)构造Hessian逆矩阵的正定近似矩阵Hp+1,满足关系式:
Hp+1yp=sp
352)利用有限数目的梯度残差和模型残差信息近似构造矩阵Hp:
353)利用拟Newton法的迭代公式,进行速度模型更新,则有:
与现有技术相比,本发明具有以下优点:
本发明基于双差思想,能够实现在光滑背景模型下减弱地震子波的误差对波形反演结果的影响,与全波形反演中子波误差项不同在于双差波形反演的子波误差项进行了一次道间作差,使得反演结果对比全波形反演更准确和稳定,获得较高精度与分辨率的速度模型。
附图说明
图1为本发明的方法流程图。
图2为实施例1中的真实速度模型。
图3为实施例1中的初始速度模型。
图4为实施例1中本发明提出的双差波形反演结果。
图5为实施例1中的全波形反演结果。
图6为实施例2中的真实速度模型。
图7为实施例2中的初始速度模型。
图8为实施例2中本发明提出的双差波形反演结果。
图9为实施例2中的全波形反演结果。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
如图1所示,本发明提出一种基于双差思想的全波形反演方法,该方法包括以下步骤:
1)对原始地震数据进行预处理,包括去噪、道均衡和去多次波;
2)建立初始速度模型;
3)利用双差波形反演方法进行速度迭代更新,获得最终的速度模型,具体包括以下步骤:
31)基于二维常密度声波方程,进行数值正演模拟计算,获得模拟地震记录u(x,t),二维常密度声波方程为:
其中,u(x,z,t)表示声波波场;v(x,z)为速度场,表示质点在(x,z)位置的纵波速度;s(t)表示震源函数,(xs,zs)表示震源所处位置;
32)利用真实地震记录dobs与模拟地震记录dcal,计算每一炮的核函数,即伴随震源;
伴随震源的计算式,
33)基于二维常密度声波方程,利用伴随震源进行波场反传;
34)利用伴随状态法,计算基于双差波形反演目标函数的梯度,目标函数E为:
目标函数的梯度表达式为:
其中u代表模型下地震的正传波场,λ通过计算得知是以(该道的模拟地震记录和观测地震记录的残差与下一道的模拟地震记录和观测地震记录的残差的差作为第该道与下一道的伴随震源)作为震源按声波方程进行反传的波场;
35)利用LBFGS法更新速度模型,具体包括以下步骤:
351)拟Newton法构造Hessian逆矩阵的正定近似矩阵Hp+1,满足关系式:
Hp+1yp=sp
352)利用有限数目的梯度残差和模型残差信息近似构造矩阵Hp:
其中,I为单位矩阵,为对Hessian逆矩阵的初始近似矩阵,单位矩阵I是一种最简单的初始矩阵选取方式,r=3~20(本例中选取8),当迭代次数p大于r时,在计算出该次迭代的模型后,采用最新的一对向量{sp,yp}替换掉中最旧的一对;
353)利用拟Newton法的迭代公式,进行速度模型更新,则有:
36)开始下一次迭代,重复步骤31-35),直到获得最终的速度反演结果。
本发明基于双差思想实现在光滑背景模型下减弱地震子波的误差对波形反演结果的影响,得到较高精度与分辨率的速度模型,具体原理为:
定义地震记录为一个算子(格林函数)作用域地震子波:
u=Gf
其值,u为地震记录,G代表格林函数,f代表震源函数。
假设用于波形反演的格林函数与地震子波为:
G′=G0+ΔG,f′=f0+Δf
式中,G0为真实速度模型的格林函数,f0为真实(观测记录)地震子波,G′为波形反演中初始速度模型的格林函数,f′为用于波形反演的地震子波,ΔG为初始模型与真实速度模型格林函数差,Δf为始模型与真实速度模型地震子波差。
在全波形反演中,目标函数(忽略高阶小项)可写作:
其中,ΔGf0为由于速度模型差异导致的误差项,在全波形反演中指导更新速度模型,当速度模型等于真实速度模型,ΔG=0,目标函数为零,但是,当存在子波误差项时,引入了G0Δf,即当初始模型等于真实模型的情况下,目标函数仍不为零,会错误的更正模型,同时子波误差项的存在也让全波形反演的速度模型更新偏离实际,造成反演结果出错,不稳定。
由此,在本发明的双差波形反演中,第k、k+1地震道的目标函数(忽略高阶小项)可写作:
式中,Gk、Gk+1表示k、k+1道的格林函数,其中(ΔGk+1-ΔGk)f0为由于速度模型差异导致的误差项,指导更新速度模型,为子波误差项,与现有的全波形反演中子波误差项不同在于,本发明的双差波形反演的子波误差项进行了一次道间作差,一般模型下,由于,该项不等于零,会遇到同全波形反演一样的问题,但是当真实速度模型较简单平滑,子波误差项接近于零,在这种情况下,双差波形反演能一定程度上克服子波误差项对反演带来的影响,主要会根据速度误差项(ΔGk+1-ΔGk)f0更新速度模型,使得反演结果对比全波形反演更准确和稳定。
实施例1:
本实施例以二维Marmousi理论模型作为真实模型(如图2所示),该模型共有298×160个网格,网格间距为20米×20米,速度最大最小值分别为1500米/秒和5500米/秒,地下200米以上为1500米/秒的均匀常速水层,在该模型上进行声波正演模拟,共正演42炮。炮点均匀分布于地表,炮间距为140米,第一炮在水平位置40米处,检波点以20米间隔均匀分布于地表,初始速度模型为光滑后速度模型(如图3所示),震源函数采用主频10Hz,能量振幅为5的雷克子波,观测时间共4s,时间采样间隔2ms,在应用本发明提出的双差波形反演方法(DDWI)的同时,基于伴随状态法的全波形反演方法(FWI)也同时被应用。
利用DDWI与FWI进行波形反演的结果分别如图4和5所示,双差波形反演与全波形反演均能有效进行速度反演,甚至前者深部的反演分辨率更高,可见Marmousi模型的典型特征和细节特征,尤其是中浅部,可见,对于复杂的Marmousi模型,双差波形反演与全波形反演结果相仿,都能有效准确地进行速度反演。
实施例2
本实施例以简单高速异常球体模型作为真实速度模型(如图6所示),该模型共有101×101个网格,网格间距为5米×5米,背景速度为3000米/秒,在中心(250米,250米)位置设置一半径50米的异常球体,速度为3300米/秒,地下25米以上、地下475米以下为3000米/秒的均匀常速层,在该模型上进行声波正演模拟,共正演25炮,炮点均匀分布于地表,炮间距为20米,第一炮在水平位置10米处,检波点以5米间隔均匀分布于底部,初始速度模型为3000米/秒背景速度,无高速异常球体(如图7所示),震源函数采用主频30Hz,能量振幅为5的雷克子波,在模拟地震记录中采用主频50Hz,能量振幅为1的雷克子波,观测时间共0.3s,时间采样间隔0.5ms,在应用本发明提出的双差波形反演方法(DDWI)的同时,基于伴随状态法的全波形反演方法(FWI)也同时被应用。
利用DDWI与FWI进行波形反演的结果分别如图8和9所示,全波形反演结果效果很差,反演结果不稳定,且不可见中心异常球体特征,偏离实际,相比之下,虽然双差波形反演结果也些许变差,中心异常球体轮廓尚存但是形状变大。但是总体上依旧符合真实速度模型特征。
因而双差波形反演在光滑背景模型下确实比全波形反演对子波误差不敏感,反演表现比传统全波形反演要好,即双差波形反演在子波误差存在的情况下,能一定程度上克服子波误差的影响。
Claims (10)
1.一种基于双差思想的全波形反演方法,其特征在于,该方法包括以下步骤:
1)对原始地震数据进行预处理;
2)建立初始速度模型:
3)利用双差波形反演方法进行速度迭代更新,获得最终的速度模型。
2.根据权利要求1所述的一种基于双差思想的全波形反演方法,其特征在于,所述的步骤1)中,预处理包括去噪、道均衡和去多次波。
3.根据权利要求1所述的一种基于双差思想的全波形反演方法,其特征在于,所述的步骤2)中,根据先验信息,采用层析成像建立初始速度模型。
4.根据权利要求1所述的一种基于双差思想的全波形反演方法,其特征在于,所述的步骤3)具体包括以下步骤:
31)基于二维常密度声波方程,进行数值正演模拟计算,获得模拟地震记录dcal;
32)根据真实地震记录dobs与模拟地震记录dcal的残差,计算每一炮的核函数,即伴随震源;
33)基于二维常密度声波方程,利用伴随震源进行波场反传,得到反传波场;
34)利用伴随状态法,计算基于双差波形反演目标函数的梯度;
35)更新速度模型;
36)开始下一次迭代,重复步骤31-35),直到获得最终的速度反演结果。
9.根据权利要求4所述的一种基于双差思想的全波形反演方法,其特征在于,所述的步骤35)中,采用LBFGS法更新速度模型。
10.根据权利要求9所述的一种基于双差思想的全波形反演方法,其特征在于,所述的步骤35)具体包括以下步骤:
351)构造Hessian逆矩阵的正定近似矩阵Hp+1,满足关系式:
Hp+1yp=sp
352)利用有限数目的梯度残差和模型残差信息近似构造矩阵Hp:
353)利用拟Newton法的迭代公式,进行速度模型更新,则有:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111574138.9A CN114325829B (zh) | 2021-12-21 | 2021-12-21 | 一种基于双差思想的全波形反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111574138.9A CN114325829B (zh) | 2021-12-21 | 2021-12-21 | 一种基于双差思想的全波形反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114325829A true CN114325829A (zh) | 2022-04-12 |
CN114325829B CN114325829B (zh) | 2023-03-28 |
Family
ID=81055212
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111574138.9A Active CN114325829B (zh) | 2021-12-21 | 2021-12-21 | 一种基于双差思想的全波形反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114325829B (zh) |
Citations (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130138408A1 (en) * | 2011-11-29 | 2013-05-30 | Sunwoong Lee | Methods for Approximating Hessian Times Vector Operation in Full Wavefield Inversion |
US20170131418A1 (en) * | 2015-11-05 | 2017-05-11 | Cgg Services Sa | Device and method for full waveform inversion |
CN107422379A (zh) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | 基于局部自适应凸化方法的多尺度地震全波形反演方法 |
CN107505654A (zh) * | 2017-06-23 | 2017-12-22 | 中国海洋大学 | 基于地震记录积分的全波形反演方法 |
CN107664771A (zh) * | 2017-09-28 | 2018-02-06 | 西南石油大学 | 一种基于相似性系数的微地震全波形定位方法 |
CN107765302A (zh) * | 2017-10-20 | 2018-03-06 | 吉林大学 | 不依赖震源子波的时间域单频波形走时反演方法 |
US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
CN110058307A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于快速拟牛顿法的全波形反演方法 |
US20210080573A1 (en) * | 2017-09-01 | 2021-03-18 | The Trustees Of Princeton University | Quantitative Ultrasound Imaging Based on Seismic Full Waveform Inversion |
CN113447981A (zh) * | 2021-06-18 | 2021-09-28 | 同济大学 | 一种基于共成像点道集的反射全波形反演方法 |
WO2021252693A1 (en) * | 2020-06-11 | 2021-12-16 | Downunder Geosolutions (America) Llc | Seismic wavefield modeling honoring avo/ava with applications to full waveform inversion and least-squares imaging |
-
2021
- 2021-12-21 CN CN202111574138.9A patent/CN114325829B/zh active Active
Patent Citations (11)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20130138408A1 (en) * | 2011-11-29 | 2013-05-30 | Sunwoong Lee | Methods for Approximating Hessian Times Vector Operation in Full Wavefield Inversion |
US20170131418A1 (en) * | 2015-11-05 | 2017-05-11 | Cgg Services Sa | Device and method for full waveform inversion |
US20180356548A1 (en) * | 2017-06-12 | 2018-12-13 | Institute Of Geology And Geophysics Chinese Academy Of Sciences | Inversion velocity model, method for establishing the same and method for acquiring images of underground structure |
CN107505654A (zh) * | 2017-06-23 | 2017-12-22 | 中国海洋大学 | 基于地震记录积分的全波形反演方法 |
CN107422379A (zh) * | 2017-07-27 | 2017-12-01 | 中国海洋石油总公司 | 基于局部自适应凸化方法的多尺度地震全波形反演方法 |
US20210080573A1 (en) * | 2017-09-01 | 2021-03-18 | The Trustees Of Princeton University | Quantitative Ultrasound Imaging Based on Seismic Full Waveform Inversion |
CN107664771A (zh) * | 2017-09-28 | 2018-02-06 | 西南石油大学 | 一种基于相似性系数的微地震全波形定位方法 |
CN107765302A (zh) * | 2017-10-20 | 2018-03-06 | 吉林大学 | 不依赖震源子波的时间域单频波形走时反演方法 |
CN110058307A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于快速拟牛顿法的全波形反演方法 |
WO2021252693A1 (en) * | 2020-06-11 | 2021-12-16 | Downunder Geosolutions (America) Llc | Seismic wavefield modeling honoring avo/ava with applications to full waveform inversion and least-squares imaging |
CN113447981A (zh) * | 2021-06-18 | 2021-09-28 | 同济大学 | 一种基于共成像点道集的反射全波形反演方法 |
Non-Patent Citations (5)
Title |
---|
S.M.A. SHOJA,等: "A Comparison between Time-Domain and Frequency-Domain Full Waveform Inversion", 《EAGE》 * |
慎国强等: "基于声波波动方程的全波形反演方法及应用", 《油气地球物理》 * |
肖卓等: "利用双差成像方法反演青藏高原东北缘及其邻区地壳速度结构", 《地球物理学报》 * |
胡勇等: "基于精确震源函数的解调包络多尺度全波形反演", 《地球物理学报》 * |
赵磊,等: "波动方程初至双差走时层析反演", 《石油物探》 * |
Also Published As
Publication number | Publication date |
---|---|
CN114325829B (zh) | 2023-03-28 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107783190B (zh) | 一种最小二乘逆时偏移梯度更新方法 | |
CN106526674B (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
CN109669212B (zh) | 地震数据处理方法、地层品质因子估算方法与装置 | |
CN111158049B (zh) | 一种基于散射积分法的地震逆时偏移成像方法 | |
CN109917454B (zh) | 基于双基准面的真地表叠前深度偏移成像方法及装置 | |
CN110018517B (zh) | 一种多尺度地面微地震逆时干涉定位方法 | |
CN109946741B (zh) | 一种TTI介质中纯qP波最小二乘逆时偏移成像方法 | |
CN110187382B (zh) | 一种回折波和反射波波动方程旅行时反演方法 | |
CN111045077B (zh) | 一种陆地地震数据的全波形反演方法 | |
CN113064203B (zh) | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 | |
CN111638551A (zh) | 地震初至波走时层析方法及装置 | |
CN111665556B (zh) | 地层声波传播速度模型构建方法 | |
CN108680968B (zh) | 复杂构造区地震勘探数据采集观测系统评价方法及装置 | |
CN106353798A (zh) | 多分量联合高斯束叠前逆时偏移成像方法 | |
CN110780341B (zh) | 一种各向异性地震成像方法 | |
CN112130199A (zh) | 一种优化的最小二乘逆时偏移成像方法 | |
CN114325829B (zh) | 一种基于双差思想的全波形反演方法 | |
CN114002742B (zh) | 一种欧拉高斯束偏移成像方法及装置 | |
CN114624766B (zh) | 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法 | |
CN111175822B (zh) | 改进直接包络反演与扰动分解的强散射介质反演方法 | |
CN110888158B (zh) | 一种基于rtm约束的全波形反演方法 | |
CN111665550A (zh) | 地下介质密度信息反演方法 | |
CN111665546A (zh) | 用于可燃冰探测的声学参数获取方法 | |
CN111665549A (zh) | 地层声波衰减因子反演方法 | |
CN114924313B (zh) | 基于行波分离的声波最小二乘逆时偏移梯度求取方法 |
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 | ||
GR01 | Patent grant | ||
GR01 | Patent grant |