CN110857999B - 一种基于全波形反演的高精度波阻抗反演方法及系统 - Google Patents

一种基于全波形反演的高精度波阻抗反演方法及系统 Download PDF

Info

Publication number
CN110857999B
CN110857999B CN201810970557.6A CN201810970557A CN110857999B CN 110857999 B CN110857999 B CN 110857999B CN 201810970557 A CN201810970557 A CN 201810970557A CN 110857999 B CN110857999 B CN 110857999B
Authority
CN
China
Prior art keywords
impedance
seismic
unit
wave
wave impedance
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.)
Active
Application number
CN201810970557.6A
Other languages
English (en)
Other versions
CN110857999A (zh
Inventor
何兵红
方伍宝
胡光辉
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Geophysical Research Institute
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China Petroleum and Chemical Corp, Sinopec Geophysical Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201810970557.6A priority Critical patent/CN110857999B/zh
Publication of CN110857999A publication Critical patent/CN110857999A/zh
Application granted granted Critical
Publication of CN110857999B publication Critical patent/CN110857999B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种基于全波形反演的高精度波阻抗反演方法及系统,属于地震资料解释领域。该基于全波形反演的高精度波阻抗反演方法采用基于CPML边界条件的速度‑阻抗方程进行地震波数值正演得到地震数值模拟记录,利用所述地震数值模拟记录与实际观测数据的误差通过全波形反演获得高精度波阻抗。本发明方法通过全波形反演得到的阻抗不仅在井位处获得了准确的阻抗值,也反演得到了真实模型中的构造形态,且断层断块边界清晰。采用全波形反演获得的阻抗模型进行地震资料解释工作将会大幅提高油气预测的准确性和可靠性。

Description

一种基于全波形反演的高精度波阻抗反演方法及系统
技术领域
本发明属于地震资料解释领域,具体涉及一种基于全波形反演的高精度波阻抗反演方法及系统。
背景技术
利用地震资料进行波阻抗反演是地震油气勘探中常用的地震特殊处理解释技术。波阻抗具有明确的物理意义,是储层预测、油藏特征描述的确定性方法,在实际资料应用中取得了显著的效果。现有的地震波阻抗反演方法也存在一定的弊端。地震波的传播在空间上是三维的,现有的地震波阻抗反演方法则是基于射线理论利用褶积模型在一维深度空间进行,没有考虑地震波传播在三维空间的综合效应,缺失了三维空间来自不同方位的地震波的影响和作用。此外,褶积模型是基于地球介质的水平层状假设,忽略了实际地球介质陡构造发育、地层倾角大等问题,是一种理想状态且严重背离了真实地下构造分布,波阻抗反演可靠性程度低。同时在射线理论假设下,只能利用地震数据的运动学信息对大尺度构造进行反演,无法建立小于地震波长的小尺度油气储层波阻抗特征,因此对于当前的复杂构造精细解释必须发展新的波阻抗反演算法。
中国专利公开文献CN201610797648.8公开了深度域地震波阻抗反演方法及系统,通过支持向量机建立测井数据中的波阻抗数据与井旁深度域地震数据之间的非线性关系,旨在在数据匮乏的情况下仍然可以获得比较准确的非线性关系;中国专利公开文献CN201610393417.8公开了全波形反演方法及装置,通过对实际记录和模拟记录进行能量归一化处理,实现了对能量不均衡的陆地资料的稳定地、高精度地全波形反演;中国公开文献“井间地震约束下的高分辨率波阻抗反演方法研究”(石油物探,2010年05期),建立井间地震约束下的地面地震高分辨率波阻抗反演目标函数,实现对目标函数的有效求解;中国公开文献“基于声学全波形反演的油气藏地震成像方法”(地球物理学报,2014年2期) 通过对西部地区某气田二维地震数据处理,实现了致密砂岩气藏成像。这些公开文献中的方案都不适用于对复杂构造的精细解释。
发明内容
本发明的目的在于解决上述现有技术中存在的难题,提供一种基于全波形反演的高精度波阻抗反演方法及系统,利用包含地震波运动学和动力学信息的全波波场通过目前最为先进的全波形反演技术建立适应于复杂地质构造的高精度波阻抗反演系统,提高复杂探区油气储层预测精度,降低勘探风险。
本发明是通过以下技术方案实现的:
一种基于全波形反演的高精度波阻抗反演方法,采用基于CPML边界条件的速度-阻抗方程进行地震波数值正演得到地震数值模拟记录,利用所述地震数值模拟记录与实际观测数据的误差通过全波形反演获得高精度波阻抗。
所述方法包括:
(1)给定一个初始步长α0,初始化n=0;输入野外采集到的地震数据,通过层析反演获得初始速度vp和初始波阻抗
Figure 100002_DEST_PATH_IMAGE002
(2)对地震数据进行地震正演,得到地震数值模拟记录d;
(3)计算所述地震数值模拟记录d和实际观测数据dobs之间的差值得到地震记录误差,利用地震记录误差计算出目标泛函J的值,判断目标泛函J的值是否小于设定的阈值,如果是,则进入步骤(10),如果否,则进入步骤(4);
(4)将所述地震记录误差作为震源,进行反向传播,得到反传波场q;
(5)利用所述反传波场q计算得到波阻抗更新的梯度项gradIJ[Ip,vp];
(6)利用所述初始步长α0和所述波阻抗更新的梯度项gradIJ[Ip,vp]计算得到新的波阻抗:
Figure 636205DEST_PATH_IMAGE002
(7)利用所述新的波阻抗计算得到当前的步长α;
(8)利用所述当前的步长α计算得到当前的全波形反演阻抗:
Figure BDA0001776059310000031
(9)将当前的步长α作为初始步长α0,n=n+1,然后返回步骤(2);
(10)输出当前的全波形反演阻抗。
所述步骤(4)的操作包括:
将所述地震记录误差作为震源,利用一阶速度-阻抗波动方程的伴随状态方程进行计算,得到反传波场q。
所述步骤(5)的操作包括:利用下式计算波阻抗更新的梯度项gradIJ[Ip,vp]:
gradIJ[Ip,vp]=-〈q,DFI[Ip,vp]>
其中,DFI[Ip,vp]表示基于CPML边界条件的速度-阻抗方程的一阶Born近似。
所述步骤(7)的操作包括:
利用所述新的波阻抗进行地震正演得到地震记录,计算该地震记录与实际观测数据之间的差值,并根据该差值计算得到目标泛函J的新值Jtemp,然后利用下式计算得到当前的步长α:
Figure BDA0001776059310000032
所述步骤(2)中的地震正演、所述步骤(7)中的地震正演均是利用基于 CPML边界条件的速度-阻抗方程实现的,具体如下式:
Figure BDA0001776059310000033
其中p为压力场,vx,vy,vz为x,y,z三个方向的质点震动速度,参数χi表示在i(x,y,z)方向的吸收系数,
Figure BDA0001776059310000041
表示在i方向的递归函数。
所述步骤(3)、步骤(7)中的目标泛函J如下:
Figure BDA0001776059310000042
一种基于全波形反演的高精度波阻抗反演系统,包括:
层析反演单元:输入野外采集到的地震数据,通过层析反演获得初始速度和初始波阻抗;
正演单元:与层析反演单元连接,用于对地震数据进行地震正演,得到地震数值模拟记录d;
地震记录误差计算单元:与正演单元连接,用于计算所述地震数值模拟记录d和实际观测数据dobs之间的差值得到地震记录误差,同时利用地震记录误差计算出目标泛函J的值;
判断单元:与地震记录误差计算单元连接,用于判断目标泛函J的值是否小于设定的阈值,如果是,则发送输出信号给输出单元,如果否,则将地震记录误差发送给反向传播单元连接;
反向传播单元:将判断单元发送来的所述地震记录误差作为震源,进行反向传播,得到反传波场q;
梯度更新单元:与反向传播单元连接,利用所述反传波场q计算得到波阻抗更新的梯度项gradIJ[Ip,vp];
波阻抗更新单元:分别与正演单元、地震记录误差计算单元、梯度更新单元、参数更新单元连接,利用参数更新单元发送来的初始步长α0、n和梯度更新单元发送来的波阻抗更新的梯度项gradIJ[Ip,vp]计算得到新的波阻抗:
Figure 969097DEST_PATH_IMAGE003
; 利用所述新的波阻抗计算得到当前的步长α,并将当前的步长α发送给参数更新单元,利用所述当前的步长α计算得到当前的全波形反演阻抗:
Figure BDA0001776059310000051
并将当前的全波形反演阻抗发送给正演单元;
参数更新单元:与波阻抗更新单元连接,给定一个初始步长α0,初始化n=0,将初始步长α0和n的值发送给波阻抗更新单元;当接收到波阻抗更新单元发送来的当前的步长α后,将当前的步长α作为初始步长α0,并计算n=n+1,然后将初始步长α0和n的值发送给波阻抗更新单元;
输出单元:接收到判断单元发送来的输出信号后,从波阻抗更新单元获得当前的全波形反演阻抗,然后将当前的全波形反演阻抗输出。
所述正演单元采用基于CPML边界条件的速度-阻抗方程对地震数据进行地震正演:
Figure BDA0001776059310000052
其中,p为压力场,vx,vy,vz为x,y,z三个方向的质点震动速度,参数χi表示在i(x,y,z)方向的吸收系数,
Figure BDA0001776059310000053
表示在i方向的递归函数;
所述反向传播单元将所述地震记录误差作为震源,利用一阶速度-阻抗波动方程的伴随状态方程进行计算,得到反传波场q;
所述梯度更新单元利用下式计算波阻抗更新的梯度项gradIJ[Ip,vp]:
gradIJ[Ip,vp]=-<q,DFI[Ip,vp]>
其中,DFI[Ip,vp]表示基于CPML边界条件的速度-阻抗方程的一阶Born近似;
所述地震记录误差计算单元采用下式计算目标泛函J的值:
Figure BDA0001776059310000061
所述波阻抗更新单元利用所述新的波阻抗计算得到当前的步长α是这样实现的:
波阻抗更新单元将所述新的波阻抗发送给正演单元,正演单元进行地震正演得到地震记录,正演单元将地震记录发送给地震记录误差计算单元,地震记录误差计算单元根据地震记录计算得到目标泛函J的新值Jtemp,并将Jtemp发送给波阻抗更新单元,波阻抗更新单元利用下式计算得到当前的步长α:
Figure BDA0001776059310000062
本发明还提供一种计算机可读存储介质,所述计算机可读存储介质存储有计算机可执行的至少一个程序,所述至少一个程序被所述计算机执行时使所述计算机执行所述的基于全波形反演的高精度波阻抗反演方法中的步骤。
与现有技术相比,本发明的有益效果是:
本发明方法通过全波形反演得到的阻抗不仅在井位处获得了准确的阻抗值,也反演得到了真实模型中的构造形态,且断层断块边界清晰。采用全波形反演获得的阻抗模型进行地震资料解释工作将会大幅提高油气预测的准确性和可靠性。
附图说明
图1Marmousi2模型真实阻抗;
图2基于稀疏脉冲反演的波阻抗;
图3基于全波形反演的波阻抗;
图4本发明方法的步骤框图。
具体实施方式
下面结合附图对本发明作进一步详细描述:
(1)方法原理
①基于速度-阻抗波动方程的地震正演
地震正演是地震反演的基础。常用的声波方程地震正演通常采用时间域速度-密度声波方程来表示:
Figure BDA0001776059310000071
其中v(x)表示速度,ρ(x)为密度,P(x,t)为压力场,f(x,t)为时间域震源,xs为震源位置坐标,δ(x-xs)表示脉冲函数。为了利用波动方程进行波阻抗反演,本发明根据波阻抗与速度密度的关系,建立了基于声波方程的一阶速度-阻抗波动方程:
Figure BDA0001776059310000072
其中p为压力场,Ip为纵波阻抗,vp为地震波传播的纵波速度。vx,vy, vz为x,y,z三个方向的质点震动速度。为了消除地震波数值模拟中由于模型在空间有限引起的边界反射,本发明采用了复频移完全匹配层吸收边界条件(CPML),
从而得到基于CPML边界条件的速度-阻抗方程:
Figure BDA0001776059310000081
其中参数χi表示在i(x,y,z)方向的吸收系数,
Figure BDA0001776059310000082
表示在i方向的递归函数。在两个参数共同作用下消除人为边界反射。
②基于全波形反演的波阻抗反演
地震反演经典的目标泛函J为波场误差的二范数形式:
Figure BDA0001776059310000083
其中d为模拟数据,dobs为观测数据。定义正演算法d=RF[m],R表示从模拟波场F[m]中提取检波点处接收到的波场中提取地震记录。
对于线性正演问题,模型参数转化为反射系数,正演F能够用显式的矩阵形式表达出来,可通过直接求解法得到最小二乘解估计。对于波动方程非线性正演模拟方法,F是算子无法实施直接求解法。对于非线性反演问题,可以通过线性化方法将目标函数局部线性化,构成一种迭代形式,用逐步逼近的方法求解。本发明根据介质弱扰动假设,利用Born近似理论得到扰动波场和阻抗扰动之间的关系:
δp=-Θ-1(DFI[Ip,vp]δIp)
式中,δp表示扰动波场,δIp表示阻抗扰动,DFI[Ip,vp]表示基于CPML边界条件的速度-阻抗方程对阻抗的一阶Born近似,Θ-1表示基于Born近似的波场反传算子。在此引入q,并且定义:
Figure BDA0001776059310000084
该式即为伴随状态方程。其中q为伴随波场,是剩余波场的反传波场。则目标泛函的扰动δJ可表示为:
δJ[Ip,vp]=-<q,(DFI[Ip,vp]δIp)> (1-5)
因此可以得到目标泛函对阻抗的梯度为:
gradIJ[Ip,vp]=-<q,DFI[Ip,vp]> (1-6)
通过迭代反演可以得到新的阻抗:
Figure BDA0001776059310000091
其中α表示迭代步长。上式就是最终波阻抗反演公式。
常规的波阻抗反演中采用的正演算法是基于射线理论利用褶积模型只能获得单道地震数据,且无法表达地震波的动力学信息,在地震数据匹配中只能对运动学信息进行一维匹配。本发明中采用基于速度-阻抗的波动方程进行地震波数值模拟和波场匹配,不仅能够得到运动学信息,同时能够充分利用地震数据的振幅、频率、相位等动力学信息。
(2)本发明方法的步骤如图4所示:
a:对野外采集到的地震数据通过层析反演建立初始速度和初始波阻抗
Figure DEST_PATH_IMAGE004
; 给定一个初始步长α0
b:利用基于CPML边界条件的速度-阻抗方程进行地震波数值正演,得到包含运动学和动力学信息的地震记录,将该地震记录称为地震数值模拟记录(即图4中的炮域数据);
c:将所述地震数值模拟记录和实际观测数据进行运动学和动力学信息的匹配,得到地震记录误差(即利用d-dobs计算得到地震记录误差),同时计算出目标泛函J的值,即利用公式(1-3)计算J的值;
d:将地震记录误差(即图4中的炮域数据误差)作为震源,利用一阶速度 -阻抗波动方程的伴随状态方程(即公式(1-4))进行反向传播,得到反传波场,即伴随波场q;
e:根据公式(1-6)计算阻抗更新的梯度项gradIJ[Ip,vp];
f:利用初始步长α0得到一个新的阻抗
Figure BDA0001776059310000101
g:利用该新的阻抗进行波动方程正演(即采用公式(1-2)进行正演),得到的地震记录与观测记录做差就得到一个新的目标泛函值Jtemp,利用三点抛物法可以得到当前的步长α:
Figure BDA0001776059310000102
h:进而得到本次迭代的全波形反演阻抗:
Figure BDA0001776059310000103
i:然后利用新的阻抗通过b-e的过程得到新的梯度,并利用当前的α作为初始步长,利用f-h的过程得到新的步长并进行阻抗更新;
j:依次循环i的过程,当目标泛函J的值足够小(是将J的值与一个设定的阈值进行比较来判断其是否足够小,该阈值的大小是根据实际情况确定的) 时,终止迭代,得到最终的全波形反演阻抗。
本发明还提供一种基于全波形反演的高精度波阻抗反演系统,包括:
层析反演单元:输入野外采集到的地震数据,通过层析反演获得初始速度和初始波阻抗;
正演单元:与层析反演单元连接,用于对地震数据进行地震正演,得到地震数值模拟记录d;
地震记录误差计算单元:与正演单元连接,用于计算所述地震数值模拟记录d和实际观测数据dobs之间的差值得到地震记录误差,同时利用地震记录误差计算出目标泛函J的值;
判断单元:与地震记录误差计算单元连接,用于判断目标泛函J的值是否小于设定的阈值,如果是,则发送输出信号给输出单元,如果否,则将地震记录误差发送给反向传播单元连接;
反向传播单元:将判断单元发送来的所述地震记录误差作为震源,进行反向传播,得到反传波场q;
梯度更新单元:与反向传播单元连接,利用所述反传波场q计算得到波阻抗更新的梯度项gradIJ[Ip,vp];
波阻抗更新单元:分别与正演单元、地震记录误差计算单元、梯度更新单元、参数更新单元连接,利用参数更新单元发送来的初始步长α0、n和梯度更新单元发送来的波阻抗更新的梯度项gradIJ[Ip,vp]计算得到新的波阻抗:
Figure 783470DEST_PATH_IMAGE003
; 利用所述新的波阻抗计算得到当前的步长α,并将当前的步长α发送给参数更新单元,利用所述当前的步长α计算得到当前的全波形反演阻抗:
Figure BDA0001776059310000112
并将当前的全波形反演阻抗发送给正演单元;
参数更新单元:与波阻抗更新单元连接,给定一个初始步长α0,初始化n=0,将初始步长α0和n的值发送给波阻抗更新单元;当接收到波阻抗更新单元发送来的当前的步长α后,将当前的步长α作为初始步长α0,并计算n=n+1,然后将初始步长α0和n的值发送给波阻抗更新单元;
输出单元:接收到判断单元发送来的输出信号后,从波阻抗更新单元获得当前的全波形反演阻抗,然后将当前的全波形反演阻抗输出。
所述正演单元采用基于CPML边界条件的速度-阻抗方程对地震数据进行地震正演;
所述反向传播单元将所述地震记录误差作为震源,利用一阶速度-阻抗波动方程的伴随状态方程进行计算,得到反传波场q;
所述梯度更新单元利用下式计算波阻抗更新的梯度项gradIJ[Ip,vp]:
gradIJ[Ip,vp]=-〈q,DFI[Ip,vp]>
其中,DFI[Ip,vp]表示基于CPML边界条件的速度-阻抗方程的一阶Born近似;
所述地震记录误差计算单元采用下式计算目标泛函J的值:
Figure BDA0001776059310000121
所述波阻抗更新单元利用所述新的波阻抗计算得到当前的步长α是这样实现的:
波阻抗更新单元将所述新的波阻抗发送给正演单元,正演单元进行地震正演得到地震记录,正演单元将地震记录发送给地震记录误差计算单元,地震记录误差计算单元根据地震记录计算得到目标泛函J的新值Jtemp,并将Jtemp发送给波阻抗更新单元,波阻抗更新单元利用下式计算得到当前的步长α:
Figure BDA0001776059310000122
本发明的一个实施例如下:
通过Marmousi2模型中的一部分对本发明进行测试。对原始的模型进行重采样,其中深度方向5米采样,共301个采样点,测线方向1500个cdp,cdp 间距定义为5米。图1为真实的阻抗模型。常规波阻抗反演方法是基于稀疏脉冲反演的波阻抗反演方法,该方法需要利用测井上的速度建立低频趋势模型,因此抽取cdp=750处的真实速度作为测井速度(如图1中黑线的位置)。图2为稀疏脉冲反演得到的波阻抗,图3为利用本发明全波形反演得到的波阻抗。通过对比常规波阻抗反演得到的波阻抗和本发明全波形反演得到的波阻抗可以发现利用稀疏脉冲反演只能在井位置处得到准确的阻抗值,而在整个二维空间中不能反演出真实模型的构造形态,横向连续性强,严重缺失横向分辨率。
本发明推导了基于速度-阻抗的波动方程,利用观测数据和实际数据的误差建立目标泛函,建立了一种基于全波形反演的波阻抗反演系统。该发明充分利用地震数据中的运动学和动力学信息,实现了复杂地质构造高精度波阻抗反演。
上述技术方案只是本发明的一种实施方式,对于本领域内的技术人员而言,在本发明公开了应用方法和原理的基础上,很容易做出各种类型的改进或变形,而不仅限于本发明上述具体实施方式所描述的方法,因此前面描述的方式只是优选的,而并不具有限制性的意义。

Claims (8)

1.一种基于全波形反演的高精度波阻抗反演方法,其特征在于:所述方法采用基于CPML边界条件的速度-阻抗方程进行地震波数值正演得到地震数值模拟记录,利用所述地震数值模拟记录与实际观测数据的误差通过全波形反演获得高精度波阻抗;
所述方法包括:
(1)给定一个初始步长α0,初始化n=0;输入野外采集到的地震数据,通过层析反演获得初始速度vp和初始波阻抗
Figure DEST_PATH_IMAGE002
(2)对地震数据进行地震正演,得到地震数值模拟记录d;
(3)计算所述地震数值模拟记录d和实际观测数据dobs之间的差值得到地震记录误差,利用地震记录误差计算出目标泛函J的值,判断目标泛函J的值是否小于设定的阈值,如果是,则进入步骤(10),如果否,则进入步骤(4);
(4)将所述地震记录误差作为震源,进行反向传播,得到反传波场q;
(5)利用所述反传波场q计算得到波阻抗更新的梯度项gradIJ[Ip,vp];
(6)利用所述初始步长α0和所述波阻抗更新的梯度项gradIJ[Ip,vp]计算得到新的波阻抗:
Figure 466458DEST_PATH_IMAGE002
Figure FDA0003114451410000013
为第n次迭代后的全波形反演阻抗;
(7)利用所述新的波阻抗计算得到当前的步长α;
(8)利用所述当前的步长α计算得到当前的全波形反演阻抗:
Figure FDA0003114451410000014
(9)将当前的步长α作为初始步长α0,n=n+1,然后返回步骤(2);
(10)输出当前的全波形反演阻抗;
所述步骤(2)中的地震正演、所述步骤(7)中的地震正演均是利用基于CPML边界条件的速度-阻抗方程实现的,具体如下式:
Figure FDA0003114451410000021
其中,p为压力场,vx,vy,vz为x,y,z三个方向的质点震动速度,参数χi表示在i(x,y,z)方向的吸收系数,
Figure FDA0003114451410000022
表示在i方向的递归函数,Ip为纵波阻抗,vp为地震波传播的纵波速度,t表示时间,
Figure FDA0003114451410000023
分别表示x、y、z三个方向上CPML层内与时间有关的记忆变量。
2.根据权利要求1所述的基于全波形反演的高精度波阻抗反演方法,其特征在于:所述步骤(4)的操作包括:
将所述地震记录误差作为震源,利用一阶速度-阻抗波动方程的伴随状态方程进行计算,得到反传波场q。
3.根据权利要求1所述的基于全波形反演的高精度波阻抗反演方法,其特征在于:所述步骤(5)的操作包括:利用下式计算波阻抗更新的梯度项gradIJ[Ip,vp]:
gradIJ[Ip,vp]=-<q,DFI[Ip,vp]>
其中,DFI[Ip,vp]表示基于CPML边界条件的速度-阻抗方程的一阶Born近似。
4.根据权利要求1所述的基于全波形反演的高精度波阻抗反演方法,其特征在于:所述步骤(7)的操作包括:
利用所述新的波阻抗进行地震正演得到地震记录,计算该地震记录与实际观测数据之间的差值,并根据该差值计算得到目标泛函J的新值Jtemp,然后利用下式计算得到当前的步长α:
Figure FDA0003114451410000031
5.根据权利要求4所述的基于全波形反演的高精度波阻抗反演方法,其特征在于:所述步骤(3)、步骤(7)中的目标泛函J如下:
Figure FDA0003114451410000032
6.一种实现权利要求1-5任一所述方法的基于全波形反演的高精度波阻抗反演系统,其特征在于:所述系统包括:
层析反演单元:输入野外采集到的地震数据,通过层析反演获得初始速度和初始波阻抗;
正演单元:与层析反演单元连接,用于对地震数据进行地震正演,得到地震数值模拟记录d;
地震记录误差计算单元:与正演单元连接,用于计算所述地震数值模拟记录d和实际观测数据dobs之间的差值得到地震记录误差,同时利用地震记录误差计算出目标泛函J的值;
判断单元:与地震记录误差计算单元连接,用于判断目标泛函J的值是否小于设定的阈值,如果是,则发送输出信号给输出单元,如果否,则将地震记录误差发送给反向传播单元连接;
反向传播单元:将判断单元发送来的所述地震记录误差作为震源,进行反向传播,得到反传波场q;
梯度更新单元:与反向传播单元连接,利用所述反传波场q计算得到波阻抗更新的梯度项gradIJ[Ip,vp];
波阻抗更新单元:分别与正演单元、地震记录误差计算单元、梯度更新单元、参数更新单元连接,利用参数更新单元发送来的初始步长α0、n和梯度更新单元发送来的波阻抗更新的梯度项gradjJ[Ip,vp]计算得到新的波阻抗:
Figure 383598DEST_PATH_IMAGE002
; 利用所述新的波阻抗计算得到当前的步长α,并将当前的步长α发送给参数更新单元,利用所述当前的步长α计算得到当前的全波形反演阻抗:
Figure FDA0003114451410000042
并将当前的全波形反演阻抗发送给正演单元;
Figure FDA0003114451410000043
为第n次迭代的全波形反演阻抗;
参数更新单元:与波阻抗更新单元连接,给定一个初始步长α0,初始化n=0,将初始步长α0和n的值发送给波阻抗更新单元;当接收到波阻抗更新单元发送来的当前的步长α后,将当前的步长α作为初始步长α0,并计算n=n+1,然后将初始步长α0和n的值发送给波阻抗更新单元;
输出单元:接收到判断单元发送来的输出信号后,从波阻抗更新单元获得当前的全波形反演阻抗,然后将当前的全波形反演阻抗输出。
7.根据权利要求6所述的基于全波形反演的高精度波阻抗反演系统,其特征在于:所述正演单元采用基于CPML边界条件的速度-阻抗方程对地震数据进行地震正演:
Figure FDA0003114451410000044
其中,p为压力场,vx,vy,vz为x,y,z三个方向的质点震动速度,参数χi表示在i(x,y,z)方向的吸收系数,
Figure FDA0003114451410000045
表示在i方向的递归函数;Ip为纵波阻抗,vp为地震波传播的纵波速度,t表示时间,
Figure FDA0003114451410000051
分别表示x、y、z三个方向上CPML层内与时间有关的记忆变量;
所述反向传播单元将所述地震记录误差作为震源,利用一阶速度-阻抗波动方程的伴随状态方程进行计算,得到反传波场q;
所述梯度更新单元利用下式计算波阻抗更新的梯度项gradIJ[Ip,vp]:
gradIJ[Ip,vp]=-<q,DFI[Ip,vp]>
其中,DFI[Ip,vp]表示基于CPML边界条件的速度-阻抗方程的一阶Born近似;
所述地震记录误差计算单元采用下式计算目标泛函J的值:
Figure FDA0003114451410000052
所述波阻抗更新单元利用所述新的波阻抗计算得到当前的步长α是这样实现的:
波阻抗更新单元将所述新的波阻抗发送给正演单元,正演单元进行地震正演得到地震记录,正演单元将地震记录发送给地震记录误差计算单元,地震记录误差计算单元根据地震记录计算得到目标泛函J的新值Jtemp,并将Jtemp发送给波阻抗更新单元,波阻抗更新单元利用下式计算得到当前的步长α:
Figure FDA0003114451410000053
8.一种计算机可读存储介质,其特征在于:所述计算机可读存储介质存储有计算机可执行的至少一个程序,所述至少一个程序被所述计算机执行时使所述计算机执行所述权利要求1-5任一所述的基于全波形反演的高精度波阻抗反演方法中的步骤。
CN201810970557.6A 2018-08-24 2018-08-24 一种基于全波形反演的高精度波阻抗反演方法及系统 Active CN110857999B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810970557.6A CN110857999B (zh) 2018-08-24 2018-08-24 一种基于全波形反演的高精度波阻抗反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810970557.6A CN110857999B (zh) 2018-08-24 2018-08-24 一种基于全波形反演的高精度波阻抗反演方法及系统

Publications (2)

Publication Number Publication Date
CN110857999A CN110857999A (zh) 2020-03-03
CN110857999B true CN110857999B (zh) 2021-12-31

Family

ID=69636158

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810970557.6A Active CN110857999B (zh) 2018-08-24 2018-08-24 一种基于全波形反演的高精度波阻抗反演方法及系统

Country Status (1)

Country Link
CN (1) CN110857999B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113970789B (zh) * 2020-07-24 2024-04-09 中国石油化工股份有限公司 全波形反演方法、装置、存储介质及电子设备

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1635193A1 (en) * 2004-09-13 2006-03-15 Services Petroliers Schlumberger Determining impedance of material behind a casing in a borehole
CN104965222A (zh) * 2015-05-29 2015-10-07 中国石油天然气股份有限公司 三维纵波阻抗全波形反演方法及装置
CN105308479A (zh) * 2013-05-24 2016-02-03 埃克森美孚上游研究公司 通过与偏移距相关的弹性fwi的多参数反演
CN105467444A (zh) * 2015-12-10 2016-04-06 中国石油天然气集团公司 一种弹性波全波形反演方法及装置

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
AU2015382333B2 (en) * 2015-02-13 2018-01-04 Exxonmobil Upstream Research Company Efficient and stable absorbing boundary condition in finite-difference calculations

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1635193A1 (en) * 2004-09-13 2006-03-15 Services Petroliers Schlumberger Determining impedance of material behind a casing in a borehole
CN105308479A (zh) * 2013-05-24 2016-02-03 埃克森美孚上游研究公司 通过与偏移距相关的弹性fwi的多参数反演
CN104965222A (zh) * 2015-05-29 2015-10-07 中国石油天然气股份有限公司 三维纵波阻抗全波形反演方法及装置
CN105467444A (zh) * 2015-12-10 2016-04-06 中国石油天然气集团公司 一种弹性波全波形反演方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Full waveform inversion of diving & reflected waves for velocity model building with impedance inversion based on scale separation;Wei Zhou et al.;《Geophysical Journal International》;20151231;第1541-1542页 *
基于褶积完全匹配吸收边界的声波方程数值模拟;柯璇等;《石油物探》;20170930;第56卷(第5期);第638-639页 *

Also Published As

Publication number Publication date
CN110857999A (zh) 2020-03-03

Similar Documents

Publication Publication Date Title
RU2693495C1 (ru) Полная инверсия волнового поля с компенсацией показателя качества
EP2491428B1 (en) Full-waveform inversion in the traveltime domain
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
CN109188519B (zh) 一种极坐标下的弹性波纵横波速度反演系统及方法
CN104570082B (zh) 一种基于格林函数表征的全波形反演梯度算子的提取方法
CN113740901B (zh) 基于复杂起伏地表的陆上地震数据全波形反演方法及装置
EP2567063A1 (en) Artifact reduction in method of iterative inversion of geophysical data
CN109557582B (zh) 一种二维多分量地震资料偏移成像方法及系统
US10345466B2 (en) Memory efficient Q-RTM computer method and apparatus for imaging seismic data
CN105319589A (zh) 一种利用局部同相轴斜率的全自动立体层析反演方法
CN109459787B (zh) 基于地震槽波全波形反演的煤矿井下构造成像方法及系统
Lamert et al. Imaging disturbance zones ahead of a tunnel by elastic full-waveform inversion: Adjoint gradient based inversion vs. parameter space reduction using a level-set method
Nguyen et al. Reconstructing disturbance zones ahead of the tunnel face by elastic waveform inversion supported by a parametric level-set representation
Qu et al. Topographic elastic least‐squares reverse time migration based on vector P‐and S‐wave equations in the curvilinear coordinates
CN104597489B (zh) 一种震源子波优化设置方法和装置
CN111665556B (zh) 地层声波传播速度模型构建方法
Trapp et al. Non-gradient full waveform inversion approaches for exploration during mechanized tunneling applied to surrogate laboratory measurements
CN110857999B (zh) 一种基于全波形反演的高精度波阻抗反演方法及系统
CN112305595B (zh) 基于折射波分析地质体结构的方法及存储介质
CN114966831A (zh) 一种基于速度-衰减解耦的粘声全波形反演方法
CN115373022A (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
CN111665546B (zh) 用于可燃冰探测的声学参数获取方法
CN111175822A (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
CN111665551B (zh) 用于桥梁基底探测的声学参数获取方法
CN113589362B (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