CN104614767A - 基于分段延拓的时变地震子波相位校正方法 - Google Patents

Abstract

本发明涉及一种基于分段延拓的时变地震子波相位校正方法，属于地震勘探信号处理领域，其特征在于：针对时变子波振幅影响已经完全消除后相位时变的地震记录，首先将其分成若干相互重叠的等长片段，采用基于纯相位滤波器寻优的方法在每段中提取子波相位残余；然后通过多项式拟合法将子波相位残余延拓至每一个时间点；最后通过基于线性时频变换的时变相位反褶积方法恢复反射系数序列，以实现时变子波相位校正。经数据仿真验证，本发明提出的方法可以有效实现非平稳地震记录的时变子波相位校正，从而提高地震记录的分辨率。

Description

基于分段延拓的时变地震子波相位校正方法

技术领域：

本发明属于地震勘探信号处理领域。

背景技术：

如今油气田勘探和开发正向小规模、薄储层的方向发展，对地震勘探的精度要求越来越高。为适应油气藏的动态预测和寻找复杂构造和岩性油气藏，要求地震处理后的剖面具有高信噪比、高分辨率和高保真度的“三高”特点。

地震勘探包括野外采集、室内处理和地震资料解释三大环节。要提高地震记录的分辨率就必须从采集、处理及解释各个环节上下功夫，挖掘各个环节的能力，使各个环节上都达到高质量，最后才能获得具有“三高”特点的地震剖面。而如果已知地震子波，将该子波与地震记录作反褶积，就可求得地层的反射系数，从而可得一个高分辨率的地震剖面。地震子波估计(子波提取)作为地震资料反褶积处理、地震波阻抗反演和地震正演模型的基础，其提取精度直接影响后续地震资料处理和地震资料解释的可靠性和准确性。

在地震信号处理的发展过程中，子波相位经历了零相位、最小相位和混合相位的假设，目前通常认为叠后地震子波相位应当为接近零相位的混合相位。传统的子波提取方法和反褶积方法中，为了方便处理，往往将子波相位假设为最小相位。虽然相关学者研究出了多种针对混合相位子波的提取方法，但由于方法本身假设过于严苛，且提取过程中受算法误差和计算工具精度影响，提取结果往往与实际不符。由于子波相位提取的不准确，在后续反射系数序列反演(反褶积)结果中残留的子波相位易造成误差的放大和假象的产生，严重影响后续地震解释工作的进行。

相位的提取和控制在地震信号处理领域扮演着十分重要的角色，传统的地震数据处理和解释方法一般基于地震记录平稳假设，即地震子波为平稳的、不随时间而改变的。然而由于地层的吸收衰减作用，地震波在实际传播过程中是动态衰减的，则采集到的原始地震记录应当是非平稳的动态衰减过程。由于时变子波相位处理的复杂性，目前时变地震记录处理技术通常针对时变振幅进行处理，往往忽略时变相位影响或对其简化处理，如时变常相位校正、时变线性相位滤波等。目前暂未检索到时变相位处理的有效方法，但传统的叠加、偏移、补偿等处理后的地震记录中由于子波时变相位的存在，严重影响了各类处理结果的精度，故亟需一种有效的时变子波相位校正方法，以提高实际地震资料的分辨率。

本发明提出一种基于分段延拓的时变地震子波相位校正方法，以解决非平稳地震记录中子波相位校正问题，实现提高地震记录分辨率的目的。

发明内容

本发明的目的在于提出一种基于分段延拓的时变地震子波相位校正方法，其中包括地震记录的时变子波模型；子波相位校正分段提取方法；相位延拓方法；基于线性时频变换的时变反褶积方法。

本发明的特征之一在于建立了地震记录的时变子波模型。

传统的平稳地震记录模型可以用Robinson褶积模型表示：

$x\left(t\right)={\∫}_{-\∞}^{\∞}w(t-\mathrm{\τ})r\left(\mathrm{\τ}\right)\mathrm{d\τ}=w\left(t\right)*r\left(t\right)---\left(1\right)$

该模型基于若干基本假设，其中最重要的假设即为地震子波在传播过程中具有时不变性质。但实际子波在传播过程中是随时间改变的，则式(1)不能够准确描述实际地震记录。

为了更好地描述实际地震记录的特性，本发明将在式(1)的基础上给出地震记录时变子波模型。此处仅考虑平面波垂直入射，并在水平层状粘弹性介质中传播的情况。地震记录的时间采样间隔用Δt表示，则t_k＝kΔt表示第k个采样点对应的时刻，反射系数序列在第k个采样点的值用c_k表示，其符号代表反射系数的极性。将每一个采样点对应的反射系数视为独立的，即当且仅当t＝t_k时反射系数序列的值为c_k，其它采样点的数值均为零，则第k个采样点的反射系数可视为以c_k为权系数的脉冲，即：

r_k(n)＝c_kδ(n-k)      (2)其中δ(n-k)为t_k点的具有单位能量的脉冲函数。令w_k表示t_k点对应的地震子波，假设检波器仅接收到t_k点反射的地震子波，则地震记录可以表示为：

x_k(n)＝w_k(n)*r_k(n)＝c_kw_k(n)*δ(n-k)      (3)将上式通过傅里叶变换映射到频域可得：

X_k(e^jω)＝c_kW_k(e^jω)e^-jωk      (4)通过傅里叶逆变化可以获得式(3)的另一种表达形式，即：

$x_k\left(n\right)=\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}{w}_k\left(e^{\mathrm{j\ω}}\right)c_k{e}^{-\mathrm{j\ωk}}{e}^{\mathrm{j\ωn}}\mathrm{d\ω}---\left(5\right)$

当检波器接收到来自所有时间点的反射波时，其时间域表示为：

$x\left(n\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{x}_k\left(n\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{c}_k{w}_k\left(n\right)*\mathrm{\δ}(n-k)---\left(6\right)$

其中N为数据长度。地震记录的频率域表示为：

$X\left(e^{\mathrm{j\ω}}\right)=\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{W}_k\left(e^{\mathrm{j\ω}}\right)c_k{e}^{-\mathrm{j\ωk}}---\left(7\right)$

式(6)和式(7)即为地震记录的时变子波模型。该模型提供了一种合成时变地震记录的方法，即将反射系数每个点与对应的地震子波褶积后在时域进行移位叠加，或者褶积后在频域叠加后再反变换回时域。将式(7)进行逆傅里叶变换可得：

$x\left(n\right)=\frac{1}{2\mathrm{\π}}{\∫}_{-\mathrm{\π}}^{\mathrm{\π}}\underset{k=1}{\overset{N}{\mathrm{\Σ}}}{W}_k\left(e^{\mathrm{j\ω}}\right)c_k{e}^{-\mathrm{j\ωk}}{e}^{\mathrm{j\ωn}}\mathrm{d\ω}---\left(8\right)$

假设时变子波振幅谱通过时频分析、谱模拟等方法已经准确提取，并且进行了反褶积处理，完全消除了子波振幅影响。若用代表第k个采样点的子波残余相位谱，类似于时变子波模型，则用以相位校正的时变子波纯相位模型可以表示为：

同样，在构造出时变纯相位子波的条件下可以由上式合成相位时变的地震记录。以上公式满足以下假设：

AS1)子波残余相位的各个频率成分均为光滑连续变化的。

AS2)地层反射系数序列r(n)为非高斯过程，满足零均值、独立同分布假设，方差其四阶累积量|γ_4r|＜∞。

本发明的特征之二在于将相位时变的地震记录分为若干重叠的等长片段，并采用相位寻优法提取每段的子波相位残余。

本发明针对消除时变子波振幅影响后的时变子波相位进行校正，虽然基于自适应分子分解或者其它自适应分段方法在时变振幅处理方面有着诸多优点，但以上方法一般采用分段后振幅的相似度作为基准进行分段，但本发明应用在完全消除振幅影响的条件下；若采用相位基准进行分段则有可能受反射系数序列相位细节的影响，故自适应分段方法不适用于本发明。

平均分段无需考虑振幅的影响，在相位连续变化的情况下可以提取出每段的平均相位。但本发明应用统计性方法提取每段中的子波相位残余需要足够的数据长度保证提取精度，分段长度的增加必然导致分段数目的减少，而较少的分段数目则会影响后续的相位延拓处理。综合各类影响因素，本发明采用平均重叠分段方法实现子波相位校正的分段，在保证分段长度的前提下提高分段数目。

在每段中提取子波相位残余，方法采用张亚南等人提出的基于纯相位滤波器的地震子波相位校正方法。(张亚南,戴永寿,丁进杰等.基于纯相位滤波器的子波相位校正方法研究[J].石油物探,2014,53(1):26-31.)该方法通过ARMA模型和振幅谱单位化的方式构造纯相位滤波器以描述地震子波相位残余，在最大方差模准则约束下采用改进的粒子群算法针对残余的子波相位进行准确寻优，以实现每段中的子波相位残余的提取，提取结果作为该段中间点所对应的相位残余。

本发明的特征之三在于采用多项式拟合的方式将每段中的子波相位残余延拓至每一个时间点。

由于分段提取的结果仅为该段中一条平均相位谱，显然该相位谱无法准确描述整段记录的相位变化，故需要将其延拓至每一时间点，从而产生一系列随时间变化的相位残余，以近似描述每一点的时变相位残余，进而实现时变子波相位的校正。本发明采用多项式拟合方式实现子波相位残余的延拓。

假设N段子波提取出的相位谱中相同频率成分的数据点为(x_i,y_i)，x_i为时间点，y_i为该时间点对应的相位谱数值，其取值范围是(-π,π)。Φ为所有阶次不超过n(n≤N)的多项式的集合，求取拟合函数 $f_n\left(x\right)=\underset{k=0}{\overset{n}{\mathrm{\Σ}}}{a}_k{x}^k\∈\mathrm{\Φ},$ 使之满足：

当拟合函数为多项式时，称之为多项式拟合，满足式(11)的f_n(x)称为最小二乘拟合多项式。当n＝1时，称之为线性拟合。

上式中为a₀,a₁,K,a_n的多元函数，因此多项式拟合的问题可以转化为求取I＝I(a₀,a₁,K,a_n)极值的问题。由多元函数极值求取的必要条件可知：

$\frac{\∂I}{{\∂a}_l}=2\underset{i=0}{\overset{N}{\mathrm{\Σ}}}(\underset{i=0}{\overset{n}{\mathrm{\Σ}}}{a}_k{x}_i^k-y_i)x_i^l=0,l=\mathrm{0,1},K,n---\left(12\right)$

上式同样可写为：

$\underset{k=0}{\overset{n}{\mathrm{\Σ}}}\left(\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^{l+k}\right)a_k=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^l{y}_i,l=\mathrm{0,1},K,n---\left(13\right)$

式(13)为关于a₀,a₁,K,a_n的线性方程组，可以写成矩阵形式：

$\left[\begin{array}{cccc}N+1& \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i& L& \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^n\\ \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i& \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^2& L& \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^{n+1}\\ M& M& O& M\\ \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^n& \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^{n+1}& L& \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^{2n}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ M\\ a_n\end{array}\right]=\left[\begin{array}{c}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{y}_i\\ \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i{y}_i\\ M\\ \underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^n{y}_i\end{array}\right]---\left(14\right)$

式(13)和(14)被称为正规方程组或法方程组，其系数矩阵为对称正定阵，存在唯一解。

求解式(14)中的a_k(k＝0,1,K,n)即可获得多项式f_n(x)，其均方误差记为：

${\left|\right|R_{\mathrm{es}}\left|\right|}_2^2=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{[f_n\left(x_i\right)-y_i-]}^2---\left(15\right)$

类似于式(12)，可以写为：

${\left|\right|R_{\mathrm{es}}\left|\right|}_2^2=\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{y}_i^2-\underset{k=0}{\overset{n}{\mathrm{\Σ}}}{a}_k\left(\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^k{y}_i\right)---\left(16\right)$

综上所述，多项式拟合方法可以总结为以下步骤：

Step1：通过已知数据点画出散点图，确定拟合多项式阶次n；

Step2：列表计算 $\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^l(l=\mathrm{0,1},K,2n)$ 和 $\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{x}_i^l{y}_i(l=\mathrm{0,1},K,2n);$

Step3：列正规方程组，求取参数a₀,a₁,K,a_n；

Step4：写出拟合多项式

在应用过程中，阶次n一般为n＜m，当n＝m时所获得的拟合多项式即为拉格朗日插值多项式或牛顿插值多项式。

通过以上方式即可得到单一频率成分的拟合多项式，针对频域每一个频率成分相位谱均进行多项式拟合，即可实现子波相位残余的延拓，得到全时间范围内的相位曲面。

本发明的特征之四在于采用基于线性时频变换的时变相位反褶积实现反射系数序列的恢复。

本发明由待校正的地震记录的时变特性出发，在提取出分段相位残余并延拓至每一时间点的基础上，采用基于线性时频变换的时变相位反褶积方法恢复原始反射系数序列。

线性时频分析可以有效表征每一时间点处的频率特性，是研究时变信号的有力工具。常用的时频变换方法主要有短时傅里叶变换、Gabor变换、小波变换和S变换等，其中S变换综合了短时傅里叶变换和小波变换的优点，并且克服了短时傅里叶变换时频分辨率固定和小波变换尺度因子与频率无关的缺点，具有良好的时域和频域分辨率，应用前景广泛。

S变换采用的可变高斯窗函数与频率有关，定义为：

$w\left(t\right)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{f^2{t}^2}{2}}---\left(17\right)$

由式(17)可以看出，S变换的窗函数变化规律固定，不能根据实际需求进行调整，故相关学者对S变换的窗函数进行改进，提出了一系列广义S变换方法。常用的广义S变换采用时窗宽度随频率f呈反比例变化的高斯窗函数，在高频段时窗较窄，以获得较高的时间分辨率；在低频段时窗较宽，以获得较高的频率分辨率。以上广义S变换对于大多数非平稳信号能获得较好的时频变换结果，但在时变地震记录处理中分辨率要求较高，以上方法不能满足要求。为此，齐春燕等在广义S变换的基础上对其进行改进，提出了一种改进的广义S变换，其窗函数定义为：(齐春艳,李彦鹏,彭继新等.一种改进的广义S变换[J].石油地球物理勘探,2010,45(2):215-218)

$G(t,f)=\frac{1}{\sqrt{2\mathrm{\π}}{\left|f\right|}^p}\exp (-\frac{t^2}{{2q}^2{f}^{2p}})---\left(18\right)$

其中p,q为调节因子，且p＞0,q＞0。由式(5-29)可以看出，该窗函数的宽度与频率成正比，即在低频处时窗较窄以获得较高的时间分辨率；在高频处时窗较宽以获得较高的频率分辨率。对信号x(t)采用上述窗函数实现改进的广义S变换：

$\mathrm{NGST}(\mathrm{\τ},f)={\∫}_{-\∞}^{\∞}x\left(t\right)\left\{\frac{1}{\sqrt{2\mathrm{\π}}q{\left|f\right|}^p}\exp \left(\frac{-{(\mathrm{\τ}-t)}^2}{{2q}^2{f}^{2p}}\right)\exp (-i2\mathrm{\πft})\right\}\mathrm{dt},f\≠0---\left(19\right)$

采用该方法对地震记录进行时频变换，可以获得更高的时频分辨率，并且时频聚焦性更好，能够有效分辨出地震记录的频率成分随时间变化的规律。

在提取分段子波相位残余并延拓至每一时间点的基础上，对待校正记录采用改进的广义S变换进行时频变换，在变换后的每一个时间点的频谱中减去相对应的子波相位残余，并将处理后的结果反变换回时域，从而实现时变相位反褶积。

时变相位反褶积处理后，通过改进的广义S反变换将处理结果转换到时域。改进的广义S反变换的结果与正变换参数无关，是正变换的逆过程。改进的广义S反变换的表达式定义为：

$x^{\′}\left(t\right)={\∫}_{-\∞}^{\∞}\left\{{\∫}_{-\∞}^{\∞}\mathrm{NGST}(\mathrm{\τ},f)\exp (-i2\mathrm{\πft})\mathrm{d\τ}\right\}\exp \left(i2\mathrm{\πft}\right)\mathrm{df}---\left(20\right)$

改进的广义S反变换的计算方法为：

Step1针对NGST(τ,f)计算τ→f的傅里叶变换，得到NGST(f,f)；

Step2计算频谱H(f)，令其满足H(f)＝NGST(f,f)；

Step3对H(f)做f→τ的反傅里叶变换，得到相位反褶积后的信号x′(t)。

通过以上方法得到的信号x′(t)与原始信号x(t)相比，仅去除了时变子波相位残余的成分，改变了信号的相位谱；对信号的振幅未作处理，故未改变信号的能量，即时变相位反褶积后的信号x′(t)相对于原始信号x(t)而言是无损的。至此，完成了基于线性时频变换的时变相位反褶积处理。

综上所述，经过对非平稳地震记录平均重叠分段、基于纯相位滤波器寻优的方法提取每段子波相位残余、采用多项式拟合法将残余相位延拓至每一时间点和基于线性时频变换的时变相位反褶积等处理，可实现非平稳地震记录中时变子波相位的校正，最终形成一种基于分段延拓的时变地震子波相位校正方法。

附图说明

图1、原始的反射系数序列

图2、合成的相位时变的待校正地震记录

图3、通过本发明方法校正后的反射系数序列

图4、子波相位残余延拓结果

图5、理论的子波相位残余

图6、提取相位与理论相位之差

图7、基于分段延拓的时变地震子波相位校正方法实施流程

具体实施方式

本发明提出了基于分段延拓的时变地震子波相位校正方法。这一方法的特点是：通过分段提取、逐点延拓的思想实现地震记录中时变子波相位的校正。这一方法首先将待校正的相位时变地震记录分成若干个等长的重叠片段，通过基于纯相位滤波器寻优的方法提取出每段的平均子波相位残余作为该段中间点对应的子波相位残余；采用多项式拟合的方法将分段提取的子波相位残余延拓至每一个时间点；最后通过基于线性时频变换的时变相位反褶积方法恢复反射系数序列，以实现时变地震子波相位的校正。本发明的实施流程如图7所示。

本发明按照以下步骤实施：

1.针对时变振幅影响已经完全消除的相位时变地震记录，将其分成若干相互重叠的等长片段。

2.采用基于纯相位滤波器寻优的方法提取每段中的相位残余，并作为该段中间点对应的相位残余。

3.通过多项式拟合法将相位残余延拓至每一个时间点。

4.采用改进的广义S变换方法将待校正的地震记录变换到时频域。

5.通过谱除法在时频域的每一个时间点中去除子波相位残余。

6.采用改进的广义S反变换将处理后的时频域地震记录反变换回时域，实现子波相位校正。

Claims (2)

1.基于分段延拓的时变地震子波相位校正方法，其特征在于将地震记录分成若干等长的重叠片段，在每段中提取子波相位残余；将分段提取的子波相位残余延拓至每一个时间点；将地震记录变换到时频域，在时频域去除每一时间点的子波相位残余，并反变换回时域，以实现时变子波相位校正。

2.基于分段延拓的时变地震子波相位校正方法，其特征在于该方法依次含有以下步骤：

步骤(1)针对时变振幅影响已经完全消除的相位时变地震记录，将其分成若干相互重叠的等长片段；

步骤(2)采用基于纯相位滤波器寻优的方法提取每段中的相位残余，并作为该段中间点对应的相位残余；

步骤(3)通过多项式拟合法将相位残余延拓至每一个时间点；

步骤(4)采用改进的广义S变换方法将待校正的地震记录变换到时频域；

步骤(5)通过谱除法在时频域的每一个时间点中去除子波相位残余；

步骤(6)采用改进的广义S反变换将处理后的时频域地震记录反变换回时域，实现子波相位校正。