CN102393532A - 地震信号反演方法 - Google Patents

Abstract

本发明公开了一种地震信号反演方法，通过改进原共轭梯度算法中的梯度更新方式，在梯度更新中引入阻尼系数，通过调节阻尼系数来调节梯度更新，从而提高求解的稳定性。本发明的积极效果是：1)速度快，适合大规模地震数据反演；2)储层预测准确，与实际测井资料吻合度高；3)稳定性高，对含噪声的地震信号抗噪能力强，无异常解现象。

Description

地震信号反演方法

技术领域

本发明涉及一种地震信号反演方法，具体涉及一种地震信号反演方法。

背景技术

地震反演是利用地表观测地震资料，以已知地质规律和钻井、测井资料为约束，对地下岩层空间结构和物理性质进行成像(求解)的过程。反演与正演相对存在，正演是利用已知地质模型求解地震响应(即地震波形)的过程。反演目的是利用地震波在地下介质中的传播规律，通过数据采集、处理与解释等流程，来推测地下岩层结构和物件参数的空间分布。这些物件参数包括：速度，密度，泊松比等，其中泊松比能够为探明油气储量提供重要依据。目前，是地震解释领域最有实用价值也是研究最热门的一门技术。

地震反演根据地震资料形式分为叠前反演和叠后反演。叠后反演经过几十年发展，其理论和应用相对成熟，然而，随着全球油气储量的减少，人们对储层预测的要求逐渐提高，叠后反演存在的问题越来越明显。

1)所要预测的储层性质越来越复杂，有些储层没有明显的声学特征，不同的岩性之间波阻抗没有显著的差别，所以无法利用传统的叠后地震反演区别不同的的岩性并划分储集单元；

2)叠后地震数据在叠加过程中损失了很多重要的地震原始信息，因而降低了其解决地质问题的精度和能力；

3)叠后反演在假设条件、子波提取及振幅、频率、保真等方面均存在着诸多问题，如垂直入射、零偏移距的假设使地震数据失掉了可贵的AVO信息，以及预测地层弹性参数和含油气性的能力较差等等。

而叠前反演能够解决上述叠后反演存在的问题，具有更高的实用价值。但是叠前反演也具有很高的困难性，其困难性主要有：

a)叠前反演的数据量巨大，耗时费力，导致成本高(包括设备成本，人力成本，时间成本等)。

b)反演问题的不适定性(即多解性)，造成多解性的原因不是反演方法或技巧上的缺陷造成的，而是地球物理数据的有限观测、带限性和存在噪声所致。解决这个问题现在已有在利用测井信息进行约束以减少多解性，但正确添加测井信息同时就增加了反演的复杂度。

这种情况下，设计一种可靠性高，稳定性强，速度快，预测准确的可以用于叠前和叠后反演的方法势在必行。

现有技术中，存在如下两大类反演方法：

1)基于非线性假设的反演

当前这类反演方法又有模拟退火法，遗传算法，人工神经网络算法，禁忌搜索算法，蚁群算法，粒子群算法等。总体来说，这类方法假设的地球模型是非线性的，虽然比较符合实际，而且部分方法已取得良好效果，然而他们都有以下致命缺陷：①运算效率很低；②同时也存在着解的稳定性问题，因为遗传算法和模拟退火等算法属于随机搜索算法，对这样病态性问题求解的时候，每次解都可能差异很大，因此无法应用到大规模地震数据反演中。

2)基于线性假设的反演

当前这类方法又有梯度法，牛顿法，共轭梯度法，最小二乘法，广义逆法等。总体来说，这类方法将非线性问题线性化，通过迭代，逐次逼近，求得反演近似解。这类方法相对于非线性方法来说更简便，易行，多数情况下效果可以接受，然而在迭代中容易陷入局部最优，出现多解或无解现象，而且反演结果很大程度上依赖于初始模型。

发明内容

为了克服现有技术的上述缺点，本发明提供了一种地震信号反演方法。可靠性高，稳定性强，速度快，预测准确，既可用于叠前反演，也可用于叠后反演。

本发明解决其技术问题所采用的技术方案是：一种地震信号反演方法，包括如下步骤：

第一步、数据预处理；

第二步、建立正演模型；

第三步、以正演模型作为反演基础，对反问题进行求解：

1)建立目标函数；

2)在N维空间内，采用迭代构造法构造n个相互共轭的向量，做一维线性搜索，通过n次线性搜索得到目标函数的极小值点；

3)取目标函数的梯度构造出解空间的基；

4)对解空间的基进行迭代构造，引入阻尼系数调整解的稳定性，控制梯度向量的方向；

5)为每次迭代结果添加硬门限，在硬门限的限制范围内求解得到唯一的稳定解；

第四步、反演后处理：通过第三步求出的唯一的稳定解来求解纵波速度、横波速度和密度，从而求得目标速度和密度结果，最后求得泊松比。

与现有技术相比，本发明的积极效果是：1)速度快，适合大规模地震数据反演；2)储层预测准确，与实际测井资料吻合度高；3)稳定性高，对含噪声的地震信号抗燥能力强，无异常解现象。

具体实施方式

首先对本发明中的相关术语做一下定义：

反演：利用地震资料通过反演方法同时获得纵波阻抗(叠后)、纵波、横波、密度和泊松比(叠前)等参数。

褶积模型：褶积模型是一种制作合成(理论)地震记录的模型，它假设每道地震记录是由地震子波与地下模型各层的反射系数之褶积所构成，必要时还可以加上随机噪声。

病态反问题：反问题是相对于正问题而言的，这里主要指从观测得到的角道集中求解得到纵波、横波和密度。所谓病态性不满足存在性条件，唯一性问题和稳定性条件三个条件中任何一个条件的反问题：

共轭梯度：是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hessen矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。

一种地震信号反演方法，包括如下步骤：第一步、数据预处理；第二步、建立正演模型；第三步、以正演模型作为反演基础，对反问题进行求解；第四步、反演后处理。其中：

第一步、数据预处理：

1)原始地震数据首先经过去噪，校正等常规处理；

2)对道集进行部分叠加，以减少道集复杂度和数据量；

3)对道集全叠加来生成叠后数据；

4)标定进行反演的目标层位；

5)从测井资料中获取速度，密度等所需信息。

第二步、建立正演模型：

采用如下单层Gidlow近似方程进行正演：

$R\left(\stackrel{\‾}{\mathrm{\θ}}\right)={\sec }^2\stackrel{\‾}{\mathrm{\θ}}\×\frac{1}{2}(\frac{\mathrm{\Δ\α}}{\stackrel{\‾}{\mathrm{\α}}}+\frac{\mathrm{\Δ\ρ}}{\stackrel{\‾}{\mathrm{\ρ}}})-8{\stackrel{\‾}{\mathrm{\γ}}}^2{\sin }^2\stackrel{\‾}{\mathrm{\θ}}\×\frac{1}{2}(\frac{\mathrm{\Δ\β}}{\stackrel{\‾}{\mathrm{\β}}}+\frac{\mathrm{\Δ\ρ}}{\stackrel{\‾}{\mathrm{\ρ}}})+({4\stackrel{\‾}{\mathrm{\γ}}}^2{\sin }^2\stackrel{\‾}{\mathrm{\θ}}-{\tan }^2\stackrel{\‾}{\mathrm{\θ}})\×\frac{1}{2}\frac{\mathrm{\Δ\ρ}}{\stackrel{\‾}{\mathrm{\ρ}}}$

$={\sec }^2\stackrel{\‾}{\mathrm{\θ}}\×R_p-8{\stackrel{\‾}{\mathrm{\γ}}}^2{\sin }^2\stackrel{\‾}{\mathrm{\θ}}\×R_s+(4{\stackrel{\‾}{\mathrm{\γ}}}^2{\sin }^2\stackrel{\‾}{\mathrm{\θ}}-{\tan }^2\stackrel{\‾}{\mathrm{\θ}})\×R_d---\left(1\right)$

其中，

表示随角度变化的PP波反射系数，θ表示PP波的入射角，α表示纵波速度，β表示横波速度，γ表示横波纵波比，ρ表示PP波的密度，R_p表示纵波波阻抗反射系数，R_s表示横波波阻抗反射系数，R_d表示密度反射系数。

为了避免动校正引起的拉伸畸变及调谐效应对反演参数的影响，反演在动校正前的地震数据上进行，同时采用子波W描述褶积模型，得到基于Gidlow方程的褶积模型：

$\left[\begin{array}{c}{d}_1\\ M\\ d_k\end{array}\right]=\left[\begin{array}{ccc}{\mathrm{WE}}_{1}D& {\mathrm{WF}}_{1}D& {\mathrm{WH}}_{1}D\\ M& M& M\\ {\mathrm{WE}}_{k}D& {\mathrm{WF}}_{k}D& {\mathrm{WHD}}_k\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_d\end{array}\right]---\left(2\right)$

其中，D为差分矩阵；W为子波序列；E、F和H均为对角矩阵，各自的对角元分别为 $e_j={\sec }^2{\stackrel{\‾}{\mathrm{\θ}}}_j,$ $f_j=-8\stackrel{\‾}{\mathrm{\γ}}\mathrm{si}{n}^2{\stackrel{\‾}{\mathrm{\θ}}}_j,$ $h_j-1-4\stackrel{\‾}{\mathrm{\γ}}{\sin }^2{\stackrel{\‾}{\mathrm{\θ}}}_j,$ i表示中心道集的第i道；k为中心道集数量；L_p，L_s，L_d分别代表纵波阻抗，横波阻抗，密度阻抗。令L_pi＝ln(Z_pi)，Z_pi＝α*ρ，(i＝1，2，……，k)，则

式(2)可写为：

d_i＝c₁L_p+c₂L_s+c₃L_d         (3)

其中，c₁＝WE_iD、c₂＝WF_iD、c₃＝WH_iD。

由于L_P与L_d，L_P与L_S满足近似的线性关系：

L_S＝kL_P+k_C+ΔL_S

                       (4)

L_d＝mL_P+m_C+ΔL_d

其中，k，m，k_c，m_c均为常数，

因此，将式(4)代入式(3)，对褶积模型进行约束，得到带岩石物理约束关系的AVO褶积模型如下式：

$\left[\begin{array}{c}{d}_1\left({\mathrm{\θ}}_1\right)\\ d_2\left({\mathrm{\θ}}_2\right)\\ d_N\left({\mathrm{\θ}}_N\right)\end{array}\right]=\left[\begin{array}{ccc}{\tilde{c}}_1\left({\mathrm{\θ}}_1\right)W\left({\mathrm{\θ}}_1\right)D& {\tilde{c}}_2\left({\mathrm{\θ}}_1\right)W\left({\mathrm{\θ}}_1\right)D& {\tilde{c}}_3\left({\mathrm{\θ}}_1\right)W\left({\mathrm{\θ}}_1\right)D\\ {\tilde{c}}_1\left({\mathrm{\θ}}_2\right)W\left({\mathrm{\θ}}_2\right)D& {\tilde{c}}_2\left({\mathrm{\θ}}_2\right)W\left({\mathrm{\θ}}_2\right)D& {\tilde{c}}_3\left({\mathrm{\θ}}_2\right)W\left({\mathrm{\θ}}_2\right)D\\ M& M& M\\ {\tilde{c}}_1\left({\mathrm{\θ}}_N\right)W\left({\mathrm{\θ}}_N\right)D& {\tilde{c}}_2\left({\mathrm{\θ}}_N\right)W\left({\mathrm{\θ}}_N\right)D& {\tilde{c}}_3\left({\mathrm{\θ}}_N\right)W\left({\mathrm{\θ}}_N\right)D\end{array}\right]\left[\begin{array}{c}{L}_p\\ {\mathrm{\ΔL}}_s\\ {\mathrm{\ΔL}}_d\end{array}\right]---\left(5\right)$

其中 ${\tilde{c}}_1=(1/2)c_1+(1/2){\mathrm{kc}}_2+{\mathrm{mc}}_3,$ ${\tilde{c}}_2=(1/2)c_2.$

式(5)可简写为：

d＝Gx                   (6)

d＝Gx即为正演模型。

第三步、以正演模型作为反演基础，对反问题进行求解：

对反问题进行求解，就是对式(6)中的x进行求解。

式(6)为良态方程的条件有以下三个：

①、存在性条件：在定义域内，对应每一个y，都有解x存在；

②、唯一性问题：解x是唯一的；

③、稳定性条件：当d有微小的扰动时，在不另加限制条件的情况下，解x只产生微小的变化。

如果上面三个条件任意一个不满足，式(6)就是病态(il1posed)的。以上三个解的存在性、唯一性和稳定性条件总称为反问题的适定性。

由于地震道d也就是角道集的数量(一般为几十道)，一般要远大于估计量的数目(纵波、横波增量和密度增量)，因此式(6)的方程为超定方程。对超定方程来说，是有解的。其解是无限趋近于最小二乘解，所有满足存在性的条件。但是由于传感器的原因，地震道d只有中频信息(一般为10-80Hz范围)，而没有低频和高频信息，因此测量误差是很大的，导致方程的解并不唯一。还有一个原因是矩阵G的条件数非常大，导致了式(6)的稳定性非常差。因此式(6)问题求解就是高度病态的反问题求解。

求解这样高度病态问题，需要发展出具有高稳定性、高抗测量误差的求解方法。为此，我们建立以下目标函数：

f(x)＝||d-Gx||₂          (7)

其中，||·||₂为二次泛函，该目标函数也可以应用于其它地震反演，比如叠后反演，只要采用其对应正演模型。

设在N维空间内，对于目标函数(7)，为了求取目标函数的极小值解，构造n个相互共轭的向量：P₀，P₁K，P_n-1做一维线性搜索，得到一极小值点序列x⁽¹⁾，x⁽²⁾Λx⁽ⁿ⁾，x^*；其中x^*即为通过n次线性搜索得到的目标函数的极小值点。

共轭梯度法具有二次截止的性质。共轭梯度法进行n次搜索就可求出最小解，因此它是一种快速，高精度的解法。若目标函数高于二次并为单峰值函数时，可以在n次搜索的基础上，再构造一组n个共轭向量，继续搜索。按照这种方式反复进行，直至达到要求的精度。一般的共轭梯度算法对庞大的系统是非常有效的。它仅仅要求计算梯度向量，而不要求计算所谓的海森矩阵，也不要求求解庞大的线性系统方程。

在共轭梯度法中n个共轭向量的构造十分关键，它既影响算法的速度，也影响存储量。目前，共轭向量多采用迭代构造法，即在计算过程中构造，在构造共轭向量前需要先构造出一组解空间的基g₀，g₁K，g_n-1。通常取目标函数的梯度：

$g_k=\▿f\left(x^k\right),$ (k＝0，1，...n-1)                   (8)

作为解空间的基。

采用Fletcher-Reeves算法对解空间的基进行迭代构造。具体如下：

x_k＝-g_k+α_kp_k-1，(k＝0，1，...n-1)         (9)

其中：

${\mathrm{\α}}_k=\frac{g_k^T{g}_k}{g_{k-1}^T{g}_{k-1}}---\left(10\right)$

为了保证反演的稳定性，我们引入下式

(G^TG+ε²I)Δx＝G^TΔd                    (11)

ΔX＝(G^TG)^-1G^TΔd                       (12)

其中，Δx为前后两次迭代更新差，ε为阻尼系数，Δd为地震数据真实值与d’之差，d’为第K次迭代结果x_k的正演值。

使用式(12)替换式(8)，将能够调整解的稳定性的阻尼系数ε引入到了共轭梯度算法中，利用式(8)来稳定地控制梯度向量的方向。

g_k＝(G^TG)^-1G^TΔd                       (13)

通过调整ε值能够提高解的稳定性，而式(9)则没有这样调整解的稳定性的机制。但是从式(11)来说，对x来说，则不容易加门限限制，因为式(11)是直接通过广义逆计算得到的，也就是式(11)得到的稳定解有可能并不是理想的解。而式(9)可以得到限制范围内的解，但是解不稳定。所以在此添加一个硬的门限：

假如m_L为下限模型向量和m_U为上限模型向量。那么式(9)使其满足：

m_L≤x_k≤m_U                            (14)

因此为了获得求解的唯一性和稳定性，该发明结合提出这样一种新的共轭梯度算法，能够在限制范围内求解得到一个稳定解。

第四步、反演后处理：

反演的目标是结合测井资料中的速度、密度信息，求解整个工区目标地层的速度，密度，泊松比。泊松比对储层预测起到很好的启示作用。第三步反演的解x提供了求解纵波速度，横波速度，密度的信息，再利用式(4)即可求得目标速度，密度结果，最后再利用式(15)求得泊松比，到此反演方法步骤结束。

$\mathrm{poisson}=\frac{L_p^2-2*L_s^2}{2(L_p^2-L_s^2)}---\left(15\right)$

Claims (3)

1.一种地震信号反演方法，其特征在于：包括如下步骤：

第一步、数据预处理；

第二步、建立正演模型；

第三步、以正演模型作为反演基础，对反问题进行求解：

1)建立目标函数；

2)在N维空间内，采用迭代构造法构造n个相互共轭的向量，做一维线性搜索，通过n次线性搜索得到目标函数的极小值点；

3)取目标函数的梯度构造出解空间的基；

4)对解空间的基进行迭代构造，引入阻尼系数调整解的稳定性，控制梯度向量的方向；

5)为每次迭代结果添加硬门限，在硬门限的限制范围内求解得到唯一的稳定解；

第四步、反演后处理：通过第三步求出的唯一的稳定解来求解纵波速度、横波速度和密度，从而求得目标速度和密度结果，最后求得泊松比。

2.根据权利要求1所述的地震信号反演方法，其特征在于：所述数据预处理包括如下步骤：

1)对原始地震数据进行常规处理；

2)对道集进行部分叠加，以减少道集复杂度和数据量；

3)对道集进行全叠加来生成叠后数据；

4)标定进行反演的目标层位；

5)从测井资料中获取所需信息。

3.根据权利要求1所述的地震信号反演方法，其特征在于：所述建立正演模型包括如下步骤：

1)采用单层Gidlow近似方程进行正演；

2)在动校正前的地震数据上进行反演，同时采用子波W描述褶积模型，得到基于Gidlow方程的褶积模型；

3)对褶积模型进行约束，得到带岩石物理约束关系的AVO褶积模型；

4)完成建立正演模型。