CN104280767A - 一种基于柯西分布的稀疏脉冲反演方法 - Google Patents

Abstract

一种基于柯西分布的稀疏脉冲反演方法，拾取得到层位数据，校正并内插和平滑；根据叠后地震数据、层位数据和已知的钻井分层数据，把深度域的声波时差曲线和密度曲线标定为时间域的曲线，生成波阻抗曲线，并在井旁地震道上提取地震子波；利用层位数据和时间域波阻抗曲线，生成初始波阻抗模型；对地震数据和子波数据分别进行归一化，反演得到道反射系数序列，推导得到波阻抗与反射系数的关系，最终得到反演的波阻抗剖面。本发明反演结果的准确性和稳定性好，能够正确反映岩层的物理特征，不易出现假象且运算效率高。

Description

一种基于柯西分布的稀疏脉冲反演方法

技术领域

本发明涉及地球物理勘探技术，属于储层预测反演技术类，是一种基于柯西分布的稀疏脉冲反演方法。

背景技术

地震勘探是通过人工激发地震波，在地表或地下用单分量或多分量传感器记录地震波的地层响应，研究它们在地层中的传播规律，通过地震数据处理及反演等方法以查明地下的地质构造岩性特征，进而寻找油气等矿产资源的地球物理勘探方法。地震勘探从认识地下的构造形态开始，发展到直接应用地震信息判断岩性、分析岩相、定量计算岩层的物性参数及寻找油气显示等。地震反演技术正是这一发展过程的产物。

地震反演的基本目的是利用地震波在地下介质中的传播规律，通过数据采集、处理与解释等流程，来推测地下岩层结构和物性参数的空间分布，为勘探开发提供重要依据。在地震反演研究中，有多种参数的反演方法，如波阻抗、速度、密度、孔隙度、渗透率、泊松比等。由于波阻抗信息是联系地质和地球物理的一座桥梁，在叠后计算数据量相对要小，在实际生产中应用方便而且效果明显，因此波阻抗反演在地震反演中具有特殊的地位，地震反演通常是指波阻抗反演。

常规的地震波阻抗反演就是指利用地震贤料反演地层/岩层波阻抗的地震特殊处理技术。与地震多参数模式识别预测储层油气、神经网络预测地层参数、振幅拟合预测储层厚度等统计性方法相比，波阻抗反演具有明确的物理意义，它是储层岩性预测、油藏特征描述的确定性方法，在实际应用中取得了显著的地质效果。

波阻抗反演技术在20世纪70年代开始出现，当时对地震反演的研究只是以基于褶积模型的叠后一维波阻抗反演为主，80年代得到了蓬勃发展。1983年，Cooke介绍了地震资料的广义线性反演方法，从而揭开了波阻抗反演技术的新篇章。周竹生等人在90年代初期提出了综合利用地质、地震和测井资料进行约束反演，克服单一的线性反演方法的缺陷。90年代中期，李宏兵提出了递推反演与宽带约束反演结合的方法。该方法的推出，解决了从单道出发的反演方法不能在根本上消除噪音的困惑。在此基础之上，有人进行了无井多道反演和有井多道反演的研究，使波阻抗反演方法更加完善。同时，非线性反演方法也得到了极大的发展，自从Rothmans利用模拟退火法成功地解决了地震勘探中大静校正量估算问题之后，模拟退火法受到地球物理上作者们的重视，人们开始将其应用到波阻抗反演中。遗传算法和模拟退火算法一样，利用转移概率来指导全局搜索过程。90年代至今，围绕一维波阻抗反演的各类算法以及应用成果层出不穷，随着研究的升温，在1997年左右开始出现了一些反思的文章，指出了波阻抗反演中存在的一些缺陷，并提出了一些解决方案。但是随着这类方法的引进和在生产单位的推广应用，我们逐渐发现，这些方法在实用中很难获得好的效果，很难获得可信的反演结果。利用这些反演结果进行的储层横向预测，在己知的井点还比较可靠，离开井点一定距离后，往往与实际钻井结果相差甚远，预测出的薄层、薄互层及储层岩性的可信度很低。

目前的反演方法多数是以模型为基础的方法，这些方法一般都依据测井及地质资料建立初始模型，通过广义线性反演方法进行迭代求取岩性参数。由于该问题的非线性，所以除了要求精细的子波外，还要求初始模型接近真实模型，才能达到可靠的结果，即反演结果强烈依赖于初始模型的选择。除此类方法外，全局优化的反演方法（如遗传算法和模拟退火算法等）虽然克服了基于模型方法对初始模型依赖性强的缺陷，但是由于其得到全局最优的反演结果，所以反演速度很慢。

发明内容

本发明目的是提供一种运算效率非常高的基于柯西分布的叠后波阻抗反演方法。

本发明通过如下技术方案实现：

1）采用常规的地震勘探方法采集地震数据，对地震数据进行常规处理得到叠后地震数据；

2）对叠后地震数据进行层位拾取得到层位数据，对确定的目的层层位进行检验和校正以及内插和平滑；

3）采用常规的测井方法得到测井数据，得到测井声波时差曲线和密度曲线；

4）根据叠后地震数据、层位数据和已知的钻井分层数据，把深度域的声波时差曲线和密度曲线标定为时间域的曲线，同时生成井中的时间域波阻抗曲线数据，并在井旁地震道上提取地震子波；

所述的标定为利用测井曲线和地震子波模拟井旁地震记录，实现测井分层到地震层位的标定和映射，由此得到时深关系曲线，由此时深关系可以将深度域的测井曲线转换为时间域曲线。

5）利用步骤2）的层位数据和步骤4）得到的时间域波阻抗曲线，生成初始波阻抗模型；

6）对地震数据和步骤4）得到的子波数据分别进行归一化，归一到范围[-1,1]之间；

7）把步骤1）读入的一道地震数据、步骤5）生成的初始波阻抗模型的一道数据和步骤4）提取的子波数据，代入以下公式中，通过反演得到该道的反射系数序列：

r=(G^TG+μQ+ρC^TC)^-1(G^Td+ρC^Tξ)   （1）

其中d=[d₁,d₂,...,d_N]^T是叠后地震数据，N为地震数据的总采样点数；r=[r₁,r₂,...,r_N]^T是反射系数序列；G是N×N维子波褶积矩阵，上角标T代表矩阵的转置；μ是稀疏约束因子，控制反射系数的稀疏程度；矩阵Q的对角元素为：Q_nn表示矩阵Q的第n行第n列元素的值，n为Q矩阵的行列号，矩阵Q除对角元素外其余元素均为零，σ_r代表噪声的标准方差，R_n是由初始波阻抗模型计算得到的第n个采样点位置处的初始反射系数；ρ为模型约束因子，控制反演结果对初始模型的依赖程度；C为积分算子矩阵，其离散形式表示为：

公式（1）中上角标-1为对矩阵求逆；为第n个采样点处的相对波阻抗值，I₀为反演时窗内第一个采样点对应的初始波阻抗值，I_n为反演时窗内第n个采样点处的初始波阻抗值，ln为自然对数符号，r_i为反演时窗内第i个采样点的反射系数值，表示对从第1个采样点到第n个采样点的r_i进行求和计算。

8）通过步骤7）中第n个采样点处的相对波阻抗的定义，通过推导可以得到波阻抗与反射系数的关系：

$I_n=I_0{e}^{2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r}_i}---\left(2\right)$

其中：I_n为反演时窗内第n个采样点的波阻抗值，I₀为反演时窗内第一个采样点对应的初始波阻抗值，r_i为步骤7）中反演得到的第i个采样点反射系数值，e表示自然对数的底，表示对从第1个采样点到第n个采样点的r_i进行求和计算。

通过公式（2）将步骤7）反演得到的一道的反射系数系列转换为波阻抗序列，便可得到该道的波阻抗反演结果；

9）对所有的地震道重复步骤7）至8）过程，得到所有道的波阻抗反演结果。

本发明具有如下特点，主要表现为：

（1）通过假设反射系数的先验概率服从柯西分布来实现，柯西分布属于长尾巴分布，相对于高斯分布来说其峰值处分布更窄，且逼近零的速度也较缓，从而得到少量的非零值和大量的零值来实现稀疏脉冲反演；

（2）通过添加波阻抗模型约束可以控制反演结果的准确性和稳定性；

（3）稀疏脉冲反演是以地震道为主的反演方法，反演结果的分辨率、信噪比以及可靠程度主要依赖于地震资料本身的品质，地震噪音对反演结果敏感，且影响大，因此用于稀疏脉冲反演的地震资料应具有较宽的频带、较低的噪声、相对振幅保持和成像准确等特征。测井资料，尤其是声波测井和密度测井资料，是地震横向预测的对比标准和解释依据，在反演处理之前应进行仔细的编辑和校正，使其能够正确反映岩层的物理特征；

（4）稀疏脉冲反演得到的波阻抗结果在井位处与井曲线吻合度高，且运算效率高。

本发明的反演方法不仅具有一般递推反演方法的特点，即反演结果忠实于地震资料，能反映储层的横向变化。而且，在迭代过程中引入地质和测井资料参与反演约束，增加了部分低频和高频成分，一定程度拓宽了反演频带。该方法对初始模型依赖较小，反演结果的唯一性较好，不易出现假象。

附图说明

图1是本发明实例1叠后地震剖面图；

图2是本发明实例1波阻抗剖面；

图3是本发明实例2某地区反演波阻抗结果与Jason软件反演结果对比，a为本发明反演结果，b为Jason反演结果。

具体实施方式

本发明结合附图反演效果，具体详述说明。

本发明所提供的方法，是一种基于柯西分布的稀疏脉冲反演方法，而且运算效率非常高。

本发明通过如下技术方案实现：

1）采用常规的地震勘探方法采集地震数据，对地震数据进行常规处理得到叠后地震数据；

2）对叠后地震数据进行层位拾取得到层位数据，对确定的目的层层位进行检验和校正以及内插和平滑；

3）采用常规的测井方法得到测井数据，得到测井声波时差曲线和密度曲线；

4）根据叠后地震数据、层位数据和已知的钻井分层数据，把深度域的声波时差曲线和密度曲线标定为时间域的曲线，同时生成井中的时间域波阻抗曲线数据，并在井旁地震道上提取地震子波；

所述的标定为利用测井曲线和地震子波模拟井旁地震记录，实现测井分层到地震层位的标定和映射，由此得到时深关系曲线，由此时深关系可以将深度域的测井曲线转换为时间域曲线。

5）利用步骤2）的层位数据和步骤4）得到的时间域波阻抗曲线，生成初始波阻抗模型；

6）对地震数据和步骤4）得到的子波数据分别进行归一化，归一到范围[-1,1]之间；

7）把步骤1）读入的一道地震数据、步骤5）生成的初始波阻抗模型的一道数据和步骤4）提取的子波数据，代入以下公式中，通过反演得到该道的反射系数序列：

r=(G^TG+μQ+ρC^TC)^-1(G^Td+ρC^Tξ)   （1）

其中d=[d₁,d₂,...,d_N]^T是叠后地震数据，N为地震数据的总采样点数；r=[r₁,r₂,...,r_N]^T是反射系数序列；G是N×N维子波褶积矩阵，上角标T代表矩阵的转置；μ是稀疏约束因子，控制反射系数的稀疏程度；矩阵Q的对角元素为：Q_nn表示矩阵Q的第n行第n列元素的值，n为Q矩阵的行列号，矩阵Q除对角元素外其余元素均为零，σ_r代表噪声的标准方差，R_n是由初始波阻抗模型计算得到的第n个采样点位置处的初始反射系数；ρ为模型约束因子，控制反演结果对初始模型的依赖程度；C为积分算子矩阵，其离散形式表示为：

公式（1）中上角标-1为对矩阵求逆；为第n个采样点处的相对波阻抗值，I₀为反演时窗内第一个采样点对应的初始波阻抗值，I_n为反演时窗内第n个采样点处的初始波阻抗值，ln为自然对数符号，r_i为反演时窗内第i个采样点的反射系数值，表示对从第1个采样点到第n个采样点的r_i进行求和计算。

8）通过步骤7）中第n个采样点处的相对波阻抗的定义，通过推导可以得到波阻抗与反射系数的关系：

$I_n=I_0{e}^{2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r}_i}---\left(2\right)$

其中：I_n为反演时窗内第n个采样点的波阻抗值，I₀为反演时窗内第一个采样点对应的初始波阻抗值，r_i为步骤7）中反演得到的第i个采样点反射系数值，e表示自然对数的底，表示对从第1个采样点到第n个采样点的r_i进行求和计算。

通过公式（2）将步骤7）反演得到的一道的反射系数系列转换为波阻抗序列，便可得到该道的波阻抗反演结果；

9）对所有的地震道重复步骤7）至8）过程，得到所有道的波阻抗反演结果。

以下通过具体实例说明本发明的效果。

实例1：本发明对一个二维模型数据（来自Jason公司）进行试验，该模型垂直方向是三个菱形砂体，每个砂体内部从上到下所含流体分别为气、油和水。图1为该模型的无噪音叠后地震剖面，依据前面的步骤1）-9）进行基于柯西分布的稀疏脉冲反演，得到图2的波阻抗剖面。由图2可以看出反演结果在井位处与井曲线吻合的很好，并且从图上可以清晰的反映出砂体内部的阻抗差异，说明反演结果很准确。

反演时窗为1250毫秒，采样间隔为4毫秒，共200道，用本发明在HP2单机上进行反演耗时约110秒。

实例2：本发明对某地区进行反演，依据前面的步骤1）-9）进行基于柯西分布的稀疏脉冲反演，图3(a)为本发明得到的波阻抗剖面，图3(b)为使用国外软件Jason系统进行反演得到的波阻抗剖面，从图中可看出总体上本发明反演结果与Jason系统反演结果相当。

反演目的层时窗约900毫秒，采样间隔为4毫秒，共300道，用本发明在HP2单机上进行反演耗时约340秒。

理论模型和实际资料的结果表明，柯西分布不仅能较好地恢复反射系数序列的稀疏性，同时能在提高地震资料的分辨率和减小对弱反射信息的压制之间达到一个平衡。由该方法得到的反演结果不仅准确，精度高，而且运算效率很高，利用此结果可以进行准确的储层预测。

Claims (2)

1.一种基于柯西分布的叠后波阻抗反演方法，特点是通过以下步骤实现：

1）采用常规的地震勘探方法采集地震数据，对地震数据进行常规处理得到叠后地震数据；

2）对叠后地震数据进行层位拾取得到层位数据，对确定的目的层层位进行检验和校正以及内插和平滑；

3）采用常规的测井方法得到测井数据，得到测井声波时差曲线和密度曲线；

4）根据叠后地震数据、层位数据和已知的钻井分层数据，把深度域的声波时差曲线和密度曲线标定为时间域的曲线，同时生成井中的时间域波阻抗曲线，并在井旁地震道上提取地震子波；

5）利用步骤2）的层位数据和步骤4）得到的时间域波阻抗曲线，生成初始波阻抗模型；

6）对地震数据和步骤4）得到的子波数据分别进行归一化，归一到范围[-1,1]之间；

7）把步骤1）读入的一道地震数据、步骤5）生成的初始波阻抗模型的一道数据和步骤4）提取的子波数据，代入以下公式中，通过反演得到该道的反射系数序列：

r=(G^TG+μQ+ρC^TC)^-1(G^Td+ρC^Tξ)   （1）

其中d=[d₁,d₂,...,d_N]^T是叠后地震数据，N为地震数据的总采样点数；r=[r₁,r₂,...,r_N]^T是反射系数序列；G是N×N维子波褶积矩阵，上角标T代表矩阵的转置；μ是稀疏约束因子，控制反射系数的稀疏程度；矩阵Q的对角元素为：Q_nn表示矩阵Q的第n行第n列元素的值，n为Q矩阵的行列号，矩阵Q除对角元素外其余元素均为零，σ_r代表噪声的标准方差，R_n是由初始波阻抗模型计算得到的第n个采样点位置处的初始反射系数；ρ为模型约束因子，控制反演结果对初始模型的依赖程度；C为积分算子矩阵，其离散形式表示为：

公式（1）中上角标-1为对矩阵求逆；为第n个采样点处的相对波阻抗值，I₀为反演时窗内第一个采样点对应的初始波阻抗值，I_n为反演时窗内第n个采样点处的初始波阻抗值，ln为自然对数符号，r_i为反演时窗内第i个采样点的反射系数值，表示对从第1个采样点到第n个采样点的r_i进行求和计算；

8）通过步骤7）中第n个采样点处的相对波阻抗的定义，通过推导得到波阻抗与反射系数的关系：

$I_n=I_0{e}^{2\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{r}_i}---\left(2\right)$

其中：I_n为反演时窗内第n个采样点的波阻抗值，I₀为反演时窗内第一个采样点对应的初始波阻抗值，r_i为步骤7）中反演得到的第i个采样点反射系数值，e表示自然对数的底，表示对从第1个采样点到第n个采样点的r_i进行求和计算；

通过公式（2）将步骤7）反演得到的一道的反射系数系列转换为波阻抗序列，得到该道的波阻抗反演结果；

9）对所有的地震道重复步骤7）至8）过程，得到所有道的波阻抗反演结果。

2.根据权利要求1的方法，特点是步骤4）所述的标定为利用测井曲线和地震子波模拟井旁地震记录，实现测井分层到地震层位的标定和映射，由此得到时深关系曲线，由此时深关系可以将深度域的测井曲线转换为时间域曲线。