CN104237937A - 叠前地震反演方法及其系统 - Google Patents

Abstract

本发明公开了一种叠前地震反演方法及其系统，该叠前地震反演方法包括以下步骤：分别基于纵波阻抗、横波阻抗以及密度的初始模型求得各自对应的纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数；利用所得各反射系数，根据Fatti近似方程及摺积原理正演得到模拟叠前角度道集；基于模拟叠前角度道集和实际叠前角度道集来构建误差函数；根据实际叠前角度道集建立叠前反演目标模型；利用误差函数对叠前反演目标模型进行处理，反演得到纵波阻抗、横波阻抗及密度。本发明可以直接通过反演得到岩层密度，具有计算速度快、稳定性好的优点。

Description

叠前地震反演方法及其系统

技术领域

本发明涉及石油地球物理勘探领域，尤其涉及一种叠前地震反演方法及其系统。

背景技术

在石油勘探开发领域，常用的地震岩性分析方法有两种：一是建立在叠后地震资料基础之上的地震属性和地震振幅分析方法，称为叠后地震分析或叠后地震反演；二是建立在叠前地震资料基础上的地震属性和地震振幅分析方法，称为叠前地震分析或叠前地震(Amplitude Versus Offset，AVO)反演。其中，叠前地震分析由于其所使用的资料包括地下介质中纵、横波速度和密度的信息，以及丰富的振幅和旅行时信息，因此在地震岩性分析中获得了广泛的应用。

纵波速度、横波速度以及密度是岩石的基本弹性参数，在流体检测及估算储层气体含量时需要同时提供纵、横波速度和密度的信息。岩石物理实验表明，纵波速度对岩石的孔隙度及孔隙流体的变化较为敏感，而横波速度主要与岩石骨架有关，流体的变化对横波速度的变化影响较小。传统的叠前AVO反演方法中为了提高反演的稳定性，通常采用两参数方程进行反演，或在三参数方程中以经验关系式去除方程中的密度项，之后再通过岩石物理经验公式变换得到密度来进行反演。密度项的缺失虽然提高了反演的稳定性，但是也引入了较大的误差。

综上，亟需提供一种叠前地震反演方法，以解决传统叠前AVO反演中存在的问题，减小反演结果的误差。

发明内容

本发明所要解决的技术问题之一是需要提供一种叠前地震反演方法，该方法能够减小密度反演所引入的误差。

为了解决上述技术问题，本发明提供了一种叠前地震反演方法，包括以下步骤：步骤一，分别基于纵波阻抗初始模型、横波阻抗初始模型以及密度初始模型求得各自对应的纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数；步骤二，利用纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数，根据Fatti近似方程及摺积原理正演得到模拟叠前角度道集；步骤三，基于模拟叠前角度道集和实际叠前角度道集来构建误差函数；步骤四，根据实际叠前角度道集建立叠前反演目标模型；步骤五，利用误差函数对叠前反演目标模型进行处理，反演得到纵波阻抗、横波阻抗及密度。

在一个实施例中，在步骤二中，利用以下表达式来得到模拟叠前角度道集s：

s_i＝c_piWr_p+c_siWr_s+c_ρiWr_ρ

式中，s_i为模拟叠前角度道集s在入射角为i时的分量；c_pi，c_si，c_ρi为由入射角和速度比决定的系数；W为输入的由地震子波构建的N行N列的子波矩阵；r_p，r_s，r_ρ分别为纵波阻抗反射系数，横波阻抗反射系数，密度反射系数，且r_p＝0.5DL_p，r_s＝0.5DL_s，r_ρ＝DL_ρ，其中，L_p＝ln(Z_p)，L_s＝ln(Z_s)，L_ρ＝ln(ρ)，Z_p为纵波阻抗初始模型，Z_s为横波阻抗初始模型，ρ为密度初始模型，L_p为纵波阻抗初始模型的自然对数，L_s为横波阻抗初始模型的自然对数，L_ρ为密度初始模型的自然对数，D为N行N+1列的差分矩阵，其形式为：

$D={\left[\begin{array}{cccccc}-1& 1& 0& 0& ...& 0\\ 0& -1& 1& 0& ...& 0\\ 0& 0& -1& 1& ...& 0\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& 0& ...& -1& 1\end{array}\right]}_{N\×(N+1)}.$

在一个实施例中，在步骤三中，根据以下表达式构建误差函数f(X)：

$f\left(X\right)={\left|\right|s-d\left|\right|}_2=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[s_{i,j}-d_{i,j}]}^2$

式中，X为依次由L_p，L_s，L_ρ组成的复合矩阵；d为实际叠前角度道集；s_i,j为模拟叠前角度道集中入射角为i的地震记录中第j个样点值；d_i,j为实际叠前角度道集中入射角为i的地震记录中第j个样点值；M为入射角个数；N为地震记录样点个数。

在一个实施例中，误差函数f(X)所需满足的终止条件以及反演结果Xinv的表达式为如下形式：

f(X)≤ε

Xinv＝e^X

式中，X是满足终止条件f(X)≤ε的输出，ε是实数。

在一个实施例中，在步骤四中，根据实际叠前角度道集所建立的叠前反演目标模型如下：

$\left[\begin{array}{ccc}{0.5c}_{p1}\mathrm{WD}& {0.5c}_{s1}\mathrm{WD}& c_{\mathrm{\ρ}1}\mathrm{WD}\\ {0.5c}_{p2}\mathrm{WD}& {0.5c}_{s2}\mathrm{WD}& c_{\mathrm{\ρ}2}\mathrm{WD}\\ .& .& .\\ .& .& .\\ .& .& .\\ {0.5c}_{\mathrm{pM}}\mathrm{WD}& {0.5c}_{\mathrm{sM}}\mathrm{WD}& c_{\mathrm{\ρM}}\mathrm{WD}\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_{\mathrm{\ρ}}\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_M\end{array}\right]$

式中，c_pi，c_si，c_ρi为由入射角和速度比决定的系数；L_p，L_s，L_ρ是待反演的地层参数；d_i为实际叠前角度道集中入射角为i的地震记录。

在一个实施例中，采用LU分解法、追赶法、奇异值分解法、共轭梯度法、松弛迭代法中的任意一种方法对叠前反演目标模型进行处理。

另一方面，还提供了一种叠前地震反演系统，包括以下模块：第一模块，其分别基于纵波阻抗初始模型、横波阻抗初始模型以及密度初始模型求得各自对应的纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数；第二模块，其利用纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数，根据Fatti近似方程及摺积原理正演得到模拟叠前角度道集；第三模块，其基于模拟叠前角度道集和实际叠前角度道集来构建误差函数；第四模块，其根据实际叠前角度道集建立叠前反演目标模型；第五模块，其利用误差函数对叠前反演目标模型进行处理，反演得到纵波阻抗、横波阻抗及密度。

在一个实施例中，所述第二模块利用以下表达式来得到模拟叠前角度道集s：

s_i＝c_piWr_p+c_siWr_s+c_ρiWr_ρ

式中，s_i为模拟叠前角度道集s在入射角为i时的分量；c_pi，c_si，c_ρi为由入射角和速度比决定的系数；W为输入的由地震子波构建的N行N列的子波矩阵；r_p，r_s，r_ρ分别为纵波阻抗反射系数，横波阻抗反射系数，密度反射系数，且r_p＝0.5DL_p，r_s＝0.5DL_s，r_ρ＝DL_ρ，其中，L_p＝ln(Z_p)，L_s＝ln(Z_s)，L_ρ＝ln(ρ)，Z_p为纵波阻抗初始模型，Z_s为横波阻抗初始模型，ρ为密度初始模型，L_p为纵波阻抗初始模型的自然对数，L_s为横波阻抗初始模型的自然对数，L_ρ为密度初始模型的自然对数；D为N行N+1列的差分矩阵，其形式为：

$D={\left[\begin{array}{cccccc}-1& 1& 0& 0& ...& 0\\ 0& -1& 1& 0& ...& 0\\ 0& 0& -1& 1& ...& 0\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& 0& ...& -1& 1\end{array}\right]}_{N\×(N+1)}.$

在一个实施例中，所述第三模块根据以下表达式构建误差函数f(X)：

$f\left(X\right)={\left|\right|s-d\left|\right|}_2=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[s_{i,j}-d_{i,j}]}^2$

式中，X为依次由L_p，L_s，L_ρ组成的复合矩阵；d为实际叠前角度道集；s_i,j为模拟叠前角度道集中入射角为i的地震记录中第j个样点值；d_i,j为实际叠前角度道集中入射角为i的地震记录中第j个样点值；M为入射角个数；N为地震记录样点个数。

在一个实施例中，误差函数f(X)所需满足的终止条件以及反演结果Xinv的表达式为如下形式：

f(X)≤ε

Xinv＝e^X

式中，X是满足终止条件f(X)≤ε的输出，ε是实数。

在一个实施例中，所述第四模块根据实际叠前角度道集所建立的叠前反演目标模型如下：

$\left[\begin{array}{ccc}{0.5c}_{p1}\mathrm{WD}& {0.5c}_{s1}\mathrm{WD}& c_{\mathrm{\ρ}1}\mathrm{WD}\\ {0.5c}_{p2}\mathrm{WD}& {0.5c}_{s2}\mathrm{WD}& c_{\mathrm{\ρ}2}\mathrm{WD}\\ .& .& .\\ .& .& .\\ .& .& .\\ {0.5c}_{\mathrm{pM}}\mathrm{WD}& {0.5c}_{\mathrm{sM}}\mathrm{WD}& c_{\mathrm{\ρM}}\mathrm{WD}\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_{\mathrm{\ρ}}\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_M\end{array}\right]$

式中，c_pi，c_si，c_ρi为由入射角和速度比决定的系数；L_p，L_s，L_ρ是待反演的地层参数；d_i为实际叠前角度道集中入射角为i的地震记录。

在一个实施例中，所述第五模块采用LU分解法、追赶法、奇异值分解法、共轭梯度法、松弛迭代法中的任意一种方法对叠前反演目标模型进行处理。

与现有技术相比，上述方案中的一个或多个实施例可以具有如下优点或有益效果：

通过利用不同角度上的模拟叠前角度道集和实际叠前角度道集来构建误差函数，并求解满足终止条件的叠前反演目标方程，将叠前正演与叠前反演相结合，实现了从叠前角度道集到纵波阻抗、横波阻抗及密度等岩性参数的快速、稳定的反演过程。该方法可以直接通过反演得到岩层密度，具有计算速度快、精度高、稳定性好的优点。

本发明的其它特征和优点将在随后的说明书中阐述，并且部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1是根据本发明第一实施例的叠前地震反演方法的流程图；

图2(a)、(b)分别是具有不同地质结构的单层界面模型示意图；

图3是某勘探区域的实际测井数据图；

图4是对某勘探区域正演得到的模拟叠前角度道集的示意图；

图5是某勘探区域的实际叠前角度道集的示意图；

图6(a)、(b)、(c)分别是某勘探区域的叠前反演结果与对应的初始模型的对比图，图6(d)是某勘探区域的叠前反演速度比与实际测井速度比的对比图；

图7是叠前反演得到的密度与使用现有技术反演得到的密度和实际测井得到的密度的对比图；

图8是根据本发明第二实施例的叠前地震反演系统的结构示意图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

第一实施例

图1是根据本发明第一实施例的叠前地震反演方法的流程图。以下参考图1来详细说明本方法的各个步骤。

步骤S101、分别基于纵波阻抗初始模型、横波阻抗初始模型以及密度初始模型求得各自对应的纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数。

首先，将纵波阻抗初始模型、横波阻抗初始模型以及密度初始模型分别记为Z_p，Z_s，ρ，并分别对这三个初始模型求取各自的自然对数，如表达式(1)所示：

L_p＝ln(Z_p),L_s＝ln(Z_s),L_ρ＝ln(ρ)   (1)

然后，通过表达式(2)分别计算得到纵波阻抗反射系数、横波阻抗反射系数和密度反射系数：

r_pj＝0.5(L_pj+1-L_pj),r_sj＝0.5(L_sj+1-L_sj),r_ρj＝(L_ρj+1-L_ρj),(j＝1…N)   (2)

式中，r_pj，r_sj，r_ρj分别为第j层纵波阻抗反射系数，第j层横波阻抗反射系数，第j层密度反射系数，L_pj和L_pj+1为第j层和第j+1层纵波阻抗的自然对数，L_sj和L_sj+1为第j层和第j+1层横波阻抗的自然对数，L_ρj和L_ρj+1为第j层和第j+1层密度的自然对数，N为纵波阻抗、横波阻抗及密度的样点个数。

进一步的，表达式(2)可以被改写成矩阵的形式，如表达式(3)所示：

r_p＝0.5DL_p,r_s＝0.5DL_s,r_ρ＝DL_ρ   (3)

式中，r_p，r_s，r_ρ分别为纵波阻抗反射系数，横波阻抗反射系数，密度反射系数，L_p，L_s，L_ρ分别为纵波阻抗初始模型的自然对数，横波阻抗初始模型的自然对数，密度初始模型的自然对数，各向量的形式如下：

$\left\{\begin{array}{c}{r}_p={\left[\begin{array}{cccc}{r}_{p1}& r_{p2}& ...& r_{\mathrm{pN}}\end{array}\right]}^T\\ r_s={\left[\begin{array}{cccc}{r}_{s1}& r_{s2}& ...& r_{\mathrm{sN}}\end{array}\right]}^T\\ r_{\mathrm{\ρ}}={\left[\begin{array}{cccc}{r}_{\mathrm{\ρ}1}& r_{\mathrm{\ρ}2}& ...& r_{\mathrm{\ρN}}\end{array}\right]}^T\end{array}\right.,\left\{\begin{array}{c}{L}_p={\left[\begin{array}{ccccc}{L}_{p1}& L_{p2}& ...& L_{\mathrm{pN}}& L_{\mathrm{pN}+1}\end{array}\right]}^T\\ L_s={\left[\begin{array}{ccccc}{L}_{s1}& L_{s2}& ...& L_{\mathrm{sN}}& L_{\mathrm{sN}+1}\end{array}\right]}^T\\ L_{\mathrm{\ρ}}={\left[\begin{array}{ccccc}{L}_{\mathrm{\ρ}1}& {\mathrm{L\ρ}}_2& ...& L_{\mathrm{\ρN}}& L_{\mathrm{\ρN}+1}\end{array}\right]}^T\end{array}\right.$

其中，D是一个N行N+1列的差分矩阵，具体如下表示：

$D={\left[\begin{array}{cccccc}-1& 1& 0& 0& ...& 0\\ 0& -1& 1& 0& ...& 0\\ 0& 0& -1& 1& ...& 0\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& 0& ...& -1& 1\end{array}\right]}_{N\×(N+1)}.$

另外，需要说明的是，上述纵波阻抗初始模型、横波阻抗初始模型及密度初始模型均是以实际测井资料为基础，通过平滑或插值得到的数据，是事先给定的。如表1和图3所示。

表1 单层界面模型的岩性参数

表1中给出的是泥岩与含水砂岩或含气砂岩的岩性参数，是一组实验室数据。与其对应的理论单界面地质模型如图2(a)和图2(b)所示。图2(a)显示的是上层储层为泥岩，下层储层为含水砂岩的情况，图2(b)显示的是上层储层为泥岩，下层储层为含气砂岩的情况。

由于纵波阻抗为纵波速度与密度的乘积，横波阻抗为横波速度与密度的乘积，因此可根据表1中的参数分别得到纵波阻抗初始模型Z_p、横波阻抗初始模型Z_s以及密度初始模型ρ。

图3显示的是某勘探区域的实际测井数据图，这是一种多层地层界面、多个偏移距的复杂情形。每条曲线都各自具有546个样点，上下每两个样点就构成了一个单界面模型，据此可以分别得到纵波阻抗初始模型Z_p、横波阻抗初始模型Z_s以及密度初始模型ρ。

步骤S102、利用纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数，根据Fatti近似方程及摺积原理正演得到模拟叠前角度道集。

具体地，Fatti近似方程是Zoeppritz方程的一种近似形式，用于计算入射角为i时的反射系数r_i,j，其形式如表达式(4)所示：

r_i,j＝c_pir_pj+c_sir_sj+c_ρir_ρj,(i＝1...M,j＝1...N)   (4)

式中，r_i,j为入射角为i时的反射系数r的第j个分量，c_pi，c_si，c_ρi为由入射角和速度比决定的系数，M为入射角个数。各系数的具体计算如表达式(5)所示：

c_pi＝1+tan²i,c_si＝-8(αsini)²,c_ρi＝0.5tan²i-2(αsini)²   (5)

式中，i为入射角，α为速度比，α一般取0.5。

进一步地，可以将表达式(4)改写成矩阵的形式r_i＝c_pi·r_p+c_si·r_s+c_ρi·r_ρ，再根据褶积原理，得到表达式(6)：

s_i＝w*r_i＝Wr_i   (6)

式中，s_i为正演得到的模拟叠前角度道集s在入射角为i时的分量，w为地震子波，是事先给定的输入数据，r_i为入射角为i时的反射系数，是一个M行1列的矩阵，即r_i(M×1)。W是由地震子波w构建的N行N列的子波矩阵，即W_N×N。

还可以将矩阵形式的r_i带入表达式(6)中，得到表达式(7)：

s_i＝c_piWr_p+c_siWr_s+c_ρiWr_ρ   (7)

依次对不同角度的叠前记录进行正演，便得到模拟叠前角度道集s：

$\left[\begin{array}{ccc}{0.5c}_{p1}\mathrm{WD}& {0.5c}_{s1}\mathrm{WD}& c_{\mathrm{\ρ}1}\mathrm{WD}\\ {0.5c}_{p2}\mathrm{WD}& {0.5c}_{s2}\mathrm{WD}& c_{\mathrm{\ρ}2}\mathrm{WD}\\ .& .& .\\ .& .& .\\ .& .& .\\ {0.5c}_{\mathrm{pM}}\mathrm{WD}& {0.5c}_{\mathrm{sM}}\mathrm{WD}& c_{\mathrm{\ρM}}\mathrm{WD}\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_{\mathrm{\ρ}}\end{array}\right]=\left[\begin{array}{c}{s}_1\\ s_2\\ .\\ .\\ .\\ s_M\end{array}\right]---\left(8\right)$

图4所示即为对图3所示测井数据利用表达式(7)进行正演所得到的模拟叠前角度道集。将其与图5所示的实际叠前角度道集相比较，可见两者非常接近。于是，验证了如下事实：利用表达式(7)正演得到的结果是符合实际情况的，说明正演建模是正确的，也说明井震标定是合理的。同时也反映出如果我们能反演得到一组参数(纵波阻抗、横波阻抗、密度)，使得这组参数正演出的模拟叠前角度道集与实际叠前角度道集匹配良好，那么这组参数就可以作为反演结果输出。基于上述考虑，进一步实施步骤S103。

步骤S103、基于模拟叠前角度道集和实际叠前角度道集来构建误差函数。

具体地，可以利用如下表达式来构建误差函数f(X)：

$f\left(X\right)={\left|\right|s-d\left|\right|}_2=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[s_{i,j}-d_{i,j}]}^2---\left(9\right)$

式中，X为依次由L_p，L_s，L_ρ组成的复合矩阵，即X(L_p，L_s，L_ρ)，s为模拟叠前角度道集，d为实际叠前角度道集，s_i,j为模拟叠前角度道集中入射角为i的地震记录中第j个样点值，d_i,j为实际叠前角度道集中入射角为i的地震记录中第j个样点值，M为入射角个数，N为地震记录样点个数。

步骤S104、根据实际叠前角度道集建立叠前反演目标模型。

将实际叠前角度道集带入表达式(8)，可以得到叠前反演目标模型，具体如表达式(10)所示：

$\left[\begin{array}{ccc}{0.5c}_{p1}\mathrm{WD}& {0.5c}_{s1}\mathrm{WD}& c_{\mathrm{\ρ}1}\mathrm{WD}\\ {0.5c}_{p2}\mathrm{WD}& {0.5c}_{s2}\mathrm{WD}& c_{\mathrm{\ρ}2}\mathrm{WD}\\ .& .& .\\ .& .& .\\ .& .& .\\ {0.5c}_{\mathrm{pM}}\mathrm{WD}& {0.5c}_{\mathrm{sM}}\mathrm{WD}& c_{\mathrm{\ρM}}\mathrm{WD}\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_{\mathrm{\ρ}}\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_M\end{array}\right]---\left(10\right)$

需要注意的是，L_p，L_s，L_ρ虽然与步骤S101中输入的纵波阻抗初始模型的自然对数，横波阻抗初始模型的自然对数以及密度初始模型的自然对数在形式上相同，但此处它们是未知量，是待反演的地层参数。d_i为实际叠前角度道集中入射角为i的地震记录。

步骤S105、利用误差函数对叠前反演目标模型进行处理，反演得到纵波阻抗、横波阻抗及密度。

为了得到最终的反演结果，求解当误差函数满足终止条件时的叠前反演目标方程组(10)，输出X(L_p，L_s，L_ρ)，对X取指数，得到反演结果Xinv。终止条件如表达式(11)所示，反演结果由表达式(12)得到。

f(X)≤ε   (11)

Xinv＝e^X   (12)

式中，ε是一个较小的实数，如0.1；Xinv是反演结果，依次包括了纵波阻抗、横波阻抗及密度。

叠前反演目标模型是线性方程组，可以采用常用线性代数方程组的解法进行求解，例如LU分解法、追赶法、奇异值分解法、共轭梯度法、松弛迭代法等等。

根据上述步骤S104及步骤S105对图5所示的实际叠前角度道集进行叠前反演。图5是实际中某勘探区域的井旁角度道集，45道，每道1度，入射角范围是1度至45度。反演得到纵波阻抗、横波阻抗和密度，其结果示于图6(a)、图6(b)和图6(c)中。

图6(a)显示的是纵波阻抗反演结果，图中曲线1是纵波阻抗初始模型曲线，曲线2是纵波阻抗反演结果曲线，从图中可以看出，反演结果在初始模型附近呈波动状分布，进一步分析表明，波动反映的是底层中的薄层信息，可将其直接用于岩层分析。

图6(b)和图6(c)分别显示了横波阻抗和密度的反演结果，其显示的信息同图6(a)，不再赘述。

图6(d)是根据反演结果计算的纵横波速度比与实际测井得到的纵横波速度比的对比图。其中曲线1是由实际测井数据得到的纵横波速度比，可以根据图3中的纵波阻抗除以横波阻抗得到。曲线2是根据反演结果计算的纵横波速度比，由反演得到的纵波阻抗除以横波阻抗得到。从对比图中可以看出，两条纵横波速度比曲线匹配良好，说明本发明方法用于勘探预测的有效性。例如在样点100、150、500处，由于这几处都属于含气储层，纵波速度会降低，但横波速度不会发生明显变化，所以纵横波速度比降低。

速度比在储层流体识别中具有重要应用，岩石物理实验表明，纵波速度对岩石的孔隙度及孔隙流体的变化较为敏感，而横波速度主要与岩石骨架有关，流体的变化对横波速度的影响较小。通常饱含流体的储层对应较小的速度比，因此利用速度比可以对储层流体进行检测。

图7是根据第一实施例的方法反演得到的密度与根据常规方法反演得到的密度，以及根据实际测井数据得到的密度的对比图。图中曲线1是常规反演得到的密度，曲线2是实际测井数据的密度，曲线3是利用第一实施例的叠前反演方法得到密度。图中椭圆标注部分中可以看出，曲线1常规反演结果与曲线2实际测量的密度相差较大，而叠前反演得到的密度曲线3与实际测量的密度更为接近。

第二实施例

图8是根据本发明第二实施例的叠前地震反演系统的结构示意图。下面参考图8来说明本实施例的各部分组成。

具体在图8中，本实施例的第一模块81执行第一实施例的步骤S101的操作，第二模块82执行第一实施例的步骤S102的操作，第三模块83执行第一实施例的步骤S103的操作，第四模块84执行第一实施例的步骤S104的操作，第五模块85执行第一实施例的步骤S105的操作。在此不再详细展开。

本领域的技术人员应该明白，上述的本发明的各模块或各步骤可以用通用的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。这样，本发明不限制于任何特定的硬件和软件结合。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (12)

1.一种叠前地震反演方法，其特征在于，包括以下步骤：

步骤一，分别基于纵波阻抗初始模型、横波阻抗初始模型以及密度初始模型求得各自对应的纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数；

步骤二，利用纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数，根据Fatti近似方程及摺积原理正演得到模拟叠前角度道集；

步骤三，基于模拟叠前角度道集和实际叠前角度道集来构建误差函数；

步骤四，根据实际叠前角度道集建立叠前反演目标模型；

步骤五，利用误差函数对叠前反演目标模型进行处理，反演得到纵波阻抗、横波阻抗及密度。

2.如权利要求1所述的方法，其特征在于，在步骤二中，利用以下表达式来得到模拟叠前角度道集s：

s_i＝c_piWr_p+c_siWr_s+c_ρiWr_ρ

式中，s_i为模拟叠前角度道集s在入射角为i时的分量；c_pi，c_si，c_ρi为由入射角和速度比决定的系数；W为输入的由地震子波构建的N行N列的子波矩阵；r_p，r_s，r_ρ分别为纵波阻抗反射系数，横波阻抗反射系数，密度反射系数，且r_p＝0.5DL_p，r_s＝0.5DL_s，r_ρ＝DL_ρ，其中，L_p＝ln(Z_p)，L_s＝ln(Z_s)，L_ρ＝ln(ρ)，Z_p为纵波阻抗初始模型，Z_s为横波阻抗初始模型，ρ为密度初始模型，L_p为纵波阻抗初始模型的自然对数，L_s为横波阻抗初始模型的自然对数，L_ρ为密度初始模型的自然对数，D为N行N+1列的差分矩阵，其形式为：

$D={\left[\begin{array}{cccccc}-1& 1& 0& 0& ...& 0\\ 0& -1& 1& 0& ...& 0\\ 0& 0& -1& 1& ...& 0\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& 0& ...& -1& 1\end{array}\right]}_{N\×(N+1)}.$

3.如权利要求2所述的方法，其特征在于，在步骤三中，根据以下表达式构建误差函数f(X)：

$f\left(X\right)={\left|\right|s-d\left|\right|}_2=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[s_{i,j}-d_{i,j}]}^2$

式中，X为依次由L_p，L_s，L_ρ组成的复合矩阵；d为实际叠前角度道集；s_i,j为模拟叠前角度道集中入射角为i的地震记录中第j个样点值；d_i,j为实际叠前角度道集中入射角为i的地震记录中第j个样点值；M为入射角个数；N为地震记录样点个数。

4.如权利要求3所述的方法，其特征在于，误差函数f(X)所需满足的终止条件以及反演结果Xinv的表达式为如下形式：

f(X)≤ε

Xinv＝e^X

式中，X是满足终止条件f(X)≤ε的输出，ε是实数。

5.如权利要求2所述的方法，其特征在于，在步骤四中，根据实际叠前角度道集所建立的叠前反演目标模型如下：

$\left[\begin{array}{ccc}{0.5c}_{p1}\mathrm{WD}& {0.5c}_{s1}\mathrm{WD}& c_{\mathrm{\ρ}1}\mathrm{WD}\\ {0.5c}_{p2}\mathrm{WD}& {0.5c}_{s2}\mathrm{WD}& c_{\mathrm{\ρ}2}\mathrm{WD}\\ .& .& .\\ .& .& .\\ .& .& .\\ {0.5c}_{\mathrm{pM}}\mathrm{WD}& {0.5c}_{\mathrm{sM}}\mathrm{WD}& c_{\mathrm{\ρM}}\mathrm{WD}\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_{\mathrm{\ρ}}\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_M\end{array}\right]$

式中，c_pi，c_si，c_ρi为由入射角和速度比决定的系数；L_p，L_s，L_ρ是待反演的地层参数；d_i为实际叠前角度道集中入射角为i的地震记录。

6.如权利要求5所述的方法，其特征在于，采用LU分解法、追赶法、奇异值分解法、共轭梯度法、松弛迭代法中的任意一种方法对叠前反演目标模型进行处理。

7.一种叠前地震反演系统，其特征在于，包括以下模块：

第一模块，其分别基于纵波阻抗初始模型、横波阻抗初始模型以及密度初始模型求得各自对应的纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数；

第二模块，其利用纵波阻抗反射系数、横波阻抗反射系数以及密度反射系数，根据Fatti近似方程及摺积原理正演得到模拟叠前角度道集；

第三模块，其基于模拟叠前角度道集和实际叠前角度道集来构建误差函数；

第四模块，其根据实际叠前角度道集建立叠前反演目标模型；

第五模块，其利用误差函数对叠前反演目标模型进行处理，反演得到纵波阻抗、横波阻抗及密度。

8.如权利要求7所述的系统，其特征在于，所述第二模块利用以下表达式来得到模拟叠前角度道集s：

s_i＝c_piWr_p+c_siWr_s+c_ρiWr_ρ

式中，s_i为模拟叠前角度道集s在入射角为i时的分量；c_pi，c_si，c_ρi为由入射角和速度比决定的系数；W为输入的由地震子波构建的N行N列的子波矩阵；r_p，r_s，r_ρ分别为纵波阻抗反射系数，横波阻抗反射系数，密度反射系数，且r_p＝0.5DL_p，r_s＝0.5DL_s，r_ρ＝DL_ρ，其中，L_p＝ln(Z_p)，L_s＝ln(Z_s)，L_ρ＝ln(ρ)，Z_p为纵波阻抗初始模型，Z_s为横波阻抗初始模型，ρ为密度初始模型，L_p为纵波阻抗初始模型的自然对数，L_s为横波阻抗初始模型的自然对数，L_ρ为密度初始模型的自然对数，D为N行N+1列的差分矩阵，其形式为：

$D={\left[\begin{array}{cccccc}-1& 1& 0& 0& ...& 0\\ 0& -1& 1& 0& ...& 0\\ 0& 0& -1& 1& ...& 0\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& 0& ...& -1& 1\end{array}\right]}_{N\×(N+1)}.$

9.如权利要求8所述的系统，其特征在于，所述第三模块根据以下表达式构建误差函数f(X)：

$f\left(X\right)={\left|\right|s-d\left|\right|}_2=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{[s_{i,j}-d_{i,j}]}^2$

式中，X为依次由L_p，L_s，L_ρ组成的复合矩阵；d为实际叠前角度道集；s_i,j为模拟叠前角度道集中入射角为i的地震记录中第j个样点值；d_i,j为实际叠前角度道集中入射角为i的地震记录中第j个样点值；M为入射角个数；N为地震记录样点个数。

10.如权利要求9所述的系统，其特征在于，误差函数f(X)所需满足的终止条件以及反演结果Xinv的表达式为如下形式：

f(X)≤ε

Xinv＝e^X

式中，X是满足终止条件f(X)≤ε的输出，ε是实数。

11.如权利要求8所述的系统，其特征在于，所述第四模块根据实际叠前角度道集所建立的叠前反演目标模型如下：

$\left[\begin{array}{ccc}{0.5c}_{p1}\mathrm{WD}& {0.5c}_{s1}\mathrm{WD}& c_{\mathrm{\ρ}1}\mathrm{WD}\\ {0.5c}_{p2}\mathrm{WD}& {0.5c}_{s2}\mathrm{WD}& c_{\mathrm{\ρ}2}\mathrm{WD}\\ .& .& .\\ .& .& .\\ .& .& .\\ {0.5c}_{\mathrm{pM}}\mathrm{WD}& {0.5c}_{\mathrm{sM}}\mathrm{WD}& c_{\mathrm{\ρM}}\mathrm{WD}\end{array}\right]\left[\begin{array}{c}{L}_p\\ L_s\\ L_{\mathrm{\ρ}}\end{array}\right]=\left[\begin{array}{c}{d}_1\\ d_2\\ .\\ .\\ .\\ d_M\end{array}\right]$

式中，c_pi，c_si，c_ρi为由入射角和速度比决定的系数；L_p，L_s，L_ρ是待反演的地层参数；d_i为实际叠前角度道集中入射角为i的地震记录。

12.如权利要求11所述的系统，其特征在于，所述第五模块采用LU分解法、追赶法、奇异值分解法、共轭梯度法、松弛迭代法中的任意一种方法对叠前反演目标模型进行处理。