CN103217716B - 一种地震数据快速相位校正处理的方法 - Google Patents

Abstract

本发明的一种地震数据快速相位校正处理的方法，使用凯泽窗希尔伯特因子序列，确定希尔伯特地震道序列，通过估算剩余子波相位信息，直接计算确定出最佳相位校正因子，并对地震数据记录进行相位校正处理，达到消除剩余子波相位影响，使剩余子波相位达到或者接近零相位，以有效提高地震数据分辨率，本发明仅对地震数据子波的相位进行处理，而不改变数据的振幅谱，具有计算量小、计算速度快、稳定性好和相位校正处理计算精度高的特点，且具有一定的抗噪能力。

Description

一种地震数据快速相位校正处理的方法

技术领域

本发明涉及油田的勘探、开发、开采技术，具体是为反映地下地层层位、油藏描述提供高分辨率的地震图形和数据的一种地震数据快速相位校正处理的方法，特别适用于地震数据处理反褶积之后的实际地震数据，通过剩余子波相位校正处理，达到消除地震数据剩余子波相位影响、提高地震数据分辨率目的。

背景技术

随着地震勘探技术的发展，油气勘探的难度和深度也越来越大，对地震资料分辨率的要求也越来越高。反褶积处理是提高地震数据分辨率最有效的途径。地震褶积模型是地震勘探数据处理中最基本的模型之一。常规的地震数据反褶积处理都是基于地震褶积模型。反褶积处理的一个根本假设就是地震子波是最小相位的。这个假设使得我们可以由子波振幅谱得到子波的相位谱并继而得到完整的子波。而实际数据的地震子波是混合相位的，这便使得以最小相位子波假设为前提的反褶积处理后地震子波不是一个脉冲，而存在剩余子波。实际地震数据处理时，即使将混合相位地震子波转换为最小相位地震子波，由于子波有限长度，反褶积处理后地震子波也还会存在剩余子波。地震数据中剩余子波的存在，降低了地震数据的分辨率。

在地震数据子波中，零相位子波具有最大分辨率，因此实际数据处理为了提高分辨率，希望处理后的地震子波是零相位的。“相同频带范围的子波中，以零相位子波的分辨率为最高”(李庆忠，走向精确勘探的道路，P14，石油工业出版社，1994)。“具有相同振幅谱的诸子波中，零相位子波的分辨率最高”(俞寿朋，高分辨率地震勘探，P17，石油工业出版社，1993)。“子波相位校正的目的是使子波零相位化”(俞寿朋，俞寿朋文集，P206，石油工业出版社，2001)。“子波波形由振幅谱和相位谱二者决定。振幅很小的频率成分对子波波形的影响很小，也就是说这些频率成分的相位谱重要性不大，有些误差也关系不大。而振幅谱比较大的频率成分是起主要作用的，这些频率成分的相位谱是重要的。但这个频带的相位谱一般很接近于直线。通过相位谱坐标原点做一条平行于主频带相位谱的直线，则主频带相位谱与此直线的差近于常数。因此可近似认为子波是常相位的。”(俞寿朋，高分辨率地震勘探，P167，石油工业出版社，1993)。

常规的消除地震数据剩余子波相位处理，采用相位校正方法。相位校正就是将剩余子波相位作为一个常数，显然这样的假设仅仅是一阶近似。当子波的振幅谱较窄时，此法会有一定效果。在相位校正处理方法中，通常相位校正因子的估算采用对目标函数方差模进行扫描的方法来确定(Levy S.and OldenburgD.W.，Automatic phase correction of common-midpoint stacked data，Geophysics，1987，VOL.52，NO.1，P51-P59；周兴元，常相位校正，石油地球物理勘探，1989，VOL.24，NO.2；李合群、周兴元，时差、常相位校正及加权叠加，石油地球物理勘探，2000，VOL.35，NO.4，P415-418)。扫描方法估算相位校正因子非常费时，计算效率低。

发明内容

本发明目的在于提供一种不需要进行扫描，快速消除地震数据剩余子波相位影响，提高地震数据分辨率的地震数据快速相位校正处理的方法。

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

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

步骤1)所述的预处理是指对地震数据置标签、定义观测系统、速度分析、动校正、反褶积、叠加处理。

2)计算凯泽窗希尔伯特因子序列；

步骤2)所述的计算凯泽窗希尔伯特因子序列包括计算第一类零阶修正贝赛尔函数和凯泽窗希尔伯特因子序列：

按照以下公式计算第一类零阶修正贝赛尔函数I₀(x)：

I₀(x)＝a₀+a₁+a₂+...+a_k       (1)

式中，I₀(x)为第一类零阶修正贝赛尔函数；x为第一类零阶修正贝赛尔函数自变量，其取值范围为0≤x≤1；a₀是第一类零阶修正贝赛尔函数初始值，且a₀＝1；k表示第一类零阶修正贝赛尔函数的第k项；a_k是第一类零阶修正贝赛尔函数第k项值；

按照以下公式计算第一类零阶修正贝赛尔函数第k项值a_k：

$a_k=\frac{1}{k^2}{\left(\frac{x}{2}\right)}^2{a}_{k-1},k=\mathrm{1,2,3},...---\left(2\right)$

式中，x为第一类零阶修正贝赛尔函数自变量，其取值范围为0≤x≤1；a₀是第一类零阶修正贝赛尔函数初始值，a_k是第一类零阶修正贝赛尔函数第k项值；

在第一类零阶修正贝赛尔函数中，变量x的取值范围为0≤x≤1，在有限个k值之后，a_k之值越来越小并且无限接近于零，因此可以取a_k＜10^-8作为k的最终结束标志；

按照以下公式计算凯泽窗希尔伯特因子序列h_w[n]：

$h_w\left[n\right]=\left\{\begin{array}{cc}\frac{I_0\left[\mathrm{\&beta;}{(1-{[(n-\mathrm{\&alpha;})/\mathrm{\&alpha;}]}^2)}^{1/2}\right]}{I_0\left(\mathrm{\&beta;}\right)}\left\{\frac{{\sin }^2[\mathrm{\&pi;}(n-\mathrm{\&alpha;})/2}{\mathrm{\&pi;}(n-\mathrm{\&alpha;})/2}\right\}& 0\&le;n\&le;M\\ 0& n<0,n>M\end{array}\right.---\left(3\right)$

式中，I₀(·)为第一类零阶修正贝赛尔函数；β为凯泽窗函数参数，M为凯泽窗函数样点个数参数，且α＝M/2，n为凯泽窗希尔伯特因子序列样点顺序号，n＝1，2，3，...，M。

3)计算希尔伯特地震道序列；

按照以下公式计算地震数据对应的希尔伯特地震道序列x_H[n]：

x_H[n]＝x[n]*h_w[n]       (4)

式中，x[n]为地震数据序列，x_H[n]为希尔伯特地震道序列，h_w[n]为凯泽窗希尔伯特因子序列，由方程(3)计算，符号“*”表示褶积运算；n为地震数据样点顺序号，n＝1，2，3，...，N，N为地震数据样点序列个数；

4)确定相位校正因子特征方程的系数值；

按照以下公式分别计算相位校正因子特征方程的系数值A、B、C、D和E：

A＝Pb

B＝2(Pc-2Pa)

C＝3(Pd-Pb)

D＝2(2Pe-Pc)

E＝-Pd            (5)

式中，Pa、Pb、Pc、Pd和Pe是五个常数值，仅仅与地震数据序列x[n]和希尔伯特地震道序列x_H[n]有关；

按照以下公式计算Pa、Pb、Pc、Pd和Pe五个常数值：

$\mathrm{Pa}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^4\left[i\right]-3{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\right[i\left]\right)}^2$

$\mathrm{Pb}=12\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H\left[i\right]-4\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^3\left[i\right]x_H\left[i\right]$

$\mathrm{Pc}=6\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]x_H^2\left[i\right]-6\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\left[i\right]-12{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\right[i\left]x_H\right[i\left]\right)}^2$

$\mathrm{Pd}=12\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H\left[i\right]-4\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H^3\left[i\right]$

$\mathrm{Pe}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^4\left[i\right]-3{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\right[i\left]\right)}^2---\left(6\right)$

式中，x[i]为地震数据序列，x_H[i]为地震数据序列对应的希尔伯特地震道序列，由方程(4)计算确定，i为地震数据样点顺序号，N表示地震数据样点个数。

5)计算相位校正因子特征方程的根；

按照以下公式构造相位校正因子特征方程：

A+Bx+Cx²+Dx³+Ex⁴＝0      (7)

式中，A、B、C、D和E为根据相位校正因子特征方程的五个系数值。x表示相位校正因子特征方程自变量，是相位校正因子的正切函数值。求解方程(7)，可以得到相位校正因子特征方程的对应四个根x₁、x₂、x₃和x₄；

按照以下公式计算相位校正因子φ：

φ＝arctgx，x＝tgφ     (8)

式中，x为相位校正因子特征方程自变量的根。

根据相位校正因子特征方程对应的四个根x₁、x₂、x₃和x₄，由方程(8)可以得到四个相位校正因子φ₁、φ₂、φ₃和φ₄；

6)计算相位校正因子峰度函数值并确定最佳相位校正因子；

按照以下公式计算峰度函数值kurt(φ)

kurt(φ)＝E{y⁴[n]}-3(E{y²[n]})²     (9)

式中，E{·}表示期望值运算符，y[n]是地震数据相位校正序列。已知地震数据序列x[n]，相位校正因子φ，由方程(11)计算地震数据相位校正序列。

按照以下公式确定最佳相位校正因子φ_best：

${\mathrm{\&phi;}}_{\mathrm{best}}=\underset{{\mathrm{\&phi;}}_{\mathrm{best}}\&Element;[{\mathrm{\&phi;}}_1,{\mathrm{\&phi;}}_2,{\mathrm{\&phi;}}_3,{\mathrm{\&phi;}}_4]}{\mathrm{Max}}\{\mathrm{kurt}\left({\mathrm{\&phi;}}_1\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_2\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_3\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_4\right)\}---\left(10\right)$

式中，kurt(φ₁)、kurt(φ₂)、kurt(φ₃)和kurt(φ₄)分别是相位校正因子特征方程计算的四个相位校正因子φ₁、φ₂、φ₃和φ₄对应的峰度函数值，其中四个峰度函数值中最大值对应的相位校正因子就是最佳相位校正因子；

7)实现相位校正；

y[n]＝x[n]cosφ_best-x_H[n]sinφ_best    (11)

式中，φ_best是最佳相位校正因子，由方程(10)确定，单位是弧度，x[n]是地震数据序列，x_H[n]是希尔伯特地震道序列，由方程(4)计算确定，y[n]是地震数据相位校正序列，n是地震数据样点顺序号，n＝1，2，3，...，N，N是地震数据样点个数；

由步骤2)----步骤6)确定出最佳相位校正因子φ_best，根据相位校正因子φ_best，由方程(11)进行相位校正处理；

8)实现串联快速相位校正；

重复步骤2)----步骤7)，就实现了串联快速相位校正因子确定和快速相位校正处理，以实现消除地震数据剩余子波的快速相位校正处理；

9)绘制消除剩余子波后的地震数据剖面和存储消除剩余子波后的地震数据。

本发明不是采用相位扫描方法确定最佳相位校正因子，而是使用凯泽窗希尔伯特因子序列，确定希尔伯特地震道序列，直接计算确定出最佳相位校正因子φ。本发明仅对地震数据子波的相位进行处理，而不改变数据的振幅谱，具有计算量小、计算速度快、稳定性好和相位校正处理计算精度高的特点，且具有一定的抗噪能力。

附图说明

图1第一类零阶修正Bessel函数；

图2不同窗参数β值时凯泽窗函数及频谱对比；

(a)窗参数，

(b)振幅谱，

(c)对数谱；

图3不同窗参数M值时凯泽窗函数及频谱对比；

(a)窗参数，

(b)振幅谱，

(c)对数谱；

图4不同子波不同相位校正因子校正处理结果对比；

(a)零相位子波，

(b)最小相位子波；

图5不同子波不同相位校正因子校正处理结果频谱对比；

(a)零相位子波，

(b)最小相位子波；

图6CMP305-CMP685反褶积处理后数据相位校正结果对比；

(a)反褶积处理后数据，

(b)现有相位校正，

(c)本发明；

图7CMP705-CMP1085反褶积处理后数据相位校正结果对比；

(a)反褶积处理后数据，

(b)现有相位校正，

(c)本发明；

图8反褶积处理后数据相位校正结果对比局部放大显示；

(a)反褶积处理后数据，

(b)现有相位校正，

(c)本发明；

图9CMP305-CMP685蓝色滤波处理后数据相位校正结果对比；

(a)反褶积处理后数据，

(b)现有相位校正，

(c)本发明；

图10CMP705-CMP1085蓝色滤波处理后数据相位校正结果对比；

(a)反褶积处理后数据，

(b)现有相位校正，

(c)本发明；

图11蓝色滤波处理后数据相位校正结果对比局部放大显示；

(a)反褶积处理后数据，

(b)现有相位校正，

(c)本发明。

具体实施方式

本发明的一种地震数据快速相位校正处理的方法，就是通过估算剩余子波相位信息，并对地震数据记录进行相位校正处理，达到消除剩余子波相位影响，使剩余子波相位达到或者接近零相位，以有效提高地震数据分辨率。

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

(1)用地震人工震源激发和采集地震数据并对地震数据置标签、定义观测系统、速度分析、动校正、反褶积、叠加等预处理工作；

(2)计算凯泽窗希尔伯特因子序列；

计算凯泽窗希尔伯特因子包括计算凯泽(Kaiser)窗函数、希尔伯特因子序列和凯泽窗希尔伯特因子序列。凯泽(Kaiser)窗函数就是由第一类零阶修正贝赛尔函数构成的最佳窗函数。

第一类零阶修正贝赛尔函数I₀(x)表达式为

$I_0\left(x\right)=\underset{k=0}{\overset{\&infin;}{\mathrm{\&Sigma;}}}\frac{1}{{(k!)}^2}{\left(\frac{x}{2}\right)}^{2k}---\left(1\right)$

式中，I₀(x)为第一类零阶修正贝赛尔函数；k表示第一类零阶修正贝赛尔函数的第k项；x为第一类零阶修正贝赛尔函数自变量，其取值范围为0≤x≤1。第一类零阶修正贝塞尔函数I₀(x)是指数增长型的，函数值在x＝0点趋于有限值1。将方程(1)展开，就构成了计算第一类零阶修正贝赛尔函数I₀(x)的算法：

按照以下公式计算第一类零阶修正贝赛尔函数I₀(x)：

I₀(x)＝a₀+a₁+a₂+...+a_k     (2)

式中，I₀(x)为第一类零阶修正贝赛尔函数；x为第一类零阶修正贝赛尔函数自变量，其取值范围为0≤x≤1；a₀是第一类零阶修正贝赛尔函数初始值，且a₀＝1；k表示第一类零阶修正贝赛尔函数的第k项；a_k是第一类零阶修正贝赛尔函数第k项值；

按照以下公式计算第一类零阶修正贝赛尔函数第k项值a_k：

$\begin{array}{cc}{a}_k=\frac{1}{k^2}{\left(\frac{x}{2}\right)}^2{a}_{k-1},& k=\mathrm{1,2,3}...\end{array}---\left(3\right)$

式中，x为第一类零阶修正贝赛尔函数的自变量，其取值范围为0≤x≤1；a₀是第一类零阶修正贝赛尔函数初始值，a_k是第一类零阶修正贝赛尔函数第k项值；

在第一类零阶修正贝赛尔函数中，变量x的取值范围为0≤x≤1，在有限个k值之后，a_k之值越来越小并且无限接近于零，因此可以取a_k＜10^-8作为k的最终结束标志；

按照以下公式计算凯泽窗函数w[n]：

$w\left[n\right]=\left\{\begin{array}{cc}\frac{I_0\left[\mathrm{\&beta;}{(1-{[(n-\mathrm{\&alpha;})/\mathrm{\&alpha;}]}^2)}^{1/2}\right]}{I_0\left(\mathrm{\&beta;}\right)}& 0\&le;n\&le;M\\ 0& n<0,n>M\end{array}\right.---\left(4\right)$

式中，I₀(·)为第一类零阶修正贝赛尔函数；β为凯泽窗函数参数，M为凯泽窗函数样点个数参数，且α＝M/2，n为凯泽窗希尔伯特因子序列样点顺序号，n＝1，2，3，...，M；

实现希尔伯特变换时，需要知道希尔伯特因子序列。时间域希尔伯特因子序列h[n]表达式是

$h\left[n\right]=\left\{\begin{array}{cc}\frac{{2\sin }^2(\mathrm{\&pi;n}/2)}{\mathrm{\&pi;n}}& n\&NotEqual;0\\ 0& n=0\end{array}\right.---\left(5\right)$

其希尔伯特因子序列的频率域表达式为

$H\left(e^{\mathrm{j\&omega;}}\right)=\begin{array}{}\end{array}\left\{\begin{array}{cc}-j& 0<\mathrm{\&omega;}<\mathrm{\&pi;}\\ j& -\mathrm{\&pi;}<\mathrm{\&omega;}<0\end{array}\right.---\left(6\right)$

式中，H(e^jω)表示希尔伯特因子序列的频谱，ω表示圆频率。显然时间域希尔伯特因子序列h[n]的长度是无穷的，因此实现时必须施以窗函数进行截断。常用矩形窗进行截断处理。矩形窗的旁瓣幅度大，会造成相位校正精度低，引起相位校正处理误差。

按照以下公式计算凯泽窗希尔伯特因子序列h_w[n]：

$h_w\left[n\right]=\left\{\begin{array}{cc}\frac{I_0\left[\mathrm{\&beta;}{(1-{[(n-\mathrm{\&alpha;})/\mathrm{\&alpha;}]}^2)}^{1/2}\right]}{I_0\left(\mathrm{\&beta;}\right)}\left\{\frac{{\sin }^2[\mathrm{\&pi;}(n-\mathrm{\&alpha;})/2}{\mathrm{\&pi;}(n-\mathrm{\&alpha;})/2}\right\}& 0\&le;n\&le;M\\ 0& n<0,n>M\end{array}\right.---\left(7\right)$

式中，I₀(·)为第一类零阶修正贝赛尔函数；β为凯泽窗函数参数，M为凯泽窗函数样点个数参数，且α＝M/2，n为凯泽窗希尔伯特因子序列样点顺序号，n＝1，2，3，...，M；

凯泽窗希尔伯特因子序列由两个参数β和M来描述。β参数与相对旁瓣幅度有关。相对旁瓣幅度A_sl定义为以dB计的主瓣幅度A_m与最大旁瓣幅度A_p之比，即

$A_{\text{sl}}\text{=}{20\log }_{10}\frac{A_m}{A_p}---\left(8\right)$

相对旁瓣幅度A_sl基本与窗长度参数M无关，只取决于参数β。凯泽通过大量研究，得到相对旁瓣幅度A_sl与参数β的关系：

$\mathrm{\&beta;}=\left\{\begin{array}{cc}0& A_{\mathrm{sl}}\&le;13.26\\ 0.76609{(A_{\mathrm{sl}}-13.26)}^{0.4}+0.09834(A_{\mathrm{sl}}-13.26)& 13.26<A_{\mathrm{sl}}<60\\ 0.12438(A_{\mathrm{sl}}+6.3)& 60<A_{\mathrm{sl}}<120\end{array}\right.---\left(9\right)$

当β＝0时，凯泽窗函数成为矩形窗函数，而矩形窗函数A_sl＝13.26。而主瓣宽度与窗长度M成反比。主瓣宽度Δ_ml定义为在中央的两个过零交点之间的对称距离。主瓣宽度Δ_ml、相对旁瓣幅度A_sl和窗长度M之间的近似关系为

$M=\frac{24\mathrm{\&pi;}(A_{\mathrm{sl}}+12)}{{155\mathrm{\&Delta;}}_{\mathrm{ml}}}---\left(10\right)$

加窗会降低频谱的分辨率，并且频谱会发生泄漏现象。对于凯泽窗希尔伯特因子，可以通过选择适当的参数β和M，以尽量克服和减少这两种影响。

对于凯泽窗希尔伯特因子序列h_w[n]，显然有

h_w[n]＝-h_w[M-n]，0≤n≤M    (11)

因此，凯泽窗希尔伯特因子序列是一个因果广义线形相位系统，其傅立叶变换为

$\begin{array}{cc}{H}_w\left[k\right]=\underset{n=0}{\overset{M}{\mathrm{\&Sigma;}}}{h}_w\left[n\right]e^{-j2\mathrm{\&pi;kn}/N},& 0\&le;k\&le;N-1\end{array}---\left(12\right)$

N表示傅立叶变换点数，且N≥M。经过简单运算，有

H_w[k]＝A[k]e^jφ[k]              (13)

其中，相位φ[k]为

$\mathrm{\&phi;}\left[k\right]=\frac{\mathrm{\&pi;}}{2}-\frac{\mathrm{\&pi;M}}{N}k---\left(14\right)$

显然是一个线性相位。振幅A[k]为

$A\left[k\right]=\underset{n=1}{\overset{M/2}{\mathrm{\&Sigma;}}}{2h}_w[M/2-n]\sin \frac{2\mathrm{\&pi;n}}{N}k---\left(15\right)$

其中，M为偶数。或者

$A\left[k\right]=\underset{n=1}{\overset{(M+1)/2}{\mathrm{\&Sigma;}}}{2h}_w[(M+1)/2-n]\sin \frac{2\mathrm{\&pi;n}}{N}(k-\frac{1}{2})---\left(16\right)$

其中，M为奇数。

(3)计算希尔伯特地震道序列；

按照以下公式计算地震数据对应的希尔伯特地震道序列x_H[n]：

x_H[n]＝x[n]*h_w[n]     (17)

式中，x[n]为地震数据序列，x_H[n]为希尔伯特地震道序列，h_w[n]为凯泽窗希尔伯特因子序列，由方程(7)计算，符号“*”表示褶积运算；n为地震数据样点顺序号，n＝1，2，3，...，N，N为地震数据样点序列个数；

(4)确定相位校正因子特征方程的系数值；

地震数据序列x[n]的相位校正为

y[n]＝x[n]cosφ-x_H[n]sinφ      (18)

这里，φ是相位校正因子，单位是弧度，x[n]是地震数据序列，x_H[n]是希尔伯特地震道序列，y[n]是地震数据相位校正序列。通过方程(18)，很容易进行相位校正处理。

相位校正处理的关键是确定相位校正因子，常规的相位校正因子是采用扫描方法确定。方差模为目标函数：

$\mathrm{kurt}\left(\mathrm{\&phi;}\right)=\frac{E\left\{y^4\right[n\left]\right\}}{{\left(E\right\{y^2\left[n\right]\left\}\right)}^2}---\left(19\right)$

这里，E{·}表示期望运算。

为了更有效确定相位校正因子φ，使用峰度函数(kurtosis function)kurt(φ)作为目标函数，来确定相位校正因子φ。kurt(φ)定义为

kurt(φ)＝E{y⁴[n]}-3(E{y²[n]})²      (20)

将(18)代入(20)，并简化，有

kurt(φ)＝Pa cos⁴φ+Pb cos³φsinφ+Pc cos²φsin²φ+Pd cosφsin³φ+Pe sin⁴φ  (21)

其中，

$\mathrm{Pa}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^4\left[i\right]-3{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\right[i\left]\right)}^2$

$\mathrm{Pb}=12\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H\left[i\right]-4\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^3\left[i\right]x_H\left[i\right]$

$\mathrm{Pc}=6\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]x_H^2\left[i\right]-6\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\left[i\right]-12{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\right[i\left]x_H\right[i\left]\right)}^2$

$\mathrm{Pd}=12\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H\left[i\right]-4\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H^3\left[i\right]$

$\mathrm{Pe}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^4\left[i\right]-3{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\right[i\left]\right)}^2---\left(22\right)$

式中，x[i]为地震数据序列，x_H[i]为地震数据序列对应的希尔伯特地震道序列，由方程(17)计算确定，i为地震数据样点顺序号，N表示地震数据样点个数。Pa、Pb、Pc、Pd和Pe是五个常数值，仅仅与地震数据序列x[n]和希尔伯特地震道序列x_H[n]有关。方程(21)对φ求导数，并令其导数为零，经过简单推导，并令

A＝Pb

B＝2(Pc-2Pa)

C＝3(Pd-Pb)

D＝2(2Pe-Pc)

E＝-Pd                   (23)

这里，A、B、C、D和E就是相位校正因子特征方程的系数值。

使用方程(22)计算Pa、Pb、Pc、Pd和Pe五个常数值，使用方程(23)计算相位校正因子特征方程的系数值A、B、C、D和E。

(5)构成相位校正因子特征方程并求解特征方程根；

按照以下公式构造相位校正因子特征方程：

A+Bx+Cx²+Dx³+Ex⁴＝0       (24)

式中，A、B、C、D和E为根据相位校正因子特征方程的五个系数值。x表示相位校正因子特征方程自变量，是相位校正因子的正切函数值。求解方程(24)，可以得到相位校正因子特征方程的对应四个根x₁、x₂、x₃和x₄；

按照以下公式计算相位校正因子φ：

φ＝arctgx，x＝tgφ            (25)

式中，x为相位校正因子特征方程自变量的根。

根据相位校正因子特征方程对应的四个根x₁、x₂、x₃和x₄，由方程(25)可以得到四个相位校正因子φ₁、φ₂、φ₃和φ₄；

(6)计算相位峰度函数值并确定最佳相位校正因子；

按照以下公式计算峰度函数值kurt(φ)

kurt(φ)＝E{y⁴[n]}-3(E{y²[n]})²    (26)

式中，E{·}表示期望值运算符，y[n]是地震数据相位校正序列。已知地震数据序列x[n]，相位校正因子φ，由方程(18)计算地震数据相位校正序列。

按照以下公式确定最佳相位校正因子φ_best：

${\mathrm{\&phi;}}_{\mathrm{best}}=\underset{{\mathrm{\&phi;}}_{\mathrm{best}}\&Element;[{\mathrm{\&phi;}}_1,{\mathrm{\&phi;}}_2,{\mathrm{\&phi;}}_3,{\mathrm{\&phi;}}_4]}{\mathrm{Max}}\{\mathrm{kurt}\left({\mathrm{\&phi;}}_1\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_2\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_3\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_4\right)\}---\left(27\right)$

式中，kurt(φ₁)、kurt(φ₂)、kurt(φ₃)和kurt(φ₄)分别是相位校正因子特征方程计算的四个相位校正因子φ₁、φ₂、φ₃和φ₄对应的峰度函数值，其中四个峰度函数值中最大值对应的相位校正因子就是最佳相位校正因子；

(7)实现相位校正；

y[n]＝x[n]cosφ_best-x_H[n]sinφ_best    (28)

式中，φ_best是最佳相位校正因子，由方程(27)确定，单位是弧度，x[n]是地震数据序列，x_H[n]是希尔伯特地震道序列，由方程(4)计算确定，y[n]是地震数据相位校正序列，n是地震数据样点顺序号，n＝1，2，3，...，N，N是地震数据样点个数；

由步骤2)----步骤6)确定出最佳相位校正因子φ_best，根据相位校正因子φ_best，由方程(28)进行相位校正处理；

(8)实现串联快速相位校正；

重复步骤(2)----步骤(7)，就实现了串联快速相位校正因子确定和快速相位校正处理，以实现消除地震数据剩余子波的快速相位校正处理。

(9)绘制消除剩余子波后的地震数据剖面和存储消除剩余子波后的地震数据。

本发明进行了相位校正处理验证。

图1是第一类零阶修正Bessel函数，图2是不同窗参数β值时凯泽窗函数及频谱对比，图3是不同参数M值时凯泽窗函数及频谱对比。图4是已知零相位子波和最小相位子波不同相位校正因子校正处理结果对比，图5是已知零相位子波和最小相位子波不同相位校正因子校正处理结果频谱对比，其中，(a)是零相位子波，(b)是最小相位子波。

反褶积和蓝色滤波是地震数据处理中两种最有效的提高分辨率处理技术，因此我们使用反褶积和蓝色滤波处理后的数据进行相位校正处理以说明本发明的有效性。图6、图7和图8是反褶积后数据相位校正结果对比，其中，(a)是反褶积处理后数据显示；(b)是现有相位校正结果；(c)是本发明结果。图6和图7是同一条地震测线不同部位CMP剖面，图6是CMP305-CMP685，图7是CMP705-CMP1085，图8是图7局部放大显示。显然本发明的快速相位校正结果整体效果以及局部细节要好于现有相位校正结果。

图9、图10和图11是蓝色滤波后数据相位校正结果对比，其中，(a)是蓝色滤波处理后数据显示；(b)是现有相位校正结果；(c)是本发明结果。图9和图10是同一条地震测线不同部位CMP剖面，图9是CMP305-CMP685，图10是CMP705-CMP1085，图11是图10局部放大显示。显然本发明快速相位校正结果整体效果以及局部细节要好于现有相位校正结果。

Claims (4)

1.一种地震数据快速相位校正处理的方法，特点是采用以下步骤：

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

2)计算凯泽窗希尔伯特因子序列；

3)计算希尔伯特地震道序列；

4)确定相位校正因子特征方程的系数值；

按照以下公式分别计算相位校正因子特征方程的系数值A、B、C、D和E：

A＝Pb

B＝2(Pc-2Pa)

C＝3(Pd-Pb)

D＝2(2Pe-Pc)

E＝-Pd      (5)

式中，Pa、Pb、Pc、Pd和Pe是五个常数值，仅仅与地震数据序列x[n]和希尔伯特地震道序列x_H[n]有关；

按照以下公式计算Pa、Pb、Pc、Pd和Pe五个常数值：

$\mathrm{Pa}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^4\left[i\right]-3{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\right[i\left]\right)}^2$

$\mathrm{Pb}=12\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H\left[i\right]-4\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^3\left[i\right]x_H\left[i\right]$

$\mathrm{Pc}=6\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]x_H^2\left[i\right]-6\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\left[i\right]-12{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\right[i\left]x_H\right[i\left]\right)}^2$

$\mathrm{Pd}=12\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\left[i\right]\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H\left[i\right]-4\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}x\left[i\right]x_H^3\left[i\right]$

$\mathrm{Pe}=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^4\left[i\right]-3{\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{x}_H^2\right[i\left]\right)}^2---\left(6\right)$

式中，x[i]为地震数据序列，x_H[i]为地震数据序列对应的希尔伯特地震道序列，由方程(4)计算确定，i为地震数据样点顺序号，N表示地震数据样点个数；

5)计算相位校正因子特征方程的根；

按照以下公式构造相位校正因子特征方程：

A+Bx+Cx²+Dx³+Ex⁴＝0      (7)

式中，A、B、C、D和E为相位校正因子特征方程的五个系数值；x表示相位校正因子特征方程自变量，是相位校正因子的正切函数值；求解方程(7)，可以得到相位校正因子特征方程的对应四个根x₁、x₂、x₃和x₄；

按照以下公式计算相位校正因子φ：

φ＝arctgx，x＝tgφ      (8)

式中，x为相位校正因子特征方程自变量的根；

根据相位校正因子特征方程对应的四个根x₁、x₂、x₃和x₄，由方程(8)可以得到四个相位校正因子φ₁、φ₂、φ₃和φ₄；

6)计算相位校正因子峰度函数值并确定最佳相位校正因子；

按照以下公式计算峰度函数值kurt(φ)

kurt(φ)＝E{y⁴[n]}-3(E{y²[n]})²      (9)

式中，E{·}表示期望值运算符，y[n]是地震数据相位校正序列,y[n]＝x[n]cosφ-x_H[n]sinφ；已知地震数据序列x[n]，相位校正因子φ；

按照以下公式确定最佳相位校正因子φ_best：

${\mathrm{\&phi;}}_{\mathrm{best}}=\underset{{\mathrm{\&phi;}}_{\mathrm{best}}\&Element;[{\mathrm{\&phi;}}_1,{\mathrm{\&phi;}}_2,{\mathrm{\&phi;}}_3,{\mathrm{\&phi;}}_4]}{\mathrm{Max}}\{\mathrm{kurt}\left({\mathrm{\&phi;}}_1\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_2\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_3\right),\mathrm{kurt}\left({\mathrm{\&phi;}}_4\right)\}---\left(10\right)$

式中，kurt(φ₁)、kurt(φ₂)、kurt(φ₃)和kurt(φ₄)分别是相位校正因子特征方程计算的四个相位校正因子φ₁、φ₂、φ₃和φ₄对应的峰度函数值，其中四个峰度函数值中最大值对应的相位校正因子就是最佳相位校正因子；

7)实现相位校正；

y[n]_best＝x[n]cosφ_best-x_H[n]sinφ_best      (11)

式中，φ_best是最佳相位校正因子，由方程(10)确定，单位是弧度，x[n]是地震数据序列，x_H[n]是希尔伯特地震道序列，由方程(4)计算确定，y[n]_best是最佳地震数据相位校正序列，n是地震数据样点顺序号，n＝1,2,3,...,N，N是地震数据样点个数；

由步骤2)----步骤6)确定出最佳相位校正因子φ_best，根据相位校正因子φ_best，由方程(11)进行相位校正处理；

8)实现串联快速相位校正；

重复步骤2)----步骤7)，就实现了串联快速相位校正因子确定和快速相位校正处理，以实现消除地震数据剩余子波的快速相位校正处理；

9)绘制消除剩余子波后的地震数据剖面和存储消除剩余子波后的地震数据。

2.根据权利要求1的方法，特点是步骤1)所述的预处理是指对地震数据置标签、定义观测系统、速度分析、动校正、反褶积、叠加处理。

3.根据权利要求1的方法，特点是步骤2)所述的计算凯泽窗希尔伯特因子序列包括计算第一类零阶修正贝赛尔函数和凯泽窗希尔伯特因子序列：

按照以下公式计算第一类零阶修正贝赛尔函数I₀(x)：

I₀(x)＝a₀+a₁+a₂+...+a_k      (1)

式中，I₀(x)为第一类零阶修正贝赛尔函数；x为第一类零阶修正贝赛尔函数自变量，其取值范围为0≤x≤1；a₀是第一类零阶修正贝赛尔函数初始值，且a₀＝1；k表示第一类零阶修正贝赛尔函数的第k项；a_k是第一类零阶修正贝赛尔函数第k项值；

按照以下公式计算第一类零阶修正贝赛尔函数第k项值a_k：

$a_k=\frac{1}{k^2}{\left(\frac{x}{2}\right)}^2{a}_{k-1},k=\mathrm{1,2,3},...---\left(2\right)$

式中，x为第一类零阶修正贝赛尔函数自变量，其取值范围为0≤x≤1；a₀是第一类零阶修正贝赛尔函数初始值，a_k是第一类零阶修正贝赛尔函数第k项值；

在第一类零阶修正贝赛尔函数中，变量x的取值范围为0≤x≤1，在有限个k值之后，a_k之值越来越小并且无限接近于零，因此可以取a_k＜10^-8作为k的最终结束标志；

按照以下公式计算凯泽窗希尔伯特因子序列h_w[n]：

$h_w\left[n\right]=\left\{\begin{array}{cc}\frac{I_0\left[\mathrm{\&beta;}{(1-{[(n-\mathrm{\&alpha;})/\mathrm{\&alpha;}]}^2)}^{1/2}\right]}{I_0\left(\mathrm{\&beta;}\right)}\left\{\frac{{\sin }^2[\mathrm{\&pi;}(n-\mathrm{\&alpha;})/2]}{\mathrm{\&pi;}(n-\mathrm{\&alpha;})/2}\right\}& 0\&le;n\&le;M\\ 0& n<0,n>M\end{array}\right.---\left(3\right)$

式中，I₀(·)为第一类零阶修正贝赛尔函数；β为凯泽窗函数参数，M为凯泽窗函数样点个数参数，且α＝M/2，_n为凯泽窗希尔伯特因子序列样点顺序号，n＝1,2,3,...,M。

4.根据权利要求1的方法，特点是步骤3)所述的地震数据对应的希尔伯特地震道序列x_H[n]按照以下公式计算：

x_H[n]＝x[n]*h_w[n]      (4)

式中，x[n]为地震数据序列，x_H[n]为希尔伯特地震道序列，h_w[n]为凯泽窗希尔伯特因子序列，由方程(3)计算，符号“*”表示褶积运算；n为地震数据样点顺序号，n＝1,2,3,...,N，N为地震数据样点序列个数。