CN102707314A - 一种多路径双谱域混合相位子波反褶积方法 - Google Patents

Abstract

公开了一种多路径双谱域混合相位子波反褶积方法。所述方法包括如下步骤：从采集到的反射地震记录中选取n道地震记录；根据选取的n道地震记录确定子波的振幅谱；根据选取的n道地震记录确定子波的相位谱；根据确定的子波的振幅谱和相位谱，进行反傅里叶变换来获得时域的混合相位子波；根据用户基于确定的子波的振幅谱输入的期望的输出子波的振幅谱，确定希望输出的零相位子波；以混合相位子波为输入以及零相位子波为输出，进行匹配滤波来确定反褶积滤波器；基于反褶积滤波器和地震记录，获得分辨率提高的地震记录。

Description

一种多路径双谱域混合相位子波反褶积方法

技术领域

本发明涉及本发明属于石油勘探地震资料处理领域，尤其涉及一种能够提高地震资料分辨率的多路径双谱域混合相位子波反褶积方法。

背景技术

地震勘探是一种利用人工地震技术探测地下构造和岩性的勘探方法。它按照一定的方式在地表附近激发地震波，产生称之为地震子波的振动信号，地震子波由震源开始向地下传播，遇到地质界面之后，在界面处发生透射和反射，透射的地震子波继续向下传播，而反射子波在界面的反射点位置向上传播，其强度和极性取决于界面反射系数的大小和符号。来自不同反深度界面的反射子波以不同的时间到达地表，通过在地表布设一种称之为检波器的接收装置，接收来自不同深度地质界面的反射波，接收到的数字信号称之为地震记录，因此，地震记录是不同深度地质界面反射子波的叠加。

地下的每一个地质界面都对应于地震记录中的一个反射子波，通过研究地震记录中不同时间的地震子波，可以研究地质界面的空间展布情况，进而预测和分析地下构造特征。但是，由于地震子波是具有一定延续时间的振动信号，当地层较薄时，来自不同界面的地震子波叠合在一起，使得从地震记录上无法分辨地质界面，这种利用地震信号分辨薄层的能力称之为地震信号的分辨率。

目前，基于地震记录褶积模型的脉冲(预测)反褶积技术是工业界提高地震勘探分辨率的主要方法。忽略噪声的影响，地震记录x(t)、地震子波w(t)和界面反射系数r(t)的关系可用褶积模型表示为：

x(t)＝w(t)*r(t)    (2-1)

其中，“*”表示褶积。

反褶积的目的就是从地震记录x(t)中恢复界面反射系数r(t)。但是，由于地震记录的子波是未知的，在方程(2-1)中同时存在两个未知量w(t)和r(t)，因此直接利用方程(2-1)无法从地震记录中恢复地下界面的反射系数r(t)。

为了能够利用方程(2-1)从地震记录x(t)中恢复地下界面的反射系数r(t)，脉冲(预测)反褶积技术对方程(2-1)进行了两个假设，首先假设地震子波w(t)是最小相位子波，其次假设反射系数r(t)具有高斯白噪的特点。在这两个假设条件下，脉冲(预测)反褶积技术利用方程(2-1)得到反射系数r(t)在最小二乘意义下的解，并以此来近似实际地下反射系数。

勘探实践证明：脉冲(预测)反褶积技术对反射系数所做的高斯白噪假设是基本合理的，而实际地震勘探的子波绝大多数为混合相位子波，因此脉冲(预测)反褶积技术对地震子波所做的最小相位假设与地震勘探的实践相矛盾，由此降低了脉冲(预测)反褶积技术提高地震勘探分辨率的能力和保真性。

如果能够从地震记录中估算出混合相位子波w(t)，那么方程(2-1)中只有反射系数r(t)是未知的，此时就能够利用匹配滤波技术由方程(2-1)计算界面的反射系数r(t)，从而达到提高地震勘探分辨率的目的。因此，从地震记录中估计混合相位子波成为利用反褶积技术提高地震勘探分辨率的关键。

利用同态反褶积方法，在复赛谱域做低通滤波来提取子波，可以避免对子波相位的要求，得到混合相位子波。但是，相位展开的困难、子波相位对噪声的敏感、反射系数的相位和子波的相位在复赛谱中的重叠这三方面原因使得该方法在子波估计时存在较大的误差。

Porsani和Ursin(1998，2000)给出了一种混合相位子波估计方法，该方法的基本做法是：假设子波的Z变换在单位圆上无零点，在单位圆内有k个零点，通过求解扩展的Yule-Walker方程首先获得最小延迟分量的反滤波器，再通过解两次YW方程来获得混合相位子波的最大延迟分量。由最大、最小分量确定出混合相位子波，然后确定混合相位子波的反滤波器，将反滤波器应用到地震记录上并计算输出的L_p模(L_p＞2)，由L_p模最大来确定出最佳的最大、最小相位组合，进而得到最终的混合相位子波。由于该方法要求解三个EYW和YW方程，为增强稳定性，每一步都要加入一定的白噪，在达到稳定性的同时也降低了子波求取的精度。另外，以L_p模最大做为确定子波相位的准则，也缺乏地震勘探理论基础和实践基础的支持。

一般来说，地震子波的振幅谱相对容易估计准确，但相位谱很难估计准确。尽管基于反射系数振幅谱，在一些情况下可以分辨地层的细节，但是子波的相位信息估计不准对反射系数或波阻抗反演结果会有一定影响。因此，估计地震子波的核心是估计子波的相位谱。20世纪80年代后期以来，许多学者开始应用高阶统计方法来估计子波，大多是基于数学领域发展起来的高阶累积量和高阶谱理论。不论是从时间域还是从频率域，为基于高阶统计量的混合相位子波估算提供了一条新的思路。如Lazear(1993)率先将高阶统计理论用于子波估计，相比二阶统计量，高阶统计量在能够保留更多的地震道信息，例如相位信息等，更符合现在日益复杂的勘探目标的需求。

尽管脉冲(预测)反褶积技术存在着明显缺陷，且混合相位反褶积方法的研究取得了一些成果，但是目前还没有一种成熟、实用的混合相位反褶积方法在地震勘探实践中取得实质性应用，脉冲(预测)反褶积技术仍然是提高地震勘探分辨率的主要技术，其局限性对地震勘探分辨率的影响依然存在。

20世纪90年代后期，出现了非线性反褶积方法。例如：最小熵反褶积、L1模的反褶积、稀疏反褶积等。如基于L1模的反褶积方法主要利用反演前后波形的最佳匹配(线性问题)和反演出的反射系数的L1模最小(非线性问题)来相互约束实现的。本质上，该方法是通过这两项的加权之和最小来反演反射系数的。同时也可以从贝叶斯的观点来对该方法进行解释，反演前后波形的最佳匹配相当于似然函数，而反演反射系数的L1模相当于假设反射系数满足拉普拉斯分布而导出的先验概率分布函数。贝叶斯准则可以把这两者联系起来，最后导出使这两项的加权之和最小的反射系数反演形式。这一非线性反演问题可以采用线性规划、同伦算法、匹配追踪和反复加权迭代等算法来求解。然而这些求解方法运算量较大，对于海量地震数据处理来说，是不可取的。

本发明在对以往的反褶积方法进行了系统分析的基础上，提出了一种新的混合相位子波反褶积方法。

发明内容

为了解决现有技术中存在的上述问题，根据本发明的实施例的一种多路径双谱域混合相位子波反褶积方法，包括如下步骤：

a)从采集到的反射地震记录中选取n道地震记录x(ix，t)，其中，n为大于1的整数且ix＝1，2，…，n；

b)根据选取的n道地震记录确定子波的振幅谱w(f)；

c)根据选取的n道地震记录确定子波的相位谱φ(ω)；

d)根据确定的子波的振幅谱w(f)和相位谱φ(ω)，按如下的数学式进行反傅里叶变换来获得时域的混合相位子波w(t)，

$w\left(t\right)=\underset{-\∞}{\overset{+\∞}{\∫}}w\left(f\right)e^{\mathrm{i\φt}}\mathrm{df};$

e)根据用户基于确定的子波的振幅谱w(f)输入的期望的输出子波的振幅谱w′(f)，利用如下的数学式来确定希望输出的零相位子波w′(t)

$w^{\′}\left(t\right)=\underset{-\∞}{\overset{+\∞}{\∫}}{w}^{\′}\left(f\right)e^{i2\mathrm{\πft}}\mathrm{df};$

f)以混合相位子波w(t)为输入以及零相位子波w′(t)为输出，进行匹配滤波来确定反褶积滤波器d(t)；

g)基于反褶积滤波器d(t)和地震记录x(ix，t)，按如下的数学式来获得分辨率提高的地震记录y(ix，t)，

y(ix，t)＝d(t)*x(ix，t)，其中，*表示褶积。

另外，在所述步骤b)中，

对选取的n道地震记录进行傅里叶变换来获得每个地震记录的振幅谱x(ix，f)；

计算n道地震记录的的振幅谱的平均值

确定在最小二乘意义下使下面的数学式具有最小值时的a_m，从而确定子波的振幅谱w(f)，

$Q=\underset{f_1}{\overset{f_2}{\mathrm{\Σ}}}{\left(\right|\stackrel{\‾}{x}\left(f\right)|-|w\left(f\right)\left|\right)}^2,$

其中，m的数值范围为，L为大于1的整数。

另外，在所述步骤c)中，

利用选取的n道地震记录和三阶累积量方式，确定双谱的相位谱Ψ(ω₁，ω₂)，

根据如下的数学式确定子波的相位谱φ(ω)，

Ψ(ω₁，ω₂)＝φ(ω₁)+φ(ω₂)-φ(ω₁+ω₂)。

另外，在所述步骤c)中

设置l＝1，且

c-1)按如下数学式计算

其中，k＝1，2，...2N-l，这里N是大于1的整数，

这里，PV[x]是相位主算子，表示将[]中的值限定在模为π的-π，π]之内，；

c-2)按如下数学式计算

${\mathrm{\φ}}_p^{l+1}(l+1)=\frac{1}{2N-2l+1}\underset{j=l}{\overset{2N-2l+1}{\mathrm{\Σ}}}{\mathrm{\φ}}_p^l\left(j\right);$

c-3)对l进行加1处理并判断l是否等于N；

c-4)若l不等于N，则重复上述步骤c-1)、c-2)和c-3)；

c-5)若l等于N，则按下述数学式计算φ(ω)，其中，ω为不小于0的整数，

若ω大于1且不大于N，则

若ω大于N且不大于2N-1，则

在本发明中，从地震资料中提取混合相位子波并进行确定性反褶积。其中，子波的振幅谱相对容易估计准确，但相位谱相对来说较难估计。本发明利用高阶统计量中的双谱，提取出地震子波的真实相位，结合谱模拟技术进一步重构地震子波并进行确定性反褶积，同时还提出了一种高阶谱中求取相位时相位解卷的方法，通过求取主值算子，并从多条路径对其平均，不仅使得求取的相位是非卷绕的，而且较好的压制了噪声。本发明提出的方法有效地提高了地震子波的估计精度，从而较大程度上提高地震资料的分辨率。

附图说明

通过下面结合附图对实施例进行的描述，本发明的这些和/或其他方面和优点将会变得清楚和更易于理解，其中：

图1是示出根据本发明的多路径双谱域混合相位子波反褶积方法的流程图。

图2是根据本发明的实施例的基于多路径迭代估计子波的相位的方法的示意图。

图3为反褶积处理前后的的地震数据。

图4为从地震资料中提取的混合相位子波。

图5为利用该子波进行确定性反褶积后的振幅谱比较图。

具体实施方式

现在对本发明实施例进行详细的描述，其示例表示在附图中，其中，相同的标号始终表示相同部件。下面通过参照附图对实施例进行描述以解释本发明。

如图1所示，在步骤S110，从采集到的反射地震记录中选取n道地震记录x(ix，t)，其中，n为大于1的整数且ix＝1，2，…，n。

在步骤S120中，根据选取的n道地震记录确定子波的振幅谱w(f)。稍后将详细描述确定子波的振幅谱w(f)的方法。

在步骤S130中，根据选取的n道地震记录确定子波的相位谱φ(ω)。稍后将详细描述确定子波的相位谱φ(ω)的方法。

在步骤S140中，根据确定的子波的振幅谱w(f)和相位谱φ(ω)，按如下的数学式进行反傅里叶变换来获得时域的混合相位子波w(t)，

$w\left(t\right)=\underset{-\∞}{\overset{+\∞}{\∫}}w\left(f\right)e^{\mathrm{i\φt}}\mathrm{df}.$

在步骤S150中，根据用户基于确定的子波的振幅谱w(f)输入的期望的输出子波的振幅谱w′(f)，利用如下的数学式来确定希望输出的零相位子波w′(t)

$w^{\′}\left(t\right)=\underset{-\∞}{\overset{+\∞}{\∫}}{w}^{\′}\left(f\right)e^{i2\mathrm{\πft}}\mathrm{df}.$

在步骤S160中，以混合相位子波w(t)为输入以及零相位子波w′(t)为输出，进行匹配滤波来确定反褶积滤波器d(t)。

在步骤S170中，基于反褶积滤波器d(t)和地震记录x(ix，t)，按如下的数学式来获得分辨率提高的地震记录y(ix，t)，

y(ix，t)＝d(t)*x(ix，t)。

下面，详细介绍在步骤S120中，确定子波的振幅谱w(f)的方法。

对选取的n道地震记录按如下的等式进行傅里叶变换来获得每个地震记录的振幅谱x(ix，f)。

$x(\mathrm{ix},f)=\left|\frac{1}{2\mathrm{\π}}\underset{-\∞}{\overset{+\∞}{\∫}}x(\mathrm{ix},t)e^{i2\mathrm{\πft}}\mathrm{dt}\right|---(3-1)$

按如下的数学式计算n到地震记录的的振幅谱的平均值

$\stackrel{\‾}{x}\left(f\right)=\frac{\underset{\mathrm{ix}=1}{\overset{n}{\mathrm{\Σ}}}x(\mathrm{ix},f)}{n}---(3-2)$

确定在最小二乘意义下使下面的数学式具有最小值时的a_m，从而确定子波的振幅谱w(f)，

$Q=\underset{f_1}{\overset{f_2}{\mathrm{\Σ}}}{\left(\right|\stackrel{\‾}{x}\left(f\right)|-|w\left(f\right)\left|\right)}^2,---(3-3)$

其中，L为大于1的整数。

下面，详细介绍在步骤S120中，确定子波的相位谱φ(ω)的方法。

首先确定选择的地震记录的三阶累积量，并获得确定双谱的相位谱Ψ(ω₁，ω₂)。

根据如下的数学式确定子波的相位谱φ(ω)。

Ψ(ω₁，ω₂)＝φ(ω₁)+φ(ω₂)-φ(ω₁+ω₂)     (3-4)

通常，如本领善于所公知的对于地震记录s(t)，由Robinson褶积模型，有

s(t)＝w(t)*r(t)+n(t)        ( 3-5)

其中，w(t)、r(t)和n(t)分别表示地震子波、反射系数序列和加性噪声，*表示褶积算子。

对于上述地震记录两边同时取三阶累积量，由累积量的Bartlett-Brillinger-Rosenblatt公式(Mendel，1991)，可以得到上述地震记录的三阶累积量表达式：

c_3s(τ₁，τ₂)＝c_3r(τ₁，τ₂)*m_3w(τ₁，τ₂)+c_3n(τ₁，τ₂)   (3-6)

其中，c_3s(τ₁，τ₂)、c_3r(τ₁，τ₂)和c_3n(τ₁，τ₂)分别表示地震记录、反射系数序列和噪声的三阶累积量，m_3w(τ₁，τ₂)表示地震子波的三阶矩。假设反射系数序列r(t)服从非高斯独立同分布(Sacchi等，2000)，则

$c_{3r}({\mathrm{\τ}}_1,{\mathrm{\τ}}_2)=\left\{\begin{array}{cc}{\mathrm{\γ}}_{3r}& {\mathrm{\τ}}_1={\mathrm{\τ}}_2=0\\ 0& \mathrm{else}\end{array}\right.---(3-7)$

其中，γ_3r为常数，称为反射系数序列r(t)的斜度。假设噪声n(t)为加性高斯噪声(Nikias和Petropulu，1993)，则

c_3n(τ₁，τ₂)＝0(3-8)

将(3-7)式和(3-8)式代入(3-6)式，有

c_3s(τ₁，τ₂)＝γ_3rm_3w(τ₁，τ₂)    (3-9)

将(3-9)式变换到频率域，即可得到基于双谱的地震子波估计表达式

B_s(ω₁，ω₂)＝γ_3rW(ω₁)W(ω₂)W^H(ω₁+ω₂)，    (3-10)

式中，B_s(ω₁，ω₂)为地震记录s(t)的双谱，W(ω)为地震子波w(t)的傅里叶谱，H为共轭算子。将上式的复信号表示成振幅谱和相位谱形式，可得基于双谱的地震子波相位谱估计表达式

Ψ(ω₁，ω₂)＝φ(ω₁)+φ(ω₂)-φ(ω₁+ω₂)            (3-4)

在此，可利用本领域中常用的诸如BMU递归算法、MU算法等方式来从数学式(3-4)获得子波的相位谱φ(ω)。由于相位谱Ψ(ω₁，ω₂)和φ(ω)的取值范围均为(-π，π]，因此，方程左端的取值范围为(-π，π]，而方程右端的取值范围为(-3π，3π]，这样，当方程右端的取值未在(-π，π]时，由

恢复φ(ω)，就会造成相位谱的卷绕问题，即：一值多解的问题。针对此问题可以按图2所述的方式获得子波的相位谱φ(ω)。

图2是根据本发明的实施例的基于多路径迭代估计子波的相位的方法的示意图。

在第l条路径(即，与横轴平行的线段)上，该路径上的相位可通过以下公式以及该路径的初始值(即，图2中的由粗实现表示的等腰三角形的腰上的坐标为(l，l)的点的相位)获得。

而每条路径的初始值为其前一路径(即，与其平行且纵坐标小于1的线条)上的相位值的平均，通过数学式表示为

${\mathrm{\φ}}_p^{l+1}(l+1)=\frac{1}{2N-2l+1}\underset{j=l}{\overset{2N-2l+1}{\mathrm{\Σ}}}{\mathrm{\φ}}_p^l\left(j\right).$

最后，估计的相位为φ(ω)在图2中的等腰三角形内的与其具有相同横坐标的点的相位的平均。

整理上述过程如下，首先，将变量l设置为1，并将

设置为0或1。

然后，按如下数学式计算

其中，k＝1，2，...2N-l，这里N是大于1的整数，

这里，PV[x]是相位主算子，表示将[]中的值限定在模为π的(-π，π]之内，

表示计算相位主值。

然后，按如下数学式计算

${\mathrm{\φ}}_p^{l+1}(l+1)=\frac{1}{2N-2l+1}\underset{j=l}{\overset{2N-2l+1}{\mathrm{\Σ}}}{\mathrm{\φ}}_p^l\left(j\right).$

然后，对l进行加1处理并判断l是否等于N。

若l不等于N，则重复上述步骤c-1)、c-2)和c-3)；若l等于N，则按下述数学式计算φ(ω)，其中，ω为不小于0的整数，

若ω大于1且不大于N，则

若ω大于N且不大于2N-1，则

如图3所示，实施例来源于四川盆地某地区地震勘探资料，该资料已经利用传统的方法进行了提高分辨率处理，但分辨率仍然不能满足确定薄砂体油气储层的目标，为此，利用本发明进行了提高分辨率的进一步处理。图4为从地震资料中提取的混合相位子波，该子波主瓣能量高，旁瓣能量较低，是一典型的混合相位子波。图3为利用该子波进行确定性反褶积后的结果比较(其中，左图为反褶积处理前，而右图为反褶积处理后)。可以看出，反褶积之后的分辨率得到明显改善，沉积内幕和接触关系清晰自然。尤为难得的是，由于消除了反褶积之后剩余相位对弱反射信号的干涉影响，2100-2300ms之间的波阻抗差异较小的弱反射信号得到了有效恢复，为精细地质解释和储层描述提供更加可靠的基础数据。而且图5为利用该子波进行确定性反褶积后的振幅谱比较图(其中，左图为反褶积处理前的振幅谱，而右图为反褶积处理后的振幅谱)，从图中可以看出处理后有效频带的带宽得到很好的拓宽。

在本发明中，从地震资料中提取混合相位子波并进行确定性反褶积。其中，子波的振幅谱相对容易估计准确，但相位谱相对来说较难估计。本发明利用高阶统计量中的双谱，提取出地震子波的真实相位，结合谱模拟技术进一步重构地震子波并进行确定性反褶积，同时还提出了一种高阶谱中求取相位时相位解卷的方法，通过求取主值算子，并从多条路径对其平均，不仅使得求取的相位是非卷绕的，而且较好的压制了噪声。本发明提出的方法有效地提高了地震子波的估计精度，从而较大程度上提高地震资料的分辨率。

虽然已表示和描述了本发明的一些实施例，但本领域技术人员应该理解，在不脱离由权利要求及其等同物限定其范围的本发明的原理和精神的情况下，可以对这些实施例进行修改。

Claims (4)

1.一种多路径双谱域混合相位子波反褶积方法，其特征在于，包括如下步骤：

a)从采集到的反射地震记录中选取n道地震记录x(ix，t)，其中，n为大于1的整数且ix＝1，2，…，n；

b)根据选取的n道地震记录确定子波的振幅谱w(f)；

c)根据选取的n道地震记录确定子波的相位谱φ(ω)；

d)根据确定的子波的振幅谱w(f)和相位谱φ(ω)，按如下的数学式进行反傅里叶变换来获得时域的混合相位子波w(t)，

$w\left(t\right)=\underset{-\∞}{\overset{+\∞}{\∫}}w\left(f\right)e^{\mathrm{i\φt}}\mathrm{df};$

e)根据用户基于确定的子波的振幅谱w(f)输入的期望的输出子波的振幅谱w′(f)，利用如下的数学式来确定希望输出的零相位子波w′(t)

$w^{\′}\left(t\right)=\underset{-\∞}{\overset{+\∞}{\∫}}{w}^{\′}\left(f\right)e^{i2\mathrm{\πft}}\mathrm{df};$

f)以混合相位子波w(t)为输入以及零相位子波w′(t)为输出，进行匹配滤波来确定反褶积滤波器d(t)；

g)基于反褶积滤波器d(t)和地震记录x(ix，t)，按如下的数学式来获得分辨率提高的地震记录y(ix，t)，

y(ix，t)＝d(t)*x(ix，t)，其中，*表示褶积。

2.如权利要求1所述的多路径双谱域混合相位子波反褶积方法，其特征在于，在所述步骤b)中，

对选取的n道地震记录进行傅里叶变换来获得每个地震记录的振幅谱x(ix，f)；

计算n道地震记录的的振幅谱的平均值

确定在最小二乘意义下使下面的数学式具有最小值时的a_m，从而确定子波的振幅谱w(f)，

$Q=\underset{f_1}{\overset{f_2}{\mathrm{\Σ}}}{\left(\right|\stackrel{\‾}{x}\left(f\right)|-|w\left(f\right)\left|\right)}^2,$

其中，

m的数值范围为，L为大于1的整数。

3.如权利要求1所述的多路径双谱域混合相位子波反褶积方法，其特征在于，在所述步骤c)中，

利用选取的n道地震记录和三阶累积量方式，确定双谱的相位谱Ψ(ω₁，ω₂)，

根据如下的数学式确定子波的相位谱φ(ω)，

Ψ(ω₁，ω₂)＝φ(ω₁)+φ(ω₂)-φ(ω₁+ω₂)。

4.如权利要求3所述的多路径双谱域混合相位子波反褶积方法，其特征在于，

设置l＝1，且

c-1)按如下数学式计算

其中，k＝1，2，...2N-l，这里N是大于1的整数，

这里，PV[x]是相位主算子，表示将[]中的值限定在模为π的(-π，π]之内，；

c-2)按如下数学式计算

${\mathrm{\φ}}_p^{l+1}(l+1)=\frac{1}{2N-2l+1}\underset{j=l}{\overset{2N-2l+1}{\mathrm{\Σ}}}{\mathrm{\φ}}_p^l\left(j\right);$

c-3)对l进行加1处理并判断l是否等于N；

c-4)若l不等于N，则重复上述步骤c-1)、c-2)和c-3)；

c-5)若l等于N，则按下述数学式计算φ(ω)，其中，ω为不小于0的整数，

若ω大于1且不大于N，则

若ω大于N且不大于2N-1，则

Images