HTI介质地震波各向异性AVO反演

李 勤,王 玮,王瀚林

(西安科技大学 地质与环境学院,陕西 西安 710054)

摘 要:常规AVO反演主要针对各向同性介质进行,含有裂隙的地层实际上是各向异性的。地震反演问题是典型的不适定问题,各向异性参数的加入,使得地震反演问题的非线性程度显著增加,比传统反演难度加大。SA-PSO算法将粒子群算法与模拟退火法进行混合优化,具有全局寻优能力更强、克服早熟现象、不易陷入局部极小值等优点。将含有垂向裂隙的介质近似为HTI介质,基于HTI介质近似反射系数方程,建立关于各向异性参数的目标函数,通过SA-PSO算法实现HTI介质各向异性参数反演,进一步,利用流体因子、裂隙弱度与各向异性参数的关联,对裂隙相关参数及流体因子进行预测。建立含有HTI型煤层的理论模型,对原始数据同时使用SA-PSO算法与PSO算法进行反演,对算法的优越性进行测试,SA-PSO算法反演精度更高,反演过程更稳定;对多层理论模型数据分别加上信噪比为5,10,20的高斯随机噪声,对反演方法的抗噪性进行测试,反演结果误差较小,反演过程能快速收敛,反演抗噪性较好;根据HTI型含水裂隙中各向异性参数的特点,对Marmousi2模型的含水砂岩部分进行改造,引入各向异性参数,通过反演获取速度、密度及各向异性参数剖面,并实现流体因子预测,反演及预测结果与理论模型吻合,论证了方法的有效性。

关键词:HTI介质;各向异性;AVO反演;流体因子;裂隙

HTI(Transverse Isotropy with a Horizontal axis of symmetry)介质是一种具有近似水平对称轴的横向各向同性介质,对于含有垂直裂隙的模型,可将其近似等效为HTI介质模型[1]。HTI介质是一种典型的各向异性模型,THOMSEN[2]提出通过各向异性参数来描述介质的各向异性程度,并推导了各向异性参数与弹性参数之间的量化关系;董守华[3]针对煤层进行了各向异性参数测试,讨论煤层各向异性参数与孔隙率的关联,并指出通过煤层各向异性的大小和方向能够预测煤层的裂隙密度和方位;陈同俊等[4]利用介质等效模型理论分析了构造煤的各向异性AVO(Amplitude Variation with Offset)特征,对不同的AVO近似方程应用于煤层的效果进行比较,对AVO近似方程进行改进,使其能适用于各向异性煤层。彭苏萍和毕银丽[5]指出煤层在开采过程中会产生裂隙。许多专家研究了各向异性参数与裂隙之间的关联,SCHOENBERG和DOUMA[6]针对含平行裂隙的介质,将其等效为无限薄无限柔顺的薄层,忽略裂隙的形状及其空间分布细节,假设裂隙面处位移张量不连续,但牵引力张量连续,且2者线性相关,提出了线性滑动理论,建立了裂隙介质等效模型;BAKULIN等[7]提出用裂隙弱度来描述裂隙,分析裂隙流体因子与裂隙充填物的关联;QUINTAL[8]研究了饱和流体研究中流体饱和度对反射系数的影响;陈怀震等[9]通过分析流体因子对不同流体的响应关系,指出HTI介质中裂隙饱和水时,裂隙流体因子值偏小,趋近于0;郎玉泉等[10]利用Gassmann方程进行流体替代,通过AVO正演,估计顶板砂岩孔隙度和干湿性。AVO反演的基础是反射系数近似式,RÜGER[11]首先给出不同的横向各向同性介质的纵波反射系数近似公式,近似式精度较高,是目前AVO理论的基础[12-13]。利用反射系数近似公式,可以实现地震叠前AVO反演[14-16]。MA[17]利用模拟退火法,通过叠前AVO反演得到了岩石特性;MISRA和SACCHI[18]利用快速模拟退火法解决了基于边界保护的AVO反演问题,粒子群算法最早由Eberhart提出,源于对鸟群捕食行为的拟态[19],MIKKI[20]将这种方法用于电磁优化,使该算法进入工程物探领域并得到发展,SUN等[21]基于粒子群算法进行储层预测,收到一定成效;SUN和LIU[22]利用混沌量子粒子群算法求解了叠前AVO弹性参数;MEHRGINI等[23]利用粒子群算法对具有中子孔隙度、体积密度、水饱和度时的地层横波速度VS进行了反演并论证;WU等[24]利用粒子群算法对地层的纵横波速度、密度进行了反演;严哲等[25]利用量子粒子群算法对Marmousi2模型的纵横波速度和密度进行了反演,反演结果稳定。对反演算法进行改进和优化是提高反演精度和反演效率的主要手段,对于叠前AVO反演来说,粒子群算法(PSO,Particle Swarm Optimization)虽然收敛速度快、原理简单,但却容易早熟,陷入局部最优;模拟退火法(SA,Simulated Annealing)虽然全局寻优能力强但其稳定性又不能得到保证。因此,在粒子群算法的基础上结合模拟退火的思想,对粒子群算法进行改进,使新算法有强行跳出局部最优解的能力,这种优化后的混合算法在反演时具有明显的优越性和有效性。笔者基于HTI介质反射系数近似公式,建立关于各向异性参数的目标函数,采用基于模拟退火的优化粒子群算法(SA-PSO)实现HTI介质叠前AVO反演,并通过裂隙流体因子与各向异性参数的关联,实现流体因子预测。

1 HTI介质各向异性AVO反演理论及步骤

1.1 基于模拟退火的优化粒子群算法反演原理

粒子群算法是通过模拟鸟类觅食的过程,达到解决优化问题的方法,算法收敛速度快。在粒子群算法中,鸟类被抽象为没有质量和体积的微粒,并延伸到n维空间,粒子in维空间的位置表示为向量Xi=(x1,x2,x3,…,xn),飞行速度表示为向量Vi=(v1,v2,v3,…,vn),粒子不仅知道自己发现的最好位置Pb,也知道整个群体中的其他粒子发现的最优位置Pb,g,粒子们根据这些信息来决定接下来的行动方向。根据目标函数的构建,每个粒子都有其适应值范围。粒子群优化算法首先将粒子初始化为一群随机粒子,然后通过迭代找到最优解,在每一次迭代中,粒子通过跟踪两个“最优解”(即Pb,Pb,g)来更新自己,在找到“最优解”(Pb,Pb,g)后,粒子通过下式更新自己的位置:

(1)

式中,Vi(t)为粒子it时刻的速度;ω为惯性因子;c1c2为学习因子,分别为粒子的局部学习能力和全局学习能力;r1r2为[0,1]的随机数;Pb,i(t),Pb,gi(t)分别为粒子i和群体中其他粒子在t时刻为止发现的最好位置;xi(t)为粒子it时刻的位置。

在搜索过程中,粒子的寻优范围局限在Pb附近,因此不能保证全局收敛得到全局最优解,而惯性因子ω表示的就是继承先前粒子速度的能力,当ω较大时,继承粒子前一时刻速度较多,全局寻优能力会有所提高。ω的值如果过大或过小就会造成算法陷入局部寻优或者无法找出最优解。因此,对惯性因子进行实时调整很有必要。针对ω在算法运行各个阶段的不同作用,加入随迭代次数非线性减小的惯性因子。减小的惯性因子公式为

ω(k)=ωs-(ωs-ωe)(k/Tmax)2

(2)

式中,ωs为初始惯性因子;ωe为迭代至最大次数时的惯性因子;k为当前迭代次;Tmax为最大迭代次数。

学习因子c1c2影响算法的全局和局部搜索能力,c1为粒子的个体认知,因此影响局部搜索能力;c2为群体认知,影响全局搜索能力。在粒子群优化算法迭代过程中,迭代前期需要粒子群优化算法的全局搜索能力强,随着迭代进行至后期,算法需要快速收敛,此时需要局部搜索能力强。因此,在粒子群优化算法加入自适应的学习因子使得算法找到更优解乃至最优解。学习因子自适应公式为

(3)

式中,R1,R2,R3,R4为学习因子的自适应度参数。

在粒子群算法中,粒子在发现的当前最优位置靠近局部最优解的时候,所有的粒子就会聚集到最小的区域,全局搜索能力就会削弱。而模拟退火法的优势正在于通过调整温度,控制概率性跳出特性,进行新一次的随机搜索,将模拟退火的这种特性引入到粒子群算法,可以有效避免陷入局部最优解,这种修正主要是针对全局最优位置Pb,g,利用模拟退火的跳出机制来选择每个粒子的当前最优位置,即通过控制模拟退火法对“变差解”的接受概率,可以提高粒子群算法的灵活性,增加粒子的多样性。模拟退火法通过在求解区间随机游走,利用Metropolis抽样准则使随机搜索最新位置逐渐收敛于最优解。Metropolis准则定义了某一温度下系统状态从状态i到状态j 的能量概率为

(4)

式中,为能量概率;T为温度;EkEj为固体在状态k和状态j下的能量;K为玻尔兹曼常数。

在模拟退火的跳出机制下,认为Pb是和Pb,g悬殊很大的特殊解,可以计算出温度为TPb相对Pb,g 的能量概率,即exp(-fPb-fPb,g)T,其中f为目标函数解值。若将所求概率作为Pb的适配值,则用Pb替代Pb,g

1.2 HTI介质反射系数与裂隙流体因子

HTI介质通常是由平行排列的垂直裂隙形成,即对称轴接近于水平方向时的TI介质。HTI介质弹性矩阵:

(5)

由于HTI介质在垂直面弹性参数相同,根据矩阵对称性,其弹性参数有如下关系:

(6)

根据Thomson各向异性介质理论,HTI介质的各向异性参数和弹性参数可由式(7)进行计算

(7)

其中,VP为纵波速度;ρ为密度;VS为横波速度;ε(v),δ(v),γ(v)为HTI介质的Thomsen各向异性参数,ε(v)为纵波各向异性程度,δ(v)为纵波在横向和垂向之间各向异性变化的快慢程度,γ(v)为快、慢横波速度的差异程度。

基于Thomsen各向异性参数,RÜGER给出了HTI介质的PP波反射系数公式:

(8)

式中,Δ′为上下介质参数的差值;Z=ρVP,为垂向纵波阻抗;上划线为上、下介质参数的平均值;为横波剪切模量;θ为入射角;φ为方位角;下标1,2分别为上层介质和下层介质。

根据线性滑动模型理论中应力与裂隙应变的关系,裂隙介质弹性系数矩阵C可以表示成各向同性背景系数矩阵Ciso加上各向异性扰动Cani之和,则描述TI介质的裂隙等效弹性矩阵为

C=Ciso+Cani

(9)

其中,令λ+2μ=L,则Ciso,Cani分别为

(10)

(11)

g=μ/(λ+2μ)=(VS/VP)2为各向同性岩石骨架横波速度和纵波速度比值的平方,经推导得到裂隙刚度矩阵中的ΔNΔT表达式:

(12)

式中,ΔNΔT为Schoenberg线性滑动理论中的裂隙法向弱度和切向弱度,其绝对值在0~1;e为裂隙密度;g为横纵波速比的平方;Kfμf为裂隙充填流体体积模量和剪切模量;λμ为不含裂隙岩石的拉梅系数;a/c为裂隙纵横比,a,c分别为裂隙纵向长度,横向宽度。

结合Thomsen各向异性参数和线性滑动模型等效理论,经推导可以得到HTI介质各向异性参数与裂隙柔度参数之间的表达式为

(13)

所以,ΔNΔT可由各向异性参数表示为

(14)

假设裂隙之间互不连通、互不影响,HTI型裂隙介质中流体识别可以用KN/KT来指示:

(15)

其中,KN为裂隙法向柔度;KT为裂隙切向柔度。对于液体充填的裂隙介质,剪切模量μf=0,体积模量Kfμ大致相差一个数量级,裂隙纵横比通常很小,a/c≪1,[(Kf+4/3μf)/ μ]/(a/c)≫1,此时ΔN≈0,ΔT与裂隙为干燥状态时相同,流体因子KN/KT接近于0。

1.3 HTI介质各向异性反演实现过程

采用SA-PSO算法对HTI介质进行叠前AVO反演,基本思想可概括为:将待反演模型的上、下层的弹性参数作为种群中寻找最优解的粒子,通过位置公式进行迭代计算,选用不同的惯性因子和学习因子,使目标函数求解达到全局最小。具体流程如下:

(1)将模型确立为HTI介质模型;

(2)通过地震波反射系数近似公式计算出模型的PP波反射系数;

(3)确定HTI介质模型的输入参数及待反演参数,构建目标问题的反演目标函数;

(4)将先验信息与待反演参数输入反演目标函数;

(5)输入先验地震数据信息或实测地震数据,利用SA-PSO算法寻找反演目标函数的最佳拟合;

(6)输出待反演参数的最优解,预测储层的相关流体因子。

假设待反演的各向异性参数符合粒子群分布,将反演的约束条件设为

(16)

其中,X=(δ(v),ε(v),γ(v))为待反演的各向异性参数;i为地震记录的道序号;dobs(θi)=[R(Sobs(θi)),I(Sobs(θi))],dsyn(θi)=[R(Ssyn(θi)),I(Ssyn(θi))],R和I分别为反射系数与子波褶积后复数的实部与虚部;θi为该道地震波对应的PP波入射角;Sobs(θi)

Ssyn(θi)分别为HTI介质模型在PP波入射角为θi时的PP波反射系数序列与子波序列的褶积。

利用上述方法进行反演的技术流程如图1所示。

图1 基于PSO-SA算法的反演流程
Fig.1 Schematic diagram of inversion flow based on PSO-SA algorithm

2 模型试算

2.1 多层模型试算

建立一个4层的模型,从地表往下分别为第1,2,3,4层,横向为400 m。其中,第1层为泥岩,厚度约为300 m;第2层为砂岩,含有垂向饱水裂隙,可近似为HTI型介质,层厚约为100 m;第3层为煤层,厚度约为10 m;第4层为泥岩,厚度约为300 m。模型对应的纵、横波速度体如图2(a),3(a)所示,模型弹性参数和各向异性参数见表1。根据模型中第2层为基于HTI型饱水裂隙的假设,其各向异性参数ε(v)的值无限接近于0,δ(v)γ(v) 取值在0.1~0.3。

表1 理论模型的物性参数和各向异性参数
Table 1 Physical and anisotropy parameters of the theoretical model

层号VP/(m·s-1)VS/(m·s-1)ρ/(g·cm-3)ε(v)δ(v)γ(v)13 5001 9002.6000023 0701 9802.400~0.010.10~0.230.20~0.3032 1001 1001.3500043 5001 9002.60000

注:表中纵横波速度以及密度参数来自文献[26-27]。

图2 基于SA-PSO算法的模型VP反演剖面对比
Fig.2 Comparison of model VP inversion profiles based on SA-PSO algorithm

图3 基于SA-PSO算法的多层模型VS反演剖面对比
Fig.3 Model values and inversion results of S-wave velocity based on SA-PSO in the model

对多层模型设计目标函数并进行弹性参数和各向异性参数反演。首先,为证明本研究所提出的基于SA-PSO算法的科学性和有效性,对多层模型的VP进行不同迭代次数的SA-PSO反演和单一的PSO反演进行对比,结果如图2所示。图2(a)为模型的VP,对比图2(b)可以看出,利用SA-PSO混合算法迭代次数较多(100次)时,寻优效果更稳定,得到的地层分异明显,反演速度也与理论速度误差较小。在图2(c)中,使用SA-PSO迭代50次时,此时全局寻优尚不稳定,精度不高,所以出现细小噪点;图2(d)为利用PSO反演,迭代次数达到100次时的结果,可以看出地层分异相较于图2(c)效果稍好一些,参数区间也更精确一些,但是相比图2(b)来说,反演结果不够理想。综上所述,利用SA-PSO反演时,由于SA-PSO反演克服了PSO的早熟现象,在迭代次数达到一定要求后反演精度高于单一的PSO反演,本模型中,当迭代次数达到100次时反演结果最好,得到的结果更稳定。

同时,对多层模型中的VS,ρ进行迭代次数为100次的SA-PSO反演,反演结果如图3,4所示。分析图3,4可以看出,从VPVS的反演剖面来看,模型的地层能被清晰地分辨,反演结果较好。对于ε(v)δ(v)γ(v),在初始模型中,由于除了砂岩层之外的其他地层均等效为各向同性介质,其各向异性参数均设置为0。在对砂岩层的各向异性参数进行反演时,抽取了横向距离200 m处的数据进行试算,试算结果如图5所示。在图5中,对各向异性参数的反演值与理论各向异性参数对比,可以看出,反演值始终围绕理论值波动,反演结果稳定,同时,ε(v)始终在0附近,这符合对砂岩层含有饱水裂隙的预设。多层模型的弹性参数和各向异性参数反演剖面均能指示地层,各向异性参数ε(v)出现在0附近,对裂隙流体的指示性强。

图4 基于SA-PSO算法的多层模型ρ反演剖面对比
Fig.4 Model values and inversion results of ρ based on SA-PSO in the model

针对抽取的横向距离200 m处的试算结果,对其中纵向深度为350 m、横向距离为200 m处的反演结果进行误差定量分析,见表2。由表2可知,对于多层模型的多参数反演,ε(v)VP的反演误差最小,ρδ(v)的反演误差略大。但整体来讲,反演误差都较小。

图5 基于SA-PSO算法的多层模型各向异性参数 反演结果对比
Fig.5 Model values and inversion results of anisotropic parameters based on SA-PSO in the model

同时,为讨论该反演方法的抗噪性,依照正演流程计算模型的反射系数并采用40 Hz的雷克子波与反射系数进行褶积,生成不同入射角的角道集,并分别添加了信噪比(SNR)为5,10,20的高斯随机噪声,并针对各向异性层,分析不同信噪比数据的反演误差(表3)。从表3来看,大致可以认为信噪比越高的数据反演误差越小,在反演过程中也能够越快收敛,但是在模型试算的结果中,信噪比与反演误差并不完全负相关;即使在SNR=5时,角道集仍然能够反映岩石物理模型的层位特征,在信噪比较低的地震记录中也可以读取有效信息,也就是说,信噪比较低的情况下,即使反演误差稍大,但是对VPε(v)δ(v)的反演结果仍然能满足反演要求,反演算法具有一定的抗噪性。

表2 理论模型的反演误差结果分析
Table 2 Analysis of inversion error results of theoretical model

反演对比VP/(m·s-1)VS/(m·s-1)ρ/(g·cm-3)ε(v)δ(v)γ(v)理论值3 0701 9802.4000.2300.300反演值3 0852 0302.530.0010.2570.289误差/%0.482.535.420.1011.743.67

表3 有噪声的理论模型反演结果误差分析
Table 3 Analysis of inversion error results of theoretical model with noise

信噪比误差/%VPVSρε(v)δ(v)γ(v)200.190.450.830.13.043.22100.070.201.670.18.0011.7650.331.113.270.12.8012.4

利用模型的物性参数和各向异性参数反演值,根据多层模型纵向深度为301~399 m、横向距离为100 m处的反演结果,利用SA-PSO算法对模型中砂岩层的流体指示因子KN/KT进行预测,结果如图6所示。由图6可知,利用抽取的部分模型物性参数和各向异性参数作为先验值,得到的流体指示因子KN/KT均在0.4以下,最优值在0.12左右,符合裂隙饱和水时流体因子较小,接近于0的规律。

图6 基于SA-PSO算法预测得到的流体因子
Fig.6 Fluid factors predicted based on SA-PSO algorithm

2.2 改进的Marmousi2部分模型试算

原始Marmousi2模型如图7所示,是复杂二维构造的地质-地球物理模型,属于各向同性介质,本文选取红框圈定的部分模型,并根据不同岩层的特性赋予各向异性参数:将含水砂岩之外的薄层作为各向同性介质,假定含水砂岩垂向裂隙发育,将其近似等效为HTI介质,且由于HTI介质在饱水裂隙发育情形下,各向异性参数ε(v)接近0,用于各向异性参数反演和裂隙流体因子的预测。采用主频为50 Hz的雷克子波序列与HTI介质的PP波反射系数序列进行褶积,并对横向距离为2 000 m、深度为 500~800 m的剖面,生成角道集如图8所示。由图8可知,反射同相轴大致能够反映所选模型区域反射界面物性差异特征,反射振幅在裂隙附近呈现正极性,能量较弱。

图7 Marmousi2模型及抽取区域示意[28]
Fig.7 Marmousi2 model and the selected area[28]

图8 Marmousi2改进模型抽取的角道集(2 000 m处)
Fig.8 Synthetic angle gather in the Marmousi2 model (at the distance of 2 000 m)

采用本文研究的反演方法对这部分数据进行迭代次数为100次的SA-PSO反演,反演结果如图9~12所示。分析图9~11可以看出,对改进的Marmousi2部分模型反演得到VPVSρ剖面,虽然反演值与理论值有一定误差,但仍然能清晰反映所抽取的Marmousi2部分模型的地层变化规律,分层效果明显,与实际模型基本吻合,反演结果较好。如图12所示的各向异性参数反演结果曲线对比,各向异性参数ε(v)γ(v)反演精度略高于各向异性参数δ(v),对反演目标方程的敏感度相对较好。从多组反演结果与理论值的对比剖面可以看出,反演结果能够较准确反映原始剖面形态特征,且与理论模型的拟合度高,误差比较小。选取模型中横向距离为2 000 m、深度为 630~660 m的区域,对流体指示因子KN/KT进行预测,所得结果如图13所示,由于抽取的区域均为饱水裂隙区域,此时进行间隔采样得到的流体因子就是PSO-SA算法中的“粒子”,“粒子们”最终迭代出最优结果(即红色较大圆点),因子粒子寻优时其初始值是随机的,粒子的总数量与KN/KT的总数量一致,但粒子的取值与KN/KT的值不具有一一对应的关系,图中所有粒子的值都小于1,与饱水裂隙中KN/KT的取值范围一致。

图9 基于SA-PSO算法的 Marmousi2改进模型VP反演剖面对比
Fig.9 Model values and inversion results of P-wave velocity based on SA-PSO in the Marmousi2 model

图10 基于SA-PSO算法的 Marmousi2改进模型VS反演剖面对比
Fig.10 Model values and inversion results of S-wave velocity based on SA-PSO in the Marmousi2 model

图11 基于SA-PSO算法的 Marmousi2改进模型ρ反演剖面对比
Fig.11 Model values and inversion results of ρ based on SA-PSO in the improved Marmousi2 model

图12 基于SA-PSO算法的Marmousi2改进模型各向异性 参数反演结果对比
Fig.12 Model values and inversion results of anisotropic para- meters based on SA-PSO in the improved Marmousi2 model

图13 基于SA-PSO算法的Marmousi2改进模型 流体因子预测
Fig.13 Fluid factor prediction of Marmousi2 improved model based on SA-PSO algorithm

3 结 论

(1)综合考虑粒子群算法快速收敛的特点和模拟退火的概率性跳出特性,通过对两种算法的混合使用,实现了利用基于模拟退火的混合粒子群算法对HTI介质的AVO反演。该算法具有全局寻优能力较强,克服早熟现象,不易陷入局部极小值等优点。

(2)将基于SA-PSO算法的HTI介质叠前AVO反演应用于多层模型和复杂模型的反演计算中,反演结果与理论模型拟合较好,计算结果稳定,抗噪性较好,反演精度较高。

(3)根据裂隙含水时流体指示因子KN/KT较小的特点,并借助HTI介质各向异性反演,通过建立裂隙流体因子KN/KT与各向异性参数的关系,可实现流体因子预测。

参考文献(References):

[1] 戴世鑫.基于物理模型的煤田地震属性响应特征的关键技术研究[D].北京:中国矿业大学(北京),2012.

DAI Shixin.Research on key technology of response characteristics of seismic attributes based on the physical model[D].Beijing:China University of Mining and Technology(Beijing),2012.

[2] THOMSEN L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966.

[3] 董守华.气煤弹性各向异性系数实验测试[J].地球物理学报,2008(3):947-952.

DONG Shouhua.Test on elastic anisotropic coefficients of gas coal[J].Chinese Journal of Geophysics,2008(3):947-952.

[4] 陈同俊,崔若飞,刘恩儒.构造煤AVO特征及正演模拟研究[J].地球物理学进展,2008,23(5):1610-1615.

CHEN Tongjun,CUI Ruofei,LIU Enru.AVO characters of tectonic coal and its modeling research[J].Progress in Geophysics,2008,23(5):1610-1615.

[5] 彭苏萍,毕银丽.黄河流域煤矿区生态环境修复关键技术与战略思考[J].煤炭学报,2020,45(4):1211-1221.

PENG Suping,BI Yinli.Strategic consideration and core technology about environmental ecological restoration in coal mine areas in the Yellow River basin of China[J].Journal of China Coal Society,2020,45(4):1211-1221.

[6] SCHOENBERG M,DOUMA J.Elastic wave propagation in media with parallel fractures and aligned cracks[J].Geophysical Prospecting,1988,36(6):571-590.

[7] BAKULIN A,GRECHKA V,TSVANKIN I.Estimation of fracture parameters from reflection seismic data—Part I:HTI model due to a single fracture set[J].Geophysics,2000,65:1788-1802.

[8] QUINTAL B,SCHMALHOLZ S M,PODLADCHIKOV Y.Impact of fluid saturation on the reflection coefficient of a poroelastic layer[J].Geophysics,2011,76:1-12.

[9] 陈怀震,印兴耀,高成国,等.基于各向异性岩石物理的缝隙流体因子AVAZ反演[J].地球物理学报,2014,57(3):968-978.

CHEN Huaizhen,YIN Xingyao,GAO Chengguo,et al.AVAZ inversion for fluid factor based on fracture anisotropic rock physics theory[J].Chinese Journal of Geophysics,2014,57(3):968-978.

[10] 郎玉泉,陈同俊,马丽,等.煤层顶板砂岩富水性AVO预测技术[J].煤田地质与勘探,2018,46(3):127-132.

LANG Yuquan,CHEN Tongjun,MA Li,et al.Water content prediction of roof sandstone using AVO technology[J].Coal Geology & Exploration,2018,46(3):127-132.

[11] RÜGER A.P-wave reflection coefficients for transversely isotropic models with vertical and horizontal axis of symmetry[J].Geophysics,1997,62(3):713-722.

[12] 王宏伟,彭苏萍,杜文凤,等.HTI 构造煤方位 AVO 正演[J].石油地球物理勘探,2014,49(6):1122-1130.

WANG Hongwei,PENG Suping,DU Wenfeng,et al.AVO forward modeling of HTI structural coal orientation[J].Oil Geophysical Prospecting,2014,49(6):1122-1130.

[13] 梁锴,印兴耀,吴国忱.TTI介质qP波入射精确和近似反射透射系数[J].地球物理学报,2011,54(1):208-217.

LIANG Kai,YIN Xingyao,WU Guochen.Exact and approximate reflection and transmission coefficient for incident qP wave in TTI media[J].Chinese Journal of Geophysics,2011,54(1):208-217.

[14] 何樵登,陶春辉.用遗传算法反演裂隙各向异性介质[J].石油物探,1995,34(3):46-50.

HE Qiaodeng,TAO Chunhui.Inversion of fractured anisotropic media by genetic algorithm[J].Petroleum Geophysical Prospecting,1995,34(3):46-50.

[15] 轩义华.倾斜横向各向同性介质(TTI)参数反演方法研究[D].长春:吉林大学,2007.

XUAN Yihua.Study on parameter inversion method of tilted transversely isotropic medium (TTI)[D].Changchun:Jilin University,2007.

[16] 张广智,陈怀震,印兴耀,等.基于各向异性AVO的裂缝弹性参数叠前反演方法[J].吉林大学学报(地球科学版),2012,42(3):845-851.

ZHANG Guangzhi,CHEN Huaizhen,YIN Xingyao,et al.Method of fracture elastic parameter inversion based on anisotropic AVO[J].Journal of Jilin University (Earth Science Edition),2012,42(3):845-851.

[17] MA X Q.Simultaneous inversion of prestack seismic data for rock properties using simulated annealing[J].Geophysics,2002,67(6):1877-1885.

[18] MISRA S,SACCHI M D.Global optimization with model-space preconditioning:Application to AVO inversion[J].Geophysics,2008,73(5):R71-R82.

[19] KENNEDY J,EBERHART R.Particle Swarm optimization[J].Proceedings of IEEE International Conference on Neural Network,1995,4:1942-1948.

[20] MIKKI S M.Quantum particle swarm optimization for electromagnetics[J].IEEE Transactions on Antenna and Propagation,2006,54(10):2764-2775.

[21] SUN S Z,LEI C,BAI Y,et al.PSO non-linear pre-stack inversion method and the application in reservoir prediction[A].SEG Technical Program Expanded Abstracts[C].Las Vegas,2012:1949-1953.

[22] SUN S Z,LIU L.A numerical study on non-linear AVO inversion using chaotic quantum particle swarm optimization[J].Journal of Seismic Exploration,2014,23(4):379-392.

[23] MEHRGINI B,IZADI H,MEMARIAN H.Shear wave velocity prediction using Elman artificial neural network[J].Carbonates and Evaporites,2017,34(4):1281-1291.

[24] WU Q,ZHU Z,YAN X,et al.An improved particle swarm optimization algorithm for AVO elastic parameter inversion problem[J].Concurrency and Computation:Practice and Experience,2019,31(9):1987-2003.

[25] 严哲,顾汉明.量子行为的粒子群算法在叠前AVO反演中的应用[J].石油地球物理勘探,2010,45(4):516-519.

YAN Zhe,GU Hanming.Application of quantum-behaved particle swarm optimization algorithm in AVO Inversion[J].Oil Geophysical Prrospecting,2010,45(4):516-519.

[26] 汤红伟,程建远,王世东.深层煤矿床的煤岩样物性测试结果与分析[J].中国煤炭,2009,35(9):75-78.

TANG Hongwei,CHENG Jianyuan,WANG Shidong.The test results and its analysis of deep coal seam and rock sample[J].China Coal,2009,35(9):75-78.

[27] 李东会,董守华,赵小翠,等.煤储层各向异性地震波场模拟[J].物探与化探,2010(5):604-609.

LI Donghui,DONG Shouhua,ZHAO Xiaocui,et al.Seismic wave simulation of anisotropy in coal-bed media[J].Geophysical and Geochemical Exploration,2010(5):604-609.

[28] MARTIN G S,WILEY R,MARFURT K J.Marmousi2:An elastic upgrade for Marmousi[J].The Leading Edge,2006,25(2):156-166.

Anisotropic AVO inversion of seismic wave in HTI media

LI Qin,WANG Wei,WANG Hanlin

(School of Geology and EnvironmentXian University of Science and TechnologyXian 710054,China)

Abstract:Conventional AVO inversion is mainly carried out for isotropic media,but the strata with fractures are actually anisotropic.Seismic inversion is a typical ill-posed problem.With the addition of anisotropic parameters,the nonlinear degree of seismic inversion is increased significantly,which makes it more difficult than traditional inversion.The SA-PSO algorithm is a hybrid optimization of particle swarm optimization and simulated annealing algorithm,which has the advantages of stronger global optimization ability,overcoming precocity phenomenon and not easy to fall into local minimum value.In this paper,the medium containing vertical cracks is approximated as HTI medium.Based on the HTI medium reflection coefficient equation,the objective function of anisotropic parameters is established.Anisotropic parameter inversion is realized by SA-PSO algorithm.Furthermore,the fracture related parameters and fluid factors are predicted based on the correlation between the fluid factor,fracture weakness and anisotropy parameters.The theoretical model of coal seam containing HTI type is established,and the original data are inverted by SA-PSO algorithm and PSO algorithm at the same time.The superiority of the algorithm is tested.The inversion accuracy of SA-PSO algorithm is higher and the inversion process is more stable.By adding Gaussian random noise with SNR of 5,10,20 to the data of multi-layer theoretical model,the anti-noise performance of the inversion method is tested.The inversion results show that the error of the inversion results is small,the inversion process can converge quickly,and the anti-noise performance of the inversion is good.According to the characteristics of anisotropic parameters in the HTI type water-bearing fracture,the water-bearing sandstone part of Marmousi2 model is modified,and anisotropic parameters are introduced.The velocity,density and anisotropic parameter profiles are obtained through inversion,and the fluid factor prediction is realized.The inversion and prediction results are consistent with the theoretical model,which demonstrates the effectiveness of the method.

Key words:HTI media;anisotropy;AVO inversion;fluid factor;fracture

中图分类号:P631.4

文献标志码:A

文章编号:0253-9993(2021)06-1925-11

移动阅读

收稿日期:2021-02-03

修回日期:2021-04-20

责任编辑:黄小雨

DOI:10.13225/j.cnki.jccs.ST21.0246

基金项目:国家自然科学基金资助项目(41674135)

作者简介:李 勤(1979—),女,湖南岳阳人,副教授,博士。E-mail:eriliqin@126.com

引用格式:李勤,王玮,王瀚林.HTI介质地震波各向异性AVO反演[J].煤炭学报,2021,46(6):1925-1935.

LI Qin,WANG Wei,WANG Hanlin.Anisotropic AVO inversion of seismic wave in HTI media[J].Journal of China Coal Society,2021,46(6):1925-1935.