CN102269824A - 一种地震数据子波相位转换处理的方法 - Google Patents

Abstract

本发明涉及石油勘探开发的一种地震数据子波相位转换处理的方法，计算地震数据零相位子波，Z变换后分解成N个相乘的因式，将特征因子按照绝对值从小到大依次排列，按绝对值的大小将特征因子分成单位圆内和单位圆外两部分，构成最小相位子波的纯相位转换因子，将零相位子波转换为最小相位子波，再波转换为零相位子波，将零相位子波地震数据转换为最小相位子波地震数据，再转换为零相位子波地震数据，实现地震子波和地震数据的纯相位转换。本发明通过反褶积处理可以有效提高地震数据分辨率，可以有效提高地震数据视分辨率，且计算简单。

Description

一种地震数据子波相位转换处理的方法

技术领域

本发明涉及油田的勘探、开发、开采技术，具体是为反映地下地层层位、油藏描述提供高分辨率的地震图形和数据的一种地震数据子波相位转换处理的方法，特别适用于地震数据处理反褶积之前的实际地震数据。

背景技术

地面地震勘探的过程，就是在地面上的一系列点上，利用人工激发地震波，地震波向地下传播，当遇到波阻抗(地震波在地层介质中向地下传播的速度与介质密度的乘积)界面(即上下地层波阻抗不相等面)时，在波阻抗界面上地震波产生反射现象，地震波传播方向发生改变，地震波开始向上传播，在地面上的一系列接收点上安置着接收器，接收向上传播的地震波数据，完成野外地震勘探。地面接收点上接收器接收到的信号就是地震数据记录，它表示地下地层介质的反射系数序列与地震子波的褶积。然而，实际接收到的地震数据还包含着激发点和接收点空间位置和排列位置的信息和各种噪声干扰等。地震数据处理就是对野外勘探过程中向上传播的地震数据记录进行处理，保留反映地下地层波阻抗界面的信息，而消除其它的信息，这种信息就是叠后地震数据。这种叠后地震数据仅反映地下地层的结构和构造。

随着地震勘探的发展，油气勘探的难度和深度也越来越大，对地震资料的分辨率要求也越来越高。反褶积是提高地震数据分辨率最有效的途径。地震褶积模型是地震勘探数据处理中最基本的模型之一。常规的地震数据反褶积处理都是基于地震褶积模型。反褶积处理的一个根本假设就是地震子波是最小相位子波。这个假设使得我们可以由子波振幅谱得到子波的相位谱并继而得到完整的子波。而实际数据的地震子波是混合相位的，这便使以该假设为前提的处理产生误差，为降低地震数据的处理误差，有必要研究把零相位子波转换为最小相位子波的方法和技术。而在所有地震子波中，零相位地震子波具有最高的视分辨率，而其它相位子波的分辨率都没有零相位子波高，因此为了提高地震数据的视分辨率，必须把最小相位子波的地震数据转换为零相位子波的地震数据。

当前我们在地震数据反褶积处理时，并不考虑地震数据中地震子波相位问题，而是直接使用地震反褶积处理，并不考虑反褶积的最小相位假设。本发明就是一种把零相位子波地震数据转换为最小相位子波地震数据，通过地震反褶积处理提高地震数据分辨率；或者把最小相位子波地震数据转换为零相位子波地震数据，以提高地震数据视分辨率。

发明内容

本发明目的在于通过反褶积处理以提高地震数据处理分辨率，以有效提高地震数据视分辨率，且计算简单的地震数据子波相位转换处理的方法。

本发明采用如下技术方案，包括以下步骤：

1)用地震人工震源激发和采集地震数据并做预处理；

步骤1)所述的预处理是指对地震数据置标签、定义观测系统。

2)确定地震数据子波参数和计算地震数据零相位子波；

步骤2)所述的确定地震数据子波参数，是指对地震数据进行振幅谱分析，根据地震数据振幅谱形状，确定原始数据子波的参数。然后根据地震数据子波参数，计算地震数据零相位子波。

3)对子波进行Z变换；

步骤3)所述的对子波进行Z变换是根据地震数据零相位子波，按照下式计算：

$B\left(z\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{b}_n{z}^{-n}---\left(1\right)$

式中，

是单一的地震零相位子波序列{b_i}，其中(N+1)表示子波长度(序列点数)。

4)对子波Z变换进行因式分解并将特征因子按照从小到大排列；

按下式分解成N个相乘的因式：

$B\left(z\right)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}(1-{\hat{b}}_i{z}^{-1})=(1-{\hat{b}}_1{z}^{-1})(1-{\hat{b}}_2{z}^{-1})\mathrm{\Λ}\left({1-\hat{b}}_N{z}^{-1}\right)---\left(3\right)$

且每个因式的特征因子

按照绝对值从小到大依次排列，

$\left|{\hat{b}}_1\right|\≤\left|{\hat{b}}_2\right|\≤\mathrm{\Λ}\≤\left|{\hat{b}}_N\right|---\left(4\right)$

5)将特征因子分成单位圆内和单位圆外两部分，按特征因子绝对值的大小，将特征因子分成两部分：

第一部分是特征因子的绝对值小于或者等于1，

共有K(K≥1)个，其中K是Z平面上单位圆内和单位圆上特征因子的总个数；

第二部分是特征因子的绝对值大于1，

共有M(M≥1)个，其中M是Z平面上单位圆外特征因子的总个数，且M+K＝N；

6)根据单位圆外的特征因子，构成将零相位子波转换为最小相位子波的纯相位转换因子H₁(z)：

$H_1\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}\frac{(1-\frac{1}{{\hat{b}}_{K+i}^{*}}{z}^{-1})}{(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})}---\left(6\right)$

根据特征方程单位圆外的特征因子，构成将最小相位子波转换为零相位子波的纯相位转换因子H₂(z)：

$H_2\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}\frac{(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})}{(1-\frac{1}{{\hat{b}}_{K+i}^{*}}{z}^{-1})}---\left(7\right)$

式中：“*”表示取共轭；

7)根据最小相位子波的纯相位转换因子H₁(z)，将零相位子波转换为最小相位子波：

W_N(z)＝H₁(z)B_N(z)                 (8)

根据零相位子波的纯相位转换因子H₂(z)，将最小相位子波转换为零相位子波：

B_N(z)＝H₂(z)W_N(z)                 (9)

式中：W_N(z)表示最小相位子波的Z变换，B_N(z)表示零相位子波的Z变换；

8)根据纯相位转换因子H₁(z)，将零相位子波地震数据转换为最小相位子波地震数据：

Y_W(z)＝H₁(z)X_B(z)                 (10)

式中：Y_W(z)表示最小相位子波地震数据的Z变换，X_B(z)表示零相位子波地震数据的Z变换；

根据零相位子波的纯相位转换因子H₂(z)，将最小相位子波地震数据转换为零相位子波地震数据：

Y_B(z)＝H₂(z)X_W(z)                   (11)

式中：Y_B(z)表示零相位子波地震数据的Z变换，X_W(z)表示最小相位子波地震数据的Z变换；

9)根据纯相位转换因子实现地震子波和地震数据的纯相位转换。

步骤9)所述的纯相位转换是在时间域内，通过递推和递归运算实现纯相位转换。

步骤9)所述的纯相位转换是将零相位子波地震数据转换为最小相位子波地震数据：

$y_i=\left\{\begin{array}{cc}\frac{u_0}{q_0}& i=0\\ \frac{1}{q_0}(u_i-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{q}_j{y}_{i-j})& i=\mathrm{1,2},K,M\\ \frac{1}{q_0}(u_i-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{q}_j{y}_{i-j})& i=M+1,M+2,K,M+L\end{array}\right.---\left(12\right)$

式中：序列{y_i}为最小相位地震数据或最小相位地震子波，长度为(L+M+1)，且

$u_i=\underset{j=0}{\overset{L}{\mathrm{\Σ}}}{x}_j{p}_{i-j}---\left(13\right)$

$p_0=1,p_i={(-1)}^i\underset{k_1=1}{\overset{M-i+1}{\mathrm{\Σ}}}\underset{k_2=2}{\overset{M-i+2}{\mathrm{\Σ}}}\mathrm{\Λ}\underset{k_i=i}{\overset{M}{\mathrm{\Σ}}}\frac{1}{{\hat{b}}_{k_1}^{*}{\hat{b}}_{k_2}^{*}\mathrm{\Λ}{\hat{b}}_{k_i}^{*}}---\left(14\right)$

$q_M=1,q_{M-i}={(-1)}^i\underset{k_1=1}{\overset{M-i+1}{\mathrm{\Σ}}}\underset{k_2=2}{\overset{M-i+2}{\mathrm{\Σ}}}\mathrm{\Λ}\underset{k_i=i}{\overset{M}{\mathrm{\Σ}}}\frac{1}{{\hat{b}}_{k_1}{\hat{b}}_{k_2}\mathrm{\Λ}{\hat{b}}_{k_i}}---\left(15\right)$

式中：i＝1，2，Λ，M，

式中：序列{x_i}为零相位地震数据或零相位地震子波，长度为(L+1)。

在以上步骤中序列{x_i}为零相位地震数据，则式(12)中相位转换后的序列{y_i}为最小相位地震数据；如果在式(13)中，序列{x_i}为零相位地震子波，则方程(12)中相位转换后的序列{y_i}为最小相位地震子波。

步骤9)所述的纯相位转换是将最小相位子波地震数据转换为零相位子波地震数据：

$y_i=\left\{\begin{array}{cc}\frac{v_0}{p_0}& i=0\\ \frac{1}{p_0}(v_i-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{p}_j{y}_{i-j})& i=\mathrm{1,2},K,M\\ \frac{1}{p_0}(v_i-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{p}_j{y}_{i-j})& i=M+1,M+2,K,M+L\end{array}\right.---\left(16\right)$

式中：序列{y_i}为零相位地震数据或零相位地震子波，长度为(L+M+1)，且

$v_i=\underset{j=0}{\overset{L}{\mathrm{\Σ}}}{x}_j{q}_{i-j}---\left(17\right)$

式中：i＝1，2，Λ，M，

在以上步骤中序列{x_i}为最小相位地震数据或最小相位地震子波，其长度为(L+1)。在式(17)中，序列{x_i}为最小相位地震数据，则式(16)中相位转换后的序列{y_i}为零相位地震数据；如果在式(17)中，序列{x_i}为最小相位地震子波，则式(16)中相位转换后的序列{y_i}为零相位地震子波。

10)绘制相位转换后的地震数据剖面和存储相位转换后的地震数据。

本发明把零相位地震子波转换为最小相位地震子波，把零相位子波的地震数据转换为最小相位子波的地震数据，并且处理后的数据子波一定是最小相位，通过反褶积处理可以有效提高地震数据分辨率。本发明把最小相位地震子波转换为零相位地震子波，把最小相位子波的地震数据转换为零相位子波的地震数据，并且处理后的数据子波一定是零相位，可以有效提高地震数据视分辨率。

本发明仅对地震数据子波的相位进行处理，而不改变数据的振幅谱，计算简单。

附图说明

图1子波对比

图1-1低通子波LP(25，35)对比

图1-2高通子波HP(35，40)对比

图1-3带通子波BP(15，20，50，60)对比

图1-4陷频子波NP(35，50，65)对比

图1-5Ricker子波RP(25)对比

图1-6俞氏子波YP(15，50)对比

(a)零相位子波；

(b)最小相位转换因子；

(c)最小相位子波

图2子波频谱对比

图2-1低通子波LP(25，35)频谱对比

图2-2高通子波HP(35，40)频谱对比

图2-3带通子波BP(15，20，50，60)频谱对比

图2-4陷频子波NP(35，50，65)频谱对比

图2-5Ricker子波RP(25)频谱对比

图2-6俞氏子波YP(15，50)频谱对比

(a)零相位子波；

(b)最小相位转换因子；

(c)最小相位子波

图3子波反褶积结果对比

(a)零相位子波；

(b)(a)的反褶积结果；

(c)最小相位子波；

(d)(c)的反褶积结果

图4本发明与OMEGA地震数据系统相位转换处理后叠加数据对比

(a)零相位子波叠加数据

(b)OMEGA最小相位子波转换后叠加数据

(c)本发明最小相位子波转换后叠加数据

图5本发明相位转换处理前后叠加数据对比

(a)本发明最小相位子波转换前叠加数据

(b)本发明最小相位子波转换后叠加数据

图6本发明最小相位相位转换处理前后反褶积处理后叠加数据对比

(a)本发明最小相位转换处理前数据反褶积叠加数据

(b)本发明最小相位转换处理后数据反褶积叠加数据

具体实施方式

本发明采用如下技术方案实现，包括以下步骤：

(1)用地震人工震源激发和采集地震数据并对地震数据置标签、定义观测系统等预处理工作；

(2)确定地震数据子波参数和计算地震数据零相位子波，即对地震数据进行振幅谱分析，根据地震数据振幅谱形状，确定原始地震数据地震子波的参数。然后根据地震数据子波参数，计算地震数据零相位子波。

理想地震子波参数和子波计算公式如下：

理想低通零相位子波表达式为

$h_{\mathrm{lp}}\left(t\right)=\frac{\sin \mathrm{\π}(f_2+f_1)t\sin \mathrm{\π}(f_2-f_1)t}{(f_2-f_1){\mathrm{\π}}^2{t}^2}---\left(1\right)$

其中，f₁和f₂是理想低通零相位子波频率参数，时间范围为-∞＜t＜∞。

理想高通零相位子波表达式为

$h_{\mathrm{hp}}\left(t\right)=\frac{\sin 2\mathrm{\π}{f}_{N}t}{\mathrm{\πt}}-\frac{\sin \mathrm{\π}(f_2+f_1)t\sin \mathrm{\π}(f_2-f_1)t}{(f_2-f_1){\mathrm{\π}}^2{t}^2}---\left(2\right)$

其中时间范围为-∞＜t＜∞，f₁和f₂是理想高通零相位子波频率参数，f_N是奈奎斯特频率。

理想带通零相位子波表达式为

$h_{\mathrm{bp}}\left(t\right)=\frac{\sin \mathrm{\π}(f_4+f_3)t\sin \mathrm{\π}(f_4-f_3)t}{(f_4-f_3){\mathrm{\π}}^2{t}^2}-\frac{\sin \mathrm{\π}(f_2+f_1)t\sin \mathrm{\π}(f_2-f_1)t}{(f_2-f_1){\mathrm{\π}}^2{t}^2}---\left(3\right)$

其中，f₁、f₂、f₃和f₄是理想带通零相位子波频率参数，时间范围为-∞＜t＜∞。

理想陷频零相位子波表达式为

$h_{\mathrm{bp}}\left(t\right)=\frac{\sin \mathrm{\π}(f_3+f_2)t\sin \mathrm{\π}(f_3-f_2)t}{(f_3-f_2){\mathrm{\π}}^2{t}^2}-\frac{\sin \mathrm{\π}(f_2+f_1)t\sin \mathrm{\π}(f_2-f_1)t}{(f_2-f_1){\mathrm{\π}}^2{t}^2}---\left(4\right)$

其中，f₁、f₂和f₃是理想陷频零相位子波频率参数，时间范围为-∞＜t＜∞。

理想Ricker零相位子波表达式为

h_r(t)＝[1-2(πft)²]exp[-(πft)²]            (5)

其中，f是理想Ricker零相位子波中心频率参数，时间范围为-∞＜t＜∞。

理想俞氏零相位子波表达式为

$h_y\left(t\right)=\frac{1}{f_2-f_1}\left\{f_2\exp \right[-{\left(\mathrm{\π}{f}_{2}t\right)}^2]-f_1\exp [-{\left(\mathrm{\π}{f}_{1}t\right)}^2\left]\right\}---\left(6\right)$

其中，f₁和f₂是理想俞氏零相位子波频率参数，时间范围为-∞＜t＜∞。

通过方程(1)----(6)，根据子波频率参数，可以计算出零相位子波序列：

它表示一个简单的、有限长的、单一的地震零相位子波序列{b_i}，其中(N+1)表示子波长度(序列点数)。

(3)对子波进行Z变换；

地震子波序列{b_i}的Z变换为

$B_N\left(z\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{b}_n{z}^{-n}---\left(8\right)$

(4)对子波Z变换进行因式分解并将特征因子按照从小到大排列；

步骤(4)所述的对子波Z变换进行因式分解并将特征因子按照从小到大排列，就是将方程(1)分解成N个相乘的因式，即

$B\left(z\right)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}(1-{\hat{b}}_i{z}^{-1})=(1-{\hat{b}}_1{z}^{-1})(1-{\hat{b}}_2{z}^{-1})\mathrm{\Λ}(1-{\hat{b}}_N{z}^{-1})---\left(9\right)$

且每个因式的特征因子

是按照绝对值从小到大依次排列，即

$\left|{\hat{b}}_1\right|\≤\left|{\hat{b}}_2\right|\≤\mathrm{\Λ}\≤\left|{\hat{b}}_N\right|---\left(10\right)$

(5)将特征因子分成单位圆内和单位圆外两部分；

步骤(5)所述的将特征因子分成单位圆内和单位圆外两部分，就是根据特征因子绝对值的大小，将特征因子分成两部分：第一部分是特征因子的绝对值小于或者等于1，即

$\left|{\hat{b}}_i\right|\≤\mathrm{1,1}\≤i\≤K---\left(4\right)$

共有K(K≥1)个，其中K是Z平面上单位圆内和单位圆上特征因子的总个数。第二部分是特征因子的绝对值大于1，即

$\left|{\hat{b}}_{K+i}\right|>\mathrm{1,1}\≤i\≤M---\left(5\right)$

共有M(M≥1)个，其中M是Z平面上单位圆外特征因子的总个数，且M+K＝N。

(6)构成纯相位转换因子；

根据特征方程单位圆外的特征因子，构成将零相位子波转换为最小相位子波的纯相位转换因子H₁(z)：

$H_1\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}\frac{(1-\frac{1}{{\hat{b}}_{K+i}^{*}}{z}^{-1})}{(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})}---\left(13\right)$

根据特征方程单位圆外的特征因子，构成将最小相位子波转换为零相位子波的纯相位转换因子H₂(z)：

$H_2\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}\frac{(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})}{(1-\frac{1}{{\hat{b}}_{K+i}^{*}}{z}^{-1})}---\left(14\right)$

这里，“*”表示取共轭。

(7)子波纯相位转换；

根据纯相位转换因子H₁(z)，将零相位子波转换为最小相位子波：

W_N(z)＝H₁(z)B_N(z)              (15)

根据纯相位转换因子H₂(z)，将最小相位子波转换为零相位子波：

B_N(z)＝H₂(z)W_N(z)              (16)

这里，W_N(z)表示最小相位子波的Z变换，B_N(z)表示零相位子波的Z变换。

(8)地震数据纯相位转换；

根据纯相位转换因子H₁(z)，将零相位子波地震数据转换为最小相位子波地震数据：

Y_W(z)＝H₁(z)X_B(z)              (17)

这里，Y_W(z)表示最小相位子波地震数据的Z变换，X_B(z)表示零相位子波地震数据的Z变换。

根据纯相位转换因子H₂(z)，将最小相位子波地震数据转换为零相位子波地震数据：

Y_B(z)＝H₂(z)X_W(z)              (18)

这里，Y_B(z)表示零相位子波地震数据的Z变换，X_W(z)表示最小相位子波地震数据的Z变换。

(9)纯相位转换的实际实现；

可以通过时间域的褶积与递归来实现纯相位因子转换。假设零相位地震数据序列为{x_i}，其长度为(L+1)，Z变换为

$X_B\left(z\right)=\underset{i=0}{\overset{L}{\mathrm{\Σ}}}{x}_i{z}^{-i}---\left(19\right)$

根据纯相位转换因子H₁(z)，将零相位子波地震数据转换为最小相位子波地震数据

$Y_W\left(z\right)=H_1\left(z\right)X_B\left(z\right)=\frac{P_M\left(z\right)X_B\left(z\right)}{Q_M\left(z\right)}---\left(20\right)$

这里，

$P_M\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}(1-\frac{1}{{\hat{b}}_{K+i}^{*}})z^{-1}---\left(21\right)$

和

$Q_M\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})---\left(22\right)$

将(21)和(22)展开，有

$P_M\left(z\right)=\underset{i=0}{\overset{M}{\mathrm{\Σ}}}{p}_i{z}^{-i}---\left(23\right)$

和

$Q_M\left(z\right)=\underset{i=0}{\overset{M}{\mathrm{\Σ}}}{q}_i{z}^{-i}---\left(24\right)$

其中，

p₀＝1，

$p_i={(-1)}^i\underset{k_1=1}{\overset{M-i+1}{\mathrm{\Σ}}}\underset{k_2=2}{\overset{M-i+2}{\mathrm{\Σ}}}\mathrm{\Λ}\underset{k_i=i}{\overset{M}{\mathrm{\Σ}}}\frac{1}{{\hat{b}}_{k_1}^{*}{\hat{b}}_{k_2}^{*}\mathrm{\Λ}{\hat{b}}_{k_i}^{*}}---\left(25\right)$

这里，i＝1，2，Λ，M。

q_M＝1，

$q_{M-i}={(-1)}^i\underset{k_1=1}{\overset{M-i+1}{\mathrm{\Σ}}}\underset{k_2=2}{\overset{M-i+2}{\mathrm{\Σ}}}\mathrm{\Λ}\underset{k_i=i}{\overset{M}{\mathrm{\Σ}}}\frac{1}{{\hat{b}}_{k_1}{\hat{b}}_{k_2}\mathrm{\Λ}{\hat{b}}_{k_i}}---\left(26\right)$

这里，i＝1，2，Λ，M。令

$U_{L+M}\left(z\right)=X_B\left(z\right)P_M\left(z\right)=\underset{i=0}{\overset{L+M}{\mathrm{\Σ}}}{u}_i{z}^{-1}---\left(27\right)$

表示时间域一个褶积运算。有

$u_i=\underset{j=0}{\overset{L}{\mathrm{\Σ}}}{x}_j{p}_{i-j}---\left(28\right)$

这里，i＝1，2，Λ，M+L。令

U_L+M(z)＝Y_W(z)Q_M(z)             (29)

表示时间域一个递归运算。有

$y_i=\left\{\begin{array}{cc}\frac{u_0}{q_0}& i=0\\ \frac{1}{q_0}(u_i-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{q}_j{y}_{i-j})& i=\mathrm{1,2},K,M\\ \frac{1}{q_0}(u_i-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{q}_j{y}_{i-j})& i=M+1,M+2,K,M+L\end{array}\right.---\left(30\right)$

这样，通过(28)和(30)式，可以实现纯相位因子转换过程，即将零相位子波地震数据转换为最小相位子波地震数据。如果将方程(19)的零相位地震数据序列换成零相位地震子波序列，同样通过(28)和(30)式，可以实现纯相位因子转换过程，即将零相位地震子波转换为最小相位地震子波。

同理，假设最小相位地震数据序列为{x_i}，其长度为(L+1)，Z变换为

$X_W\left(z\right)=\underset{i=0}{\overset{L}{\mathrm{\Σ}}}{x}_i{z}^{-i}---\left(31\right)$

根据纯相位转换因子H₂(z)，将最小相位子波地震数据转换为零相位子波地震数据

$Y_B\left(z\right)=H_2\left(z\right)X_W\left(z\right)=\frac{Q_M\left(z\right)X_B\left(z\right)}{P_M\left(z\right)}---\left(32\right)$

令

$V_{L+M}\left(z\right)=X_B\left(z\right)Q_M\left(z\right)=\underset{i=0}{\overset{L+M}{\mathrm{\Σ}}}{v}_i{z}^{-i}---\left(33\right)$

表示时间域一个褶积运算。有

$v_i=\underset{j=0}{\overset{L}{\mathrm{\Σ}}}{x}_j{q}_{i-j}---\left(34\right)$

这里，i＝1，2，Λ，M+L。令

V_L+M(z)＝Y_B(z)Q_M(z)           (35)

表示时间域一个递归运算。有

$y_i=\left\{\begin{array}{cc}\frac{v_0}{p_0}& i=0\\ \frac{1}{p_0}(v_i-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{p}_j{y}_{i-j})& i=\mathrm{1,2},K,M\\ \frac{1}{p_0}(v_i-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{p}_j{y}_{i-j})& i=M+1,M+2,K,M+L\end{array}\right.---\left(36\right)$

这样，通过(34)和(36)式，可以实现纯相位因子转换过程，即将最小相位子波地震数据转换为零相位子波地震数据。如果将方程(31)的最小相位地震数据序列换成最小相位地震子波序列，同样通过(34)和(36)式，可以实现纯相位因子转换过程，即将最小相位地震子波转换为零相位地震子波。

10)绘制相位转换后的地震数据剖面和存储相位转换后的地震数据。

本发明进行了相位之间的相位转换处理验证。

图1-1是低通子波相位转换处理对比，其中，(a)是零相位子波，子波频率参数分别是LP(25，35)；(b)是计算得到的最小相位转换因子；(c)是由(b)最小相位转换因子计算得到的最小相位子波。图2-1是低通子波相位转换处理振幅谱对比，显然零相位子波振幅谱(a)与最小相位子波振幅谱(c)完全相同，说明低通子波相位转换处理方法是正确的。图1-2是高通子波相位转换处理对比，其中，(a)是零相位子波，子波频率参数分别是HP(35，40)；(b)是计算得到的最小相位转换因子；(c)是由(b)最小相位转换因子计算得到的最小相位子波。图2-2是高通子波相位转换处理振幅谱对比，显然零相位子波振幅谱(a)与最小相位子波振幅谱(c)完全相同，说明高通子波相位转换处理方法是正确的。图1-3是带通子波相位转换处理对比，其中，(a)是零相位子波，子波频率参数分别是BP(15，20，50，60)；(b)是计算得到的最小相位转换因子；(c)是由(b)最小相位转换因子计算得到的最小相位子波。图2-3是带通子波相位转换处理振幅谱对比，显然零相位子波振幅谱(a)与最小相位子波振幅谱(c)完全相同，说明带通子波相位转换处理方法是正确的。图1-4是陷频子波相位转换处理对比，其中，(a)是零相位子波，子波频率参数分别是NP(35，50，65)；(b)是计算得到的最小相位转换因子；(c)是由(b)最小相位转换因子计算得到的最小相位子波。图2-4是陷频子波相位转换处理振幅谱对比，显然零相位子波振幅谱(a)与最小相位子波振幅谱(c)完全相同，说明陷频子波相位转换处理方法是正确的。图1-5是Ricker子波相位转换处理对比，其中，(a)是零相位子波，子波中心频率参数是RP(25)；(b)是计算得到的最小相位转换因子；(c)是由(b)最小相位转换因子计算得到的最小相位子波。图2-5是Ricker子波相位转换处理振幅谱对比，显然零相位子波振幅谱(a)与最小相位子波振幅谱(c)完全相同，说明Ricker子波相位转换处理方法是正确的。图1-6是俞氏子波相位转换处理对比，其中，(a)是零相位子波，子波中心频率参数是YP(15，50)；(b)是计算得到的最小相位转换因子；(c)是由(b)最小相位转换因子计算得到的最小相位子波。图2-6是俞氏子波相位转换处理振幅谱对比，显然零相位子波振幅谱(a)与最小相位子波振幅谱(c)完全相同，说明俞氏子波相位转换处理方法是正确的。图3是子波反褶积结果对比，期中(a)是带通零相位子波；(b)是(a)的反褶积结果，由于零相位子波并不满足地震数据反褶积处理的最小相位假设条件，因此反褶积后存在很强的剩余子波，并且存在相位移；这也充分说明在地震数据反褶积处理之前，必须对地震数据进行处理以使地震数据子波相位是最小相位。否则不会得到理想的反褶积处理效果。(c)是(a)经过最小相位转换后得到最小相位子波；(d)是(c)的反褶积结果，显然由于最小相位子波满足地震数据反褶积最小相位假设条件，反褶积后将子波压缩为单位脉冲，是地震数据子波反褶积的最理想情况。

本发明对实际地震数据地震子波应用验证。

图4是本发明与OMEGA地震数据系统相位转换处理后叠加数据对比，(a)是零相位子波叠加数据，(b)是OMEGA系统对叠前数据进行最小相位转换处理后的最小相位子波叠加数据，(c)是本发明对叠前数据进行最小相位转换处理后的最小相位子波叠加数据。OMEGA系统对最小相位转换处理算子进行滤波、平滑处理，本方法对转换算子并未做任何处理。

从处理的叠加剖面效果上看，本发明与OMEGA系统处理结果相当，但在局部细节上优于OMEGA系统。图5是本发明最小相位相位转换处理前后叠加数据对比，(a)是零相位子波叠加数据，(b)是本发明最小相位转换处理后的最小相位子波叠加数据。从处理前后的叠加剖面效果上看，总体上本发明相位转换处理前后结果变化不大，但在局部细节上优于相位转换处理前的地震数据。图6是本发明最小相位相位转换处理前后反褶积处理后叠加数据对比，(a)是零相位子波反褶积叠加数据，(b)是本发明最小相位转换处理后的最小相位子波反褶积叠加数据。从处理前后的反褶积叠加剖面效果上看，最小相位转换处理之后，最小相位子波反褶积处理明显优于相位转换处理前反褶积处理的地震数据。

Claims (9)

1.一种地震数据子波相位转换处理的方法，其特征是包括以下步骤：

1)用地震人工震源激发和采集地震数据并做预处理；

2)确定地震数据子波参数和计算地震数据零相位子波；

3)对子波进行Z变换；

4)对子波Z变换按下式分解成N个相乘的因式：

$B\left(z\right)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}(1-{\hat{b}}_i{z}^{-1})=(1-{\hat{b}}_1{z}^{-1})(1-{\hat{b}}_2{z}^{-1})\mathrm{\Λ}(1-{\hat{b}}_N{z}^{-1})---\left(3\right)$

且每个因式的特征因子

按照绝对值从小到大依次排列，

$\left|{\hat{b}}_1\right|\≤\left|{\hat{b}}_2\right|\≤\mathrm{\Λ}\≤\left|{\hat{b}}_N\right|---\left(4\right)$

5)将特征因子分成单位圆内和单位圆外两部分，按特征因子绝对值的大小，将特征因子分成两部分：

第一部分是特征因子的绝对值小于或者等于1，

共有K(K≥1)个，其中K是Z平面上单位圆内和单位圆上特征因子的总个数；

第二部分是特征因子的绝对值大于1，共有M(M≥1)个，其中M是Z平面上单位圆外特征因子的总个数，且M+K＝N；

6)根据单位圆外的特征因子，按下式构成将零相位子波转换为最小相位子波的纯相位转换因子H₁(z)：

$H_1\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}\frac{(1-{\hat{b}}_{K+i}^{*}{z}^{-1})}{(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})}---\left(6\right)$

根据特征方程单位圆外的特征因子，按下式构成将最小相位子波转换为零相位子波的纯相位转换因子H₂(z)：

$H_2\left(z\right)=\underset{i=1}{\overset{M}{\mathrm{\Π}}}\frac{(z^{-1}-\frac{1}{{\hat{b}}_{K+i}})}{(1-\frac{1}{{\hat{b}}_{K+i}^{*}}{z}^{-1})}---\left(7\right)$

式中：“*”表示取共轭；

7)根据最小相位子波的纯相位转换因子H₁(z)，按下式将零相位子波转换为最小相位子波：

W_N(z)＝H₁(z)B_N(z)            (8)

根据零相位子波的纯相位转换因子H₂(z)，按下式将最小相位子波转换为零相位子波：

B_N(z)＝H₂(z)W_N(z)            (9)

式中：W_N(z)表示最小相位子波的Z变换，B_N(z)表示零相位子波的Z变换；

8)根据纯相位转换因子H₁(z)，按下式将零相位子波地震数据转换为最小相位子波地震数据：

Y_W(z)＝H₁(z)X_B(z)            (10)

式中：Y_W(z)表示最小相位子波地震数据的Z变换，X_B(z)表示零相位子波地震数据的Z变换；

根据零相位子波的纯相位转换因子H₂(z)，按下式将最小相位子波地震数据转换为零相位子波地震数据：

Y_B(z)＝H₂(z)X_W(z)            (11)

式中：Y_B(z)表示零相位子波地震数据的Z变换，X_W(z)表示最小相位子波地震数据的Z变换；

9)根据纯相位转换因子实现地震子波和地震数据的纯相位转换，绘制相位转换后的地震数据剖面和存储相位转换后的地震数据。

2.根据权利要求1的方法，步骤1)所述的预处理是指对地震数据置标签、定义观测系统。

3.根据权利要求1的方法，步骤2)所述的确定地震数据子波参数，是指对地震数据进行振幅谱分析，根据地震数据振幅谱形状，确定原始数据子波的参数。然后根据地震数据子波参数，计算地震数据零相位子波。

4.根据权利要求1的方法，步骤3)所述的对子波进行Z变换是根据地震数据零相位子波，按照下式计算：

$B\left(z\right)=\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{b}_n{z}^{-n}---\left(1\right)$

式中，

是单一的地震零相位子波序列{b_i}，其中(N+1)表示子波长度(序列点数)。

5.根据权利要求1的方法，步骤9)所述的纯相位转换是在时间域内，通过递推和递归运算实现纯相位转换。

6.根据权利要求1的方法，步骤9)所述的纯相位转换是将零相位子波地震数据转换为最小相位子波地震数据：

$y_i=\left\{\begin{array}{cc}\frac{u_0}{q_0}& i=0\\ \frac{1}{q_0}(u_i-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{q}_j{y}_{i-j})& i=\mathrm{1,2},K,M\\ \frac{1}{q_0}(u_i-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{q}_j{y}_{i-j})& i=M+1,M+2,K,M+L\end{array}\right.---\left(12\right)$

式中：序列{y_i}为最小相位地震数据或最小相位地震子波，长度为(L+M+1)，且

$u_i=\underset{j=0}{\overset{L}{\mathrm{\Σ}}}{x}_j{p}_{i-j}---\left(13\right)$

$p_0=1,p_i={(-1)}^i\underset{k_1=1}{\overset{M-i+1}{\mathrm{\Σ}}}\underset{k_2=2}{\overset{N-i+2}{\mathrm{\Σ}}}\mathrm{\Λ}\underset{k_i=i}{\overset{M}{\mathrm{\Σ}}}\frac{1}{{\hat{b}}_{k_1}^{*}{\hat{b}}_{k_2}^{*}\mathrm{\Λ}{\hat{b}}_{k_i}^{*}}---\left(14\right)$

$p_0=1,p_i={(-1)}^i\underset{k_1=1}{\overset{M-i+1}{\mathrm{\Σ}}}\underset{k_2=2}{\overset{N-i+2}{\mathrm{\Σ}}}\mathrm{\Λ}\underset{k_i=i}{\overset{M}{\mathrm{\Σ}}}\frac{1}{{\hat{b}}_{k_1}^{*}{\hat{b}}_{k_2}^{*}\mathrm{\Λ}{\hat{b}}_{k_i}^{*}}---\left(14\right)$

式中：i＝1，2，Λ，M，

式中：序列{x_i}为零相位地震数据或零相位地震子波，长度为(L+1)。

7.根据权利要求1或6的方法，在以上步骤中序列{x_i}为零相位地震数据，则式(12)中相位转换后的序列{y_i}为最小相位地震数据；如果在式(13)中，序列{x_i}为零相位地震子波，则方程(12)中相位转换后的序列{y_i}为最小相位地震子波。

8.根据权利要求1的方法，步骤9)所述的纯相位转换是将最小相位子波地震数据转换为零相位子波地震数据：

$y_i=\left\{\begin{array}{cc}\frac{v_0}{p_0}& i=0\\ \frac{1}{p_0}(v_i-\underset{j=1}{\overset{i}{\mathrm{\Σ}}}{p}_j{y}_{i-j})& i=\mathrm{1,2},K,M\\ \frac{1}{p_0}(v_i-\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{p}_j{y}_{i-j})& i=M+1,M+2,K,M+L\end{array}\right.---\left(16\right)$

式中：序列{y_i}为零相位地震数据或零相位地震子波，长度为(L+M+1)，且

$v_i=\underset{j=0}{\overset{L}{\mathrm{\Σ}}}{x}_j{q}_{i-j}---\left(17\right)$

式中：i＝1，2，A，M，

9.根据权利要求1或8的方法，在以上步骤中序列{x_i}为最小相位地震数据或最小相位地震子波，其长度为(L+1)。在式(17)中，序列{x_i}为最小相位地震数据，则式(16)中相位转换后的序列{y_i}为零相位地震数据；如果在式(17)中，序列{x_i}为最小相位地震子波，则式(16)中相位转换后的序列{y_i}为零相位地震子波。

Images