CN112505772B - 一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法 - Google Patents

Abstract

本发明公开了一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，包括如下处理步骤：步骤一、测量岩心在不同压力下的饱和、干燥速度，确定岩石密度等性质；步骤二、用多形态裂隙的孔裂隙理论模拟岩石的弹性波速度；步骤三、对于岩石中存在多种形态的裂隙，不同压力点下的孔隙纵横比谱可以由0有效压力的孔隙纵横比谱计算；步骤四、建立反演目标函数；步骤五、设置多种纵横比的裂隙，包含孔隙；步骤六、反复调节岩石在0有效应力下的各形态裂隙的纵横比与裂隙密度，使目标函数达到最小值，就得到了岩石在各个压力点下的孔隙纵横比谱。本发明可以解决更准确地获取岩石孔隙结构特征，分析岩石的力学、声学和流体渗透性质等技术问题。

Description

一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征 的方法

技术领域

本发明属于岩石物理领域，具体涉及一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法。

背景技术

岩石中的孔隙分布特征对岩石的力学、声学和流体渗透性质有着十分重要的影响，是储层岩石声学关注的重点。用纵横比这一参数来度量，岩石孔隙分布的显著特征是包含纵横比～1的孔隙和纵横比<<1的裂隙，而不同形态的裂隙分布可以用纵横比谱来很好地描述。多年来，确定岩石孔隙的纵横比谱一直是岩石物理学研究的一个方向；实验压力加载条件下的弹性波速测量为此研究提供了一条有效途径。许多学者利用不同纵横比的裂隙对压力的不同响应来反演岩石的孔隙纵横比谱。Cheng(1978)和Cheng和

等(1979)使用Kuster和

(1974)的裂缝模型(K-T模型)，结合孔隙体积随压力的变化(

等，1976)，提出了从实验测量的纵、横波速度反演孔隙纵横比谱的方法。David和Zimmerman(2012)提出了另一种用干燥岩石的速度-压力曲线计算孔隙纵横比谱的方法。近年来又有很多学者在这方面做了大量的工作(Yan等，2014；Yan等，2015；邓继新等，2015；Duan等，2018；Han等，2019；李闯等，2020)。

事实上，唐等人(2013)应用该理论成功地模拟和反演了实验室数据，得到了岩石裂隙密度和纵横比随压力的变化曲线。不足的是，该反演模拟的是单一纵横比裂隙体系在压力作用下的变化，不能得到反映岩石孔隙形态分布的纵横比谱。

迄今为止，岩石孔隙结构反演的研究已有相当的进展(如上述的Cheng和

(1979)与David和Zimmerman(2012)等的工作)，但改进理论方法来准确获取岩石的孔隙结构一直是岩石物理专业领域的努力方向。

发明内容

本发明的目的是提供一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，以解决进一步扩展传统理论，使之包括孔隙与多形态裂隙体系的相互作用，将原来单一纵横比裂隙体系随压力的变化描述为不同形态(纵横比)裂隙在不同压力下的闭合，以达到更准确地获取岩石孔隙结构特征，确定岩石的力学、声学和流体渗透性质等技术问题。

为实现上述发明目的，本发明采用如下处理步骤：

步骤一、测量岩心在不同压力下干燥与饱含流体的速度，并确定岩石密度、岩石渗透率、岩石孔隙度、流体密度、流体模量；

步骤二、用多形态裂隙的孔裂隙理论隙模拟和计算岩石的弹性波速度；

步骤三、基于压力对孔隙结构影响的计算公式，对于岩石中存在多种形态的裂隙，不同压力点下的孔隙纵横比谱可以由0有效压力的孔隙纵横比谱计算；

步骤四、将实验测量的不同压力点的岩石速度(可以仅为饱和条件下的纵、横波速度；或联合饱和、干燥条件下的纵、横波速度)与由步骤二理论计算的各压力点的速度做最小二乘拟合，其方差作为目标函数，并利用步骤三的压力对孔隙结构影响的计算公式，化简目标函数；

步骤五、设置多种形态的裂隙(包含孔隙)，其中孔隙的纵横比大于0.01(通常不超过2种)，裂隙的纵横比小于0.01；

步骤六、反复调节岩石在0有效应力下的各形态裂隙的纵横比与裂隙密度，使目标函数达到最小值，就得到了岩石在各个压力点下的孔隙纵横比谱。孔隙结构的反演可以看作是搜寻目标函数最小值的过程，为了降低反演结果的非唯一性，采取全局优化的GA算法(Goldberg，1988)对目标函数进行求解。将反演得到的孔隙纵横比谱代入到步骤二的多形态裂隙并存的孔裂隙理论,计算和预测岩石样品的纵横波波速随压力的变化关系。通过对比理论计算得到的纵横波速度与实验超声测量得到的纵横波速度的吻合程度,可以对反演得到的孔隙纵横比谱进行质量监控和验证。

所述步骤二具体为：

(1)唐(2011)和Tang等(2012)的孔、裂隙并存的双孔介质弹性波理论给出饱和岩石的体积模量、剪切模量的表达式为：

其中K_d为干燥岩石的体积模量，α＝1-K_d/K_s，β＝(α-φ)/K_s+φ/K_f，K_s为岩石基质的体积模量，φ为岩石的孔隙度，S(ω)为描述孔隙与裂隙相互作用的挤喷流函数，包含了裂隙密度和裂隙纵横比这两个描述裂隙的重要参数，K₀、μ₀分别为S(ω)＝0时的饱和岩石体积模量与剪切模量。

针对硬币型的裂隙，唐(2011)推导出S(ω)的表达式为：

其中

ω为角频率，η为孔隙流体黏度，ε为裂隙密度，γ为裂隙纵横比，K₀、μ₀为S(ω)＝0时背景介质的体积模量、剪切模量，ν₀为干燥介质的泊松比，K_d、μ₀和ν₀可由Biot相洽理论计算(Thomsen，1985)，K₀可由Gassmann方程计算，J_n(n＝0,1)为第一类n阶贝塞尔函数。

上述硬币模型中孔隙与裂隙的流体交换在硬币的边缘，但力学上裂隙在此是闭合的，作为模型的改进，Tang等(2012)提出了钹状的孔、裂隙模型，将流体交换放到硬币模型的中部，由此推导出的S(ω)表达式为：

其中

值得指出的是，通过调节与模型有关的流体挤喷驰豫频率，可以使二模型的计算结果一致，因此这两种模型都可以用作以下反演模型的正演计算，但钹状模型的计算相对简单快捷，本发明的叙述采用钹状模型及其相关的理论公式。

孔、裂隙介质中的快纵波、慢纵波和横波的波数由下式计算(Tang等，2012)：

其中下标p和s分别代表纵波和横波，+和-分别代表快纵波和慢纵波，上述公式的符号表达式为：

b₀＝-[β+S(ω)](K^Tang+4μ^Tang/3)/α

c＝(α-b_sρ/(ρ_fb₀))/(α+b_s)

b_s＝ρ_fθω²

ρ＝ρ_s(1-φ)+ρ_fφ

其中ρ_s与ρ_f分别为岩石固体基质和孔隙流体的密度，θ＝iκ(ω)/(ηω)，其中κ(ω)为Johnson等(1987)推导出的动态渗透率，κ(ω)的具体表达式为：

其中κ₀为达西渗透率，τ为孔隙内流体的弯曲度。

在得到三种弹性波的波数后，可计算弹性波的速度频散和衰减，具体表达式为(Tang and Cheng，2004)：

其中v与Q分别为速度与品质因子，Re{k}与Im{k}分别取k的实部和虚部。

(2)多形态裂隙并存的孔裂隙理论：现在考虑岩石中存在多个形态裂隙的情况，在Tang等(2012)提出的仅含有单一形态裂隙的孔裂隙介质中，孔隙流体压力本构方程为：

其中w＝φ(U-u)，U和u分别为流体和固体的位移，q为单位体积岩石由于裂隙挤喷流入孔隙空间的流体体积。当M个形态(用纵横比衡量)不一的裂隙并存时，将各个裂隙挤喷的贡献相加，(10)式的压力本构方程变为：

其中下标m代表第m种纵横比(形态)的裂隙，用S_m(ω)＝φq_m/p表示该种裂隙的挤喷函数，(11)式可化简为：

其中第m种裂隙的S_m(ω)由(4)式给出，其具体表达式为：

其中

ε_m和γ_m分别代表第m种形态裂隙的裂隙密度与纵横比。

基于单一裂隙的公式(1)和(2)，多裂隙体系岩石的体积模量、剪切模量的计算公式变为：

在孔裂隙理论计算公式(5)～(9)式中，将K^new与μ^new代替原来的K^Tang与μ^Tang即可计算多裂隙体系的情况。

所述步骤三具体为：

压力是影响孔隙结构的重要因素，岩石中裂隙的体积随着压力的增大而减小，形态变得扁平，随着压力继续增大，裂隙逐渐闭合；而压力的改变对硬孔的形态没有太大影响。为了分析这一现象，

等(1976)将硬孔和裂隙视为具有一定纵横比的扁球形包含体，并给出了包含体的体积变化率和压力的关系，对于纵横比为γ的扁球形裂隙，其体积变化率与有效压力的关系为：

其中P_e为有效压力，K^*为干岩石的静态有效模量，E_i(i＝1,2,3,4)为纵横比γ与岩石基质模量K和μ的函数，E_i的表达式为：

其中

等(1976)建议K^*、K和μ用干燥岩石动态有效模量代替，本专利采用实验测量的干燥岩石纵、横波速度计算干燥岩石动态有效模量。

Cheng(1978)和Cheng和

(1979)建立了压力与孔隙结构的函数关系，压力为P_n时的裂隙体积含量c_nm(即裂隙孔隙度)可以由0有效压力下的裂隙孔隙度c_0m表示：

其中下标n代表第n个压力点，下标m代表第m个形态的纵横比。当

时，裂隙闭合。

对于压力作用下的扁球形裂隙，其短轴的变化远大于长轴变化，在忽略后一变化的条件下，纵横比的变化率即为孔隙体积变化率：

因此，压力P_n下的裂隙纵横比γ_nm可以由0有效压力下的裂隙纵横比γ_0m表示：

结合(18)～(20)式，得到：

因此，任意压力点下的孔隙纵横比谱(纵横比与孔隙度)可以由0有效压力的孔隙纵横比谱计算得到。

扁球形裂隙的裂隙密度与裂隙孔隙度的关系为：

由(22)式，可以将裂隙孔隙度转化为裂隙密度。

所述步骤四具体为：

(1)首先介绍仅反演饱和数据的方法。采用实验室测量的饱和岩石纵、横波速度反演岩石的裂隙密度谱与纵横比谱，而后利用(22)式将裂隙密度转化为裂隙孔隙度，最终得到孔隙纵横比谱。在反演中，将理论计算速度与实验速度的误差平方和作为目标函数，调节参数使目标函数达到最小值，目标函数的具体表达式如下：

其中运算符Min表示取函数的最小值；

与

分别为实验测得的饱和岩石的纵、横波速，

与

分别为理论计算得到的饱和岩石的纵、横波速度；ε_N×M与γ_N×M分别为裂隙密度与纵横比，下标N代表压力加载N次，下标M代表0有效压力时共存在M个形态的裂隙，考虑到裂隙随着压力的增大逐渐闭合，因此每个压力点的裂隙形态总数不会大于M；φ_p为岩石除去裂隙部分的孔隙度，这部分孔隙对应的纵横比大于0.01；K_s和μ_s分别为岩石基质的体积模量和剪切模量；流体体积模量K_f在实验中是已知的。(23)式所描述的反演是多参数反演，对于N个压力测量点的纵、横波速度数据，最多需要反演2×M×N+3个未知参数。

下面分析压力对裂隙形态的影响，并简化目标函数。如果已知0有效压力下的孔隙纵横比谱，任意压力点下的孔隙纵横比谱均可由其计算。因此，在反演中计算

与

时采用的各压力点的孔隙纵横比谱均可由0有效压力下的孔隙纵横比谱计算。由此(23)式可化简为：

(24)式只需要反演0有效压力下的孔隙纵横比谱，再由反演得到的0有效压力下的孔隙纵横比谱计算各压力点下的孔隙纵横比谱，此时需要反演的参数减少为2×M+3个。因此，裂隙体积变化理论(

等，1976)建立裂隙形态与压力关系的同时也减少了反演的参数个数，降低了反演的不确定性。

(2)对于采用饱和、干燥条件的纵、横波速度数据联合反演岩石的孔隙纵横比谱，可以将目标函数(24)式变为

其中实验速度(

与

)与理论速度(

与

)均包括饱和与干燥两种状态，分别对应下标sat与dry。其中饱和条件下的纵、横波速度由(9)式计算，干燥条件下的纵、横波速度的计算公式如下：

其中K_d与μ₀由Biot相恰理论(Thomsen，1985)计算，ρ_d为干燥介质的密度。将(25)式作为目标函数代替(24)式，其他步骤不变，即可实现采用饱和、干燥条件的纵、横波速度数据联合反演岩石的孔隙纵横比谱。

所述步骤五具体为：

下面给出0有效压力下裂隙纵横比的取值及其数量M的确定方法。裂隙纵横比的取值及其数量M的取值受到裂隙闭合、测量压力范围、测量压力点数N以及岩性的影响，为了达到最优的拟合效果，最好令每个压力点下都有裂隙闭合(Cheng，1978；Yan等，2014)。对于裂隙纵横比的取值，应满足如下条件：1、由于把纵横比为0.01的孔隙作为裂隙与硬孔的分界，在孔裂隙理论中硬孔不易变形，且对硬孔纵横比没有限定，因此硬孔的纵横比可在0.01至1之间任取；2、裂隙纵横比的设置要考虑压力的影响，基于(22)式，可以计算出0有效压力下纵横比的取值范围，以此保证每个压力点下都有裂隙闭合，且设置的裂隙纵横比应较均匀地分布；3、考虑到岩性的差异，致密岩石(或高压段速度基本不变的岩石)的裂隙纵横比应该较小，通常要小于高孔隙岩石的纵横比。对于裂隙数量M的确定，应满足如下条件：1、考虑到每个压力点下都有裂隙闭合，M的数量应通常应大于N；2、不同岩性的M有所不同，致密岩石(或高压段速度基本不变的岩石)的M个数应较少，甚至可以小于N。为了提高反演结果的准确性，将0有效压力的总孔隙度作为约束条件：

在孔裂隙理论中，硬孔的孔隙度远大于裂隙的孔隙度，因此对于(27)式，可令0压下的硬孔孔隙度φ_p保持不变，反演时调节0压下的裂隙密度，使0压下的裂隙孔隙度与φ_p之和满足(27)式。

本发明具有以下优点：

(1)本方法扩展了孔裂隙岩石物理理论，将孔裂隙理论扩展至多形态裂隙并存的情况；

(2)本方法结合了有效压力对孔隙结构的影响，建立压力对孔隙结构影响的同时，也使得反演的参数个数大大减少，使反演方法更加快速实用。反演结果与利用其它反演方法和扫描电镜(SEM)技术获得的孔隙纵横比谱具有相似的规律，也与岩石的孔隙特征吻合，可以很好的表征岩石的裂隙分布形态。相比于扫描电镜技术，本方法成本较低。

附图说明

图1为本发明提供的一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法的工作流程图。

图2-a为利用本发明所述方法反演Kayenta砂岩饱和数据得到的反演结果计算的速度与实验速度的对比；

图2-b为利用本发明所述方法反演Kayenta砂岩饱和数据得到的0有效压力下的孔隙纵横比谱；

图2-c为利用本发明所述方法反演Kayenta砂岩饱和数据得到的0有效压力下的孔隙纵横比谱与利用扫描电镜及Cheng反演方法获取的孔隙纵横比谱的对比。

图3-a为利用本发明所述方法联合反演Navajo砂岩饱和、干燥数据得到的反演结果计算的速度与实验速度的对比；

图3-b为利用本发明所述方法联合反演Navajo砂岩饱和、干燥数据得到的0有效压力下的孔隙纵横比谱；

图3-c为利用本发明所述方法联合反演Navajo砂岩饱和、干燥数据得到的0有效压力下的孔隙纵横比谱与利用扫描电镜及Cheng反演方法获取的孔隙纵横比谱的对比。

具体实施方式

本发明将孔裂隙理论扩展至多形态裂隙共存的情况，再结合有效压力对孔隙结构影响(

等,1976)，提出一种考虑多形态裂隙挤喷流效应的孔隙结构反演方法。本发明相较于先前的反演方法，考虑了多形态裂隙挤喷流对速度的影响。结合有效压力对孔隙结构的影响，建立压力对孔隙结构影响的同时，也使得所需反演的参数个数大大减少。反演得到的孔隙纵横比谱可以更加详细准确地描述孔隙结构。结合反演结果，可以实现全频带的频散曲线与衰减曲线预测，更好地解释实验测量的多个弛豫频率与多个衰减峰叠加的现象。

如图1所示，本发明提出了一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，工作流程如下：

步骤一、测量岩心在不同压力下干燥与饱含流体的速度，并确定岩石密度、岩石渗透率、岩石孔隙度、流体密度、流体模量；

步骤二、用多形态裂隙的孔裂隙理论隙模拟和计算岩石的弹性波速度；

步骤二、用多形态裂隙的孔裂隙理论隙模拟和计算岩石的弹性波速度；

(1)唐(2011)和Tang等(2012)的孔、裂隙并存的双孔介质弹性波理论给出饱和岩石的体积模量、剪切模量的表达式为：

其中K_d为干燥岩石的体积模量，α＝1-K_d/K_s，β＝(α-φ)/K_s+φ/K_f，K_s为岩石基质的体积模量，φ为岩石的孔隙度，S(ω)为描述孔隙与裂隙相互作用的挤喷流函数，包含了裂隙密度和裂隙纵横比这两个描述裂隙的重要参数，K₀、μ₀分别为S(ω)＝0时的饱和岩石体积模量与剪切模量。

针对硬币型的裂隙，唐(2011)推导出S(ω)的表达式为：

其中

ω为角频率，η为孔隙流体黏度，ε为裂隙密度，γ为裂隙纵横比，K₀、μ₀为S(ω)＝0时背景介质的体积模量、剪切模量，ν₀为干燥介质的泊松比，K_d、μ₀和ν₀可由Biot相洽理论计算(Thomsen，1985)，K₀可由Gassmann方程计算，J_n(n＝0,1)为第一类n阶贝塞尔函数。

上述硬币模型中孔隙与裂隙的流体交换在硬币的边缘，但力学上裂隙在此是闭合的，作为模型的改进，Tang等(2012)提出了钹状的孔、裂隙模型，将流体交换放到硬币模型的中部，由此推导出的S(ω)表达式为：

其中

值得指出的是，通过调节与模型有关的流体挤喷驰豫频率，可以使二模型的计算结果一致，因此这两种模型都可以用作以下反演模型的正演计算；但钹状模型的计算相对简单快捷，本发明的叙述采用钹状模型及其相关的理论公式。

孔、裂隙介质中的快纵波、慢纵波和横波的波数由下式计算(Tang等，2012)：

其中下标p和s分别代表纵波和横波，+和-分别代表快纵波和慢纵波，上述公式的符号表达式为：

b₀＝-[β+S(ω)](K^Tang+4μ^Tang/3)/α

c＝(α-b_sρ/(ρ_fb₀))/(α+b_s)

b_s＝ρ_fθω²

ρ＝ρ_s(1-φ)+ρ_fφ

其中ρ_s与ρ_f分别为岩石固体基质和孔隙流体的密度，θ＝iκ(ω)/(ηω)，其中κ(ω)为Johnson等(1987)推导出的动态渗透率，κ(ω)的具体表达式为：

其中κ₀为达西渗透率，τ为孔隙内流体的弯曲度。

在得到三种弹性波的波数后，可计算弹性波的速度频散和衰减，具体表达式为(Tang and Cheng，2004)：

其中v与Q分别为速度与品质因子，Re{k}与Im{k}分别取k的实部和虚部。

(2)多形态裂隙并存的孔裂隙理论：现在考虑岩石中存在多个形态裂隙的情况，在Tang等(2012)提出的仅含有单一形态裂隙的孔裂隙介质中，孔隙流体压力本构方程为：

其中w＝φ(U-u)，U和u分别为流体和固体的位移，q为单位体积岩石由于裂隙挤喷流入孔隙空间的流体体积。当M个形态(用纵横比衡量)不一的裂隙并存时，将各个裂隙挤喷的贡献相加，(10)式的压力本构方程变为：

其中下标m代表第m种纵横比(形态)的裂隙，用S_m(ω)＝φq_m/p表示该种裂隙的挤喷函数，(11)式可化简为：

其中第m种裂隙的S_m(ω)由(4)式给出，其具体表达式为：

其中

ε_m和γ_m分别代表第m种形态裂隙的裂隙密度与纵横比。

基于单一裂隙的公式(1)和(2)，多裂隙体系岩石的体积模量、剪切模量的计算公式变为：

在孔裂隙理论计算公式(5)～(9)式中，将K^new与μ^new代替原来的K^Tang与μ^Tang即可。

所述步骤三具体为：

压力是影响孔隙结构的重要因素，岩石中裂隙的纵横比随着压力的增大而减小，形态变得扁平，随着压力继续增大，裂隙逐渐闭合；而压力的改变对硬孔的形态没有太大影响。为了分析这一现象，

等(1976)将硬孔和裂隙视为具有一定纵横比的扁球形包含体，并给出了包含体的体积变化率和压力的关系，对于纵横比为γ的扁球形裂隙，其体积变化率与有效压力的关系为：

其中P_e为有效压力，K^*为干岩石的静态有效模量，E_i(i＝1,2,3,4)为纵横比γ与岩石基质模量K和μ的函数，E_i的表达式为：

其中

等(1976)建议K^*、K和μ用干燥岩石动态有效模量代替，本专利采用实验测量的干燥岩石纵、横波速度计算干燥岩石动态有效模量。

Cheng(1978)和Cheng和

(1979)建立了压力与孔隙结构的函数关系，压力为P_n时的裂隙体积含量c_nm(即裂隙孔隙度)可以由0有效压力下的裂隙孔隙度c_0m表示：

其中下标n代表第n个压力点，下标m代表第m个形态的纵横比。当

时，裂隙闭合。

对于压力作用下的扁球形裂隙，其短轴的变化远大于长轴变化，在忽略后一变化的条件下，纵横比的变化率即为孔隙体积变化率：

因此，压力P_n下的裂隙纵横比γ_nm可以由0有效压力下的裂隙纵横比γ_0m表示：

结合(18)～(20)式，得到：

因此，任意压力点下的孔隙纵横比谱(纵横比与孔隙度)可以由0有效压力的孔隙纵横比谱计算得到。

扁球形裂隙的裂隙密度与裂隙孔隙度的关系为：

由(22)式，可以将裂隙孔隙度转化为裂隙密度。

所述步骤四具体为：

(1)首先介绍仅反演饱和数据的方法。采用实验室测量的饱和岩石纵、横波速度反演岩石的裂隙密度谱与纵横比谱，而后利用(22)式将裂隙密度转化为裂隙孔隙度，最终得到孔隙纵横比谱。在反演中，将理论计算速度与实验速度的误差平方和作为目标函数，调节参数使目标函数达到最小值，目标函数的具体表达式如下：

其中运算符Min表示取函数的最小值；

与

分别为实验测得的饱和岩石的纵、横波速，

与

分别为理论计算得到的饱和岩石的纵、横波速度；ε_N×M与γ_N×M分别为裂隙密度与纵横比，下标N代表压力加载N次，下标M代表0有效压力时共存在M个形态的裂隙，考虑到裂隙随着压力的增大逐渐闭合，因此每个压力点的裂隙形态总数不会大于M；φ_p为岩石除去裂隙部分的孔隙度，这部分孔隙对应的纵横比大于0.01；K_s和μ_s分别为岩石基质的体积模量和剪切模量；流体体积模量K_f在实验中是已知的。(23)式所描述的反演是多参数反演，对于N个压力测量点的纵、横波速度数据，最多需要反演2×M×N+3个未知参数。

下面分析压力对裂隙形态的影响，并简化目标函数。如果已知0有效压力下的孔隙纵横比谱，任意压力点下的孔隙纵横比谱均可由其计算。因此，在反演中计算

与

时采用的各压力点的孔隙纵横比谱均可由0有效压力下的孔隙纵横比谱计算。由此(23)式可化简为：

(24)式只需要反演0有效压力下的孔隙纵横比谱，再由反演得到的0有效压力下的孔隙纵横比谱计算各压力点下的孔隙纵横比谱，此时需要反演的参数减少为2×M+3个。因此，裂隙体积变化理论(

等，1976)建立裂隙形态与压力关系的同时也减少了反演的参数个数，降低了反演的不确定性。

(2)对于采用饱和、干燥条件的纵、横波速度数据联合反演岩石的孔隙纵横比谱，可以将目标函数(24)式变为：

其中实验速度(

与

)与理论速度(

与

)均包括饱和与干燥两种状态，分别对应下标sat与dry。其中饱和条件下的纵、横波速度由(9)式计算，干燥条件下的纵、横波速度的计算公式如下：

其中K_d与μ₀由Biot相恰理论(Thomsen，1985)计算，ρ_d为干燥介质的密度。将(25)式作为目标函数代替(24)式，其他步骤不变，即可实现采用饱和、干燥条件的纵、横波速度数据联合反演岩石的孔隙纵横比谱。

步骤五、设置多种形态的裂隙(包含孔隙)，其中孔隙的纵横比大于0.01(通常不超过2种)，裂隙的纵横比小于0.01。下面给出0有效压力下裂隙纵横比的取值及其数量M的确定方法。裂隙纵横比的取值及其数量M的取值受到裂隙闭合、测量压力范围、测量压力点数N以及岩性的影响，为了达到最优的拟合效果，最好令每个压力点下都有裂隙闭合(Cheng，1978；Yan等，2014)。对于裂隙纵横比的取值，应满足如下条件：1、由于把纵横比为0.01的孔隙作为裂隙与硬孔的分界，在孔裂隙理论中硬孔不易变形，且对硬孔纵横比没有限定，因此硬孔的纵横比可在0.01至1之间任取；2、裂隙纵横比的设置要考虑压力的影响，基于(22)式，可以计算出0有效压力下纵横比的取值范围，以此保证每个压力点下都有裂隙闭合，且设置的裂隙纵横比应较均匀地分布；3、考虑到岩性的差异，致密岩石(或高压段速度基本不变的岩石)的裂隙纵横比应该较小，通常要小于高孔隙岩石的纵横比。对于裂隙数量M的确定，应满足如下条件：1、考虑到每个压力点下都有裂隙闭合，M的数量应通常应大于N；2、不同岩性的M有所不同，致密岩石(或高压段速度基本不变的岩石)的M个数应较少，甚至可以小于N。为了提高反演结果的准确性，将0有效压力的总孔隙度作为约束条件：

在孔裂隙理论中，硬孔的孔隙度远大于裂隙的孔隙度，因此对于(27)式，可令0压下的硬孔孔隙度φ_p保持不变，反演时调节0压下的裂隙密度，使0压下的裂隙孔隙度与φ_p之和满足(27)式。

步骤六、反复调节岩石在0有效应力下的各形态裂隙的纵横比与裂隙密度，使目标函数达到最小值，就得到了岩石在各个压力点下的孔隙纵横比谱。孔隙结构的反演可以看作是搜寻目标函数最小值的过程，为了降低反演结果的非唯一性，采取全局优化的GA算法(Goldberg，1988)对目标函数进行求解。将反演得到的孔隙纵横比谱代入到步骤二的多形态裂隙并存的孔裂隙理论,计算和预测岩石样品的纵横波波速随压力的变化关系。通过对比理论计算得到的纵横波速度与实验超声测量得到的纵横波速度的吻合程度,可以对反演得到的孔隙纵横比谱进行质量监控和验证。

以下，结合两种经典岩石实验数据的实例处理成果，进一步说明本发明所述的利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法的应用效果。实例1代表采用饱和数据反演的结果；实例2代表采用饱和、干燥数据联合反演的结果。

实施实例1：Kayenta砂岩(饱和的纵、横波速度)

首先用反演方法处理Kayenta砂岩(Coyner，1984)的岩芯测量数据，Kayenta砂岩为高孔隙度岩石，孔隙度为22.2％，骨架密度为2017kg/m³，渗透率取为63mD，声速测量的频率为1MHz，实验的有效压力范围为1MPa至70Mpa，图2-a给出了该砂岩样品在苯饱和条件下纵、横波速随压力变化的测量数据(圆形标识符)。反演得到的岩石基质的体积模量和剪切模量分别为31.00GPa和24.04GPa。

在图2-a中，横坐标为有效压力(MPa)，纵坐标为速度(m/s)，黑色圆圈为测量的饱和条件下的纵、横波速度；黑色实、虚线代表用反演的岩石弹性模量和裂隙纵横比谱计算的饱和纵、横波速度，与实测数据吻合很好。

图2-b是反演得到的0有效压力下的孔隙纵横比谱，其中，图2-b的横坐标为纵横比，纵坐标为孔隙度；裂隙纵横比分布范围在3.5×10^-4至1×10^-2。

Burns等(1990)采用Cheng方法反演了Kayenta砂岩的孔隙纵横比谱，还利用扫描电镜技术给出了岩石的孔隙纵横比谱。

图2-c给出了本发明反演的结果与以上两种结果的对比，其中，图2-c的横坐标为纵横比，纵坐标为归一化孔隙度，小圆圈代表本发明反演的结果，大圆圈代表Cheng方法的结果，直方图代表扫描电镜的结果。相比于Cheng方法的反演结果，本发明的反演结果与扫描电镜结果更加吻合，结合图2-a中用反演结果计算的速度与实测速度的良好吻合关系，证明了利用孔隙、裂隙介质弹性波理论可以很好地反演岩石的孔隙结构。

实施实例2：Navajo砂岩(联合饱和、干燥的纵、横波速度)

下面用本发明反演方法处理Navajo砂岩的岩芯测量数据(Coyner，1984)，Navajo砂岩较为致密，孔隙度为11.8％，骨架密度为2316kg/m³，声速测量的频率为1MHz，实验的有效压力范围为1MPa至100MPa。图3-a给出了在苯饱和、干燥条件下纵、横波速随压力变化的测量数据(圆形标识符)。与实例1不同的是，我们采取饱和、干燥的纵、横波速度联合反演。特别地，对于砂岩，为了更好地拟合数据，我们在反演中针对饱和、干燥条件两种条件，分别设置两个岩石基质剪切模量，反演得到的岩石基质的体积模量、饱和(干燥)时基质剪切模量分别为36.60Gpa、35.00(32.50)Gpa。

在图3-a中的横坐标为有效压力(MPa)，纵坐标为速度(m/s)，黑色圆圈为测量的饱和条件下的纵、横波速度，白色圆圈为干燥条件下的纵、横波速度；黑色实线代表用反演的岩石弹性模量和裂隙纵横比谱计算的饱和纵、横波速度；黑色虚线代表用反演的岩石弹性模量和裂隙纵横比谱计算的干燥纵、横波速度。图3-a中用反演的岩石弹性模量和裂隙纵横比谱计算的饱和、干燥的纵横波速度，与实测数据吻合很好。

图3-b是反演得到的0有效压力下的孔隙纵横比谱，其中，图3-b的横坐标为纵横比，纵坐标为孔隙度；裂隙纵横比分布范围在8×10^-5至6×10^-3。

Burns等(1990)采用Cheng方法和扫描电镜技术给出了Navajo砂岩的孔隙纵横比谱。

图3-c给出了本发明反演的结果与以上两种结果的对比，其中，图3-c的横坐标为纵横比，纵坐标为归一化孔隙度，小圆圈代表本发明反演的结果，大圆圈代表Cheng方法的结果，直方图代表扫描电镜的结果。本发明反演的结果与扫描电镜结果更加吻合。

反演实例所显示的良好的吻合性证明了本方法的可靠性、准确性以及其广阔的应用前景。

附参考文献

Adam L,Batzle M,Brevik I.2006.Gassmann’s fluid substitution and shearmodulus variability in carbonates at laboratory seismic and ultrasonicfrequencies.Geophysics,71(6),F173–F183.

Burns D R,Cheng C H,Wilkens R H.1990.Sandstone pore aspect ratio fromdirect observations and velocity inversions.International Journal of RockMechanics and Mining Sciences and Geomechanics Abstracts,21:315–323.

Cheng C H.1978.Seismic velocities in porous rocks—Direct and inverseproblems[Ph.D.thesis].Boston:Massachusetts Institute of Technology.

Cheng C H,

M N.1979.Inversion of seismic velocities for the poreaspect ratio spectrum of a rock.J.Geophys.Res.,84,7533–7543.

Coyner K B.1984.Effects of stress,pore pressure and pore fluids onbulk strain,velocity and permeability in rocks[Ph.D.thesis].Boston:Massachusetts Institute of Technology.

David E C,Zimmerman R W,2012,Pore structure model for elastic wavevelocities in fluid-saturated sandstones.J.Geophys.Res.,117,B07210–B07224.

Duan C,Deng J,Li Y,et al.2018,Effect of pore structure on thedispersion and attenuation of fluid-saturated tight sandstone.J.Geophys.Eng.,15,449-460.

Goldberg D E.1988.Genetic Algorithms in Search,Optimization andMachine Learning.Boston:Addison-Wesley Publishing Co.

Han X,Tang G,He Y,et al.2019.Effect of pore aspect ratio spectrum onultrasonic velocities of tight sandstones saturated with differentfluids.89th Ann.Internat.Mtg.,Soc.Expi.Geophys..Expanded Abstracts:3588-3592.

Hadley K.1976.Comparison of calculated and observed crack densitiesand seismic velocities in Westerly granite.J.Geophys.Res.,81(20):3484-3494.

Johnson D L,J Koplik,R Dashen.1987.Theory of dynamic permerabilityand tortuosity in fluid saturated porous media.Journal of Fluid Mechanics,176,379-402.

King M S and Marsden J R.2002.Velocity dispersion between ultrasonicand seismic frequencies in brine-saturated reservoir sandstones.Geophysics,67:254-258.

Kuster G T,

M N.1974.Velocity and attenuation of seismic waves intwo-phase media:part I.Theoretical formulations.Geophysics,39:587-606.

Tang X M,Cheng C H.2004.Quantitative borehole acousticmethods.Elsevier Science Publishing Company Inc.

Tang X M,Chen X L,Xu X K.2012.A cracked porous medium elastic wavetheory and its application to interpreting acoustic data from tightformation.Geophysics,77(6):D245-D252.

Thomsen L.1985.Biot-consistent elastic moduli of porous rocks:Low-frequency limit.Geophysics,50,2797-2807.

M N,Cheng C H,Timur A.1976.Velocities of seismic waves inporous rocks.Geophysics,41,621-645.

Yan F,Han D-H,Yao Q,et al.2014.Prediction of seismic wave dispersionand attenuatio n from ultrasonic velocity measurements.Geophysics,79:WB1-WB8.

Yan F,Han D-H,Chen X L.2015.Pore aspect ratio spectrum inversion fromultrasonic m-easurements and its application.Journal of ComputationalAcoustics,23(04):1540009.

附中文参考文献

邓继新,周浩,王欢等.2015.基于储层砂岩微观孔隙结构特征的弹性波频散响应分析.地球物理学报,58(09):3389-3400.

李闯,赵建国,王宏斌等.2020.致密碳酸盐岩跨频段岩石物理实验及频散分析.地球物理学报,63(02):627-637.

唐晓明.含孔隙、裂隙介质弹性波动的统一理论—Biot理论的推广.中国科学:地球科学,2011,41(06):784-795.

唐晓明,钱玉萍,陈雪莲.孔隙、裂隙介质弹性波理论的实验研究.地球物理学报,2013,56(12):4226-4233.

Claims (5)

1.一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，包括如下处理步骤：

步骤一、测量岩心在不同压力下干燥与饱含流体的速度，并确定包括岩石密度、岩石渗透率、岩石孔隙度、流体密度和流体模量的性质；

步骤二、用多形态裂隙的孔裂隙理论模拟和计算岩石的弹性波速度；具体如下：

(1)孔、裂隙并存的双孔介质弹性波理论给出饱和岩石的体积模量、剪切模量的表达式为：

其中K_d为干燥岩石的体积模量，α＝1-K_d/K_s，β＝(α-φ)/K_s+φ/K_f，K_s为岩石基质的体积模量，φ为岩石的孔隙度，S(ω)为描述孔隙与裂隙相互作用的挤喷流函数，包含了裂隙密度和裂隙纵横比这两个描述裂隙的重要参数，K₀、μ₀分别为S(ω)＝0时的饱和岩石体积模量与剪切模量；

针对硬币型的裂隙，推导出S(ω)的表达式为：

其中

ω为角频率，η为孔隙流体黏度，ε为裂隙密度，γ为裂隙纵横比，K₀、μ₀为s(ω)＝0时背景介质的体积模量、剪切模量，v₀为干燥介质的泊松比，K_d、μ₀和v₀可由Biot相洽理论计算，K₀可由Gassmann方程计算，

为第一类n*阶贝塞尔函数，其中n*＝0，1；

硬币模型中孔隙与裂隙的流体交换在硬币的边缘，但力学上裂隙在此是闭合的，作为模型的改进，钹状的孔、裂隙模型将流体交换放到硬币模型的中部，由此推导出的S(ω)表达式为：

其中

通过调节与模型有关的流体挤喷驰豫频率，可以使二模型的计算结果一致，因此这两种模型都可以用作以下反演模型的正演计算；采用钹状模型及其理论公式，具体公式如(5)～(9)式；

孔、裂隙介质中的快纵波、慢纵波和横波的波数由下式计算：

其中下标p和s分别代表纵波和横波，+和-分别代表快纵波和慢纵波，上述公式的符号表达式为：

其中ρ_s与ρ_f分别为岩石固体基质和孔隙流体的密度，θ＝iκ(ω)/(ηω)，其中κ(ω)为动态渗透率，κ(ω)的具体表达式为：

其中κ₀为达西渗透率，τ为孔隙内流体的弯曲度；

在得到三种弹性波的波数后，可计算弹性波的速度频散和衰减，具体表达式为：

其中v与Q分别为速度与品质因子，Re{k}与Im{k}分别取k的实部和虚部；

(2)多形态裂隙并存的孔裂隙理论：仅含有单一形态裂隙的孔裂隙介质中，孔隙流体压力本构方程为：

其中w＝φ(U-u)，U和u分别为流体和固体的位移，q为单位体积岩石由于裂隙挤喷流入孔隙空间的流体体积；当M个形态不一的裂隙并存时，将各个裂隙挤喷的贡献相加，(10)式的压力本构方程变为：

其中下标m代表第m种纵横比的裂隙，用S_m(ω)＝φq_m/p表示该种裂隙的挤喷函数，(11)式可化简为：

其中第m种裂隙的S_m(ω)由(4)式给出，其具体表达式为：

其中

ε_m和γ_m分别代表第m种形态裂隙的裂隙密度与纵横比；

基于孔裂隙介质弹性模量的计算公式(1)和(2)，多裂隙体系岩石的体积模量、剪切模量的计算公式变为：

在孔裂隙理论计算公式(5)～(9)式中，将K^new与μ^new代替原来的K^Tang与μ^Tang即可计算多裂隙体系的情况；

步骤三、基于压力对孔隙结构影响的计算公式，对于岩石中存在多种形态的裂隙，不同压力点下的孔隙纵横比谱可以由0有效压力的孔隙纵横比谱计算；

步骤四、将实验测量的不同压力点的岩石速度与由步骤二理论计算的各压力点的速度做最小二乘拟合，其方差作为目标函数，并利用步骤三的压力对孔隙结构影响的计算公式，化简目标函数；

步骤五、设置多种形态的裂隙，包含孔隙，其中孔隙的纵横比大于0.01，并且不超过2种，裂隙的纵横比小于0.01；

步骤六、反复调节岩石在0有效应力下的各形态裂隙的纵横比与裂隙密度，使目标函数达到最小值，就得到了岩石在各个压力点下的孔隙纵横比谱；孔隙结构的反演可以看作是搜寻目标函数最小值的过程，为了降低反演结果的非唯一性，采取全局优化的GA算法对目标函数进行求解；将反演得到的孔隙纵横比谱代入到步骤二的多形态裂隙并存的孔裂隙理论，计算和预测岩石样品的纵横波波速随压力的变化关系；通过对比理论计算得到的纵横波速度与实验超声测量得到的纵横波速度的吻合程度，可以对反演得到的孔隙纵横比谱进行质量监控和验证。

2.根据权利要求1所述的一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，所述不同压力点的岩石速度可以仅为饱和条件下的纵、横波速度，或为联合饱和、干燥条件下的纵、横波速度。

3.根据权利要求1所述的一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，所述步骤三具体为：

将硬孔和裂隙视为具有一定纵横比的扁球形包含体，并给出了包含体的体积变化率和压力的关系，对于纵横比为γ的扁球形裂隙，其体积变化率与有效压力的关系为：

其中P_e为有效压力，K^*为干岩石的静态有效模量，E_i为纵横比γ与岩石基质模量K和μ的函数，其中i＝1，2，3，4，E_i的表达式为：

其中

K^*、K和μ用干燥岩石动态有效模量代替，以下反演计算中采用实验测量的干燥岩石纵、横波速度计算干燥岩石动态有效模量；

压力为P_n时的裂隙体积含量c_nm，即裂隙孔隙度，可以由0有效压力下的裂隙孔隙度c_0m表示：

其中下标n代表第n个压力点，下标m代表第m种纵横比的裂隙；当

时，裂隙闭合；

对于压力作用下的扁球形裂隙，其短轴的变化远大于长轴变化，在忽略后一变化的条件下，纵横比的变化率即为孔隙体积变化率：

因此，压力P_n下的裂隙纵横比γ_nm可以由0有效压力下的裂隙纵横比γ_0m表示：

结合(18)～(20)式，得到：

因此，任意压力点下的孔隙纵横比谱，即纵横比与孔隙度，可以由0有效压力的孔隙纵横比谱计算得到；

扁球形裂隙的裂隙密度与裂隙孔隙度的关系为：

由(22)式，可以将裂隙孔隙度转化为裂隙密度。

4.根据权利要求3所述的一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，所述步骤四具体为：

(1)对于仅反演饱和数据的方法，采用实验室测量的饱和岩石纵、横波速度反演岩石的裂隙密度谱与纵横比谱，，而后利用(22)式将裂隙密度转化为裂隙孔隙度，最终得到孔隙纵横比谱；

在反演中，将理论计算速度与实验速度的均方差作为目标函数，调节参数使目标函数达到最小值，目标函数的具体表达式如下：

其中运算符Min表示取函数的最小值；

与

分别为实验测得的饱和岩石的纵、横波速，

与

分别为理论计算得到的饱和岩石的纵、横波速度；ε_N×M与γ_N×M分别为裂隙密度与纵横比，下标N代表压力加载N次，下标M代表0有效压力时共存在M个形态的裂隙，考虑到裂隙随着压力的增大逐渐闭合，因此每个压力点的裂隙形态总数不会大于M；φ_p为岩石除去裂隙部分的孔隙度，这部分孔隙对应的纵横比大于0.01；K_s和μ_s分别为岩石基质的体积模量和剪切模量；流体体积模量K_f在实验中是已知的；(23)式所描述的反演是多参数反演，对于N个压力测量点的纵、横波速度数据，最多需要反演2×M×N+3个未知参数；

分析压力对裂隙形态的影响，并简化目标函数；如果已知0有效压力下的孔隙纵横比谱，任意压力点下的孔隙纵横比谱均可由其计算；因此，在反演中计算

与

时采用的各压力点的孔隙纵横比谱均可由0有效压力下的孔隙纵横比谱计算；由此(23)式可化简为：

(24)式只需要反演0有效压力下的孔隙纵横比谱，再由反演得到的0有效压力下的孔隙纵横比谱计算各压力点下的孔隙纵横比谱，此时需要反演的参数减少为2×M+3个；因此，裂隙体积变化理论建立裂隙形态与压力关系的同时也减少了反演的参数个数，降低了反演的不确定性；

(2)对于采用饱和、干燥条件的纵、横波速度数据联合反演岩石的孔隙纵横比谱，可以将目标函数(24)式变为：

其中实验速度

与

与理论速度

与

均包括饱和与干燥两种状态，分别对应下标sat与dry；其中饱和条件下的纵、横波速度由(9)式计算，干燥条件下的纵、横波速度的计算公式如下：

其中K_d与μ₀由Biot相恰理论计算，ρ_d为干燥介质的密度；将(25)式作为目标函数代替(24)式，其他步骤不变，即可实现采用饱和、干燥条件的纵、横波速度数据联合反演岩石的孔隙纵横比谱。

5.根据权利要求3所述的一种利用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征的方法，所述步骤五具体为：

0有效压力下裂隙纵横比的取值及其数量M的确定方法为：裂隙纵横比的取值及其数量M的取值受到裂隙闭合、测量压力范围、测量压力点数N以及岩性的影响，对于裂隙纵横比的取值，应满足如下条件：(1)由于把纵横比为0.01的孔隙作为裂隙与硬孔的分界，在孔裂隙理论中硬孔不易变形，且对硬孔纵横比没有限定，因此硬孔的纵横比在0.01至1之间任取；(2)裂隙纵横比的设置要考虑压力的影响，基于(22)式，计算出0有效压力下纵横比的取值范围，以此保证每个压力点下都有裂隙闭合，且设置的裂隙纵横比应较均匀地分布；(3)考虑到岩性的差异，致密岩石或高压段速度基本不变的岩石的裂隙纵横比应小于高孔隙岩石的纵横比；

对于裂隙数量M的确定，应满足如下条件：(1)考虑到每个压力点下都有裂隙闭合，M的数量通常应大于N；(2)不同岩性的M有所不同，致密岩石或高压段速度基本不变的岩石的M个数应小于高孔隙岩石的M；

为了提高反演结果的准确性，将0有效压力的总孔隙度作为约束条件：

在孔裂隙理论中，硬孔的孔隙度远大于裂隙的孔隙度，因此对于(27)式，可令0有效压力下的硬孔孔隙度φ_p保持不变，反演时调节0有效压力下的裂隙密度，使0有效压力下的裂隙孔隙度与φ_p之和满足(27)式。

Images