基于快速反射率法的AVA反演技术在致密砂岩薄储层勘探中的应用

2023-02-09刘俊州吴高奎

石油物探 2023年1期

杨震,刘俊州,时磊,韩磊,吴高奎

(1.页岩油气富集机理与有效开发国家重点实验室,北京100083;2.中国石化弹性波理论与探测技术重点实验室,北京100083;3.中国石油化工股份有限公司石油勘探开发研究院,北京102206)

我国致密砂岩气地质资源量为17.4×1012～25.1×1012m3,可采储量为8.8×1012～12.1×1012m3[1],目前已形成四川盆地须家河组和鄂尔多斯盆地上古生界等2个致密砂岩气大气区。自2009年以来,四川盆地须家河组天然气勘探潜力逐渐引起关注[2-5],该致密砂岩储层成层性明显且具有储层厚度薄的特点,优质储层的单层厚度在4～9m。使得基于精确Zoeppritz方程和相应近似正演模拟不再适用于致密砂岩薄储层。这是因为Zoeppritz方程只适用于单一反射界面,对于层状模型,尤其是含有薄互层的储层,地震记录中的多次波、转换波和透射损失等更为复杂,Zoeppritz方程难以适用。

为解决上述问题,需要能模拟薄储层波场响应的正演模拟算法。反射率法使用了波动方程的解析解,可模拟层内的复杂波场响应,包括反射波、透射波、转换波和多次波。MALLICK等[6]和KENNETT[7-8]指出,虽然反射率法基于一维模型假设,但由于采用了波动方程的解析解,因此该类方法比有限差分法等正演方法有更高的模拟精度且计算速度更快。FUCHS[9]首次提出将反射率法用于合成地震记录,该方法可模拟层内的复杂地震波场响应,包括反射波、透射波、转换波和多次波。在此基础上,KENNETT等[7-8,10-11]首次提出了用于计算层状地层的总反射和透射系数的递归算法,该方法的特点是对所有频率和慢度无条件稳定。为提高计算速率,PHINNEY等[12]在KENNETT研究成果的基础上提出了快速反射率法,即将原有的反射率法进行重组,实现矢量化,使计算速度提高了约20倍。YANG等[13]基于反射率法推导了针对单薄层的二阶近似公式,并利用递归算法将其推广到薄互层。

在地震叠前反演方面,国内外学者已经尝试将反射率法用于薄储层AVA反演。SEN等[14]比较了KENNETT的反射率法与快速反射率法,并建立基于Gauss-Newton的反演方法用于叠前波形反演。王建花[15]在层参数等效的情况下,推导了P波入射到单薄层时的波场能量再分配的简化公式,并进行了初步的P波反演试验。在这些研究中,反演均是在截距时间和射线参数域(即τ-p域)进行的。然而,当可用偏移距数据较少时,τ-p域变换会受到混叠现象的影响。为了避免混叠现象,MALLICK等[16]建议可生成角道集与实际道集数据相匹配。LIU等[17]针对P波AVA反演,提出基于快速反射率法的Bayesian反演方法。YANG等[18]基于快速反射率法对PP波进行反演,并在反演中考虑到稀疏约束。LIU等[19]提出基于快速反射率法的遗传算法,用于AVA叠前联合反演。尽管已有多位学者发表相关文章,但受到公式本身复杂性的制约,基于反射率法的叠前反演并没有被广泛应用,仍需要进一步的研究。

本文将基于快速反射率法的AVA反演应用于致密砂岩储层。在快速反射率法理论的基础上,利用Gauss-Newton法建立目标函数。为了估计弹性参数(纵波和横波速度以及密度),目标函数是通过角度域中模拟数据和观测数据之间的最小差异实现。最后,分别使用薄互层模型和实际数据进行基于反射率法的反演,并与基于精确Zoeppritz方程的反演方法进行比较,结果表明,基于快速反射率法的反演结果与真实值更加接近。

1 技术方法

1.1 反射系数的计算

KENNETT[7-8]推导出了反射率法的迭代公式,可以计算目标区域的总反射系数及总透射系数,包括多次波和转换波等。为提高运算速度,PHINNEY等[12]对KENNETT[7-8]反射率法进行了改进,将最内层的循环改为计算频率,并拓展到整个频率域,实现矢量化,不仅加快了计算速度,还能实现对弹性参数偏导数的解析计算,提高了后续反演工作的精度和效率。据PHINNEY等[12]的研究成果,频率域PP波的总反射系数公式可以写为:

(1)

式中:p为慢度;ω为角频率;v01、v04分别为数组v0的第1个和第4个元素;v0为PHINNEY基于Haskell-Dunkin minor matrix方法在慢度-频率域(p,ω)定义的含有6个元素的数组,则有:

v0=

(2)

式中:Δ为比例系数;Rij(i=P或S,j=P或S)为对应波类型的总反射系数;detR为反射系数行列式。

v0可以通过递归法计算得到:

v0=Q0Q1…Qn…QN-1vN

(3)

其中,vN为初始项。

(4)

Qn为地层n的传播矩阵。

Qn=EnFn

(5)

式中:En为层穿越矩阵,描述了平面波在第n层地层中传播时相位的移动;Fn为界面矩阵。En和Fn可以具体写为:

(6)

(7)

式中:qP,qS分别为纵波、横波的垂直慢度;Tn源自Dunkin矩阵,是第n层大小为6×6的矩阵,与频率无关,依赖于第n层的水平慢度p和弹性参数,具体形式可以参看PHINNEY等[12]的文献。

1.2 反演原理

本文的反演基于Gauss-Newton法,其核心思想是利用快速反射率法计算层状地层PP波反射系数,最大程度地模拟目标地震记录,以此来估算目标区域的弹性参数(纵波速度α,横波速度β和密度ρ)。为得到合适的模型参数,使得模型响应和观测数据的差异最小,反演的目标函数定义为:

Q(M)=‖Φ-Φobs‖2

(8)

式中:M=(α,β,ρ)为模型弹性参数向量;Φobs为实际地震记录;Φ为模拟地震记录。

需要注意的是,这里使用的Φobs和Φ均为时间域角道集,而本文所述的快速反射率法在慢度-频率域计算反射系数,因此需要进行域的转换。由慢度-频率域转换到时间-角度域的过程可以分成以下两个步骤:

1) 慢度转换为角度。在反射系数R(p,ω)的计算中有3个循环,即最外层循环是慢度p,中间循环为层循环(n)和最内层循环为频率循环。一般来说,慢度在每一个最外层循环上是固定的,入射角在每一层上都会发生变化,即θn=arcsin(αnp)。但将入射角设置在最外层循环,每层的慢度都会发生变化,因此,通过公式(9),R(p,ω)被重采样为R(θ,ω),最后得到角度道集。

(9)

式中:下角标n表示第n层。

2) 频率域转到时间域。类比FRYER[20],KENNETT[7-8]的处理方法,对反射系数做一次逆Fourier变换,最终得到截距时间-角度域(τ-θ)的反射系数为:

(10)

得到时间-角度域的反射系数后,地震记录Φ可以通过反射系数R与子波W的褶积得到。

Φ=W*R

(11)

给定符合假设条件的初始模型M0(能反映出地层的大致形态),利用Gauss-Newton法[21-24],模型的更新参数公式可以写为:

ΔM=[JTJ]-1JT(Φ-Φobs)

(12)

式中:J为雅可比矩阵。

(13)

通过迭代ΔM更新初始弹性参数向量M0直到目标函数Q足够小,得到反演结果。公式(12)的难点在于其不稳定性,这里应用Levenberg-Marquardt方法修改公式(12),得到:

ΔM=[H+λI]-1JT(Φ-Φobs)

(14)

式中:H=JTJ,为海森矩阵;λ为比例系数;I为单位矩阵。

根据公式(10)和公式(13),雅可比矩阵J具体可以写为:

(15)

其中,

(16)

其中,v0对Mn的偏导数为:

(17)

其中,

(18)

(19)

式中:Mn代表第n层的α,β,ρ。

对于Fn的偏导数具体可以写为:

(20)

其中,

(21)

En对M的偏导数可以求出解析解,具体可以写成以下形式:

(22)

式中:A=qP-p2/qS,B=qP+p2/qS。

(23)

(24)

在E拓展后,对应的v0也需要做同样的拓展,拓展方法可参考PHINNEY[12]。

2 模型测试

2.1 反射系数对比

鉴于实际地层中薄层不能单独存在的原因,为验证快速反射率法在AVA建模和反演中的有效性,采用多层模型测试其中薄层的AVA曲线特点,模型参数如表1所示。PP波合成角道集由快速反射率法和精确Zoeppritz方程分别进行模拟,子波采用40Hz的Ricker子波。与基于精确Zoeppritz方程的AVA道集相比(图1),基于快速反射率法的AVA道集包含更完整的信息,例如图1c中所示的多次波。这是因为使用精确Zoeppritz方程的合成地震记录没考虑多次波,可能会导致解释和反演中出现错误。反射率法考虑了波在地层中传播时的传播效应,可以计算传播过程中所有的反射系数,由于现有的地震资料处理技术难以在不影响初级反射振幅的情况下去除波在地层中的传播效应,如透射损失、层间多次波和几何扩散等,因此,采用快速反射率法可以更逼真地模拟地震记录。

为了进一步显示快速反射率法和精确Zoeppritz方程方法之间的差异,沿图1中黄线位置提取瞬时振幅随入射角的变化数值,如图2所示。受到传播效应的影响,使用快速反射率法计算两个界面的振幅随深度的增加而逐渐减小,比使用Zoeppritz方程模拟的振幅小,主要是由于该方法考虑了波在传播时的透射损失。而基于Zoeppritz方程计算的AVA曲线不会随深度的变化而变化。如果在地震数据处理程序中未正确校正传输损耗和层间多次波,则不建议基于Zoeppritz进行正演模拟,利用快速反射率法模拟的道集将与现场数据更匹配。

表1 等厚模型参数

图1 PP波AVA道集对比(红色曲线代表P波速度)a 基于快速反射率法的AVA道集; b 基于精确Zoeppritz方程的AVA道集; c 两种方法模拟的AVA道集之间的差异

图2 AVA曲线对比(Phinney 1和Phinney 2是图1a中的两个界面,Zoeppritz 1和Zoeppritz 2是图1b中的两个界面)

2.2 薄互层模型反演

将薄互层模型的单个薄层的厚度设置为8m,具体参数如表2所示。利用快速反射率法计算反射系数时,频率的计算范围为0～125Hz,合成地震记录所用子波为40Hz Ricker子波,采样率为1ms。考虑到实际数据的角道集数据入射角范围通常较窄,该薄互层模型角道集的角度范围是5°～30°,如图3a所示。图3b为对地震记录添加15%随机噪声后的地震记录。此外,该模型为深度域模型,后续反演得到的模型更新量是时间域模型,需将深度域模型转换到时间域(PP波时间域),各个界面对应于表2中的双程旅行时。

表2 薄互层模型参数

图3 薄互层模型PP波角道集记录(入射角为5°～30°,红色曲线为纵波速度)a 不加噪声的地震记录; b 添加15%噪声的地震记录

以平滑后的真实模型为初始模型,反演的目标数据为图3所示地震记录。分别采用基于快速反射率法和精确Zoeppritz的反演方法进行反演,反演结果如图4所示。由反演结果可以看出,基于快速反射率法的反演结果基本与真实值重合。而基于Zoeppritz方程的反演结果会有一些波动,与地震记录中的多次波有关。图5为反演迭代次数收敛曲线,红色为基于快速反射率法的收敛曲线;绿色为基于Zoeppritz方程的收敛曲线。若是将收敛精度定义在剩余误差绝对值小于1时反演收敛,那么基于快速反射率法的反演在迭代3次后收敛,基于Zoeppritz方程的反演在迭代4次后收敛。

图4 基于快速反射率法和精确Zoeppritz方程的反演结果

与之前的反演类似,均以平滑模型为初始模型,图6为两种反演方法的反演结果。受地震记录中随机噪声的影响,两种反演方法的反演结果中均含有随机噪声。但是在目标层位置处,特别是单独的PP波反演,反演结果仍然能较好地反演地层的真实参数。同时,基于快速反射系数法的反演结果更逼近于真实值,说明该方法具有更强的抗噪性。

图5 不添加噪声的薄互层模型AVA反演迭代次数收敛曲线

图6 基于快速反射率法和精确Zoeppritz方程的反演结果(对模型地震记录添加15%随机噪声后)

需要强调的是,这里的模拟地震记录包含层间多次波,并考虑了透射损失,因此将该方法应用到实际数据的时候,实际数据含透射损失和多次波。此外,本反演方法没有反演地层厚度参数,未讨论厚度对反演结果的影响,这要求使用的初始模型能大体反应地层厚度的变化。

3 应用实例

实际数据选自四川某工区,其目的层段为须家河组的须三段。由于龙门山造山带北段和西部山系的持续隆升,造成须三段沉积时期盆地边缘大量沉积物快速堆积,物源以近南北向为主,在盆地中央深处发育湖相沉积,目的层岩性组合为致密砂岩、泥岩互层,埋深4000～5000m。

在应用本文反演方法之前,首先,要确保数据已经被合理地处理过。由于反射率法在计算反射系数时会考虑透射损失以及多次波,因此地震数据处理时,可以不用考虑透射损失以及多次波的处理。其次,利用井数据做正演模拟,然后对地震数据进行标定,这个过程是将井数据(深度域)与地震数据(时间域)进行对应,对于该数据,标定是基于商业软件完成,所以标定的过程是在Zoeppritz方程的基础上完成的。标定完成后,模拟的地震记录与实际数据之间的相似度约0.78。基于Zoeppritz和基于反射率法的合成地震记录与实际数据的对比如图7所示,为了更细致地观察二者的区别,在图7中选3个时间点(T1,T2,T3)提取AVO曲线,AVO曲线对比如图8所示。由图8可以看出,无论是数据振幅的幅值,还是AVO曲线的走向趋势,均是基于反射率法的地震记录与实际地震记录更加接近。最后是初始模型的选择,将实际测井数据限频后的曲线作为初始模型,频率的选取范围为0～60Hz。需要注意的是,反射率法的基本假设条件是数据为局部一维数据或是水平层状。结合叠前时间偏移,该方法可以适用于地质条件简单的地区。在地质条件复杂的地区,需要进行叠前时间偏移或逆时偏移。本文数据因经过叠前偏移处理,所以是水平地层的剖面。

图7 井旁地震道及不同方法合成地震记录(红色曲线为纵波速度测井曲线)a 实际地震数据; b 基于快速反射率法的合成地震记录; c 基于Zoeppritz方程的合成地震记录

图8 图7中T1,T2,T3位置处的AVO曲线对比(散点代表振幅数据点;曲线代表回归曲线)a T1时间; b T2时间; c T3时间

为了验证基于反射率法反演的优越性,分别基于精确Zoeppritz方程和反射率法对实际资料进行AVA反演。反演的目标区是2150～2350ms,井旁角道集数据如图9所示。共有6个角道集,范围为3°～28°。该区域主要的地层特征为砂泥岩互层,图中绿色方框所圈出区域为目标层段。图10为井旁道地震数据反演结果,由反演结果可以看出,基于反射率法的反演结果与测井曲线更加贴合,能显示出更多的细节。值得注意的是,相对于纵波速度反演结果,横波速度和密度的反演结果相对较差,这是因为反射系数公式本身对纵波速度更加敏感[25]。图11和图12 分别为过井的基于快速反射率法和基于精确Zoeppritz方程的反演剖面。反演中所使用的子波需要对每个入射角分别提取。尽管受到地震分辨率的限制,反演结果与测井相比频率较低,但从反演结果可以看出,基于反射率法的反演结果要更接近于测井曲线,能更准确地反映地层变化的趋势,对识别有效的致密砂岩储层非常有帮助。