技术领域
[0001] 本
发明属于
岩石孔隙度
地震反演技术领域，具体为一种提高分辨率的孔隙度反演方法。
背景技术
[0002] 孔隙度是控制岩石或沉积体中孔隙
流体存储和运输的重要因素之一，其被定义为岩石中孔隙体积与岩石总体积的比值。储层预测中岩石孔隙对
碳氢化合物的实际存储能
力被定义为有效孔隙度。孔隙度预测在含油气储层的岩性判识、流体识别及储量估算等方面有重要作用。
[0003] 现有孔隙度预测方法多是基于岩石物理理论，以等效介质思想为
桥梁，建立孔隙度与实际地震数据之间的定量关系，实现孔隙度的间接定量反演与预测。然而，该方法对孔隙度的预测能力有限，孔隙度反演结果横向分辨率差，特别是针对复杂的含油高孔
砂岩储层或碳酸盐岩储层。因此，本发明领域需要研发一种方法，建立孔隙度与实际地震数据之间的直接定量关系，实现孔隙度的直接定量反演，提高孔隙度反演分辨率，为精确储层预测和流体识别提供理论和方法技术
支撑。
发明内容
[0004] 为解决
现有技术中孔隙度反演分辨率差的技术问题，本发明提供了一种提高分辨率的孔隙度反演方法，其实现的目的为，提供了一种基于含饱和流体孔隙介质地震平面波近似解析理论的高分辨率孔隙度反演方法。假定存在待反演地震数据、
测井反射系数及测井孔隙度，建立孔隙度与实际地震数据之间的直接定量关系，提高孔隙度反演分辨率。
[0005] 为了实现上述目的，本发明采用的技术方案为，假定存在待反演地震数据、测井反射系数及测井孔隙度，利用模拟
退火算法实现高
精度孔隙度反演的流程如下，[0006] (1)推导垂直入射至含流体孔隙介质分界面上的快纵波反射系数近似式，基于岩石物理理论建立孔隙度和地震记录之间的直接定量关系；
[0007] (2)根据孔隙度和地震记录之间的关系，以测井孔隙度和测井反射系数作为先验约束，构建非线性反演目标函数；
[0008] (3)利用测井数据信息计算孔隙度调节参数，给定约束参数、最大反演
迭代次数及孔隙度初始解等先验信息；
[0009] (4)将上述已知参数带入非线性反演目标函数，结合快速模拟退火算法，实现孔隙度的高精度直接定量反演。
[0010] 进一步的，1、孔隙度与地震数据之间的定量关系
[0011] 基于Biot弹性孔隙介质理论，快纵波垂直入射到含饱和流体孔隙介质分界面上的快纵波反射系数近似式为，
[0012]
[0013]
[0014]
[0015]
[0016]
[0017]
[0018]
[0019] 其中，r为快纵波反射系数，Kl,μl,l＝1,2分别表示岩石骨架的体积模量和
剪切模量，αl,Ml,l＝1,2为Biot参数；vplj,l＝1,2,j＝1,2为纵
波速度，γlj,l＝1,2,j＝1,2为孔隙流体相对于岩石骨架的位移位振幅和岩石骨架的位移位振幅之比；下标l＝1,2分别表示上、下层孔隙介质，j＝1,2分别表示快、慢纵波。
[0020] 岩石骨架模量、岩石基质模量和孔隙度之间存在一定联系，本发明利用 (8)和(9)式描述他们之间的定量关系。
[0021] Kd＝(1.0-f)Ks/(1.0+ck·f), (8)
[0022] μd＝(1.0-f)μs/(1.0+cμ·f)， (9)
[0023] 其中，f为孔隙度，Kd、μd分别为岩石骨架的体积模量和剪切模量，Ks、μs分别为岩石基质的体积模量和剪切模量。ck、cμ为孔隙介质调节参数。
[0024] 含饱和流体孔隙介质中存在快、慢两类纵波，慢纵波在孔隙介质中传播时会迅速衰减消失，因此地面
检波器接收到的地震记录事实上是快纵波的反射。孔隙度f是时间t的函数f(t)，而由(1)-(9)式可知快纵波反射系数r是孔隙度f(t)的函数r[f(t)]。根据褶积公式，与孔隙度直接相关的快纵波合成地震记录d[f(t)]可以表示为，
[0025] d[f(t)]＝w(t)*r[f(t)]+n(t), (10)
[0026] 将上式离散化，并写成向量的形式为，
[0027] dN×1＝GN×N·rN×1+nN×1. (11)
[0028] 其中，w(t)为时间域子波，n(t)为噪声，下标N×1表示N维列向量，N×N 表示N维矩阵，N为一道地震数据的
采样点数；符号‘*’表示褶积算子。
[0029] 进一步的，2、非线性反演目标函数
[0030] 根据(11)式中的正演过程，以测井中的先验信息为约束，含流体孔隙介质孔隙度反演目标函数F为，
[0031]
[0032] 其中，f为预测孔隙度，R为预测反射系数，d为合成地震记录；f0和R0分别为由测井数据得到的孔隙度和反射系数，ddata为实际地震数据；β为反射系数约束参数，γ为孔隙度约束参数，目标函数中的所有向量均为N×1维列向量， 表示向量的L2范数。
[0033] 上式中，第一项为观测地震数据和合成地震记录的拟合程度，第二项表示由预测孔隙度计算的反射系数与测井反射系数的相似性，第三项表示预测孔隙度与测井孔隙度的相似性。后两项是对反演结果的测井约束，为反演目标函数提供了更多的先验信息，提高反演结果的
稳定性和抗噪性。目标函数中的约束参数β和γ控制着测井约束项对反演结果的贡献程度，调节该参数可以提高反演精度。
[0034] 进一步的，3、模拟退火算法孔隙度反演
[0035] 模拟退火算法源于固体退火原理，模拟了物理学中晶体的退火过程。首先将固体加温至充分高，再让其徐徐冷却。加温时，固体内部粒子随
温度升高变为无序状态，内能增大；而冷却时，粒子渐趋有序，每个粒子温度都达到平衡态；最后在常温时达到基态，内能减小为最小。算法具有渐进收敛性，是一种以概率1收敛的全局
优化算法。
[0036] (1)Metropolis收敛准则
[0037] 反演算法通过不断迭代更新孔隙度来获得最优解。在第i次迭代时，首先计算与孔隙度fi对应的目标函数Fi，再更新孔隙度为fi+1，并计算其对应的目标函数值Fi+1。两函数值之差ΔF为，
[0038] ΔF＝Fi+1-Fi,i≥1, (13)
[0039] 若ΔF＜0说明孔隙度的更新方向使目标函数值减小，则接受孔隙度的
修改；若ΔF≥0，则进一步判断ΔF是否满足(14)式，若满足条件，则仍接受孔隙度的修改，否则拒绝修改。
[0040] 0＜p＜K, (14)
[0041]
[0042] 其中，p是概率
密度，K是一个随机数，其值域为0＜K＜1；kb为玻尔兹曼常数，实际应用中常取1，exp{·}为指数函数。T为当前退火温度，常表示成指数形式，[0043] Ti＝T0·(0.95)i-1,i≥1, (16)
[0044] 其中,T0为初始温度，Ti为第i次退火温度。
[0045] (2)解的迭代与更新
[0046] 模拟退火算法的解与初始状态无关，但研究经验表明：若给定恰当的初始解，可以减少算法迭代求解次数、提高算法运行速度。本文将测井孔隙度 f0经平滑后作为初始孔隙度f1。给定初始解后，可以采用依赖于温度的似 Cauchy分布产生更新向量Δfi，即[0047] fi+1＝fi+ξ·Δfi,i≥1, (17)
[0048] Δfi＝T·sign(q-0.5)·[(1+1/T)|2·q-1|-1]. (18)
[0049] 上式中，ξ为更新系数，fi为当前预测孔隙度，fi+1为更新后的孔隙度，其值域为fi+1∈[0,0.3]；q是服从[0,1]分布的N维随机向量，sign为符号函数。
[0050] 采用依赖于温度的似Cauchy分布产生更新向量的优点是：在高温情况下对孔隙度进行大范围的搜索，而在低温时搜索尽可能在当前预测值附近进行，似Cauchy分布有一平坦的‘尾巴’易于跳出局部极小值。
[0051] (3)非线性孔隙度反演流程
[0052] 假设存在待反演地震数据ddata、测井反射系数R0、测井孔隙度f0、初始孔隙度f1、最大迭代次数L及测井约束参数β和γ，基于模拟退火算法实现高精度孔隙度反演的流程如下，
[0053] 1)对于第i次迭代的孔隙度fi，利用(17)式计算新的孔隙度值fi+1，再基于(12)式计算孔隙度更新前后反演目标函数值Fi和Fi+1；
[0054] 2)结合(13)式计算两目标函数值之差ΔF，并判断是否满足条件ΔF＜0，若条件满足则接受孔隙度值的更新，否则继续判断是否满足(14)式，若满足则仍可接受孔隙度值的更新，否则拒绝更新；
[0055] 3)判断是否达到最大迭代次数L，若是则退出循环，获得最优孔隙度值；否则利用(16)式迭代更新当前退火温度，并返回1)继续迭代更新孔隙度值 f。
[0056] 本发明采用上述技术方案，具备以下有益效果：本发明以含流体孔隙介质垂直入射地震平面波近似解析式为理论
基础，结合快速模拟退火算法实现了高精度孔隙度反演。与已有孔隙度反演技术相比，本技术突破了传统依赖于单相介质和等效介质理论的限制，从根本上建立了孔隙度与地震记录间的直接定量关系，这是本发明的理论先进性。巧妙利用模拟退火算法，将测井数据、岩石物理及地质等先验信息引入至反演目标函数中，避免了反演陷入局部极小值，降低了反演的多解性，最终实现了高精度孔隙度反演。
附图说明
[0058] 图2为单道砂岩孔隙介质模型，(a)，(b)，(c)，(d)分别为模型对应的孔隙度、反射系数、合成地震记录及低频孔隙度；
[0059] 图3为不同初始解情况下孔隙度的
能量比；
[0060] 图4为不同初始解情况下孔隙度反演结果，(a)，(b)，(c)，(d)分别为第100、200、300、400次迭代反演；
[0061] 图5为不同初始解情况下反射系数反演结果，(a)，(b)，(c)，(d)分别为第100、200、300、400次迭代反演；
[0062] 图6为不同初始解情况下合成地震记录，(a)，(b)，(c)，(d)分别为第100、200、300、400次迭代反演结果重构的地震记录；
[0063] 图7参数β为0.3，γ分别取0.005(实线)、0.01(虚线)、0.1(点划线) 及0.5时对应的孔隙度能量比；
[0064] 图8为参数β取0.3，γ分别取0.005(a)、0.01(b)、0.1(c)及0.5 (d)时对应的孔隙度反演结果(虚线)与原始模型(实线)对比图；
[0065] 图9实际地震数据及孔隙度反演结果.(a)实际地震数据，(b)低频孔隙度模型，(c)传统方法反演的孔隙度，(d)本发明方法反演的孔隙度。
具体实施方式
[0066] 下面通过具体的
实施例对本发明做进一步的详细描述。
[0067] 实施例：参照图1所示，为本发明方法流程图，假定存在待反演地震数据、测井反射系数及测井孔隙度，利用模拟退火算法实现高精度孔隙度反演的流程如下，[0068] (1)推导垂直入射至含流体孔隙介质分界面上的快纵波反射系数近似式，基于岩石物理理论建立孔隙度和地震记录之间的直接定量关系；
[0069] (2)根据孔隙度和地震记录之间的关系，以测井孔隙度和测井反射系数作为先验约束，构建非线性反演目标函数；
[0070] (3)利用测井数据信息计算孔隙度调节参数，给定约束参数、最大反演迭代次数及孔隙度初始解等先验信息；
[0071] (4)将上述已知参数带入非线性反演目标函数，结合快速模拟退火算法，实现孔隙度的高精度直接定量反演。
[0072] 进一步的，1、孔隙度与地震数据之间的定量关系
[0073] 基于Biot弹性孔隙介质理论，快纵波垂直入射到含饱和流体孔隙介质分界面上的快纵波反射系数近似式为，
[0074]
[0075]
[0076]
[0077]
[0078]
[0079]
[0080]
[0081] 其中，r为快纵波反射系数，Kl,μl,l＝1,2分别表示岩石骨架的体积模量和剪切模量，αl,Ml,l＝1,2为Biot参数；vplj,l＝1,2,j＝1,2为纵波速度，γlj,l＝1,2,j＝1,2为孔隙流体相对于岩石骨架的位移位振幅和岩石骨架的位移位振幅之比；下标l＝1,2分别表示上、下层孔隙介质，j＝1,2分别表示快、慢纵波。
[0082] 岩石骨架模量、岩石基质模量和孔隙度之间存在一定联系，本发明利用 (8)和(9)式描述他们之间的定量关系。
[0083] Kd＝(1.0-f)Ks/(1.0+ck·f), (8)
[0084] μd＝(1.0-f)μs/(1.0+cμ·f)， (9)
[0085] 其中，f为孔隙度，Kd、μd分别为岩石骨架的体积模量和剪切模量，Ks、μs分别为岩石基质的体积模量和剪切模量。ck、cμ为孔隙介质调节参数。
[0086] 含饱和流体孔隙介质中存在快、慢两类纵波，慢纵波在孔隙介质中传播时会迅速衰减消失，因此地面检波器接收到的地震记录事实上是快纵波的反射。孔隙度f是时间t的函数f(t)，而由(1)-(9)式可知快纵波反射系数r是孔隙度f(t)的函数r[f(t)]。根据褶积公式，与孔隙度直接相关的快纵波合成地震记录d[f(t)]可以表示为，
[0087] d[f(t)]＝w(t)*r[f(t)]+n(t), (10)
[0088] 将上式离散化，并写成向量的形式为，
[0089] dN×1＝GN×N·rN×1+nN×1. (11)
[0090] 其中，w(t)为时间域子波，n(t)为噪声，下标N×1表示N维列向量，N×N 表示N维矩阵，N为一道地震数据的采样点数；符号‘*’表示褶积算子。
[0091] 进一步的，2、非线性反演目标函数
[0092] 根据(11)式中的正演过程，以测井中的先验信息为约束，含流体孔隙介质孔隙度反演目标函数F为，
[0093]
[0094] 其中，f为预测孔隙度，R为预测反射系数，d为合成地震记录；f0和R0分别为由测井数据得到的孔隙度和反射系数，ddata为实际地震数据；β为反射系数约束参数，γ为孔隙度约束参数，目标函数中的所有向量均为N×1维列向量， 表示向量的L2范数。
[0095] 上式中，第一项为观测地震数据和合成地震记录的拟合程度，第二项表示由预测孔隙度计算的反射系数与测井反射系数的相似性，第三项表示预测孔隙度与测井孔隙度的相似性。后两项是对反演结果的测井约束，为反演目标函数提供了更多的先验信息，提高了反演结果的稳定性和抗噪性。目标函数中的约束参数β和γ控制着测井约束项对反演结果的贡献程度，调节该参数可以提高反演精度。
[0096] 进一步的，3、模拟退火算法孔隙度反演
[0097] 模拟退火算法源于固体退火原理，模拟了物理学中晶体的退火过程。首先将固体加温至充分高，再让其徐徐冷却。加温时，固体内部粒子随温度升高变为无序状态，内能增大；而冷却时，粒子渐趋有序，每个粒子温度都达到平衡态；最后在常温时达到基态，内能减小为最小。算法具有渐进收敛性，是一种以概率1收敛的全局优化算法。
[0098] (1)Metropolis收敛准则
[0099] 反演算法通过不断迭代更新孔隙度来获得最优解。在第i次迭代时，首先计算与孔隙度fi对应的目标函数Fi，再更新孔隙度为fi+1，并计算其对应的目标函数值Fi+1。两函数值之差ΔF为，
[0100] ΔF＝Fi+1-Fi,i≥1, (13)
[0101] 若ΔF＜0说明孔隙度的更新方向使目标函数值减小，则接受孔隙度的修改；若ΔF≥0，则进一步判断ΔF是否满足(14)式，若满足条件，则仍接受孔隙度的修改，否则拒绝修改。
[0102] 0＜p＜K, (14)
[0103]
[0104] 其中，p是概率密度，K是一个随机数，其值域为0＜K＜1；kb为玻尔兹曼常数，实际应用中常取1，exp{·}为指数函数。T为当前退火温度，常表示成指数形式，[0105] Ti＝T0·(0.95)i-1,i≥1, (16)
[0106] 其中,T0为初始温度，Ti为第i次退火温度。
[0107] (2)解的迭代与更新
[0108] 模拟退火算法的解与初始状态无关，但研究经验表明：若给定恰当的初始解，可以减少算法迭代求解次数、提高算法运行速度。本文将测井孔隙度 f0经平滑后作为初始孔隙度f1。给定初始解后，可以采用依赖于温度的似 Cauchy分布产生更新向量Δfi，即[0109] fi+1＝fi+ξ·Δfi,i≥1, (17)
[0110] Δfi＝T·sign(q-0.5)·[(1+1/T)|2·q-1|-1]. (18)
[0111] 上式中，ξ为更新系数，fi为当前预测孔隙度，fi+1为更新后的孔隙度，其值域为fi+1∈[0,0.3]；q是服从[0,1]分布的N维随机向量，sign为符号函数。
[0112] 采用依赖于温度的似Cauchy分布产生更新向量的优点是：在高温情况下对孔隙度进行大范围的搜索，而在低温时搜索尽可能在当前预测值附近进行，似Cauchy分布有一平坦的‘尾巴’易于跳出局部极小值。
[0113] (3)非线性孔隙度反演流程
[0114] 假设存在待反演地震数据ddata、测井反射系数R0、测井孔隙度f0、初始孔隙度f1、最大迭代次数L及测井约束参数β和γ，基于模拟退火算法实现高精度孔隙度反演的流程如下，
[0115] 1)对于第i次迭代的孔隙度fi，利用(17)式计算新的孔隙度值fi+1，再基于(12)式计算孔隙度更新前后反演目标函数值Fi和Fi+1；
[0116] 2)结合(13)式计算两目标函数值之差ΔF，并判断是否满足条件ΔF＜0，若条件满足则接受孔隙度值的更新，否则继续判断是否满足(14)式，若满足则仍可接受孔隙度值的更新，否则拒绝更新；
[0117] 3)判断是否达到最大迭代次数L，若是则退出循环，获得最优孔隙度值；否则利用(16)式迭代更新当前退火温度，并返回1)继续迭代更新孔隙度值 f。
[0118] 以下列举试验用以说明采用本发明的方法进行孔隙度反演能够实现提高分辨率的技术效果。
[0119] 1、模型测试
[0120] (1)单道砂岩模型
[0121] 建立砂岩孔隙介质模型，通过分析初始解的取值和不同的约束参数，验证本发明的有效性和稳定性。图2为单道砂岩孔隙介质模型，(a)，(b)，(c)， (d)分别为模型对应的孔隙度、反射系数、合成地震记录及低频孔隙度。模型的采样率为2ms、延时为0.6s，其他参数包括：砂岩基质的体积模量、剪切模量及密度分别为37GPa、44GPa、2650kg/m3；砂岩的渗透率和孔隙因子分别为0.2darcy和2.0；孔隙介质的调节参数ck和cμ都为常数6.0。
[0122] (2)不同初始解对应的反演结果
[0123] 本发明选取低频孔隙度(图2d)和一列随机向量分别作为初始解进行孔隙度反演。同时利用预测孔隙度和测井孔隙度误差的能量与测井孔隙度的能量之比，即来描述反演结果的收敛情况。图3为不同初始解情况下孔隙度的能量比，图中黑色虚线表示初始解为低频孔隙度，黑色实线表示随机产生初始解。图4、5、6分别为不同初始解情况下反演孔隙度、反演反射系数及重构地震记录，图(a)、(b)、(c)、(d)分别表示第100、200、 300、
400次迭代反演，黑色实线表示原始模型，黑色虚线和点划线分别是初始解为低频孔隙度和随机产生时对应的反演或重构结果。
[0124] 分析图4-6可知，在不确定初始解的情况下，反演方法仍可以得到较理想的孔隙度反演结果，反演反射系数和重构地震记录分别与原始模型吻合较好，说明所提方法适用于无初始解的反演，且反演结果有效、稳定。采用随机产生的初始解进行反演时，虽然算法采用了Metropolis收敛准则避免陷入局部极小值解，但反演需要更多的迭代次数，反演速率相对较慢；而选择低频孔隙度作为初始解时，反演孔隙度不但分辨率高，而且可以快速收敛到最优解，提高了反演速度，如图3所示。
[0125] (3)不同参数对应的反演结果测试
[0126] 反演目标函数中存在两种约束参数β和γ，这些参数决定了地震记录、测井反射系数及测井孔隙度约束对反演结果的相对贡献程度。本发明分析了两种约束参数取不同值时，经过600次反演迭代后对孔隙度反演结果的影响。图8是参数β取0.3、γ分别为0.005(a)、0.01(b)、0.1(c)及0.5(d)时对应的孔隙度反演结果(虚线)，黑色实线为原始孔隙度模型。图7为与图8 中参数相对应的孔隙度能量比，其中黑色实线、虚线、点划线、短线与星号连续线分别对应参数γ的值为0.005、0.01、0.1及0.5。
[0127] 由图7可以看出，经过600次的迭代反演后，不同约束参数下对应的反演结果都收敛至孔隙度模型。但只有恰当选取测井反射系数和测井孔隙度的约束参数，才会获得最优的反演结果。随约束参数γ的逐渐增大，孔隙度反演结果中的低频成分逐渐得到补偿，最终得到最优反演孔隙度。
[0128] 2、实际地震数据反演结果
[0129] 截取过A井的部分实际地震数据，验证本发明的有效性。图9为实际地震数据及孔隙度反演结果，(a)过A井的实际地震剖面，时间域深度 3.2s-4.1s，采样间隔为4ms，CDP为400；(b)是低频孔隙度模型，由测井孔隙度、层位及地震数据构建；(c)是传统方法对应的孔隙度反演结果；(d) 为经过2100次迭代后孔隙度反演结果，反演中测井反射系数的约束参数为 0.3、测井孔隙度的约束参数为0.6。
[0130] 反演实际地震数据时，适当减小测井的约束参数可以增加实际地震数据对孔隙度反演结果的贡献程度，使反演结果更多依赖于实际地震数据、更加真实可靠。从图9可以看出：经过2100次迭代反演后，各道反演均收敛至最优解，获得了较理想的孔隙度反演结果；与传统孔隙度反演方法相比，本发明具有较高的孔隙度分辨率，且与测井孔隙度吻合程度较好，验证了本发明的有效性。孔隙度的反演对含油气储层中岩性和孔隙流体的预测提供了指导意义，为油气田的勘探开发提供了理论和数据支撑。
[0131] 以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
