0 引言

页岩储层中的非常规资源已经成为油气勘探的研究热点。页岩结构复杂，内部物质的不均匀分布和定向排列使其通常表现出各向异性；可以通过研究页岩储层内部结构属性和各向异性的联系，为识别和勘探富含油气区域提供依据。

岩石物理模型是建立页岩微观结构与各向异性联系的常用工具^[1-2]。但页岩各向异性受多种因素影响，包括黏土矿物的定向排列、非球形孔隙、微裂隙和裂缝属性^[3]。黏土矿物定向排列是由页岩的压实作用引起的，不同程度的压实会引起黏土矿物不同的排列状态，进而影响各向异性。Hornby等^[4]和Sayers^[5]使用方位分布函数(ODF)来模拟黏土矿物的方位分布状态，然后使用各向异性微分等效介质模型(DEM)来估计页岩各向异性。Sayers^[6]和Bandyopadhyay^[3]指出非球形孔隙是页岩各向异性的另一个来源。Kwon等^[7]发现页岩中裂隙的形状、大小和定向排列影响各向异性。Chapman模型^[8]使用了多尺度孔隙系统，包括球形孔隙、裂隙和裂缝，并假设各向异性只来源于水平排列的裂缝。Jiang等^[9]建立了一个双孔隙系统模型，即结合自相容等效介质理论(SCM)和Chapman模型来克服Chapman模型中不包含非球形孔隙的缺点，并将该模型应用于Haynesville页岩。

对于裂缝型页岩储层，各向异性由裂缝属性主导。本文在岩石物理建模过程中参照Jiang等^[9]的工作，使用双孔隙系统模型以准确描述页岩的孔隙系统，并针对实际研究区域对该模型进行校正。由于模型的复杂性，反演需要确定非球形孔隙度、孔隙纵横比和裂缝密度，因此反演过程使用模拟退火优化粒子群算法(SA-PSO)^[10]。作为一种智能寻优算法，它能够快速准确地寻找到各参数最优解，适用于复杂的非线性反演问题。该算法被应用到位于中国西南部的龙马溪组页岩，估计了该储层物性参数，如非球形孔隙度、孔隙纵横比、裂缝密度和各向异性参数。

1 岩石物理模型

横向各向同性(TI)页岩的各向异性属性可以由Thomsen^[11]各向异性参数表示：

式中：c_ij是刚度张量分量，由岩石内部属性决定；v_P、v_S和ρ分别是纵波速度、横波速度和密度；ε、γ、δ是各向异性参数。岩石物理模型能够描述岩石内部属性并得到等效弹性张量，然后使用实测纵、横波速度对模拟结果进行校准。

包络体理论是页岩建模的一种常用方法。该理论把岩石中的组分当作不同形状的包含物，形状由纵横比决定。SCM^[12]是一种各向同性包络体理论，能够结合多种矿物和孔隙类型构建等效介质模型，等效体积模量K_sc和等效剪切模量μ_sc由下式给出：

式中：N是组分数目；f_i是每种组分的体积分数；K_i和μ_i是每种组分的体积模量和剪切模量；U_i和V_i是几何因子，由纵横比和模量决定，具体公式由Berryman^[13]给出。图 1显示了SCM模型中不同孔隙纵横比(α)和孔隙度(φ)下v_P和v_S的变化情况。所用基质是黏土，颗粒纵横比为0.1，孔隙为水饱和。由图 1可见：随着α的增加，v_P和v_S逐渐增大；而随着φ的增加，v_P和v_S都变小。

Chapman多尺度孔隙模型^[8]可用来描述流体喷射流动现象对岩石等效模量的影响。波长小于孔隙尺寸的波在岩石中传播时会产生压力梯度，喷射流动是指由此产生的流体相互作用^[14]。该模型考虑了频率和孔隙压力的影响，适用于各向同性和各向异性介质^[9]。Chapman模型认为岩石包含球形孔隙、裂隙和裂缝。图 2是该模型的示意图^[15]，其中球形孔隙和裂隙随机分布，裂缝水平排列，各向异性由裂缝引起。除了裂缝与裂缝之间，球形孔隙、裂隙和裂缝相互连通来保证流体的流动性^[8]。Chapman模型中的孔隙度(φ_chap)可表示为球形孔隙度(φ_rp)、裂隙孔隙度(φ_c)和裂缝孔隙度(φ_f)三者之和：

裂隙和裂缝被认为是理想椭球体，具有相同的纵横比，因此裂隙密度(ε_c)和裂缝密度(ε_f)由下式得出：

Chapman模型的具体公式由Jiang等^[9]总结给出。模型中参数都有相应的物理意义，并且大部分可以直接测定或通过测定值估计得到。图 3a是孔隙纵横比对各向异性参数的影响。由图 3a可见：v_S、ε和δ随着α的增大而增大，其中δ的变化最大；而v_P和γ却随着α的增大而减小。图 3b为各向异性参数随裂缝密度的变化。由图 3b可见：v_P、v_S和δ随着ε_f的增大而减小；而ε和γ却随着ε_f的增大而增大，其中γ变化最大。模型所用基质是黏土，孔隙为水饱和。

Chapman模型中孔隙被认为是球形，但大部分页岩包含非球形孔隙。Jiang等^[9]将SCM和Chapman模型相结合构建了双孔隙系统模型来估计Haynesville页岩的储层属性。在建模过程中，SCM模型用来描述矿物和非球形孔隙，非球形孔隙间相互不连通。Chapman模型包含球形孔隙、裂隙和裂缝，其中裂缝间相互不连通，裂隙和球形孔隙至少与一个裂缝相连通，用来描述流体喷射流动现象对岩石等效模量的影响。该双孔隙系统模型将总孔隙度(φ)划分为非球形孔隙(φ_np)和Chapman模型孔隙系统，这两个孔隙系统相互不连通，可表示为

使用SCM模型时需要考虑各包络体的纵横比，本文中黏土矿物纵横比是0.10，脆性矿物，例如石英、方解石、黄铁矿等是1.00，干酪根是0.01。SCM模型中影响速度的主要因素是孔隙形状，在建模流程中，SCM模型中的非球形孔隙与Chapman模型孔隙系统中的裂隙和裂缝有相同的纵横比，且非球形孔隙尺寸与矿物颗粒以及球形孔隙、裂隙尺寸一致^[9]。裂缝尺寸对喷射流动现象影响很大，通常可通过岩心数据确定。

双孔隙系统页岩岩石物理建模流程由图 4给出。建模过程中，首先使用SCM模型构建包含矿物和非球形孔隙的各向同性基质，矿物和孔隙在该基质中均匀并且随机方位分布且孔隙互不连通。然后使用Chapman模型向该基质中加入球形孔隙、裂隙和裂缝的影响，其中裂隙的方位随机，裂缝呈水平排列，最终得到各向异性页岩模型。在该岩石物理模型中：各向异性来源于水平裂缝，主要受孔隙纵横比和裂缝密度影响；而各向同性基质的纵、横波速度受非球形孔隙度和孔隙纵横比影响；Chapman模型中其他参数可根据实际数据给定。图 5是岩石物理模型中非球形孔隙度、纵横比、裂缝密度这3个重要储层参数对各向异性参数的影响，所用基质是黏土，孔隙是水饱和。其中：φ_np对v_S和γ的影响较大；α较小时对各向异性参数的影响大，但随着α增加，影响逐渐减小；ε_f的取值对各向异性参数的影响始终较大。

2 实际数据分析

研究对象是位于中国西南四川盆地的龙马溪组页岩。测井数据在图 6中显示，包括伽马射线、纵横波速度、密度、孔隙度和组分信息。目标层段页岩致密，孔隙度较低，干酪根和黏土含量高。岩心资料显示该层段水平裂缝发育，裂缝尺寸平均值是100 μm。

目标层段是包含水平裂缝的横向各向同性(TI)页岩。邓继新等^[16]对龙马溪组页岩微观结构的研究发现：该页岩的孔隙类型包含原生孔隙、次生孔隙和裂隙，各种孔隙尺度不同，但大部分小于2 μm；孔隙大部分是非球形孔隙，且呈孤立分布，连通性差。

岩石物理模型中参数众多，部分参数可由前人研究给定。Guo等^[17]直接使用了Chapman^[8]给出的频率和弛豫时间参数值。Jiang等^[9]使用双孔隙系统模型对Haynesville页岩储层参数进行了估计。龙马溪组页岩类比于Haynesville页岩，两者主要矿物组分都是黏土、石英和方解石，孔隙度和渗透性都很低，并且都具有横向各向同性，因此认为两者具有相似性。据此，我们为岩石物理建模过程中部分参数赋予了相同值。其中：Chapman模型中频率(f)取10 kHz，接近测井频率；非球形孔隙、球形孔隙、裂隙的尺寸(ξ)相同，取1 μm；裂隙密度由裂缝密度得到，

水饱和裂隙的弛豫时间(τ_m)由Chapman^[8]给出；而裂缝的弛豫时间(τ_f)由它们的尺寸(a_f)决定，

矿物和流体的模量和密度由表 1给出。最终，岩石物理模型中的未知参数可确定为非球形孔隙度、孔隙纵横比和裂缝密度。通过使用测量得到的纵、横波速度可反演得到这3个储层参数，进而估计页岩各向异性。

3 反演方法和结果

反演目标函数通常可写作

式中：v_P-real和v_S-real是纵、横波速度的测量值；v_P-model和v_S-model是由岩石物理模型计算的纵、横波速度。

传统反演方法基于遍历搜索，即将各参数取值空间划分成n个区间，遍历各参数的各取值区间来寻找满足目标函数的值^[18]。对于本文中的三参数反演问题，传统方法需要迭代n³次。岩石物理模型复杂，若区间数过多则计算量大，过少则无法找到最优解。而且，多参数反演存在多解性，需要加入更多反演约束条件来减少多解性。地震反演方法通常在目标函数中加入约束项来减少多解性，李坤等^[19]使用平滑模型作为约束项来减少地震反演中低频趋势的多解性。本文考虑到实际储层属性存在垂向连续性、各深度点储层参数并不完全独立，在目标函数中加入平滑约束项后为

式中：a、b、c是各参数平滑项的权值；、、是各参数上一深度点处的反演值。使用模拟退火优化粒子群算法来求解储层参数φ_np、α、ε_f。该智能算法结合了粒子群优化算法和模拟退火算法，具有收敛快速、避免陷入局部极值的优点，常用于求解非线性问题，对于复杂岩石物理模型有很好适用性。

基于岩石物理模型和模拟退火优化粒子群算法的页岩储层参数反演流程(图 7)如下。

1) 初始化参数。从参数取值区间抽取M个随机点，被称作粒子：

式中：x₀ⁱ和v₀ⁱ是第i个粒子的初始化位置和速度，位置也代表参数的取值，粒子可以是多维的，每个维度对应一个参数；U表示从区间等概率取值；S是参数取值区间；S_u和S_d是S的上下界；s是速度初始化区间。初始化退火温度T₀应足够大，取T₀=10 000。

本问题中3个储层参数的取值范围分别是：

2) 计算粒子适应度。

式中：fit是粒子适应度函数，由目标函数值决定；xⁱ为粒子位置。

3) 移动粒子。

式中：学习因子c₁和c₂是非负常数，且c₁+c₂=4.1；r₁和r₂是属于[0, 1]的随机数；x_kⁱ和v_kⁱ是第k次迭代时粒子的位置和速度；v_k+1ⁱ指经过本次迭代后粒子的速度；p_g是全局最优值，指本次迭代时所有粒子中适应度值最大的粒子所在位置；p_i是个体最优值，指单个粒子移动至今适应度最大时的位置；xⁱ_temp是本次迭代后粒子位置的可能取值。若满足式(14)中的任何一个，则粒子移动到位置xⁱ_temp；否则粒子位置不变。

式中：T_k是当前温度；r₃是属于[0, 1]的随机数。

4) 降温操作。

式中，c是降温系数，取0.93。

5) 重复步骤2)至4)直到满足收敛条件。

对目标层段进行反演，使用100个粒子，迭代60次。图 8展示了一个深度点处的反演迭代收敛过程。由图 8可以看出，粒子的适应度由低变高，并且逐渐聚集到最优解处。最优解满足目标函数公式(9)，此时粒子在三维参数空间的位置也就是各参数的取值。如果传统方法中每个参数区间划分50份，则需要迭代125 000次，而模拟退火优化粒子群算法共需迭代6 000次，计算量得以大大降低。

图 9是目标层段使用模型计算出的纵、横波速度值和实测值，两者拟合很好，表明了本文方法的准确性。图 10是反演得到的储层参数曲线，可清楚地看出储层参数随深度的变化趋势，并首次提供了该研究区域的非球形孔隙度。Qian等^[20]对龙马溪组页岩的研究表明，该层段孔隙纵横比较小，平均低于0.05，与本文反演值(图 10)一致。Jiang等^[9]对Haynesville页岩的研究给出它的裂缝密度是0.04，由于本文目标层裂缝发育，反演值较高(图 10)。图 11是估计的各向异性参数，其中ε和γ具有强相关性，而γ和δ在部分区域出现负相关；该现象与Bandyopadhyay^[3]及Sayers^[6]的研究一致。目标层段各向异性参数γ总体大于ε，并且δ为负值；Qian等^[20]的研究也给出了相似的结果，但目标层各向异性强度由于裂缝发育而更大。

4 结论与展望

本文使用改进的岩石物理模型和模拟退火优化粒子群算法对页岩裂缝属性和各向异性进行了反演。岩石物理建模为使用等效介质模型和Chapman模型相结合建立的双孔隙系统模型，以准确地描述页岩微观孔隙结构，反演过程使用智能算法准确求解复杂岩石物理模型下三参数反演问题。将方法应用到龙马溪组页岩气储层，反演结果与实测数据和已有研究结果相匹配，说明了本文方法的有效性。

本文的研究为页岩储层的描述提供了更多的有效信息。页岩结构复杂，影响各向异性的因素很多，进一步的研究需要使用更加准确复杂的岩石物理模型进行模拟。而求解复杂岩石物理模型，反演方法同样需要进一步提高。