李 勤，王 玮，王瀚林

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

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维空间，粒子i在n维空间的位置表示为向量Xi=(x1,x2,x3,…,xn)，飞行速度表示为向量Vi=(v1,v2,v3,…,vn)，粒子不仅知道自己发现的最好位置Pb，也知道整个群体中的其他粒子发现的最优位置Pb,g，粒子们根据这些信息来决定接下来的行动方向。根据目标函数的构建，每个粒子都有其适应值范围。粒子群优化算法首先将粒子初始化为一群随机粒子，然后通过迭代找到最优解，在每一次迭代中，粒子通过跟踪两个“最优解”(即Pb,Pb,g)来更新自己，在找到“最优解”(Pb,Pb,g)后，粒子通过下式更新自己的位置:

(1)

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

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

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

(2)

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

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

(3)

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

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

(4)

在模拟退火的跳出机制下，认为Pb是和Pb,g悬殊很大的特殊解，可以计算出温度为T时Pb相对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)

根据线性滑动模型理论中应力与裂隙应变的关系，裂隙介质弹性系数矩阵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

图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可以看出，从VP和VS的反演剖面来看，模型的地层能被清晰地分辨，反演结果较好。对于ε(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

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

利用模型的物性参数和各向异性参数反演值，根据多层模型纵向深度为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部分模型反演得到VP，VS和ρ剖面，虽然反演值与理论值有一定误差，但仍然能清晰反映所抽取的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与各向异性参数的关系，可实现流体因子预测。