CN103954992A - 一种反褶积方法及装置 - Google Patents

Abstract

本发明提供了一种反褶积方法及装置。该方法包括以下处理步骤：将采集的地震数据处理成叠后地震数据记录；构建包括反射系数序列在内的目标函数；在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料中反射系数序列的振幅和位置；读取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置；计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。本发明提供的反褶积方法能进一步提高地震资料的分辨率。

Description

一种反褶积方法及装置

技术领域

本方法涉及地震勘探中地震资料处理领域，特别涉及一种反褶积方法及装置。

背景技术

在对地层蕴藏的石油天然气资源进行开发之前，需要对地层进行勘测以测定地下的构造，为油气资源的开采和地层识别提供依据。地震勘探是勘测地层结构的重要手段。

地震勘探是一种利用人工地震技术探测地下构造的勘探方法。它以人工方法按照一定的方式在地表附近激发地震波，产生称之为地震子波的振动信号。来自不同反射深度界面的地震子波以不同的时间到达地表，通过在地表布设一种称之为检波器的接收装置，接收来自不同深度地质界面反射的地震子波，其接收的信号的合集称之为地震记录。地震记录可以认为是地震子波与一系列反射系数序列的褶积。当地层较薄时，薄层上下界面地震子波难以区分，来自不同界面的地震子波叠合在一起，使得从地震记录上很难分辨这种很薄的地质界面。这种利用地震记录分辨薄层的能力称之为地震资料的分辨率。为了利用地震记录分辨薄层，提高地震资料的分辨率，通常需要对地震记录进行反褶积处理。

所述的反褶积处理，通常是指通过压缩地震子波，将延续几十甚至一百多毫秒的地震子波压缩成原来的震源尖脉冲形式，使地震记录变为反射系数序列的尖脉冲组合，提高分辨率的数据处理过程。现有技术中常用的反褶积方法通常是基于两种基本假设（假设反射系数序列是白噪声的，假设地震子波序列是最小相位的），其反褶积过程主要包括：假设地震记录为

$x\left(t\right)=S\left(t\right)+n\left(t\right)=\underset{t=0}{\overset{m}{\mathrm{\Σ}}}b\left(t\right)\mathrm{\ξ}(t-\mathrm{\τ})+n\left(t\right)---(1-1)$

其中S(t)为有效信号，n(t)为干扰波，b(τ)为地震子波，ξ(t)为反射系数序列，τ为反射系数序列的延时，m为反射系数个数。通常会假设不存在干扰波n(t)，即：

x(t)＝S(t)＝b(t)*ξ(t)             （1-2）

对上式（1-2）两边求傅氏变换，则得到频率域的地震记录表示式：

X(ω)＝B(ω)·ξ(ω)            （1-3）

上式（1-3）中，X(ω)、B(ω)和ξ(ω)分别为地震频谱、子波频谱和反射系数的频谱。

上式（1-3）经过变化可以得到：

$\mathrm{\ξ}\left(\mathrm{\ω}\right)\frac{1}{B\left(\mathrm{\ω}\right)}\·X\left(\mathrm{\ω}\right)---(1-4)$

如果令：

$A\left(\mathrm{\ω}\right)=\frac{1}{B\left(\mathrm{\ω}\right)}---(1-5)$

则有：

ξ(ω)＝A(ω)·X(ω)            （1-6）

再对上式（1-6）式做反傅氏变换，将其变换至时间域，可得到：

ξ(t)＝a(t)*x(t)＝a(t)*b(t)*ξ(t)          （1-7）

上式（1-7）中，a(t)为A(ω)的时间函数。根据（1-7）式可知：

a(t)*b(t)＝δ(t)           （1-8）

因为b(t)为地震子波，而a(t)和b(t)之间又存在着频谱互为倒数的关系（即A(ω)＝1/B(ω)），所以把通常a(t)称为反子波，又叫做反褶积因子。

由此可知，若已知地震子波，利用数学方法求出a(t)，再利用（1-7）式让反子波a(t)与地震记录x(t)做褶积，就可以求出反射系数序列ξ(t)，即

$\mathrm{\ξ}\left(t\right)=\underset{\mathrm{\τ}}{\mathrm{\Σ}}\mathrm{\α}\left(\mathrm{\τ}\right)x(t-\mathrm{\τ})---(1-9)$

估计出地震子波的方法很多，常用的有同态滤波方法或统计性子波估算方法。估计出地震子波后，就可以求出反射系数的序列，达到地震子波压缩成尖脉冲、提高地震记录纵向分辨能力的目的。

上述反褶积方法基于两种基本假设估计出地震子波，求出反射系数序列。但是所述的两种基本假设常常与实际地层结构存在较大的差异，尤其是地震子波对弱信号的压制，使得该方法在地震数据处理时降低了地震数据记录的分辨率。

发明内容

本发明目的在于提供一种新的反褶积方法，该方法摆脱了两种基本假设，减少了地震子波对地震数据中弱信号的干涉作用，能提高地震数据记录的分辨率。

本发明的实现方法包括以下处理步骤：

S1：将所需处理的地震数据处理成叠后地震数据记录；

S2：构建包括反射系数序列在内的目标函数；

S3：在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅和位置；

S4：读取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置；

S5：计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。

上述所述的一种反褶积方法，其优选方案为，所述构建目标函数包括：

S201：建立地震数据褶积模型；

S202：将所述地震数据褶积模型转换成矩阵形式；

S203：由所述矩阵形式构建目标函数；所述构建的目标函数为：

J_α＝||Aα-S||²

上式中，J_α为关于反射系数序列α的范数，A为一个N×M的矩阵，其表达式为其中，ω_t为地震子波，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，τ_M第M个反射系数的时延，S为t时刻所述叠后地震数据记录的值。

上述所述的一种反褶积方法，其优选方案为，所述构建的目标函数为：

J＝||Aα-S||²+μ||Cα-ε||²＝J_α+μJ_c

上式中，μ为阻抗权系数，C为左下三角全为1的矩阵，ε为波阻抗。

上述所述的一种反褶积方法，其优选方案为，所述构建的目标函数为：

J＝J_α+μJ_c+βJ_α

上式中，β为正则化的阻尼系数，其中，β＝β₀f₀，β₀为预设值，f₀＝max{(A^TA+μC^TC_ii)}，i为矩阵(A^TA+μC^TC)的下标。

上述所述的一种反褶积方法，其优选方案为，从所述构建的目标函数中计算出所述地震资料中反射系数序列的振幅包括步骤：

S310：从所述叠后地震数据记录中估计出地震子波；

S320：将所述估计出的地震子波代入所述目标函数中，计算出所述反射系数序列的振幅；

其中，所述估计地震子波的方法包括步骤：

S311：建立所述叠后地震数据记录的褶积模型；

S312：将所述建立的褶积模型变换到复赛谱域；

S313：设计低通滤波器，在复赛谱域将所述建立的褶积模型的地震子波和反射系数序列进行分离，得到复赛谱域的地震子波；

S314：将所述复赛谱域的地震子波转换到时间域，得到时间域的地震子波。

上述所述的一种反褶积方法，其优选方案为，所述地震子波的估计方法还包括步骤：

315：对所述估计出的地震子波相位进行校正。

上述所述的一种反褶积方法，其优选方案为，采用模拟退火算法计算反射系数序列中每一个反射系数的时间位置，其操作过程包括以下步骤：

S301：在所述时间窗长度N内随机给定反射系数序列a(j)的初始时间位置，由所述初始时间位置计算所述目标函数J_a当前解，j＝1,2…M；

S302：在所述时间窗长度N内对所述反射系数序列a(j)进行扰动，将所述扰动后的反射系数序列a(j)代入所述目标函数J_a中，计算出所述目标函数J_a的新解；

S303：将所述J_α的新解与当前解进行比较，根据Metropolis准则判断是否接受所述目标函数J_α的新解；

如果判断结果为接受所述新解，则接受反射系数α(j)当前的时间位置；此时，所述新解作为目标函数J_a当前解；

S304：重复S302、S303直到达到收敛条件后，停止迭代，得到所述反射系数α(j)的时间位置。

上述所述的一种反褶积方法，其优选方案为，所述模拟退火算法还包括在所述S303之后：

S3031：设置误差项σ，将所述新解与当前解的差值与所述误差项σ比较；如果所述新解与当前解的差值大于所述误差项σ，则不接受所述目标函数J_α的新解。

上述所述的一种反褶积方法，其优选方案为，在达到收敛条件后，将所述模拟退火算法得出的反射系数序列α(j)的振幅值作为地震资料中反射系数序列的振幅。

本发明还提供一种反褶积装置，该装置包括：数据预处理模块、目标函数建立模块、目标函数计算模块、资料读取模块、反射系数序列输出模块，其中：

所述数据预处理模块，用于将所需处理的地震数据处理成叠后地震数据记录；

所述目标函数建立模块，用于构建包括反射系数序列在内的目标函数；

所述目标函数计算模块，用于在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅和位置；

所述资料读取模块，用于读取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置；

所述反射系数序列输出模块，用于计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。

上述所述的一种反褶积装置，其优选方案为，所述目标函数建立模块构建的目标函数为：

J_α＝||Aα-S||²

上式中，J_α为关于反射系数序列α的范数形式的函数，A为一个N×M的矩阵，其表达式为其中，ω_t为地震子波，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，τ_M第M个反射系数的时延，S为t时刻所述叠后地震数据记录的值。

上述所述的一种反褶积装置，其优选方案为，所述目标函数建立模块构建的目标函数为：

J＝||Aα-S||²+μ||Cα-ε||²＝J_α+μJ_c

上式中，μ为阻抗权系数，C为左下三角全为1的矩阵，ε为波阻抗。

上述所述的一种反褶积装置，其优选方案为，所述目标函数建立模块构建的目标函数为：

J＝J_α+μJ_c+βJ_α

上式中，β为正则化的阻尼系数，其中，β＝β₀f₀，β₀为预设值，f₀＝max{(A^TA+μC^TC_ii)}，i为矩阵(A^TA+μC^TC)的下标。

上述所述的一种反褶积装置，其优选方案为，所述目标函数计算模块包括振幅计算模块、位置计算模块，其中：

所述振幅计算模块，用于从所述目标函数建立模块中计算出所述地震资料中反射系数序列的振幅；

所述位置计算模块，用于从所述目标函数建立模块中计算出所述地震资料中反射系数序列的时间位置。

上述所述的一种反褶积装置，其优选方案为，位置计算模块包括初始模块、扰动模块、迭代模块、收敛模块，其中：

所述初始模块1，用于在所述时间窗长度N内随机给定反射系数序列a(j)的初始时间位置，由所述初始时间位置计算所述目标函数J_a当前解，j＝1,2…M；

在所述扰动模块，用于在所述时间窗长度N内对所述反射系数序列a(j)进行扰动，将所述扰动后的反射系数序列a(j)代入所述目标函数J_a中，计算出所述目标函数J_a的新解；

所述迭代模块，用于将所述J_α的新解与当前解进行比较，根据Metropolis准则判断是否接受所述目标函数J_α的新解；

如果判断结果为接受所述新解，则接受反射系数α(j)当前的时间位置；此时，所述新解作为目标函数J_a当前解；

所述收敛模块，用于继续计算所述新解，并与所述当前解进行比较直到达到收敛条件后，停止迭代，得到所述反射系数α(j)的时间位置。

本发明提供的反褶积方法是一种稀疏脉冲非线性反褶积方法，该方法将传统的以子波压缩为核心理念的反褶积方法转移到求解反射系数位置和振幅上来，减小了地震子波对地震数据处理的影响，较大幅度地提高地震记录的分辨率。同时通过对反射系数统计特征的有效约束，减小了反褶积结果的多解性，提高了该反褶积方法数据处理的结果稳定性。

附图说明

图1是本发明提供的一种反褶积方法实施例1的方法流程图；

图2是本发明实施例1中模拟退火算法的方法流程图；

图3是本发明实施例1中模拟退火算法的另一种实施方式的方法流程图；

图4是本发明提供的一种反褶积装置模块结构示意图；

图5是本发明反褶积装置中目标函数计算模块的模块结构示意图；

图6是本发明所述目标函数中位置计算模块的模块结构示意图；

图7是采用基于两种基本假设对数据进行反褶积处理后的效果图；

图8是采用本发明反褶积方法对同一地震数据反褶积处理后的效果图。

具体实施方式

为了使本技术领域的人员更好地理解本申请中的技术方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都应当属于本发明保护的范围。

图1是本发明提供的一种反褶积方法实施例1的方法流程图。如图1所示，该方法包括以下处理步骤：

S1：采集地震数据，将所需处理的地震数据处理成叠后地震数据记录。

地震数据的采集，其主要是指在油气勘探的区域，布置二维或三维测线，使用炸药震源或可控震源激发地震波，在测线上等间距布置多个检波器来接收地震波信号，以等时间间隔离散采样地震数据，并以数字形式收集记录。在地震数据采集记录中，每个检波点上记录的地震波数据称为地震道。一个检波点上的地震道为一个单道或一个道。可以为一个地震道设置一个地震道序号。采集后的地震数据通常会包括坏道、空道、随机干扰及其他误差，因此还需要进行预处理、校正等消除或减小地震数据记录的误差。

所述的预处理通常包括对采集的原始数据进行解码、抽取地震道、剔除空道或坏道等一系列的前期处理过程。在校正过程中，可以包括静校正和动校正处理。所述的静校正可以是按静校正量所确定的移动量，把地震道的振幅离散值进行整体移动，消除地面山谷、地表风化层和低速带等地表异常的影响。所述的动校正可以是通过动校正公式计算出动校正量来消除地震波到达各不同检波点的正常时差。可以在经动、静校正处理后的地震数据记录上地震道序号相同的地震道中，取其相对应采样点的算术平均值作为所述相对应采样点叠加后的新值，这样由多个地震数据叠加组成的地震道的集合称为叠后地震数据记录。地震数据叠加的方法有很多，本发明方法地震数据叠加的目的是压制多次波和随机干扰，因此不限于本实施例中所述的地震数据叠加方法。

在地震数据处理中，对地震数据记录的反褶积可以是对整个地震数据记录反褶积处理，也可以仅对某一时间窗内的一段地震数据记录反褶积处理。本申请中所述的所需处理的地震资料，可以为上述所述的两种情况中任意一种。

采集地震数据后，将所需处理的地震数据处理成叠后地震数据记录。

S2：构建包括反射系数序列在内的目标函数。

所述构建的目标函数可以包括反射系数序列，也可以同时包括地震子波和反射系数序列。所构建的目标函数可以为多次试验或理论推导得出，并且能够求解出至少一组可接受的反射系数序列的解。本实施例1中所述的构建目标函数可以包括以下步骤：

S201：建立地震数据褶积模型。

将地震数据记录的褶积模型表示为s_t＝w_t*r_t+n_t，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度。上式中，s_t为地震数据记录，ω_t为地震子波，r_t为稀疏脉冲（即反射系数序列），n_t为随机噪音。其中，所述稀疏脉冲r_t的表达式可以为如下的形式：

$r_t=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_j{\mathrm{\δ}}_{t-{\mathrm{\τ}}_j}$

上式中，M是非零脉冲的个数，δ_t是单位脉冲函数，τ_j是脉冲的时延，α_j是脉冲的振幅，j＝1,2…M。

将所述系数脉冲r_t表达式代入所述褶积模型中，可以得到：

$s_t=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_j{\mathrm{\ω}}_{t-{\mathrm{\τ}}_j+1}+n_t,t=\mathrm{1,2}...N---(2-1)$

S202：将所述地震数据褶积模型转换成矩阵形式。

上述公式（2-1）写成矩阵形式进行展开，可以表示为：

$\left[\begin{array}{cccc}{w}_{1-{\mathrm{\τ}}_1+1}& w_{1-{\mathrm{\τ}}_2+1}& \·\·\·& w_{1-{\mathrm{\τ}}_M+1}\\ w_{2-{\mathrm{\τ}}_1+1}& w_{2-{\mathrm{\τ}}_2+1}& \·\·\·& w_{2-{\mathrm{\τ}}_M+1}\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ \·& \·& \·& \·\\ w_{N-{\mathrm{\τ}}_1+1}& w_{N-{\mathrm{\τ}}_2+1}& \·\·\·& w_{N-{\mathrm{\τ}}_m+1}\end{array}\right]\left[\begin{array}{c}{\mathrm{\α}}_1\\ {\mathrm{\α}}_2\\ \·\\ \·\\ \·\\ {\mathrm{\α}}_M\end{array}\right]+\left[\begin{array}{c}{n}_1\\ n_2\\ \·\\ \·\\ \·\\ n_N\end{array}\right]=\left[\begin{array}{c}{S}_1\\ S_2\\ \·\\ \·\\ \·\\ S_n\end{array}\right]$

上述矩阵表达式中，至可以表达为一个N×M矩阵A，因此上述矩阵表达式也可以表示为：

Aα+n=s            （2-2）

上式（2-2）中，A是一个N×M的矩阵，其中

S203：由所述矩阵形式构建目标函数。

上述公式（2-2）中的A是一个N维的非线性方程组，表示为在地震数据记录中t时刻地震子波的值。褶积模型合成的地震数据记录是由地震子波和反射系数序列褶积而成，若已知实际地震数据记录和地震子波，可以通过所述叠后地震记录与合成的地震记录的最小平方差来求解反射系数序列的振幅。基于上述思想，本实施例1中可以将所述目标函数定义如下：

J_α＝||Aα-S||²        （2-3）

上式（2-3）定义的目标函数J_α是一个关于反射系数序列α的范数。A为一个N×M的矩阵，其表达式为其中，ω_t为地震子波，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，τ_M第M个反射系数的时延，S为t时刻所述叠后地震数据记录的值。

S3：在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅和位置。

从所构建的目标函数中求解出反射系数序列的振幅和位置，可以根据所构建的目标函数不同而有不同计算方法，具体的可以根据所构建的目标函数的特性进行计算。在本实施例1中，所构建的目标函数J_α＝||Aα-S||²中，其计算反射系数序列振幅的过程可以包括以下步骤：

S310：从所述叠后地震数据记录中估计出地震子波。

在本实施例1构建的目标函数中，涉及到未知参数量地震子波。因此，在求解目标函数时可以先由上述S1所述的叠后地震数据记录估计出地震子波ω_t。地震数据子波的估计方法有很多，包括同态滤波方法和统计性子波估计方法。不同的估计方法得出地震子波会有偏差。在常规的反褶积处理中，所述地震子波ω_t的偏差对最后反褶积处理的结果影响很大。而在本发明的反褶积方法中，即使所述估计出的地震子波ω_t存在一定的误差，也可以得出较为准确反射系数序列。

本实施例1中所述地震子波的估计方法可以包括以下步骤：

S311：建立所述叠后地震数据记录的褶积模型。

所述叠后地震数据记录的褶积模型可以表示为：

s_t＝ω_t*r_t+n_t

其中，s_t是已知的叠后地震数据记录，ω_t为地震子波，r_t为反射系数序列，n_t为噪声项。

S312：将所述建立的褶积模型变换到复赛谱域。

不考虑噪声项n_t。对上述褶积模型等号两边同时做傅里叶变换，得到所述叠后地震数据记录褶积模型在频率域的表达式：

x(ω)＝w(ω)ξ(ω)          （3-1）

对上述公式（2-1）x(ω)＝w(ω)ξ(ω)等号两边取对数，转化为线性系统的公式：

lnx(ω)＝lnw(ω)+lnξ(ω)              （3-2）

对上述公式（2-2）lnx(ω)＝lnw(ω)+lnξ(ω)作反傅里叶变换，得到：

$\tilde{x}\left(t\right)=\tilde{w}\left(t\right)+\widetilde{\mathrm{\ζ}}\left(t\right)---(3-3)$

其中，分别称为x(t)，w(t)，ξ(t)的复赛谱序列。

S313：设计低通滤波器，在复赛谱域将所述建立的褶积模型的地震子波和反射系数序列进行分离，得到复赛谱域的地震子波。

将地震记录由时间域转换为复赛谱域后，地震记录可以表示为地震子波与反射系数序列之和。其中地震子波能量集中在原点附近。此时，可以通过设定一个“分割点”，对复赛谱域的地震记录进行低通滤波，可以得到复赛谱域的地震子波。

S314：将所述复赛谱域的地震子波转换到时间域，得到时间域的地震子波ω_t。

在估计出地震子波后，为了得到最佳的地震子波相位，还可以包括以下处理步骤：

S315：对所述估计出的地震子波相位进行校正。

对所述地震子波相位进行校正的方法有很多，常用的相位校正方法有常相位校正方法。所述的常相位校正方法通常是指基于地震子波的剩余相位不依赖频率变化的假设，对记录的相位谱进行一个不依赖于频率的常数的调整。本实施例1中所述的对地震子波相位进行校正不限于所述的常相位校正方法。

估计出一个地震子波后，将所述地震子波ω_t作为已知量求解所述目标函数。

S320：将所述估计出的地震子波代入所述目标函数中，计算出所述反射系数序列的振幅。

将上述所述建立的目标函数（2-3）展开为：

J_α＝（Aα-S）^T(Aα-S）

＝(α^TA^T-S^T)(Aα-S）           （3-4）

＝α^TA^TAα-α^TA^TS-S^TAα+S^TS

对上述公式（3-4）进行求导，可以得到：

$\frac{{\∂J}_{\mathrm{\α}}}{{\∂\mathrm{\α}}^T}=A^T\mathrm{A\α}-A^{T}S=0---(3-5)$

由上述公式（3-5）可以得到求解所述反射系数序列α的表达式：

α＝(A^TA)^-1A^TS           （3-6）

上式（3-6）中，A为一个N×M的矩阵，A^T为A的转置，其表达式为其中，ω_t为地震子波，t=1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，S为t时刻所述叠后地震数据记录的值。

所述M值是一个经验值，该值可以由作业人员根据采集的地震数据记录中单位时间内可能含有的反射系数序列个数自行设定。多次试验证明，所述M值的设置可以略大于正常的经验值，尽量减少遗漏的反射系数序列。例如，一个地震道的一段时间窗N内设置的M为20，实际该时间窗长度N内有15个反射系数序列，那么在本实施例1求解过程中可以求解出所述15个反射系数序列，其余的5个反射系数序列为0或近似为0。相反，如果所述时间窗N内实际有20个反射系数序列，而在实施本实施例1过程中设置为15，那么在求解过程中将会仅得出15个反射系数序列，其余5个反射系数序列将会被遗漏。

上述公式（3-6）计算得出的是一个地震道中所选取的地震资料中一组M个反射系数序列的振幅，因此还需要计算出所述M个反射系数序列中每个反射系数的时间位置。计算所述反射系数序列位置的方法可以根据所构建的目标函数选择合适的方法。针对本实施例1所述构建的目标函数，本实施例1中可以采用模拟退火算法计算所述反射系数序列中每一个反射系数的时间位置。模拟退火（Simulated Annealing，简称SA）是一种通用概率算法，用来在一个大的搜寻空间（解空间）内采用迭代算法找寻命题的最优解。

图2是本发明实施例1中模拟退火算法的方法流程图。如图2所示，所述采用模拟退火算法计算所述每一个反射系数的时间位置的过程可以包括以下步骤：

S301：在所述时间窗长度N内随机给定反射系数序列a(j)的初始时间位置，由初始时间位置计算所述目标函数J_a当前解，j＝1…M。

M个反射系数序列中的每个反射系数都在所述时间窗长度N内随机给定初始时间位置。

S302：在所述时间窗长度N内对所述反射系数序列a(j)进行扰动，将所述扰动后的反射系数序列a(j)代入所述目标函数J_a中，计算出所述目标函数J_a的新解。

每次扰动后更新所述反射系数序列a(j)的时间位置。

所述的对反射系数序列a(j)进行扰动，在模拟退火算法迭代中通常是指对反射系数序列a(j)的时间位置给定一个扰动量，使其时间位置发生变化，进而使与所述反射系数序列a(j)有关的目标函数J_α的值发生变化。所述的扰动量可以是由其他函数随机产生设定的，可以是一个增加的量，也可以是一个减小的量。当然，所述的扰动量也可以是自行设定的具有一定变化规律的扰动量。本实施例1中所述的扰动量可以使所述该模拟退火算法在所述时间窗长度N内对反射系数序列的时间位置搜索一遍即可，本发明方法对具体的扰动量的定义不做限定。

所述目标函数J_α的新解是相对于上一次扰动或初始值的解。所述目标函数J_α的上一次扰动的解或初始值的解可以统称为当前解。

S303：将所述J_α的新解与当前解进行比较，根据Metropolis准则判断是否接受所述目标函数J_α的新解。

若由Metropolis准则判断接受所述目标函数J_α的新解，则将所述目标函数J_α的新解作为当前解继续迭代计算，同时也接受所述反射系数α(j)当前的时间位置。此时，将所述新解作为目标函数J_a当前解。

S304：重复S302、S303直到达到收敛条件后，停止迭代，得到所述反射系数α(j)的时间位置。

所述的收敛条件，可以是所述模拟退火中定义的“退火计划”，所述的“退火计划”通常是指模拟退火中稳定的衰减函数，该函数是一个收敛函数。所述的收敛条件也可以是自定义的迭代次数或搜索解结束的停止条件。本实施例1中所述的收敛条件可以为所述模拟退火算法对当前所求的反射系数在所述时间窗长度N内的完成一遍时间位置搜索。

停止迭代后，可以得到的所述反射系数α(j)的最佳时间位置。

图3是所述模拟退火算法的另一种实施方式。如图3所示，上述所述的模拟退火中，为消除异常坏值，所述S303中可以包括以下步骤：

S3031：设置误差项σ，将所述新解与当前解的差值与所述误差项σ比较；如果所述新解与当前解的差值大于所述误差项σ，则不接受所述目标函数J_α的新解。

所述误差项σ为一个经验值，其目的是防止求解结果与原所述叠后地震数据记录误差太大。该值可以根据计算需求自行设定。

需要说明的是，在上述所述采用模拟退火计算得出所述反射系数序列最佳时间位置的同时，也计算得出了所述目标J_α的新解，即计算得出了所述反射系数序列α(j)的振幅。试验证明，该模拟退火算法计算得出的所述反射系数序列α(j)最佳时间位置时其所述反射系数序列α(j)的振幅也为最优的解。因此，本发明方法实施例1中可以在达到收敛条件后，将所述模拟退火算法得出的所述反射系数序列α(j)的振幅值作为所述反射系数序列的振幅。

S4：读取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置；

计算完所述选取的地震资料反射系数序列的振幅和位置后，继续选取下一道地震资料计算其反射系数序列的振幅和位置。所述下一道地震资料的选取，可以按照当前所选取的地震道序列的顺序依次选取，也可以按照其他可以选取完所有所需处理的地震道的方法。

S5：计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。

计算完所需处理的地震数据中所有地震道的反射系数序列的振幅和位置，即完成了所述叠后地震数据记录的反褶积处理。

本实施例1所述的反褶积方法将传统的以子波压缩或两种基本假设的方法转移到求解反射系数序列的位置和振幅上，提供了一种新的反褶积思路。在本发明方法中可以采用模拟退火算法进行多次迭代计算出最佳的所述反射系数序列的振幅和位置，在此过程中可以将一些由人为制造的假象所破坏的地震数据横向连续性降到最低，很大程度上提高反褶积处理结果的准确性。

本发明还提供所述反褶积方法的第二种实施方式。本发明实施例2提供的反褶积方法允许在最优化过程中可引入波阻抗约束。所述波阻抗具有阻力的含义，在反褶积过程中将波阻抗的因素考虑在内，可以提高地震数据记录反褶积的准确性和横向连续性。地震波在介质中传播时，作用于某个面积上的压力与单位时间内垂直通过此面积的质点流量之比，称为波阻抗。

在地震数据记录处理中，反射系数序列α_j与波阻抗ε_i有如下的关系：

$\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\α}}_j{u}_{T_i-{\mathrm{\τ}}_j}\text{=}{\mathrm{\ϵ}}_i---(4-1)$

将公式（4-1）经整理和变换后，写成矩阵形式为：

$\left[\begin{array}{cccc}1& 0& \·\·\·& 0\\ 1& 1& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 1& 1& \·\·\·& 1\end{array}\right]\left[\begin{array}{c}{r}_1\\ r_2\\ \·\·\·\\ r_M\end{array}\right]=\left[\begin{array}{c}{\mathrm{\ϵ}}_0\\ {\mathrm{\ϵ}}_1\\ \·\·\·\\ {\mathrm{\ϵ}}_{k-1}\end{array}\right]$

上述矩阵形式可以表示为：Cα＝ε，其中，C为左下三角全为1的矩阵。

取上述反射系数序列与波阻抗差的范数形式引入所述构建的目标函数中，则所述目标函数可以表示为：

J＝||Aα-S||²+μ||Cα-ε||²＝J_α+μJ_c          （4-2）

公式（4-2）中μ为阻抗权系数，是一个可以预先设置经验值。将公式（4-2）展开为：

J＝(Aα-S)^T(Aα-S)+μ(Cα-ε)^T(Cα-ε)

＝(α^TA^T-S^T)(Aα-S)+μ(α^TC^T-ε^T)(Cα-ε)

＝α^TA^TAα-α^TA^TS-S^TAα+S^TS+μ（α^TC^TCα-α^TC^Tε-ε^TCα+ε^Tε）

对公式（4-2）求导可得：即A^TAα-A^TS+μC^TCα-μC^Tε＝0，

由上可以得出所述反射系数序列α的表达式可以为：

α＝（A^TA+μC^TC）^-1（A^TS+μC^Tε）         （4-3）

上式（4-3）中，A为一个N×M的矩阵，其表达式为A^T为A的转置，ω_t为地震子波，t=1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，S为t时刻所述叠后地震数据记录的值，μ为阻抗权系数，C为左下三角全为1的矩阵，ε为波阻抗。

上述实施例2提供的反褶积方法中引入波阻抗因素，可以提高该反褶积方法处理的准确性。

在上述所述求解目标函数J_α＝||Aα-S||²的过程中，还可以引入正则化操作。所述的正则化是为了防止矩阵求解过程中出现“病态”，而加入的辅助性因子，可以增加方程求解的稳定性。在本实施例中1，在所述目标函数的正则化可以为加入因子βJ_α，则，所述构建的目标函数可以表达为：

J＝J_α+μJ_c+βJ_α         （4-4）

上式（3-7）中，β为正则化的阻尼系数，其中β＝β₀f₀。其中β₀可以根据经验值进行预设定，通常设置在0.001左右。本实施例1中所述的β₀可以设置为0.001。可以f₀可以定义为：

f₀＝max{(A^TA+μC^TC)_ii}           （4-5）

上式中，i为矩阵(A^TA+μC^TC)的下标。上式f₀的含义可以表示为取矩阵(A^TA+μC^TC)对角线中的最大值。

所述目标函数加入正则化处理，很大程度上减少了解的非唯一性的可能性，保障了解的稳定性，提高了本申请发明方法对地震数据记录进行反褶积的处理结果的准确性和稳定性。

本发明根据所述反褶积方法提供一种反褶积装置。图4是本发明提供的一种反褶积装置模块结构示意图。如图4所述，该装置包括：数据预处理模块1、目标函数建立模块2、目标函数计算模块3、资料读取模块4、反射系数序列输出模块5，其中：

所述数据预处理模块1，可以用于将采集的所需处理的地震数据处理成叠后地震数据记录。

所述目标函数建立模块2，可以用于构建包括反射系数序列在内的目标函数。所述构建的目标函数可以为：

J_α＝||Aα-S||²

上式中，J_α为关于反射系数序列α的范数形式的函数，A为一个N×M的矩阵，其表达式为其中，ω_t为地震子波，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，τ_M第M个反射系数的时延，S为t时刻所述叠后地震数据记录的值。

所述目标函数建立模块2构建的目标函数也可以为：

J＝||Aα-S||²+μ||Cα-ε||²＝J_α+μJ_c

上式中，μ为阻抗权系数，C为左下三角全为1的矩阵，ε为波阻抗。

所述目标函数建立模块2构建的目标函数还可以为：

J＝J_α+μJ_c+βJ_α

上式中，β为正则化的阻尼系数，其中，β＝β₀f₀，β₀为预设值，f₀＝max{(A^TA+μC^TC_ii)}，i为矩阵(A^TA+μC^TC)的下标。

所述目标函数计算模块3，可以用于在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅和位置。

所述资料读取模块4，可以用于度取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置。

所述反射系数序列输出模块5，可以用于计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。

图5是所述目标函数计算模块3的模块结构示意图。如图5所示，所述目标函数计算模块3，可以包括振幅计算模块310、位置计算模块300，其中：

所述振幅计算模块310，可以用于从所述目标函数建立模块2中计算出所述地震资料中反射系数序列的振幅。

所述位置计算模块300，可以用于从所述目标函数建立模块2中计算出所述地震资料中反射系数序列的时间位置。

图6是所述位置计算模块300的模块结构示意图。如图6所示，所述位置计算模块300可以包括初始模块301、扰动模块302、迭代模块303、收敛模块304，

其中，所述初始模块301，可以用于在所述时间窗长度N内随机给定反射系数序列a(j)的初始时间位置，由所述初始时间位置计算所述目标函数J_a当前解，j＝1,2…M。

所述扰动模块302，可以用于在所述时间窗长度N内对所述反射系数序列a(j)进行扰动，将所述扰动后的反射系数序列a(j)代入所述目标函数J_a中，计算出所述目标函数J_a的新解。

所述迭代模块303，可以用于将所述J_α的新解与当前解进行比较，根据Metropolis准则判断是否接受所述目标函数J_α的新解。如果判断结果为接受所述新解，则接受反射系数α(j)当前的时间位置。此时，所述新解作为目标函数J_a当前解。

所述收敛模块304，可以用于继续计算所述新解，并与所述当前解进行比较直到达到收敛条件后，停止迭代，得到所述反射系数α(j)的时间位置。

本发明提供的一种反褶积装置，采用了最小二乘法和正则化方法能准确预测地震数据记录中反射系数序列的振幅，并且用模拟退火算法计算反射系数序列振幅的位置，提高地震数据处理后的分辨率。

试验证明，本发明方法和装置相比常规的基于两种基本假设的方法，将传统的以子波压缩为核心理念的反褶积方法转移到求解反射系数位置和振幅上来，可以较大幅度地提高地震记录的分辨率。图7采用基于两种基本假设对数据进行反褶积处理后的效果图，图8是采用本发明反褶积方法对同一地震数据反褶积处理后的效果。由图7和图8对比可以看出，采用本发明反褶积方法处理之后的地震数据分辨率得到明显改善，沉积内幕和接触关系清晰自然。尤为难得的是，由于消除了反褶积之后剩余相位对弱反射信号的干涉影响，1700-1900ms之间的波阻抗差异较小的弱反射信号得到了有效恢复，为精细地质解释和储层描述提供更加可靠的基础数据。

Claims (15)

1.一种反褶积方法，其特征在于，包括以下处理步骤：

S1：将所需处理的地震数据处理成叠后地震数据记录；

S2：构建包括反射系数序列在内的目标函数；

S3：在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅和位置；

S4：读取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置；

S5：计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。

2.如权利要求1所述的一种反褶积方法，其特征在于，所述构建目标函数包括：

S201：建立地震数据褶积模型；

S202：将所述地震数据褶积模型转换成矩阵形式；

S203：由所述矩阵形式构建目标函数；所述构建的目标函数为：

J_α＝||Aα-S||²

上式中，J_α为关于反射系数序列α的范数，A为一个N×M的矩阵，其表达式为其中，ω_t为地震子波，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，τ_M第M个反射系数的时延，s为t时刻所述叠后地震数据记录的值。

3.如权利要求2所述的一种反褶积方法，其特征在于，所述构建的目标函数为：

J＝||Aα-S||²+μ||Cα-ε||²＝J_α+μJ_c

上式中，μ为阻抗权系数，C为左下三角全为1的矩阵，ε为波阻抗。

4.如权利要求3所述的一种反褶积方法，其特征在于，所述构建的目标函数为：

J＝J_α+μJ_c+βJ_α

上式中，β为正则化的阻尼系数，其中，β＝β₀f₀，β₀为预设值，f₀＝max{(A^TA+μC^TC_ii)}，i为矩阵(A^TA+μC^TC)的下标。

5.如权利要求1-4任意一项所述的一种反褶积方法，其特征在于，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅包括步骤：

S310：从所述叠后地震数据记录中估计出地震子波；

S320：将所述估计出的地震子波代入所述目标函数中，计算出所述反射系数序列的振幅；

其中，所述估计地震子波的方法包括步骤：

S311：建立所述叠后地震数据记录的褶积模型；

S312：将所述建立的褶积模型变换到复赛谱域；

S313：设计低通滤波器，在复赛谱域将所述建立的褶积模型的地震子波和反射系数序列进行分离，得到复赛谱域的地震子波；

S314：将所述复赛谱域的地震子波转换到时间域，得到时间域的地震子波。

6.如权利要求5所述的一种反褶积方法，其特征在于，所述地震子波的估计方法还包括步骤：

315：对所述估计出的地震子波相位进行校正。

7.如权利要求5所述的一种反褶积方法，其特征在于，采用模拟退火算法计算反射系数序列中每一个反射系数的时间位置，其操作过程包括以下步骤：

S301：在所述时间窗长度N内随机给定反射系数序列a(j)的初始时间位置，由所述初始时间位置计算所述目标函数J_a当前解，j＝1,2…M；

S302：在所述时间窗长度N内对所述反射系数序列a(j)进行扰动，将所述扰动后的反射系数序列a(j)代入所述目标函数J_a中，计算出所述目标函数J_a的新解；

S303：将所述J_α的新解与当前解进行比较，根据Metropolis准则判断是否接受所述目标函数J_α的新解；

如果判断结果为接受所述新解，则接受反射系数α(j)当前的时间位置；此时，所述新解作为目标函数J_a当前解；

S304：重复S302、S303直到达到收敛条件后，停止迭代，得到所述反射系数α(j)的时间位置。

8.如权利要求7所述的一种反褶积方法，其特征在于，所述模拟退火算法还包括在所述S303之后：

S3031：设置误差项σ，将所述新解与当前解的差值与所述误差项σ比较；如果所述新解与当前解的差值大于所述误差项σ，则不接受所述目标函数J_α的新解。

9.如权利要求7所述的一种反褶积方法，其特征在于，在达到收敛条件后，将所述模拟退火算法得出的反射系数序列α(j)的振幅值作为地震资料中反射系数序列的振幅。

10.一种反褶积装置，其特征在于，该装置包括：数据预处理模块、目标函数建立模块、目标函数计算模块、资料读取模块、反射系数序列输出模块，其中：

所述数据预处理模块，用于将所需处理的地震数据处理成叠后地震数据记录；

所述目标函数建立模块，用于构建包括反射系数序列在内的目标函数；

所述目标函数计算模块，用于在所述叠后地震数据记录中选取一道所需处理的地震资料，从所述构建的目标函数中计算出所述地震资料的反射系数序列的振幅和位置；

所述资料读取模块，用于读取下一道所需处理的地震资料，计算出所述下一道地震资料反射系数序列的振幅和位置；

所述反射系数序列输出模块，用于计算出所述叠后地震数据记录中所有地震道反射系数序列的振幅和位置，输出反射系数序列，完成所述叠后地震数据记录的反褶积处理。

11.如权利要求10所述的一种反褶积装置，其特征在于，所述目标函数建立模块构建的目标函数为：

J_α＝||Aα-S||²

上式中，J_α为关于反射系数序列α的范数形式的函数，A为一个N×M的矩阵，其表达式为其中，ω_t为地震子波，t＝1，2，…N，N为所述叠后地震数据记录的时间窗长度，M为预先设置的所述时间窗长度N内反射系数序列的个数，τ_M第M个反射系数的时延，S为t时刻所述叠后地震数据记录的值。

12.如权利要求11所述的一种反褶积装置，其特征在于，所述目标函数建立模块构建的目标函数为：

J＝||Aα-S||²+μ||Cα-ε||²＝J_α+μJ_c

上式中，μ为阻抗权系数，C为左下三角全为1的矩阵，ε为波阻抗。

13.如权利要求12所述的一种反褶积装置，其特征在于，所述目标函数建立模块构建的目标函数为：

J＝J_α+μJ_c+βJ_α

上式中，β为正则化的阻尼系数，其中，β＝β₀f₀，β₀为预设值，f₀＝max{(A^TA+μC^TC_ii)}，i为矩阵(A^TA+μC^TC)的下标。

14.如权利要求11所述的一种反褶积装置，其特征在于，所述目标函数计算模块包括振幅计算模块、位置计算模块，其中：

所述振幅计算模块，用于从所述目标函数建立模块中计算出所述地震资料中反射系数序列的振幅；

所述位置计算模块，用于从所述目标函数建立模块中计算出所述地震资料中反射系数序列的时间位置。

15.如权利要求14所述的一种反褶积装置，其特征在于，位置计算模块包括初始模块、扰动模块、迭代模块、收敛模块，其中：

所述初始模块1，用于在所述时间窗长度N内随机给定反射系数序列a(j)的初始时间位置，由所述初始时间位置计算所述目标函数J_a当前解，j＝1,2…M；

在所述扰动模块，用于在所述时间窗长度N内对所述反射系数序列a(j)进行扰动，将所述扰动后的反射系数序列a(j)代入所述目标函数J_a中，计算出所述目标函数J_a的新解；

所述迭代模块，用于将所述J_α的新解与当前解进行比较，根据Metropolis准则判断是否接受所述目标函数J_α的新解；

如果判断结果为接受所述新解，则接受反射系数α(j)当前的时间位置；此时，所述新解作为目标函数J_a当前解；

所述收敛模块，用于继续计算所述新解，并与所述当前解进行比较直到达到收敛条件后，停止迭代，得到所述反射系数α(j)的时间位置。