一种基于模型约束的多分辨率波阻抗反演方法

摘要

本发明提供了一种模型约束的多分辨率波阻抗反演方法，属于油气及煤层气地震勘探与开发领域。本发明方法利用地震记录、测井数据及地震解释层位数据建立用于波阻抗反演的初始模型，并通过引入尺度因子加速算法的计算速度和降低反演的多解性。本发明将地面地震数据、测井数据及解释层位数据进行有机结合，对断层、楔状地层及薄层的反演具有一定的适应性，具有计算速度快、稳定性好、计算精度高以及抗噪能力强的特点。

权利要求

1.一种基于模型约束的多分辨率波阻抗反演方法，其特征在于：所述方法利用地震记录、测井数据及地震解释层位数据建立用于波阻抗反演的初始模型，并通过引入尺度因子加速算法的计算速度和降低反演的多解性；

所述方法包括以下步骤：

(1)输入叠后地震数据，从所述叠后地震数据得到地震记录的采样率为Δt，地震记录的时间长度为T,则样点个数M＝T/Δt；对所述叠后地震数据进行归一化处理后得到经过归一化处理后的叠后地震数据Data；

(2)由测井数据及地震解释层位数据建立反演的初始模型L[i],i＝1,M+1；

(3)设反演参数个数为P,尺度因子为n×Δt,n为正整数，计算能够整除尺度因子且大于地震记录时间长度T的最小整数T′，则反演参数个数P＝T′/nΔt；设经尺度扩展后的波阻抗L用向量Lp[i],i＝1,T′/Δt表示，Lp中的数值由以下公式确定：式中，L[n×i]是初始模型，表示在步骤(2)得到的L[i]上以间隔n取值；

(4)根据地震记录的质量确定控制约束范围的约束系数alphadi和alphaui，计算反演算法的约束范围[Ldi,Lui]，其中Ldi＝(1-alphadi)×Li，Lui＝(1+alphaui)×Li；所述约束系数alphadi和alphaui取值范围是0和1之间；

(5)计算反射系数Ri：

(6)从地震记录提取子波Wj,j＝1,N,并对子波的数据进行归一化处理，然后利用褶积公式计算得到合成地震记录，其中N表示子波的数量：(7)由步骤(6)得到的合成地震记录与步骤(1)得到的经过归一化处理后的叠后地震数据Data残差的二范数建立用于反演的目标函数fun：所述目标函数fun的最小值是0；

(8)计算所述目标函数在初始模型处的灵敏度向量G

其中

(9)在步骤(2)得到的初始模型及第(8)步得到的灵敏度向量的基础上利用具有全局优化能力的优化算法对步骤(7)得到的目标函数进行优化，即进行模型约束的波阻抗反演，当目标函数值达到足够小或者达到满意的值时停止算法迭代，输出反演结果。

2.根据权利要求1所述的基于模型约束的多分辨率波阻抗反演方法，其特征在于：所述步骤(2)具体如下：所述测井数据包括声波时差测井数据和密度测井数据Den，首先由所述声波时差测井数据计算得到地层速度Vp，然后根据地层速度Vp和密度测井数据Den利用波阻抗公式Z＝Vp*Den计算得到井中波阻抗曲线，再对所述井中波阻抗曲线进行带通滤波得到反演的初始模型。

3.根据权利要求1所述的基于模型约束的多分辨率波阻抗反演方法，其特征在于：所述步骤(9)中，所述具有全局优化能力的优化算法采用共轭梯度法或阻尼最小二乘法。

说明书

一种基于模型约束的多分辨率波阻抗反演方法

技术领域

[0001] 本发明属于油气及煤层气地震勘探与开发领域，具体涉及一种基于模型约束的多分辨率波阻抗反演方法。

背景技术

[0002] 地震勘探技术是油气勘探中应用最为广泛的一种地球物理学方法，利用地震波在不同介质中传播的速度、振幅、频率、相位、波形等参数的变化来分析、预测油气储层分布范围及储层物性特征。随油气勘探开发的不断推进，常规油气勘探已由原来的寻找构造油气藏转向寻找岩性油气藏和隐蔽油气藏，非常规(煤层气、页岩气和致密砂岩气等)油气勘探开发逐渐受到人们的重视，而地震反演是储层预测的有效技术之一，地震波阻抗反演是实现岩性油气藏储层预测和隐蔽油气藏识别的有效技术手段。目前波阻抗反演方法较为丰富，从地震波阻抗反演所使用的反演算法类型，可以将地震反演分为确定性反演和随机反演两大类。依反演所基于的理论(物理模型)类型，地震波阻抗反演可分为道积分法、递推法、模型法、地质统计法等。目前商用软件中主要采用道积分法、递推反演及基于模型的波阻抗反演方法。这些方法均存在一定适用条件，无法适应薄储层等复杂岩性储层的储层预测工作，具体如下：

[0003] 道积分反演是在假设地层波阻抗连续可微条件下，从平面波法线入射反射系数公式出发对反射系数(或经振幅标定后的地震记录)求积分，得到地层的相对波阻抗。该方法的特点是无需钻井控制，反演结果直接反应了地层波阻抗的相对变化，但是该方法受地震数据固有频带的限制，分辨率低，无法适于薄层解释的需要。求取的相对波阻抗和绝对波阻抗不能用于储层参数的定量计算。

[0004] 递推反演是利用反摺积从地震剖面求取反射系数剖面，进而将反射系数剖面转换为波阻抗剖面的一种转换方法，由于该方法直接将地震信息从一种形式转换到另一种形式，在很大程度上保留了叠后地震剖面的基本地质特征(断层、产状等)，能够反应岩相、岩性的空间变化，在岩性相对稳定条件下能较好的反应储层的物性变化。该方法的核心是从地震记录求取反射系数，通常由反摺积求取，选择有效的反摺积方法是递推反演的关键。由于受到地震资料频带的限制，递推反演的分辨率相对较低，不能满足薄储层的研究。

[0005] 基于模型的波阻抗反演技术将测井数据与地震数据有机结合，反演的波阻抗分辨率高，适于薄层油气藏的的精细描述。主要取决于初始模型与实际地质情况的符合程度，在钻井资料较多的情况下应用效果较好，多解性是基于模型反演的固有特性，在初始模型建立不准的情况下通常导致反演结果的错误。

[0006] 在勘探开发前期，通常使用道积分和递推反演来定性的进行岩相分析确定地层的沉积体系，结合钻井资料揭示的储层特征进行横向预测，确定评价井位。在油田开发阶段，由于钻井资料较为丰富，可以结合地震数据建立用于反演的初始模型，利用基于模型的波阻抗反演来进行定量的反演，实现储层精细描述。但是常规叠后波阻抗反演技术存在以下缺点：分辨率低，反演结果具有多解性，抗噪能力差，不适于储层参数定量计算及薄储层预测，应用受到限制等缺点。

发明内容

[0007] 本发明的目的在于解决上述现有技术中存在的难题，提供一种基于模型约束的多分辨率波阻抗反演方法，提高反演方法的抗噪能力，在约束范围内反演得到具有地质意义的高分辨率波阻抗信息，用于地震资料的岩性参数解释，便于精细储层描述，为寻找岩性(和地层)圈闭油气藏和非常规油气藏(煤层气、页岩气等)提供可靠的资料。

[0008] 本发明是通过以下技术方案实现的：

[0009] 一种基于模型约束的多分辨率波阻抗反演方法，所述方法利用地震记录、测井数据及地震解释层位数据建立用于波阻抗反演的初始模型，并通过引入尺度因子加速算法的计算速度和降低反演的多解性。

[0010] 所述方法包括以下步骤：

[0011] (1)输入叠后地震数据，从所述叠后地震数据得到地震记录的采样率Δt和时间长度T，则样点个数M＝T/Δt，对所述叠后地震数据进行归一化处理后得到经过归一化处理后的叠后地震数据Data；

[0012] (2)由测井数据及地震解释层位数据建立反演的初始模型L[i]，i＝1，M+1；

[0013] (3)设反演参数个数为P，尺度因子为n×Δt，n为正整数，计算能够整除尺度因子且大于地震记录时间长度T的最小整数T′，则反演参数个数P＝T′/nΔT；设经尺度扩展后的波阻抗L用向量Lp[i]，i＝1，T′/ΔT表示，Lp中的数值由以下公式确定：

[0014]

[0015] 式中，L[n×i]是初始模型，表示在步骤(2)得到的L[i]上以间隔n取值；

[0016] (4)根据地震记录的质量确定控制约束范围的约束系数alphadi和alphaui，计算反演算法的约束范围[Ldi，Lui]，其中Ldi＝(1-alphadi)×Li，Lui＝(1+alphaui)×Li；所述约束系数alphadi和alphaui取值范围是0和1之间；

[0017] (5)计算反射系数Ri：

[0018] (6)从地震记录提取子波Wj，j＝1，N，并对子波的数据进行归一化处理，然后利用摺积公式计算得到合成地震记录：

[0019]

[0020] (7)由步骤(6)得到的合成地震记录与步骤(1)得到的经过归一化处理后的叠后地震数据Data残差的二范数建立用于反演的目标函数fun：

[0021]

[0022] 所述目标函数fun的最小值是0；

[0023] (8)计算所述目标函数在初始模型处的灵敏度向量

[0024]

[0025] 其中

[0026] (9)在步骤(2)得到的初始模型及第(8)步得到的灵敏度向量的基础上利用具有全局优化能力的优化算法对步骤(7)得到的目标函数进行优化，即进行模型约束的波阻抗反演，当目标函数值达到足够小或者达到满意的值时停止算法迭代，输出反演结果。

[0027] 所述步骤(2)具体如下：所述测井数据包括声波时差测井数据和密度测井数据Den，首先由所述声波时差测井数据计算得到地层速度Vp，然后根据地层速度Vp和密度测井数据Den利用波阻抗公式Z＝Vp*Den计算得到井中波阻抗曲线，再对所述井中波阻抗曲线进行带通滤波得到反演的初始模型。

[0028] 所述步骤(9)中，所述具有全局优化能力的优化算法采用共轭梯度法或阻尼最小二乘法。

[0029] 与现有技术相比，本发明的有益效果是：本发明是一种适宜在油田开发阶段使用的反演方法，其结合地震解释层位数据和测井数据建立区块内波阻抗反演模型，在约束范围内反演得到了具有地质意义的高分辨率波阻抗信息，提高了反演方法的抗噪能力，得到的波阻抗信息可用于地震资料的岩性参数解释，便于精细储层描述，为寻找岩性(和地层)圈闭油气藏和非常规油气藏(煤层气、页岩气等)提供了可靠的资料。

附图说明

[0030] 图1(a)是实施例1中楔状模型合成地震记录。

[0031] 图1(b)是实施例1中在图1(a)基础上建立的初始模型。

[0032] 图1(c)是实施例1中图1(a)的第50道合成地震记录。

[0033] 图1(d)是实施例1中图1(a)的第30道合成地震记录。

[0034] 图1(e)是实施例1中图1(a)的第1道合成地震记录。

[0035] 图1(f)是实施例1中图1(b)的第50道初始模型。

[0036] 图1(g)是实施例1中图1(b)的第30道初始模型。

[0037] 图1(h)是实施例1中图1(b)的第1道初始模型。

[0038] 图1(i)是实施例1中图1(b)的第50道的反演结果与初始模型曲线，黑色粗线为初始模型，黑色细线为反演结果。

[0039] 图1(j)是实施例1中图1(b)的第30道的反演结果与初始模型曲线，黑色粗线为初始模型，黑色细线为反演结果。

[0040] 图1(k)是实施例1中图1(b)的第1道的反演结果与初始模型曲线，黑色粗线为初始模型，黑色细线为反演结果。

[0041] 图1(l)是实施例1中第50道反演结果。

[0042] 图1(m)是实施例1中第30道反演结果。

[0043] 图1(n)是实施例1中第1道反演结果。

[0044] 图1(o)是实施例1中楔状模型的反演结果。

[0045] 图2(a)是实施例2中叠后地震数据。

[0046] 图2(b)是实施例2中初始模型。

[0047] 图2(c)是实施例2中反演结果。

具体实施方式

[0048] 下面结合附图对本发明作进一步详细描述：

[0049] 为了从叠后地震数据反演得到能够用于计算储层岩性参数的绝对波阻抗，本发明利用地震数据、测井数据及地震解释层位数据建立用于波阻抗反演的初始模型，然后在初始模型的基础上构建反演算法的约束范围，由地震数据与摺积合成地震记录残差的二范数建立反演的目标函数，最后采用具有全局寻优能力的优化算法优化目标函数，得到约束范围内具有地质意义的反演结果；通过引入尺度因子(指的是步骤(3)中的尺度因子为n×Δt)加速了算法的计算速度，并在一定程度上降低了反演的多解性。

[0050] 假设地震道时间长度T为6000ms，采样率Δt为2ms，那么这个地震道的样点个数为：T/Δt＝6000/2＝3000个。进行反演时就要同时反演3000个参数，然而每个参数之间的关系不是简单的加减乘除，而是耦合的、非线性的，对于这种多参数的非线性问题的反演，在数学上就被叫做超高维最优化问题。这个问题很难解决。本发明通过引入“尺度因子”(尺度因子为n×Δt，n为正整数如1，2，3，4，5...)，比如此处将n取为5，采样率Δt为2ms，则尺度因子大小为5*2＝10ms，将地震道的长度T除以尺度因子，即6000/10＝600，在反演时反演600个参数就可以了。

[0051] 本发明方法利用摺积合成地震记录与实际地震数据的二范数构建反演的目标函数，由井中波阻抗低频模型(反演时将这个低频模型作为初始模型，低频模型是地球物理学专业的叫法，初始模型是数学上的叫法。)建立多分辨率波阻抗反演算法的搜索范围，以井中低频模型为初始解，通过“尺度因子”将超高维非线性最优化问题进行降维处理，最后采用具有全局优化性能的最优化算法对降维处理后的目标函数进行最优化处理，反演得到约束范围内具有高频信息的地震波阻抗。根据不同的地质问题来设置不同的“尺度因子”，即根据某个地区的地质特点，在对这一地区的地震资料进行反演时设定这个尺度因子的大小，算法将反演出这一尺度下的结果；当调整尺度因子大小，再次进行反演时，算法即会反演出该尺度因子下的反演结果；在反演另外一个地区地震资料时也是根据这一地区的地质特点以及所要解决的问题来设定尺度因子。比如对于A地区所要解决的问题，在分辨率为10ms即可，那么在反演时就可以设定这个尺度因子，让算法就在这个尺度上进行反演，这样在一定程度上减小了反演的多解性，加速了算法的收敛速度。

[0052] 本发明采用的技术原理：

[0053] 地球物理参数反演问题的求解通常是通过优化目标函数来进行的，这个过程就是寻找使得目标函数值取得最小时所对应的参数的过程，数学上将这个过程称之为最优化。

[0054] 设向量 若函数 在欧式空间中有定义，且至少存在一阶导数，则在 邻域内可将其表示成：

[0055]

[0056] 式中，

[0057] 为 在 处的灵敏度向量，是在 的邻域内使得

[0058] 目标函数值下降最快的方向。在确定了灵敏度向量后，再求取步长因子 来构建反演算法的迭代格式

[0059]

[0060] 算法的迭代步骤可以描述为：

[0061] (1)给定目标函数的初始寻优点 允许误差为ε，k起始值为0，则可将初始寻优点 表示为

[0062] (2)计算目标函数在 处的灵敏度向量

[0063] (3)确定步长因子 使得

[0064] (4)令 若 则停止迭代，得最优解

[0065] 否则令k＝k+1，返回步骤(2)。

[0066] 本发明具体的步骤如下：

[0067] (1))输入叠后地震数据，从所述叠后地震数据得到地震记录的采样率为Δt，地震记录的时间长度为T，则样点个数M＝T/Δt；采样率Δt和时间长度T为已知的，这里提到这两个参数的目的是为了说明样点个数M；

[0068] (2)由测井数据及解释层位数据建立反演的初始模型L[i]，i＝1，M+1，所述测井数据包括声波时差测井数据和密度测井数据Den，具体如下：首先由所述声波时差测井数据计算得到地层速度Vp，然后根据地层速度Vp和密度测井数据Den利用波阻抗公式Z＝Vp*Den计算得到井中波阻抗曲线，再对所述井中波阻抗曲线进行带通滤波(即在解释层位数据的控制下将所述井中阻抗曲线进行内插和外推)就得到了反演的初始模型，也就是前面所说的低频模型；

[0069] (3)设反演参数个数为P，尺度因子为n×Δt，n为正整数，计算能够整除尺度因子且大于地震记录的时间长度T的最小整数T′，则反演参数个数P＝T′/Δt。设经尺度扩展后的波阻抗L用向量Lp[i]，i＝1，T′/Δt表示，Lp中数值由以下公式确定：

[0070]

[0071] 该公式的含义是以尺度因子为标准，在初始模型上均匀的采样，式中L[n×i]是初始模型，即步骤(2)中的初始模型L[i]，i＝1，M+1(注：Ｌ[i]，i＝1，M+1的下标取值范围是从1到M+1，L[n×i]的下标i的范围是从1到P，见Lp[n×i]＝L[n×i]，i＝1，2，…，P，L[n×i]表示在L[i]上以间隔n取值)。本发明的多分辨率也体现在该步骤中；

[0072] (4)根据地震记录的质量确定控制约束范围的约束系数alphadi和alphaui(约束系数alphadi和alphaui取值范围是0和1之间，0＜＝alphadi，alphaui＜＝1)，计算反演算法的约束范围[Ldi，Lui ]，其中Ldi＝(1-alphadi)×Li，Lui＝(1+alphaui)×Li；

[0073] (5)计算反射系数

[0074] (6)从地震记录中提取子波Wj，j＝1，N(可以采用统计法提取子波，请参考文献：刘明洋，郭韬，沈铭成等的“多道统计法提取子波及其应用”，《石油天然气学报》，2010年10月第32卷第5期)，将提取的子波进行归一化处理(即首先统计出子波数据中的最大值，然后将子波数据除以该最大值，经过这个处理的结果还是子波)，继而由摺积公式(褶积也叫卷积。设f1(t)和f2(t)是以t为自变量的函数，二者的卷积积分被定义为：其中“*”表示卷积运算符号。地震信号是离散信号，
设离散信号x(l)和h(m)，lmin≤l≤lmax，mmin≤m≤mmax，nmin≤n≤nmax，nmin＝lmin+mmin，nmax＝lmax+mmax，x(l)和h(m)的卷积公式为： )计算合成地震记录
Syn＝R*W，具体计算公式为：

[0075]

[0076] 式中，R为反射系数，W为子波数据。

[0077] (7)利用合成地震记录与经过归一化处理后的叠后地震数据残差(残差就是指公式中的Synk-Datak)的二范数建立用于反演的目标函数：

[0078]

[0079] Data是经过归一化处理后的叠后地震数据，反演的目的是从这个叠后地震数据中得到地下介质的波阻抗。

[0080] (8)计算灵敏度向量

[0081]

[0082] 其中，

[0083] (9)在第(2)步的初始模型及第(8)步的灵敏度向量的基础上利用具有全局优化能力的优化算法对第(7)步的目标函数进行优化，即在给定了初始模型以及灵敏度向量后，优化算法就在初始模型处开始，并以灵敏度向量为优化方向，对目标函数进行优化，当目标函数值达到足够小或者达到满意的值时停止算法迭代(详见前述的“算法的迭代步骤”中的4个步骤。)，输出反演结果。从目标函数的公式可以看出其最小值是0，这是目标函数的最理想情况。但是由于实际地震数据含有各种噪音以及采集仪器的误差限制，就造成了不可能反演到最理想的情况，那么越接近于0越好，考虑到计算效率问题，通常设定目标函数值达到一个认为设定的值即可，比如0.1或者0.00001。

[0084] 所述具有全局优化能力的优化算法是现有技术，比如可以采用共轭梯度法或阻尼最小二乘法等。

[0085] 下面以理论模型数据和某地区的实际地震数据的反演应用实例来说明本发明的效果。图1(a)至图1(o)为理论模型数据的反演实例，图2(a)至图2(c)为实际地震数据的反演实例。

[0086] 实施例1

[0087] 本实施例是将本发明方法用于理论模型数据反演的实施例，具体步骤如下：

[0088] 1)输入叠后地震数据，得到采样率和时间长度，并对地震数据进行预处理，提高地震数据的信噪比；

[0089] 2)根据地震数据频谱特征确定尺度因子，由尺度因子确定反演的参数个数；

[0090] 3)输入解释层位数据和测井数据，并计算井中波阻抗曲线：由声波时差测井数据计算得到地层速度Vp，另外一条测井数据是密度测井数据Den，利用波阻抗公式Z＝Vp*Den计算井中波阻抗曲线；

[0091] 4)在解释层位的控制下将井中阻抗曲线进行内插、外推，建立反演的初始模型；

[0092] 5)根据地震资料质量确定控制约束范围的约束系数alphadi和alphaui，在初始模型的基础上计算反演算法的约束范围；

[0093] 6)提取子波，并对子波数据进行归一化处理；

[0094] 7)由地震数据和摺积合成地震记录的二范数建立目标函数，计算该目标函数在初始模型处的灵敏度向量；

[0095] 8)在目标函数值和灵敏度向量的指引下进行模型约束的波阻抗反演。

[0096] 图1(a)是实施例1的理论合成地震记录剖面。地震记录从右侧向左侧分辨率依次降低，用常规波阻抗反演方法及现有商业软件难以有效的反演出该地震剖面左侧的薄层，本发明的目的是从地震数据中反演出地下介质的绝对波阻抗，得到从地震数据难以识别的薄层。

[0097] 图1(b)是实施例1的初始模型。该模型是由井中波阻抗曲线在地震解释层位的指引下内插外推计算得到。可见未经反演的波阻抗模型无法适用于薄层、岩性油气藏及隐蔽油气藏的解释，无法用于岩性参数的反演。本发明在初始模型的基础上引入多尺度因子及控制搜索范围的约束系数，结合地震数据进行模型约束的波阻抗反演，得到地层的绝对波阻抗，并反演出具有地质意义的薄层。

[0098] 图1(c)至图1(e)和图1(f)至图1(h)为该实施例剖面第50、30及第1道对应的地震记录和初始模型，图1(i)至图1(k)为该实施例剖面第50、30及第1道的反演结果(细线)和初始模型(粗线)。图1(l)至图1(n)为该实施例剖面第50、30及第1道的反演结果。图1(b)是初始模型，图1(o)为利用本发明方法以图1(b)为初始模型得到的反演结果，图1(b)未包含高频信息，而图1(o)包含了高频信息，分辨率达到采样率大小，从图中可以看出本发明将薄层反演出来了，所得的反演结果为绝对波阻抗，可用于地震资料解释工作中的岩性参数反演。

[0099] 实施例2：

[0100] 本实施例选用实际地震数据。利用本发明基于模型约束的波阻抗反演方法进行反演，步骤同实施例1。

[0101] 图2(a)是该实施例叠后地震数据，图2(b)是该实施例初始模型，图2(c)是该实施例反演结果。从图2(c)的反演结果可以显示高频信息，反演结果显示了薄层信息，反映了地层的接触关系，对精细储层描述，岩性(和地层)圈闭油气藏和非常规油气藏(煤层气、页岩气等)提供了重要的参考资料。

[0102] 上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。