CN113031058A - 一种基于反射系数精确式的页岩vti储层叠前混合反演方法 - Google Patents
一种基于反射系数精确式的页岩vti储层叠前混合反演方法 Download PDFInfo
- Publication number
- CN113031058A CN113031058A CN202110218907.5A CN202110218907A CN113031058A CN 113031058 A CN113031058 A CN 113031058A CN 202110218907 A CN202110218907 A CN 202110218907A CN 113031058 A CN113031058 A CN 113031058A
- Authority
- CN
- China
- Prior art keywords
- iteration
- inversion
- reflection coefficient
- particles
- matrix
- 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
- 238000000034 method Methods 0.000 title claims abstract description 56
- 238000004364 calculation method Methods 0.000 claims abstract description 15
- 239000002245 particle Substances 0.000 claims description 88
- 239000011159 matrix material Substances 0.000 claims description 57
- 239000013598 vector Substances 0.000 claims description 30
- 239000000126 substance Substances 0.000 claims description 8
- 238000005070 sampling Methods 0.000 claims description 6
- 230000002194 synthesizing effect Effects 0.000 claims description 5
- 150000001875 compounds Chemical class 0.000 claims description 4
- 230000017105 transposition Effects 0.000 claims description 2
- PUZPDOWCWNUUKD-UHFFFAOYSA-M sodium fluoride Chemical compound [F-].[Na+] PUZPDOWCWNUUKD-UHFFFAOYSA-M 0.000 claims 1
- 230000007547 defect Effects 0.000 abstract description 2
- 238000005457 optimization Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 2
- 238000013398 bayesian method Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 230000035945 sensitivity Effects 0.000 description 1
- 239000003079 shale oil Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
- G01V1/362—Effecting static or dynamic corrections; Stacking
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
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
本发明公开了一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,该混合反演策略将线性和非线性反演相结合,利用线性反演结果作为约束,生成较小的搜索范围和初始粒子,开展非线性反演以同时获取五参数的反演结果。该方法克服常规线性反演易落入局部极值的缺点,提高非线性反演的计算效率和准确性。采用反射系数精确式可克服近似表达式小角度、弱反射和弱各向异性的限制,提高反演的适用性。测井模型数据测试验证了方法的有效性,这为后续VTI介质各向异性研究奠定了基础。
Description
技术领域
本发明涉及非常规各向异性储层地震勘探技术领域,特别是一种基于反射系数精确式的页岩VTI储层叠前混合反演方法。
背景技术
叠前地震反演通过挖掘地震数据中的AVO/AVA(振幅随偏移距/入射角度变化)响应来获取地下介质的弹性信息。目前,叠前反演已是地震勘探的重要技术。以往,针对常规砂泥岩储层,均采用基于各向同性假设的反演方法,并在过去几十年里取得了巨大的成功。但随着油气需求不断增加,勘探的重点逐渐由常规储层向非常规油气转变。其中,页岩油气已经成为了目前勘探的重点与难点。页岩岩石表现为较强的VTI各向异性特征。VTI介质,即具有垂直对称轴的横向各向同性介质(Transverse isotropy with vertical axis ofsymmetry,VTI)是地下介质中存在最普遍的各向异性。以往基于各向同性假设的地震技术已不再适用于页岩储层的勘探。因此,为提高页岩储层的勘探精度,需开展适用于VTI介质的反演方法研究,用以有效挖掘页岩储层的弹性及各向异性信息。
目前,针对VTI介质的正演研究较为完善,但是针对VTI介质的叠前反演方法并不完善,主要是选用反射系数近似式为正演算子。侯栋甲等(2014)将反射系数近似式引入VTI五参数直接反演中,并采用贝叶斯方法建立目标函数。但由于各向异性介质反演需要同时提取五个目标参数,且各向异性参数的地震敏感性较低,这导致VTI反演的不适定性更强,很难同时获取较好的五参数结果,尤其是各向异性参数的估算一直不理想。因此,Zhang等(2019)基于近似式提出了分步反演策略,通过数学转换将五参数转换为其他三项目标参数,然后通过间接反演方法换算得到各向异性参数,但是该方法引入了额外假设,且仍采用精度较低的近似式。VTI介质反射系数精确表达式虽然已经被给出,虽然其复杂的数学表达式,导致偏导数求取较难,因此目前较少被应用于叠前反演中。
常规线性反演易落入局部极值,而非线性反演存在搜索范围大计算效率低的问题,近似式反演远角度精度低、弱各向异性限制等问题。
发明内容
本发明所要解决的技术问题是克服现有技术的不足而提供一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,针对目前页岩等VTI(具有垂直同相轴的横向各向同性)介质的反演现存的难点问题,提出了一种将线性和非线性优化算法相结合的混合反演策略,并将VTI介质的反射系数精确式引入作为正演算子,用以提高页岩等VTI介质的各向异性参数的反演精度。
本发明为解决上述技术问题采用以下技术方案:
根据本发明提出的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,包括以下步骤:
步骤一:获取PP波叠前角度道集,并从中提取地震子波,并建立用于反演的子波矩阵;
步骤二:获取工区的井孔数据和三维层位解释结果,建立反演所用的初始模型;
步骤三:预设线性反演循环迭代的最大迭代次数,开始线性反演循环迭代,利用初始模型和子波矩阵获取合成地震记录;
步骤四:利用初始模型计算反射系数精确式正演算子的偏导数;
步骤五:建立目标函数,基于正演算子的偏导数和合成地震记录,计算目标函数的雅可比矩阵和伪海森矩阵;
步骤六:利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量的更新方向,利用迭代格式更新目标参数向量;
步骤七:若线性反演循环迭代次数未达到预设的线性反演循环迭代的最大迭代次数,则开始下一次线性反演迭代,重复步骤三至六;若达到,则结束线性反演循环迭代过程,输出线性反演结果,继续执行步骤八;
步骤八:利用线性反演结果,生成搜索范围和初始粒子;
步骤十:将目标函数应用于所有粒子来计算替换适应度值,并计算不同粒子的接收概率,利用轮盘赌获取群体最优位置xgbest;
步骤十一:根据个体最优位置和群体最优位置来计算更新速度,进而采用迭代格式更新粒子;
步骤十二:若迭代次数未达到预设的非线性反演迭代的最大迭代次数,则开始下一次非线性反演循环迭代,重复步骤九至十一;若达到最大迭代次数,则结束非线性反演,输出所有粒子,选取所有粒子中的最好结果作为最终的非线性反演结果。
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤二中,井孔数据包括纵波速度α、横波速度β、密度ρ,各向异性参数ε和δ曲线,与四个刚度参数c11、c13、c33和c55之间的关系为:
c33=α2ρ,c55=β2ρ
c11=(2ε+1)α2ρ,c13=ρη-β2ρ
其中,η为中间变量;
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤三中,获取合成地震记录具体如下:
利用VTI介质反射系数精确方程计算反射系数r:
A·r=b
其中,A和b均为中间过程矩阵,具有如下表达形式:
其中,为上层介质纵波相关的中间变量,为上层介质横波相关的中间变量,为下层介质纵波相关的中间变量,为下层介质横波相关的中间变量; 为上层介质的刚度参数,为下层介质的刚度参数;分别为上、下层纵波相关的垂直慢度,分别为上、下层横波相关的垂直慢度,p为水平慢度; p这些参数均由刚度参数c11、c13、c33、c55和密度ρ计算获得,上标T表示转置运算符;
基于褶积理论,合成地震记录dsyn利用子波矩阵W和反射系数r计算求得:
dsyn=G(m)=W·r
其中,m为设定的目标参数向量:
m=(α1,...,αM,β1,...,βM,ρ1,...,ρM,ε1,...,εM,δ1,...,δM)T
αi、βi、ρi、εi和δi分别表示为第i个时间采样点所对应的纵波速度、横波速度、密度、两个各向异性参数,i=1,....,M;M为模型时间深度最大样点个数,G(m)为反射系数精确式正演算子,即合成地震记录。
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤四中,计算反射系数精确式正演算子的偏导数具体如下:
反射系数精确式正演算子G(m)对目标参数向量m的偏导数为:
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤五中,
建立目标函数S(m),表示为:
S(m)=[d-G(m)]T[d-G(m)]+λ(m-u)TC-m1(m-u)
其中,d为实际叠前地震道集数据,λ为正则化参数,Cm为m计算得到的协方差矩阵,u为m的期望;
计算目标函数S(m)的雅可比矩阵J(m),即一阶偏导数矩阵:
针对目标函数S(m)的伪海森矩阵H(m),采用了内部小循环迭代形式来计算获得:
对于线性反演循环迭代中第k次迭代、内部小循环首次迭代的伪海森矩阵H k,1(m),由下式计算得到:
对于线性反演循环迭代中第k次迭代、内部小循环第i次迭代的为海森矩阵Hk,i(m),由下式计算获得:
Hk,i(m)-1=(Qk,i-1)THk,i-1(m)-1(Qk,i-1)+Kk,i-1sk,i-1(sk,i-1)T
其中,Qk,i-1、Kk,i-1、sk,i-1为线性反演循环中第k次迭代、内部小循环第i-1次迭代的中间过程矩阵,Hk,i-1(m)为线性反演循环中第k次迭代内部小循环第i-1次迭代的伪海森矩阵。
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤六中,利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量更新方向,利用迭代格式更新目标参数向量,具体形式如下:
线性反演循环迭代的第k次迭代时,目标参数向量的更新方向-gk(m)由下式计算得到:
-gk(m)=-Hk(m)·Jk(m)
其中,Hk(m)和Jk(m)分别代表线性反演循环中第k次迭代所对应的伪海森矩阵和雅可比矩阵;
利用如下的更新迭代格式来更新mk+1:
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤八中,针对目标参数向量中某一深度样点,用mLi表示线性反演结果,则基于线性反演结果生成搜索范围[mmin mmax]:
其中,mmin为搜索范围的最小值,mmax为搜索范围的最大值,mRe表示井口处抽取到的真实参考值,b为尺度因子;
设定N个初始粒子,粒子群x表示为:
其中,j的范围是1到N,xj表示第j个粒子,也是第j个可能的模型参数;并将线性反演结果mLi赋值给N个粒子中的任意一个,最终形成初始粒子群。
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤九具体如下:
其中,l为非线性反演循环迭代次数,xj,l表示第l次迭代中的第j个粒子,S(xj,l)、S(xj,l-1)分别为第l次、第l-1次迭代中的第j个粒子的适应度。
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤十中,将目标函数应用于所有粒子来计算替换适应度值,并计算不同粒子的接收概率,利用轮盘赌获取群体最优位置xgbest包括:
计算替换适用度值Zr(xj,xgbest,Tl):
其中,xj表示第j个粒子,xgbest为群体最优位置,S(xj)和S(xgbest)分别表示将xj和xgbest代入目标函数计算得到的适用度,为xj和xgbest适用度的差;l为非线性反演循环中的迭代次数,Tl为第l次迭代所对应的温度值,τ表示为一个常数;
计算好所有粒子的替换适用度值后,计算不同粒子的接收概率Pg:
其中,∑为求和运算符;
在得到所有粒子的接收概率后,采用轮盘赌的方式获得群体最优位置xgbest。
作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案,步骤十一中,根据个体最优位置和群体最优位置来计算更新速度,进而采用迭代格式更新粒子包括:
其中,vj,l+1、vj,l表示第j个粒子在第l+1次、第l次迭代中的更新速度,rand1(·)和rand2(·)表示两个介于[0,1]的随机数;κ2和κ3为两个控制更新步长的常数;xj,l表示第l次迭代的第j个粒子;表示第l次迭代的第j个粒子的个体最优位置;
基于更新速度,更新不同的粒子:
xj,l+1=xj,l+κ1·vj,l+1
其中,κ1为用以控制更新的步长的常数;xj,l+1表示第l+1次迭代的第j个粒子。
本发明采用以上技术方案与现有技术相比,具有以下技术效果:
(1)针对目前页岩等VTI(具有垂直同相轴的横向各向同性)介质的反演现存的难点问题,提出了一种将线性和非线性优化算法相结合的混合反演策略,并将VTI介质的反射系数精确式引入作为正演算子,用以提高页岩等VTI介质的各向异性参数的反演精度;
(2)该混合反演策略将线性和非线性反演相结合,利用线性反演结果作为约束,生成较小的搜索范围和初始粒子,开展非线性反演以同时获取五参数的反演结果;该方法克服常规线性反演易落入局部极值的缺点,提高非线性反演的计算效率和准确性;采用反射系数精确式可克服近似表达式小角度、弱反射和弱各向异性的限制,提高反演的适用性。
附图说明
图1是本发明混合反演算法的流程图。
图2为基于VTI反射系数近似式的常规线性直接反演所获得的五参数结果。
图3为本发明混合反演算法的线性反演结果。
图4为基于线性反演结果所获得的搜索范围。
图5为本发明混合反演算法的非线性最终结果。
具体实施方式
下面结合附图对本发明的技术方案做进一步的详细说明:
为克服常规线性反演易落入局部极值,而非线性反演搜索范围大计算效率低的问题,本发明提出了线性及非线性混合反演策略。为克服基于近似式反演远角度精度低、弱各向异性限制等问题,提出选用反射系数精确式的作为正演算子,以提高各向异性参数的反演精度。
这里以一个测井模型测试为例进行说明:
图1是本次申请一种基于VTI介质反射系数精确式的混合反演策略的实施流程,具体包括以下步骤:
一、基于VTI介质反射系数精确式的正演
在井口处获取弹性参数曲线,包括纵波速度α、横波速度β、密度ρ、各向异性参数ε和δ。同时设定该五参数为反演的目标参数,由五参数所组成的目标参数向量m为:
m=[α1,...,αM,...1,...,βM,ρ1,...,ρM,ε1,...,εM,δ1,...,δM] (1)
其中,βi、βi、ρi、εi和δi(i=1,....,M)表示为第i个时间采样点所对应的纵波速度、横波速度、密度、两个各向异性参数。
从五个目标参数曲线中,通过Backus平均算法或者低通滤波来获取提取低频数据信息;而针对二维以及三维情况,采用井位置处的低频数据,并通过克里金插值法获得反演所用的三维初始模型。其中,地质层位解释被作为插值的横向约束。
从PP地震角度道集的不同角度范围中提取不同角度下的地震子波序列,并利用子波序列建立子波矩阵W(θ),表达式为:
其中,w(θi)表示θi角度所对应的子波矩阵:
基于褶积理论,合成地震记录dsyn可以利用子波矩阵W和反射系数序列r计算求得:
dsyn=G(m)=W·r (4)
G(m)为精确式正演算子的非线性实现过程,等同于合成地震记录。
利用VTI介质反射系数精确方程可以计算反射系数r:
A·r=b (5)
其中,A和b均为中间过程矩阵,具有如下表达形式:
其中,具有上标U和L的参数分别代表反射界面上层和下层的参数,而具有下标P和S的参数分别代表纵波和横波相关的参数。c11、c13、c33和c55为弹性刚度参数。符号p为水平慢度,s为垂直慢度。n均为中间变量参数,具体计算公式如下:
其中,
式中,sP和sS表示纵波和横波所对应的水平慢度。
上面反射系数精确式均为刚度参数c11、c13、c33、c55的函数,而纵波速度α、横波速度β、密度ρ,各向异性参数ε和δ曲线,与刚度参数c11、c13、c33、c55之间的关系如下:
c33=α2ρ,c55=β2ρ (11)
c11=(2ε+1)α2ρ,c13=ρη-β2ρ (12)
其中,η为中间变量
二、基于VTI反射系数精确式的叠前线性反演方法
建立反演目标函数S(m),表示为:
其中,d为实际叠前地震道集数据,G(m)为精确式正演算子的非线性实现过程,即合成地震记录。λ为正则化参数,Cm为目标参数向量计算得到的协方差矩阵,u为目标参数向量的期望。
式中,上标U和L表示反射界面上层介质和下层介质所对应的参数。和分别为中间变量和对上层介质模型参数mU的偏导数。和分别为中间变量和对下层介质模型参数mL的偏导数。aU、bU、dU、eU的具体计算式如下:
计算目标函数S(m)的雅可比矩阵J(m),即一阶偏导数矩阵:
针对目标函数S(m)的伪海森矩阵H(m),采用了内部小循环迭代形式来计算获得:
对于线性反演循环迭代中第k次迭代、内部小循环首次迭代的伪海森矩阵H k,1(m),可由下式计算得到:
对于线性反演循环迭代中第k次迭代、内部小循环第i次迭代的为海森矩阵Hk,i(m),可由下式计算获得:
Hk,i(m)-1=(Qk,i-1)THk,i-1(m)-1(Qk,i-1)+Kk,i-1sk,i-1(sk,i-1)T (27)
其中,Hk,i-1(m)为外部大循环第k次迭代内部小循环第i-1次迭代的伪海森矩阵;外部大循环是指线性反演循环,Qk,i-1、Kk,i-1、sk,i-1为外部大循环第k次迭代内部小循环第i-1次迭代的中间过程矩阵,计算形式如下:
Qk,i-1=I-Kk,i-1yk,i-1(sk,i-1)T (28)
Kk,i-1=1/(yk,i-1)T sk,i-1 (29)
式(27)至(29)中,
sk,i-1=mk+1-mk,yk,i-1=J(mk+1)-J(mk) (30)
上式中,mk为线性反演循环迭代中第k次迭代的目标参数向量,J(mk)为线性反演循环迭代中第k次迭代的雅可比矩阵。
线性反演循环迭代的第k次迭代时,目标参数向量的更新方向-gk(m)由下式计算得到:
-gk(m)=-Hk(m)·Jk(m) (31)
其中,Hk(m)和Jk(m)分别代表外部大循环中第k次迭代所对应的伪海森矩阵和雅可比矩阵。
利用如下的更新迭代格式来更新参数向量mk+1:
若线性反演循环迭代的迭代次数未达到最大迭代次数,则开始下一次反演迭代,重复上述内容;若达到,则结束线性反演迭代过程,输出线性反演结果
三、基于VTI反射系数精确式的叠前非线性反演方法
基于线性反演结果以生成较小范围的搜索区间。以模型中某一深度样点为例,用mLi表示线性反演结果,则基于线性反演结果可生成搜索范围[mmin mmax]:
其中,mmin为搜索范围的最小值,mmax为搜索范围的最大值,mRe表示井口处抽取到的真实参考值,b为尺度因子,可随参数不同而变化。
基于线性的反演结果设定初始粒子。设定初始粒子时,粒子群x可以表示为:
x=[x1,...,xj,...,xN] (35)
其中,xj表示第j个粒子,也是第j个可能的模型参数;并将线性反演结果mLi赋值给N个粒子中的任意一个,最终形成初始粒子群。
其中,l为非线性反演循环迭代次数,xj,l表示第l次迭代中的第j个粒子,S(xj,l)、S(xj,l-1)分别为第l次、第l-1次迭代中的第j个粒子的适应度。
计算替换适用度Zr(xj,xgbest,Tl),用以计算接收概率Pg:
其中,xj表示第j个粒子,xgbest为群体最优位置,S(xj)和S(xgbest)分别表示将xj和xgbest代入目标函数计算得到的适用度;为xj和xgbest适用度的差,l为非线性反演循环迭代的迭代次数,Tl为第l次迭代中的温度,τ表示为一个常数;
计算好所有粒子的替换适用度后,可计算获得不同粒子的接收概率Pg:
其中,vj,l+1,vj,l表示第j个粒子在第l次迭代中的更新速度,rand1(·)和rand2(·)表示两个介于[0,1]的随机数;κ2和κ3为两个控制更新步长的常数;xj,l表示第l次迭代的第j个粒子;表示第l次迭代的第j个粒子的个体最优位置。
基于更新速度,更新不同的粒子:
xj,l+1=xj,l+κ1·vj,l+1 (42)
其中,κ1为用以控制更新的步长的常数;xj,l+1表示第l+1次迭代的第j个粒子。若迭代次数未达到最大迭代次数,则开始下一次反演迭代,重复上述内容;若达到最大迭代次数,则结束非线性反演,输出所有粒子,选取其中最好结果作为最终的非线性反演结果。最好结果就是粒子所对应的适应度值最低的结果。
图2为基于VTI反射系数近似式的线性直接反演所获得五参数结果,包括纵波速度α、横波速度β、密度ρ、各向异性参数ε和各向异性参数δ;其中,黑色实线为真实测井模型,黑色虚线为初始模型,灰色点线为反演结果;可以看出近似式直接反演获得结果并不理想。图3是本发明混合反演算法线性反演部分获得结果,同样包括有α、β、ρ、ε和δ;其中,黑色实线为真实测井模型,黑色虚线为初始模型,灰色点线为反演结果;可以看出该结果比图2中近似式结果有很大提升,但是与实际模型相比仍有很大误差。图4为利用线性反演结果生成的搜索范围;其中,黑色实线为真实测井模型,灰色点线为线性反演结果,灰色虚线为生成的搜索范围;可见由于线性结果大程度贴近真实值,基于此可以生成较小范围的搜索区域,从而提升后续非线性反演的效率与精度。图5为本发明混合反演方法的非线性反演最终结果;同样地,黑色实线为真实模型,灰色点线为反演结果;相比于图3中线性反演结果,非线性结果有很大的提升,与真实模型之间的一致性更好。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围内。
Claims (10)
1.一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,其特征在于,包括以下步骤:
步骤一:获取PP波叠前角度道集,并从中提取地震子波,并建立用于反演的子波矩阵;
步骤二:获取工区的井孔数据和三维层位解释结果,建立反演所用的初始模型;
步骤三:预设线性反演循环迭代的最大迭代次数,开始线性反演循环迭代,利用初始模型和子波矩阵获取合成地震记录;
步骤四:利用初始模型计算反射系数精确式正演算子的偏导数;
步骤五:建立目标函数,基于正演算子的偏导数和合成地震记录,计算目标函数的雅可比矩阵和伪海森矩阵;
步骤六:利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量的更新方向,利用迭代格式更新目标参数向量;
步骤七:若线性反演循环迭代次数未达到预设的线性反演循环迭代的最大迭代次数,则开始下一次线性反演迭代,重复步骤三至六;若达到,则结束线性反演循环迭代过程,输出线性反演结果,继续执行步骤八;
步骤八:利用线性反演结果,生成搜索范围和初始粒子;
步骤十:将目标函数应用于所有粒子来计算替换适应度值,并计算不同粒子的接收概率,利用轮盘赌获取群体最优位置xgbest;
步骤十一:根据个体最优位置和群体最优位置来计算更新速度,进而采用迭代格式更新粒子;
步骤十二:若迭代次数未达到预设的非线性反演迭代的最大迭代次数,则开始下一次非线性反演循环迭代,重复步骤九至十一;若达到最大迭代次数,则结束非线性反演,输出所有粒子,选取所有粒子中的最好结果作为最终的非线性反演结果。
3.根据权利要求2所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,其特征在于,步骤三中,获取合成地震记录具体如下:
利用VTI介质反射系数精确方程计算反射系数r:
A·r=b
其中,A和b均为中间过程矩阵,具有如下表达形式:
其中,为上层介质纵波相关的中间变量,为上层介质横波相关的中间变量,为下层介质纵波相关的中间变量,为下层介质横波相关的中间变量; 为上层介质的刚度参数,为下层介质的刚度参数;分别为上、下层纵波相关的垂直慢度,分别为上、下层横波相关的垂直慢度,p为水平慢度; p这些参数均由刚度参数c11、c13、c33、c55和密度ρ计算获得,上标T表示转置运算符;
基于褶积理论,合成地震记录dsyn利用子波矩阵W和反射系数r计算求得:
dsyn=G(m)=W·r
其中,m为设定的目标参数向量:
m=(α1,...,αM,β1,...,βM,ρ1,...,ρM,ε1,...,εM,δ1,...,δM)T
αi、βi、ρi、εi和δi分别表示为第i个时间采样点所对应的纵波速度、横波速度、密度、两个各向异性参数,i=1,....,M;M为模型时间深度最大样点个数;G(m)为反射系数精确式正演算子,即合成地震记录。
5.根据权利要求4所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,其特征在于,步骤五中,建立目标函数S(m),表示为:
其中,d为实际叠前地震道集数据,λ为正则化参数,Cm为m计算得到的协方差矩阵,u为m的期望;
计算目标函数S(m)的雅可比矩阵J(m),即一阶偏导数矩阵:
针对目标函数S(m)的伪海森矩阵H(m),采用了内部小循环迭代形式来计算获得:
对于线性反演循环迭代中第k次迭代、内部小循环首次迭代的伪海森矩阵Hk,1(m),由下式计算得到:
对于线性反演循环迭代中第k次迭代、内部小循环第i次迭代的为海森矩阵Hk,i(m),由下式计算获得:
Hk,i(m)-1=(Qk,i-1)THk,i-1(m)-1(Qk,i-1)+Kk,i-1sk,i-1(sk,i-1)T
其中,Qk,i-1、Kk,i-1、sk,i-1为线性反演循环中第k次迭代、内部小循环第i-1次迭代的中间过程矩阵,Hk,i-1(m)为线性反演循环中第k次迭代内部小循环第i-1次迭代的伪海森矩阵。
6.根据权利要求5所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,其特征在于,步骤六中,利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量更新方向,利用迭代格式更新目标参数向量,具体形式如下:
线性反演循环迭代的第k次迭代时,目标参数向量的更新方向-gk(m)由下式计算得到:
-gk(m)=-Hk(m)·Jk(m)
其中,Hk(m)和Jk(m)分别代表线性反演循环中第k次迭代所对应的伪海森矩阵和雅可比矩阵;
利用如下的更新迭代格式来更新mk+1:
9.根据权利要求8所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,其特征在于,步骤十中,将目标函数应用于所有粒子来计算替换适应度值,并计算不同粒子的接收概率,利用轮盘赌获取群体最优位置xgbest包括:
计算替换适用度值Zr(xj,xgbest,Tl):
其中,xj表示第j个粒子,xgbest为群体最优位置,S(xj)和S(xgbest)分别表示将xj和xgbest代入目标函数计算得到的适用度,为xj和xgbest适用度的差;l为非线性反演循环中的迭代次数,Tl为第l次迭代所对应的温度值,τ表示为一个常数;
计算好所有粒子的替换适用度值后,计算不同粒子的接收概率Pg:
其中,∑为求和运算符;
在得到所有粒子的接收概率后,采用轮盘赌的方式获得群体最优位置xgbest。
10.根据权利要求9所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法,其特征在于,步骤十一中,根据个体最优位置和群体最优位置来计算更新速度,进而采用迭代格式更新粒子包括:
其中,vj,l+1、vj,l表示第j个粒子在第l+1次、第l次迭代中的更新速度,rand1(·)和rand2(·)表示两个介于[0,1]的随机数;κ2和κ3为两个控制更新步长的常数;xj,l表示第l次迭代的第j个粒子;表示第l次迭代的第j个粒子的个体最优位置;
基于更新速度,更新不同的粒子:
xj,l+1=xj,l+κ1·vj,l+1
其中,κ1为用以控制更新的步长的常数;xj,l+1表示第l+1次迭代的第j个粒子。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110218907.5A CN113031058B (zh) | 2021-02-26 | 2021-02-26 | 一种基于反射系数精确式的页岩vti储层叠前混合反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110218907.5A CN113031058B (zh) | 2021-02-26 | 2021-02-26 | 一种基于反射系数精确式的页岩vti储层叠前混合反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113031058A true CN113031058A (zh) | 2021-06-25 |
CN113031058B CN113031058B (zh) | 2022-06-24 |
Family
ID=76462400
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110218907.5A Active CN113031058B (zh) | 2021-02-26 | 2021-02-26 | 一种基于反射系数精确式的页岩vti储层叠前混合反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113031058B (zh) |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102313903A (zh) * | 2011-07-01 | 2012-01-11 | 中国海洋石油总公司 | 基于波动方程外推算子的vti介质中叠前时间偏移方法 |
CN108254780A (zh) * | 2018-01-22 | 2018-07-06 | 河海大学 | 一种微地震定位及各向异性速度结构层析成像方法 |
CN109143346A (zh) * | 2017-06-19 | 2019-01-04 | 中国石油化工股份有限公司 | 叠前混合非线性反演方法及计算机可读存储介质 |
CN111025387A (zh) * | 2019-12-19 | 2020-04-17 | 河海大学 | 一种页岩储层的叠前地震多参数反演方法 |
CN111175821A (zh) * | 2020-01-17 | 2020-05-19 | 河海大学 | 一种vti介质的各向异性参数分步反演方法 |
CN111308550A (zh) * | 2020-03-16 | 2020-06-19 | 河海大学 | 一种页岩vti储层的各向异性参数多波联合反演方法 |
CN111913218A (zh) * | 2019-05-09 | 2020-11-10 | 中国石油天然气股份有限公司 | 一种基于多尺度混合反演的薄层反演方法、设备以及系统 |
-
2021
- 2021-02-26 CN CN202110218907.5A patent/CN113031058B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102313903A (zh) * | 2011-07-01 | 2012-01-11 | 中国海洋石油总公司 | 基于波动方程外推算子的vti介质中叠前时间偏移方法 |
CN109143346A (zh) * | 2017-06-19 | 2019-01-04 | 中国石油化工股份有限公司 | 叠前混合非线性反演方法及计算机可读存储介质 |
CN108254780A (zh) * | 2018-01-22 | 2018-07-06 | 河海大学 | 一种微地震定位及各向异性速度结构层析成像方法 |
CN111913218A (zh) * | 2019-05-09 | 2020-11-10 | 中国石油天然气股份有限公司 | 一种基于多尺度混合反演的薄层反演方法、设备以及系统 |
CN111025387A (zh) * | 2019-12-19 | 2020-04-17 | 河海大学 | 一种页岩储层的叠前地震多参数反演方法 |
CN111175821A (zh) * | 2020-01-17 | 2020-05-19 | 河海大学 | 一种vti介质的各向异性参数分步反演方法 |
CN111308550A (zh) * | 2020-03-16 | 2020-06-19 | 河海大学 | 一种页岩vti储层的各向异性参数多波联合反演方法 |
Non-Patent Citations (4)
Title |
---|
CONG LUO 等: "Joint PP and PS pre-stack AVA inversion for VTI medium based on the exact Graebner equation", 《JOURNAL OF PETROLEUM SCIENCE AND ENGINEERING》 * |
刘斌 等: "基于不等式约束的三维电阻率探测混合反演方法", 《地球物理学报》 * |
方中于 等: "基于混合智能优化算法的非线性AVO反演", 《石油地球物理勘探》 * |
李术才 等: "多同性源阵列电阻率法隧道超前探测方法与物理模拟试验研究", 《地球物理学报》 * |
Also Published As
Publication number | Publication date |
---|---|
CN113031058B (zh) | 2022-06-24 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN112083482B (zh) | 基于模型驱动深度学习的地震超分辨反演方法 | |
CN108802812B (zh) | 一种井震融合的地层岩性反演方法 | |
CN111025387B (zh) | 一种页岩储层的叠前地震多参数反演方法 | |
WO2015014762A2 (en) | Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media | |
CN110456417B (zh) | 一种地震数据多次波压制方法 | |
CN110780351B (zh) | 纵波和转换波叠前联合反演方法及系统 | |
CN110031895B (zh) | 一种基于图像缝合的多点地质统计学随机反演方法及装置 | |
CN105319589A (zh) | 一种利用局部同相轴斜率的全自动立体层析反演方法 | |
CN112394414B (zh) | 两步法地震绕射波场叠前分离的方法 | |
CN111368247A (zh) | 基于快速正交字典的稀疏表征正则化叠前avo反演方法 | |
CN110687597B (zh) | 一种基于联合字典的波阻抗反演方法 | |
CN113077386A (zh) | 基于字典学习和稀疏表征的地震资料高分辨率处理方法 | |
CN115598697A (zh) | 薄层结构高分辨率地震反演方法、装置、介质和设备 | |
CN113156500B (zh) | 数据驱动的快速构造约束叠前地震多道反演方法 | |
CN110737018A (zh) | Vsp地震资料各向异性建模方法 | |
CN111175821B (zh) | 一种vti介质的各向异性参数分步反演方法 | |
CN113031058B (zh) | 一种基于反射系数精确式的页岩vti储层叠前混合反演方法 | |
WO2021127382A1 (en) | Full waveform inversion in the midpoint-offset domain | |
CN113253350B (zh) | 一种基于联合字典的孔隙度反演方法 | |
CN111308550B (zh) | 一种页岩vti储层的各向异性参数多波联合反演方法 | |
Sun et al. | Seismic AVO inversion method for viscoelastic media based on a tandem invertible neural network model | |
Li et al. | Seismic full-waveform inversion based on superresolution for high-precision prediction of reservoir parameters | |
CN113721293B (zh) | 一种基于深度学习的多波地震信号人工智能匹配方法 | |
CN110398773B (zh) | 一种针对部分缺失的地震数据的恢复重构方法 | |
Sun et al. | Convolution model theory-based intelligent AVO inversion method for VTI media |
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 |