CN104062680B - 一种计算波阻抗反演目标函数梯度的方法 - Google Patents

Abstract

本发明提供了一种计算波阻抗反演目标函数梯度的方法，属于油气及煤层气地震勘探与开发领域。所述方法根据线性褶积模型和叠后地震反射系数定义，基于叠后地震振幅数据、波阻抗初始模型和地震子波数据，构建波阻抗反演的目标函数，并推导出波阻抗反演的目标函数梯度的解析解。利用本发明能够正确地计算出目标函数的梯度，在一定程度上提高了地震波阻抗反演的精度，降低了波阻抗反演的多解性，可直接用于实际地震数据的波阻抗反演。同时，本发明的计算精度和实用性高，并且具有计算速度快、稳定性好的优点。

Description

一种计算波阻抗反演目标函数梯度的方法

技术领域

本发明属于油气及煤层气地震勘探与开发领域，具体涉及一种计算波阻抗反演目标函数梯度的方法。

背景技术

地震勘探技术是油气勘探中应用最为广泛的一种地球物理学方法，利用地震波在不同介质中传播的速度、振幅、频率、相位、波形等参数的变化来分析、预测油气储层分布范围及储层物性特征。随勘探开发程度逐渐深入，地震勘探已经从认识地下构造形态的构造勘探，逐渐发展成直接应用地震信息判断岩性、分析岩相、定量计算岩层物性参数的岩性勘探。这些复杂沉积环境下的储层通常表现为具有厚度薄、物性高度非均质、有效储层规模小、分布分散、岩石物理关系复杂、储层岩性差异小等特征。波阻抗是反映储集岩储层特征的一种重要物性参数，与振幅等反映界面性质的参数相比，波阻抗具有更高的储层横向预测能力，因而地震波阻抗反演技术是岩性油气勘探的重要技术手段。近几十年来地震波阻抗反演技术获得了长足的进展，如道积分法、模型法、随机反演、谱反演等。其中模型法和谱反演均采用经典的最优化算法进行反演，这一类算法均需要计算目标函数的梯度，目前大都采用中心差分法计算梯度。采用中心差分法计算目标函数梯度时需要计算2N次目标函数(N为参数个数)，计算量较大。并且由于差分法在计算每一个参数的偏导数时需要除以一个较小的数值，该较小值的取法没有一定的标准，使得该参数的偏导数误差较大，从而导致由各个参数的偏导数构成的梯度的误差也较大。正是由于差分法计算梯度存在误差因素，使得基于经典最优化算法的波阻抗反演迭代次数增加、计算量大，且难以得到目标函数意义下的全局最优解，不适于储层参数定量计算及薄储层预测，计算速度慢不适于实时解释，应用受到限制。

发明内容

针对波阻抗反演中常规的梯度计算方法计算速度慢，精度低，不适于储层参数定量计算及薄储层预测，应用受到限制等缺点，本发明的目的在于解决上述现有技术中存在的难题，提供一种计算波阻抗反演目标函数梯度的方法，提高地震波阻抗反演的精度，加速算法收敛速度，降低波阻抗反演的多解性。用于地震资料的岩性参数解释，为寻找岩性(和地层)圈闭油气藏和非常规油气藏(煤层气、页岩气等)提供重要的资料。

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

一种计算波阻抗反演目标函数梯度的方法，所述方法根据线性褶积模型和叠后地震反射系数定义，基于叠后地震振幅数据、波阻抗初始模型和地震子波数据，构建波阻抗反演的目标函数，并推导出波阻抗反演的目标函数梯度的解析解。

所述方法包括以下步骤：

步骤1，从测井数据中获取波阻抗初始模型，即得到波阻抗Z_j；

步骤2，从井旁地震道获取地震子波数据w和叠后地震振幅数据d；

步骤3，根据地震勘探中叠后地震反射系数定义，利用下面的反射系数公式计算反射系数r_j：

$r_j=\frac{\ln Z_{j+1}-\ln Z_j}{2};---\left(9\right)$

步骤4，采用线性褶积模型来描述地震响应，其数学表达式为：

s＝w*r (10)

(10)式中，s为摺积合成地震记录；根据反射系数公式(9)及褶积公式将所述数学表达式(10)展开得到：

$s_i=\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{r}_j\·w_{i-j}=\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\·w_{i-j},i=1,\·\·\·,n---\left(11\right)$

步骤5，求取合成地震记录s中的最大值，记录其所在序列中的序号为ix，由叠后地震振幅数据d和摺积合成地震记录s构建波阻抗反演的目标函数，如公式(12)所示：

$f={\left|\right|s-d\left|\right|}_2=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{[s_i-d_i]}^2=\underset{i=0}{\overset{n}{\mathrm{\Σ}}}[{\frac{\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\·w_{i-j}}{\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\·w_{\mathrm{ix}-j}}-d_i]}^2---\left(12\right)$

步骤6，采用公式(13)计算波阻抗反演目标函数的偏导数：

$\frac{\∂f}{{\∂Z}_j}=2\underset{i=0}{\overset{n}{\mathrm{\Σ}}}[\frac{\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\·w_{i-j}}{\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\·w_{\mathrm{ix}-j}}-d_i]\·$ (13)

$\frac{\frac{1}{{2Z}_j}\left\{\right[w_{i-j+1}-w_{i-j}\left]\right[\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{2}(\ln Z_{j+1}-\ln Z_j)\·w_{\mathrm{ix}-j}]-[w_{\mathrm{ix}-j+1}-w_{\mathrm{ix}-j}\left]\right[\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{1}{2}(\ln Z_{j+1}-\ln Z_j)\·w_{i-j}\left]\right\}}{{[\underset{j=0}{\overset{m}{\mathrm{\Σ}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\·w_{\mathrm{ix}-j}]}^2}$

步骤7，采用公式(14)由步骤6计算得到的偏导数构建波阻抗反演目标函数的梯度：

$\▿f=(\frac{\∂f}{{\∂Z}_1},\frac{\∂f}{{\∂Z}_2},\·\·\·\frac{\∂f}{{\∂Z}_j},\·\·\·,\frac{\∂f}{{\∂Z}_m})---\left(14\right).$

与现有技术相比，本发明的有益效果是：

(1)利用本发明能够正确地计算出目标函数的梯度，在一定程度上提高了地震波阻抗反演的精度，降低了波阻抗反演的多解性，可直接用于实际地震数据的波阻抗反演。

(2)本发明的计算精度和实用性高，并且具有计算速度快、稳定性好的优点。

附图说明

图1(a)是实施例1中地震子波数据。

图1(b)是实施例1中波阻抗(Imp)，反射系数(r)、地震振幅数据(Amp)、初始波阻抗模型(Model)及本发明计算出的梯度(gradient)；

图2(a)是本发明实施例2中地震子波数据。

图2(b)是实施例2中楔状模型及地震合成记录。

图2(c)是实施例2中初始模型及本发明计算的梯度。

图3是本发明方法的步骤框图。

具体实施方式

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

针对现有梯度计算方法在地震波阻抗反演中存在的不足，围绕差分法求取梯度时存在的误差大、计算速度慢、无法处理大数据体的实时解释等问题，本发明提供一种计算波阻抗反演目标函数梯度的方法，

为了从叠后地震数据反演得到能够用于计算储层岩性参数的绝对波阻抗，本发明根据线性褶积模型和叠后地震反射系数定义，基于地震数据、波阻抗初始模型和地震子波数据，由实际地震数据和摺积合成地震记录的二范数构建波阻抗反演的目标函数，根据线性褶积模型和叠后地震反射系数定义推导出波阻抗反演的目标函数梯度的解析解，即精确解，将得到的梯度信息用于波阻抗反演，提高了算法收敛速度和反演结果的精度，便于地震资料的岩性参数解释和精细储层描述，为寻找岩性(和地层)圈闭油气藏和非常规油气藏(煤层气、页岩气等)提供可靠的资料。

如图3所示，本发明方法包括：

(1)根据地震勘探中叠后地震反射系数定义，由地层的波阻抗计算反射系数，其数学表达式为：

$r_j=\frac{Z_{j+1}-Z_j}{Z_{j+1}+Z_j},j=1,\·\·\·,m---\left(1\right)$

上式中，r_j为反射系数，Z_j为波阻抗，m为样点数；对公式(1)进行变换，得到公式(2)：

$Z_{j+1}=Z_j\frac{1+r_j}{1-r_j}---\left(2\right)$

(2)在公式(2)两端取自然对数，得到公式(3)：

$\ln Z_{j+1}=\ln Z_j+\ln \frac{1+r_j}{1-r_j}=\ln Z_j+\ln (1+r_j)-\ln (1-r_j)---\left(3\right)$

由自然对数展开公式(4)，对公式(3)右端项的后两项进行展开，得到公式(5)：

$\ln (1+x)=\underset{k=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-1)}^{k-1}\frac{x^k}{k},\left|x\right|<1---\left(4\right)$

$\ln (1+r_j)-\ln (1-r_j)=\underset{k=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-1)}^{k-1}\frac{r_j^k}{k}-\underset{k=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}{(-1)}^{k-1}\frac{{\left({-r}_j\right)}^k}{k}=2\underset{k=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{r_j^{2k-1}}{2k-1}---\left(5\right)$

(3)将公式(5)右端的通项公式进行展开，得到(6)：

$\ln (1+r_j)-\ln (1-r_j)=2\underset{k=1}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{r_j^{2k-1}}{2k-1}=2[r_j+\frac{r_j^3}{3}+\frac{r_j^5}{5}+\&CenterDot;\&CenterDot;\&CenterDot;]---\left(6\right)$

(4)将公式(6)代入公式(3)，得到公式(7)：

$\ln Z_{j+1}={\ln Z}_j+2[r_j+\frac{r_j^3}{3}+\frac{r_j^5}{5}+\&CenterDot;\&CenterDot;\&CenterDot;]---\left(7\right)$

对于实际地震资料，反射系数较小，一般满足|r|≤0.2，因此可以省略三次幂及三次幂以上的项，对所述的公式(7)进行简化得到公式(8)：

1nZ_j+1≈1nZ_j+2r_j (8)

(5)由公式(8)得到反射系数公式(9)：

$r_j=\frac{\ln Z_{j+1}-\ln Z_j}{2}---\left(9\right)$

(6)根据经典地震勘探理论，采用线性褶积模型来描述地震响应，其数学表达式为：

s＝w*r (10)

(10)式中，w为地震子波，r为反射系数，s为合成地震记录。w为给定的地震子波；所述的数学表达式(10)由反射系数公式(9)及褶积公式得到如下展开形式：

$s_i=\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}{r}_j\&CenterDot;w_{i-j}=\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{i-j},i=1,\&CenterDot;\&CenterDot;\&CenterDot;,n---\left(11\right)$

(7)求取合成地震记录s中的最大值，记录其所在序列中的序号为ix，由实际地震振幅数据(即d_i，i＝1，…，n)和摺积合成地震记录(即公式(11)得到的s)构建波阻抗反演的目标函数，如公式(12)所示：

$f={\left|\right|s-d\left|\right|}_2=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{[s_i-d_i]}^2=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{[\frac{\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{i-j}}{\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{\mathrm{ix}-j}}{-d}_i]}^2---\left(12\right)$

所述公式(12)中的d_i，i＝1，…，n为地震振幅数据；

(8)计算波阻抗反演目标函数的偏导数，采用如下的公式(13)：

$\frac{\&PartialD;f}{{\&PartialD;Z}_j}=2\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}[\frac{\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{i-j}}{\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{\mathrm{ix}-j}}-d_i]\&CenterDot;$ (13)

$\frac{\frac{1}{{2Z}_j}\left\{\right[w_{i-j+1}-w_{i-j}\left]\right[\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{2}(\ln Z_{j+1}-\ln Z_j)\&CenterDot;w_{\mathrm{ix}-j}]-[w_{\mathrm{ix}-j+1}-w_{\mathrm{ix}-j}\left]\right[\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{1}{2}(\ln Z_{j+1}-\ln Z_j)\&CenterDot;w_{i-j}\left]\right\}}{{[\underset{j=0}{\overset{m}{\mathrm{\&Sigma;}}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{\mathrm{ix}-j}]}^2}$

(13)式中，由于公式太长分成了两行，第一行最后有个圆点，表示乘积。

(9)由步骤(8)计算的偏导数构建波阻抗反演目标函数的梯度，采用如下的公式(14)：

$\&dtri;f=(\frac{\&PartialD;f}{{\&PartialD;Z}_1},\frac{\&PartialD;f}{{\&PartialD;Z}_2},\&CenterDot;\&CenterDot;\&CenterDot;\frac{\&PartialD;f}{{\&PartialD;Z}_j},\&CenterDot;\&CenterDot;\&CenterDot;,\frac{\&PartialD;f}{{\&PartialD;Z}_m})---\left(14\right)$

公式(14)是目标函数的梯度表达式，它是由目标函数对各个参数的偏导数组成的，关于偏导数的解析解就是公式(13)所描述的。

下面以某地区的实际井旁地震数据和理论楔状模型的应用实例来说明本发明的效果。图1(a)和图1(b)为实际井旁地震数据的应用实例，图2(a)至图2(c)为楔状模型数据的应用实例。

实施例1

本实施例是本发明用于井旁地震数据的实施实例，一种计算波阻抗反演目标函数梯度的方法的步骤如下：

第一步，从测井数据中获取波阻抗初始模型(Model)，即步骤(1)中的波阻抗(Z)；

第二步，从井旁地震道获取地震子波(w)和地震振幅数据(Amp)，即步骤(6)中的地震子波(w)和步骤(7)中的地震数据(d)；

第三步，按照步骤(5)中的公式(9)由波阻抗数据(Z)计算反射系数(r)；

第四步，由地震子波数据w、地震振幅数据d和波阻抗初始模型数据Z由步骤(8)的公式(13)计算波阻抗反演目标函数的偏导数；

第五步，按照步骤(9)中的公式(14)构建波阻抗反演目标函数的梯度；

图1(a)为实施例1中的地震子波数据，该地震子波数据通常由统计法子波提取或井旁道子波提取，本发明中令该地震子波数据为已知。

图1(b)为实施例1中的波阻抗、反射系数、振幅数据、波阻抗初始模型及本发明计算的梯度，图1(b)中的波阻抗数据(Imp)从测井数据中得到，反射系数由步骤(5)中的公式(9)计算得到，图1(b)中的振幅数据由井旁地震道得到，图1(b)中的波阻抗初始模型(Model)是对波阻抗数据(Imp)低通滤波得到，图1(b)中的梯度(gradient)按照步骤(9)中的公式(14)计算得到。

实施例2：

本实施例选用楔状模型数据。利用本发明计算波阻抗反演目标函数的梯度，步骤同实施例1。

图2(a)是给定的理论子波，图2(b)的上部为楔状模型，图2(b)的下部为按照步骤(6)中公式(10)计算得到的地震合成记录，图2(c)的上部为波阻抗初始模型，图2(c)的下部是由图2(a)的地震子波和图2(b)中地震合成记录及图2(c)中的波阻抗初始模型按照步骤(9)中的公式(14)计算得到的梯度。

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

Claims (1)

1.一种计算波阻抗反演目标函数梯度的方法，其特征在于：所述方法根据线性褶积模型和叠后地震反射系数定义，基于叠后地震振幅数据、波阻抗初始模型和地震子波数据，构建波阻抗反演的目标函数，并推导出波阻抗反演的目标函数梯度的解析解，

所述方法包括以下步骤：

步骤1，从测井数据中获取波阻抗初始模型，即得到波阻抗Z_j；

步骤2，从井旁地震道获取地震子波数据w和叠后地震振幅数据d；

步骤3，根据地震勘探中叠后地震反射系数定义，利用下面的反射系数公式计算反射系数r_j：

$r_j=\frac{\ln Z_{j+1}-\ln Z_j}{2};---\left(9\right)$

步骤4，采用线性褶积模型来描述地震响应，其数学表达式为：

S＝w*r_j (10)

(10)式中，s为褶积合成地震记录；根据反射系数公式(9)及褶积公式将所述数学表达式(10)展开得到：

$s_i=\underset{j=0}{\overset{m}{\&Sigma;}}{r}_j\&CenterDot;w_{i-j}=\underset{j=0}{\overset{m}{\&Sigma;}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{i-j},i=1,...,n---\left(11\right)$

步骤5，求取合成地震记录s中的最大值，记录其所在序列中的序号为ix，由叠后地震振幅数据d和褶积合成地震记录s构建波阻抗反演的目标函数，如公式(12)所示：

$f=\left|\right|s-d|{|}_2=\underset{i=0}{\overset{n}{\&Sigma;}}{\&lsqb;s_i-d_i\&rsqb;}^2=\underset{i=0}{\overset{n}{\&Sigma;}}{\&lsqb;\frac{\underset{j=0}{\overset{m}{\&Sigma;}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{i-j}}{\underset{j=0}{\overset{m}{\&Sigma;}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{ix-j}}-d_i\&rsqb;}^2---\left(12\right)$

步骤6，采用公式(13)计算波阻抗反演目标函数的偏导数：

$\begin{array}{c}\frac{\&part;f}{\&part;Z_j}=2\underset{i=0}{\overset{n}{\&Sigma;}}\&lsqb;\frac{\underset{j=0}{\overset{m}{\&Sigma;}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{i-j}}{\underset{j=0}{\overset{m}{\&Sigma;}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{ix-j}}-d_i\&rsqb;\&CenterDot;\\ \frac{\frac{1}{2Z_j}\{\&lsqb;w_{i-j+1}-w_{i-j}\&rsqb;\&lsqb;\underset{j=0}{\overset{m}{\&Sigma;}}\frac{1}{2}(\ln Z_{j+1}-\ln Z_j)\&CenterDot;w_{ix-j}\&rsqb;-\&lsqb;w_{ix-j+1}-w_{ix-j}\&rsqb;\&lsqb;\underset{j=0}{\overset{m}{\&Sigma;}}\frac{1}{2}(\ln Z_{j+1}-\ln Z_j)\&CenterDot;w_{i-j}\&rsqb;\}}{{\&lsqb;\underset{j=0}{\overset{m}{\&Sigma;}}\frac{\ln Z_{j+1}-\ln Z_j}{2}\&CenterDot;w_{ix-j}\&rsqb;}^2}\end{array}---\left(13\right)$

步骤7，采用公式(14)由步骤6计算得到的偏导数构建波阻抗反演目标函数的梯度：

$\&dtri;f=(\frac{\&part;f}{\&part;Z_1},\frac{\&part;f}{\&part;Z_2},...\frac{\&part;f}{\&part;Z_j},...,\frac{\&part;f}{\&part;Z_m})---\left(14\right).$