CN104199095A - 提高地震记录分辨率的反褶积方法 - Google Patents

Abstract

本发明提供了一种提高地震记录分辨率的反褶积方法。所述方法包括将预处理后的地面地震记录和VSP地震记录分别划分为多个时窗，并对每个时窗的地震记录执行以下步骤：a、在地面地震记录的双谱的复赛谱域中提取第一混合相位子波；b、在VSP地震记录的双谱的复赛谱域中提取第二混合相位子波；c、从第一混合相位子波向第二混合相位子波进行匹配滤波，计算得到子波替换反褶积算子；d、将所述子波替换反褶积算子与所述地面地震记录进行子波反褶积，得到提高分辨率后的地面地震记录。根据本发明的方法，提高分辨率效果明显、同相轴连续性较好、噪声控制较好，算法稳定。

Description

提高地震记录分辨率的反褶积方法

技术领域

本发明涉及石油地质勘探领域，更具体地讲，涉及一种提高地震资料分辨率的方法。

背景技术

子波反褶积是提高地震资料分辨率的主要方法之一，其关键在于子波的准确性。地震子波估计的准确性严重影响了地震反褶积剖面的保真度以及波阻抗反演的可靠性。在现有技术的子波处理(子波反褶积)方法都是围绕二阶统计-自相关来进行的，利用自相关函数估计地震子波时，通常假设子波为最小相位的或零相位，但在实际地震波传播过程中，由于地层的吸收作用，使得子波不满足最小相位，因此在此假设条件下估计的子波是不准确的。

并且地面地震资料在地面地震接收时干扰严重，高频成分在中深层衰减和吸收作用较大，故地面地震得到的叠后剖面中深层频率偏低，现有技术中纯粹利用地面地震资料的信息来改善地震资料质量的能力是有限的。

因此，亟需开发一种能够克服现有技术中的不足并提高地面地震记录分辨率的子波反褶积方法。

发明内容

针对现有技术中存在的不足，本发明的目的之一在于解决上述现有技术中存在的一个或多个问题。例如，本发明的目的之一在于提供一种在高地面地震记录分辨率的子波反褶积方法，以突破纯粹地面地震记录提高分辨率能力的限制。

为了实现上述目的，本发明提供了一种提高地震记录分辨率的反褶积方法。所述方法包括将预处理后的地面地震记录和VSP地震记录分别划分为多个时窗，并对每个时窗的地震记录执行以下步骤：

a、在地面地震记录的双谱的复赛谱域中提取第一混合相位子波；b、在VSP地震记录的双谱的复赛谱域中提取第二混合相位子波；c、从第一混合相位子波向第二混合相位子波进行匹配滤波，计算得到子波替换反褶积算子；d、将所述子波替换反褶积算子与所述地面地震记录进行子波反褶积，得到提高分辨率后的地面地震记录。

根据本发明提高地震记录分辨率的反褶积方法的一个实施例，所述步骤a包括：

地面地震记录基于褶积原理表示为：

x(t)＝w(t)*r(t)+n(t)   (1)

在等式(1)中，x(t)为地面地震记录，w(t)为地震子波，r(t)为地层脉冲响应，n(t)为噪声信号。

假设地层脉冲响应r(t)是平稳随机独立非高斯过程，噪声n(t)为高斯白噪声，则根据高阶累积量的性质，地层脉冲响应和噪声的三阶累积量分别为：

${\mathrm{ACCUM}}_{3r}({\mathrm{\τ}}_1,{\mathrm{\τ}}_2)={\mathrm{\γ}}_{3r}{\mathrm{\δ}}_{{\mathrm{\τ}}_1,{\mathrm{\τ}}_2}---\left(\text{2}\right)$

ACCUM_3n(τ₁,τ₂)＝0   (3)

在等式(2)和(3)中，γ_3r为地层脉冲响应的斜度，γ_3r为常数，为多维脉冲函数，ACCUM_3r(τ₁,τ₂)为地层脉冲响应的三阶累积量，ACCUM_3n(τ₁,τ₂)为噪声的三阶累积量，τ₁、τ₂为三阶累加量中时间延迟，而地面地震记录的三阶累积量ACCUM_3x(τ₁,τ₂)为：

ACCUM_3x(τ₁,τ₂)＝γ_3r∑w(t)w(t+τ₁)w(t+τ₂)   (4)

对等式(4)做二维傅里叶变换可得地面地震记录的双谱：

$B_x({\mathrm{\ω}}_1,{\mathrm{\ω}}_2)={\mathrm{\γ}}_{3r}W\left(e^{{\mathrm{i\ω}}_1}\right)W\left(e^{{\mathrm{i\ω}}_2}\right)W\left(e^{-i({\mathrm{\ω}}_1+{\mathrm{\ω}}_2)}\right)---\left(5\right)$

在等式(5)中，ω₁和ω₂表示频率，B_x(ω₁,ω₂)为地面地震记录的双谱，为地震子波频率为ω₁的傅里叶谱，为地震子波频率为ω₂的傅里叶谱，为地震子波频率为ω₁+ω₂的傅里叶谱。

令z＝e^iω，并定义双谱域中复赛谱表达式：

C_x(z₁,z₂)＝lnB_x(z₁,z₂)   (6)

可利用Z变换将混合相位地震子波分解为最小相位分量和最大相位分量：

W(z)＝W_min(z^-1)W_max(z)   (7)

在等式(7)中，W_min(z^-1)为地震子波最小相位多项式，W_max(z)为地震子波最大相位多项式，其中：

$\left\{\begin{array}{c}{W}_{\min }\left(z^{-1}\right)=\underset{i=1}{\overset{p}{\mathrm{\&Pi;}}}(1-a_i{z}^{-1}),\left|a_i\right|<1\\ W_{\max }\left(z\right)=\underset{i=1}{\overset{q}{\mathrm{\&Pi;}}}(1-b_{i}z),\left|b_i\right|<1\end{array}\right.---\left(\text{8}\right)$

在等式(8)中，a_i和b_i为多项式系数，p和q为多项式次数。

将等式(5)、(7)和(8)代入等式(6)得：

$\begin{array}{c}{C}_x(z_1,z_2)=\ln B_x(z_1,z_2)\\ =\ln \mathrm{\&gamma;}3r+\ln W\left(z_1\right)+\ln W\left(z_2\right)+\ln W(z_1^{-1},z_2^{-1})\\ =\ln \mathrm{\&gamma;}3r+\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_1^{-1})+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_i)\\ +\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_2^{-1})+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_2)\\ +\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_1{z}_2)+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_1^{-1}{z}_2^{-1})\end{array}---\left(9\right)$

对等式(9)作逆Z变换获得时间域复赛谱：

令等式(10)中对等式(9)取关于z₁或者z₂的偏微分，并把等式(10)代入，以获得三阶累积量和复赛谱的关系：

$\begin{array}{c}\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}\left\{A^{\left(i\right)}\right[{\mathrm{ACCUM}}_x(m-i,n)-{\mathrm{ACCUM}}_x(m+i,n+i)\left]\right\}\\ +\underset{j=1}{\overset{q}{\mathrm{\&Sigma;}}}\left\{B^{\left(j\right)}\right[{\mathrm{ACCUM}}_x(m-j,n-j)-{\mathrm{ACCUM}}_x(m+j,n)\left]\right\}\\ =-{\mathrm{mACCUM}}_x(m,n)\end{array}---\left(11\right)$

对于等式(11)，令ω＝max[p,q]取整数，选择m＝-ω,...,0,...ω，n＝-ω/2,...,0,...,ω/2构建以下线性方程：

Da＝y   (12)

在等式(12)中，D是由ACCUM_3x不同时间延迟三阶累积量差构成的矩阵，a＝[A⁽¹⁾,...,A^(p),B⁽¹⁾,...,B^(q)]^T是未知向量，y是由ACCUM_3x不同时间延迟三阶累积量构成的列向量，则等式(12)的最小平方解为：

a＝[D^TD]^-1D^Ty   (13)

求出未知向量a并利用等式(7)得到所述第一混合相位子波。

根据本发明提高地震记录分辨率的反褶积方法的一个实施例，所述步骤b对时窗内的VSP地震记录采用与所述步骤a相同的操作得到所述第二混合相位子波。

根据本发明提高地震记录分辨率的反褶积方法的一个实施例，所述步骤c包括：根据维纳滤波原理，得到等式(14)并通过求解等式(14)得到所述子波替换反褶积算子：

$\left[\begin{array}{cccc}{\mathrm{\&phi;}}_{11}\left(0\right)& {\mathrm{\&phi;}}_{11}\left(1\right)& ...& {\mathrm{\&phi;}}_{11}\left(k\right)\\ {\mathrm{\&phi;}}_{11}\left(1\right)& {\mathrm{\&phi;}}_{11}\left(0\right)& ...& {\mathrm{\&phi;}}_{11}(k-1)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\&phi;}}_{11}\left(k\right)& {\mathrm{\&phi;}}_{11}(k-1)& ...& {\mathrm{\&phi;}}_{11}\left(0\right)\end{array}\right]\left[\begin{array}{c}f\left(0\right)\\ f\left(1\right)\\ .\\ .\\ .\\ f\left(k\right)\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&phi;}}_{12}\left(0\right)\\ {\mathrm{\&phi;}}_{12}\left(1\right)\\ .\\ .\\ .\\ {\mathrm{\&phi;}}_{12}\left(k\right)\end{array}\right]---\left(14\right)$

在等式(14)中，φ₁₁(0),φ₁₁(1),...,φ₁₁(k-1),φ₁₁(k)为所述第一混合相位子波的自相关序列，φ₁₂(0),φ₁₂(1),...,φ₁₂(k-1),φ₁₂(k)为所述第一混合相位子波与所述第二混合相位子波的互相关序列，f(0),f(1),...,f(k)为所述子波替换反褶积算子，k为子波长度。

与现有技术相比，本发明的有益效果包括：在高阶双谱域中进行复赛谱子波分离，充分利用高阶理论和复赛谱同态理论提取子波的优势，得到精度高、抗噪强的混合相位子波，实际资料在提高分辨率的子波反褶积后，提高了地面反射波的主频，拓宽了地震信号的频带宽度，从而提高了地面地震资料的分辨率，同时也能改善了剖面同相轴的连续性，加强了波阻特征，使叠后处理剖面质量得到较大的改进。

附图说明

通过下面结合附图进行的描述，本发明的上述和其他目的和特点将会变得更加清楚，其中：

图1a示出了根据本发明示例性实施例的提高地震记录分辨率的反褶积方法在地面地震剖面上提取的子波；

图1b示出了根据本发明示例性实施例的提高地震记录分辨率的反褶积方法在VSP地震剖面上提取的子波；

图1c示出了某地区经预处理后的原始地面地震剖面；

图1d示出了根据本发明示例性实施例的提高地震记录分辨率的反褶积方法处理之后的地面地震剖面；

图1e是图1c的振幅谱；

图1f是图1d的振幅谱。

具体实施方式

在下文中，将结合附图和示例性实施例详细地描述根据本发明的提高地震记录分辨率的反褶积方法。

本发明的主要技术思路是：在高阶双谱域中对地面地震记录和VSP地震数据进行复赛谱子波分离，得到地精度高、抗噪强的混合相位子波，然后对地面地震记录和VSP地震数据的地震子波进行有效的匹配滤波得到反褶积算子，从而进行子波反褶积得到高质量的高分辨率地震剖面。

基于以上技术思路，还需要解决以下技术难题：

(1)从复赛谱域分离子波，不需要子波最小相位假设，但是必须假设在复赛谱域内子波和反射系数的频率是分开的，而实际上是部分重合的，如何提高其子波提取精度成为难点。

(2)随着高阶统计理论的进一步发展和完善，高阶统计信号处理已经渗透到信号处理各个应用领域，并涌现出了大量的理论和应用研究成果。高阶统计(高阶谱)较之常规二阶统计(功率谱)包含信号更为完整的相位信息，但其运算效率和抗噪能力成为其在实际生产中运用的难点。

也就是说，利用本发明的地震记录分辨率的子波反褶积方法的关键在于：是否能在高阶双谱域中对地震记录和VSP数据进行复赛谱子波分离，得到地精度高、抗噪强的混合相位子波；能否对地震记录和VSP数据的地震子波进行有效的匹配滤波得到反褶积算子，从而进行子波反褶积得到高质量的高分辨率地震剖面。

根据本发明示例性实施例的提高地震记录分辨率的反褶积方法包括以下步骤：

第一步：对从野外采集到的地面地震记录(又称为地面地震资料、地面地震数据、地面地震剖面)和VSP地震数据进行预处理，得到叠加或偏移后的地面地震数据和VSP地震数据，然后，对预处理后的地面地震数据和VSP地震数据分别采用由浅层到深层划分时窗，并对每个时窗内的数据采用下面的步骤进行进行处理。以下，除另有说明之外，所处理的地面地震数据和VSP地震数据均是指的预处理后某个时窗内的地面地震数据和VSP地震数据。

这里，预处理是采用本领域的常规方法，例如将原始采集的地面地震数据和VSP地震数据经过静校正、去噪、振幅补偿、动校正、叠加以及偏移等处理，最终形成叠后或偏移后的地面地震数据和VSP地震数据，在此不再赘述。

地面地震数据是指地面激发地面接收的观测方式，即炮点设置在地面上且检波点也设在地面。

VSP(Vertical Seismic Profiling，即垂直地震剖面法)地震记录(或称为VSP地震数据、VSP地震资料、VSP数据、VSP剖面)为地面激发井中接收的观测方式，即炮点设置在地面而检波点设置在井内的不同深度，相对地面地震资料而言，VSP观测到的反射波传播路径短、受近地表低速带和环境噪声影响均较小，因此VSP资料具有高分辨率、高信噪比等优势，特别零偏移距VSP地震资料往往具有高频子波信息，所以利用VSP地震资料提高地面地震资料的分辨率，可增强地面地震资料提高分辨率的能力。

时窗：由于地震资料通常以时间作为纵向的标尺展示(地震波接收地下岩层反射信号以时间标度，乘以速度，可以转换为地下的深度)，对目的层的纵向范围可以用“时窗”来表述。

双谱：高阶累计量谱简称为高阶谱或多谱，最常用的高阶谱是三阶谱(也称为双谱)。

第二步：对时窗内的地面地震记录和VSP地震数据在双谱域中提取复赛谱子波。

本步骤按下列方式完成：

地面地震记录剖面基于褶积原理表示为：

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

在等式(2-1)中，x(t)为地面地震记录，w(t)为地震子波(简称为子波)，r(t)为地层脉冲响应(又称为地层反射序列)，n(t)为噪声信号。

假设地层脉冲响应r(t)是平稳随机独立非高斯过程(非对称)，噪声n(t)为高斯白噪声，则根据高阶累积量的性质，地层脉冲响应和噪声的三阶累积量分别为：

${\mathrm{ACCUM}}_{3r}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2)={\mathrm{\&gamma;}}_{3r}{\mathrm{\&delta;}}_{{\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2}---\left(\text{2-2}\right)$

ACCUM_3n(τ₁,τ₂)＝0   (2-3)

在等式(2-2)和(2-3)中，γ_3r为地层脉冲响应的斜度，γ_3r为常数，为多维脉冲函数，ACCUM_3r(τ₁,τ₂)为地层脉冲响应的三阶累积量，ACCUM_3n(τ₁,τ₂)为噪声的三阶累积量，τ₁、τ₂为三阶累加量中时间延迟。而地面地震记录剖面的三阶累积量ACCUM_3x(τ₁,τ₂)为：

ACCUM_3x(τ₁,τ₂)＝γ_3r∑w(t)w(t+τ₁)w(t+τ₂)   (2-4)

对等式(4)做二维傅里叶变换可得地面地震记录的双谱：

$B_x({\mathrm{\&omega;}}_1,{\mathrm{\&omega;}}_2)={\mathrm{\&gamma;}}_{3r}W\left(e^{{\mathrm{i\&omega;}}_1}\right)W\left(e^{{\mathrm{i\&omega;}}_2}\right)W\left(e^{-i({\mathrm{\&omega;}}_1+{\mathrm{\&omega;}}_2)}\right)---(2-5)$

在等式(2-5)中，ω₁和ω₂表示频率，B_x(ω₁,ω₂)为地面地震记录的双谱，为地震子波的傅里叶谱，具体地，为地震子波频率为ω₁的傅里叶谱，为地震子波频率为ω₂的傅里叶谱，为地震子波频率为ω₁+ω₂的傅里叶谱。这里，等式(2-4)中小写的w表示时间域中的子波，等式(2-5)中大写的W表示傅里叶变换后频率域中的子波。

令z＝e^iω，并定义双谱域中复赛谱表达式：

C_x(z₁,z₂)＝lnB_x(z₁,z₂)   (2-6)

地震子波为混合相位的地震子波(简称为混合相位子波)，则可利用Z变换将其分解为最小相位分量和最大相位分量：

W(z)＝W_min(z^-1)W_max(z)   (2-7)

在等式(2-7)中，W_min(z^-1)为地震子波最小相位多项式，W_max(z)为地震子波最大相位多项式，其中：

$\left\{\begin{array}{c}{W}_{\min }\left(z^{-1}\right)=\underset{i=1}{\overset{p}{\mathrm{\&Pi;}}}(1-a_i{z}^{-1}),\left|a_i\right|<1\\ W_{\max }\left(z\right)=\underset{i=1}{\overset{q}{\mathrm{\&Pi;}}}(1-b_{i}z),\left|b_i\right|<1\end{array}\right.---\left(\text{2-8}\right)$

在等式(2-8)中，a_i和b_i表示多项式系数，p和q表示多项式次数。

将等式(2-5)、(2-7)和(2-8)代入等式(2-6)得：

$\begin{array}{c}{C}_x(z_1,z_2)=\ln B_x(z_1,z_2)\\ =\ln \mathrm{\&gamma;}3r+\ln W\left(z_1\right)+\ln W\left(z_2\right)+\ln W(z_1^{-1},z_2^{-1})\\ =\ln \mathrm{\&gamma;}3r+\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_1^{-1})+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_i)\\ +\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_2^{-1})+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_2)\\ +\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_1{z}_2)+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_1^{-1}{z}_2^{-1})\end{array}---(2-9)$

对等式(2-9)作逆Z变换可以获得时间域复赛谱：

令等式(2-10)中则包含了地震子波最小相位和最大相位信息，对等式(2-9)取关于z₁或者z₂的偏微分，并把等式(2-10)代入，可获得三阶累积量和复赛谱的关系：

$\begin{array}{c}\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}\left\{A^{\left(i\right)}\right[{\mathrm{ACCUM}}_x(m-i,n)-{\mathrm{ACCUM}}_x(m+i,n+i)\left]\right\}\\ +\underset{j=1}{\overset{q}{\mathrm{\&Sigma;}}}\left\{B^{\left(j\right)}\right[{\mathrm{ACCUM}}_x(m-j,n-j)-{\mathrm{ACCUM}}_x(m+j,n)\left]\right\}\\ =-{\mathrm{mACCUM}}_x(m,n)\end{array}---(2-11)$

对于等式(2-11)，令ω＝max[p,q]取整数，选择m＝-ω,...,0,...ω，n＝-ω/2,...,0,...,ω/2可以构建一个线性方程：

Da＝y   (2-12)

在等式(2-12)中，D是由ACCUM_3x不同时间延迟三阶累积量差构成的矩阵，a＝[A⁽¹⁾,...,A^(p),B⁽¹⁾,...,B^(q)]^T是未知向量，y是由ACCUM_3x不同时间延迟三阶累积量构成的列向量，则等式(2-12)的最小平方解为：

a＝[D^TD]^-1D^Ty   (2-13)

求出未知向量a并利用等式(2-7)得到双谱域提取的第一复赛谱子波W₁(t)(又称地面子波、第一混合相位子波)。

同理，可求得VSP剖面在双谱域中提取第二复赛谱子波W₂(t)(又称VSP子波、第二混合相位子波)，在此不再赘述。

第三步：利用第二步得到的子波进行匹配滤波，得到反褶积算子。

对于提取的地面子波与VSP子波，可以从地面子波向VSP子波进行匹配，计算得到子波替换反褶积算子。

根据第二步可知，地面子波数据W₁(t)为：

W₁(0),W₁(1),...,W₁(i),...,W₁(k)   (3-1)

在等式(3-1)中，k为子波长度。

地面子波的自相关序列可表示为：

φ₁₁(0),φ₁₁(1),...,φ₁₁(i),...,φ₁₁(k)   (3-2)

并且，根据第二步，VSP子波数据W₂(t)为：

W₂(0),W₂(1),...,W₂(i),...,W₂(k)   (3-3)

那么地震子波与VSP子波的互相关序列可表示为：

φ₁₂(0),φ₁₂(1),...,φ₁₂(i),...,φ₁₂(k)   (3-4)

根据维纳滤波原理，得到TOEPLITZ方程：

$\left[\begin{array}{cccc}{\mathrm{\&phi;}}_{11}\left(0\right)& {\mathrm{\&phi;}}_{11}\left(1\right)& ...& {\mathrm{\&phi;}}_{11}\left(k\right)\\ {\mathrm{\&phi;}}_{11}\left(1\right)& {\mathrm{\&phi;}}_{11}\left(0\right)& ...& {\mathrm{\&phi;}}_{11}(k-1)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\&phi;}}_{11}\left(k\right)& {\mathrm{\&phi;}}_{11}(k-1)& ...& {\mathrm{\&phi;}}_{11}\left(0\right)\end{array}\right]\left[\begin{array}{c}f\left(0\right)\\ f\left(1\right)\\ .\\ .\\ .\\ f\left(k\right)\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&phi;}}_{12}\left(0\right)\\ {\mathrm{\&phi;}}_{12}\left(1\right)\\ .\\ .\\ .\\ {\mathrm{\&phi;}}_{12}\left(k\right)\end{array}\right]---(3-5)$

求解等式(3-5)，可得到地震子波与VSP子波匹配后的反褶积算子f(t)，t＝0,1,...,k。

第四步：利用第三步得到的反褶积算子进行子波反褶积，得到高分辨率地震剖面。即，对地面地震剖面数据与第三步得到的反褶积算子进行褶积，得到提高分辨率后的地震剖面。

图1a示出了根据本发明示例性实施例的提高地震记录分辨率的反褶积方法在第二步地面地震剖面上提取的子波；图1b示出了根据本发明示例性实施例的提高地震记录分辨率的反褶积方法在第二步VSP地震剖面上提取的子波；图1c示出了某地区经预处理后的原始地面地震剖面，其中在420道～450道嵌入了VSP叠加数据；图1d示出了根据本发明示例性实施例的提高地震记录分辨率的反褶积方法处理之后的地面地震剖面，其中在420道～450道嵌入了VSP叠加数据；图1e是图1c的振幅谱；图1f是图1d的振幅谱。

根据图1a和图1b可知：VSP地震剖面提取的子波要比地面地震剖面提取的子波瘦，即：VSP地震剖面提取的子波的分辨率要相对高一些，利用地震子波向VSP地震剖面提取的子波进行匹配，可以计算出分辨率更高的子波反褶积算子。

分析图1c和图1d可知：经本发明反褶积后，分辨率得到明显提高，剖面同相轴的连续性及波阻特征均得到加强，而且处理后的地面地震剖面与嵌入的VSP剖面的层位能一一对应，即子波反褶积使地面地震数据和VSP数据具有更高的闭合度。

从图1e和1f可知：采用本发明方法进行褶积后，地震信号的高频成分得到了加强，频带宽度得到了拓宽，主频向高频端移动。

综上所述，本发明利用高阶统计理论进行复赛谱子波反褶积来克服现有反褶积技术的不足，同时提取VSP子波与地震记录的子波进行匹配滤波得到反褶积算子，利用该反褶积算子与地震记录进行子波反褶积得到的高分辨率剖面，能够改善剖面同相轴的连续性，加强波阻特征，为后续解释提供质量更高的剖面。由实际生产应用效果也可以看出，根据本发明提高分辨率效果明显、同相轴连续性较好、噪声控制较好，算法稳定。

尽管上面已经通过结合示例性实施例描述了本发明，但是本领域技术人员应该清楚，在不脱离权利要求所限定的精神和范围的情况下，可对本发明的示例性实施例进行各种修改和改变。

Claims (4)

1.一种提高地震记录分辨率的反褶积方法，其特征在于，所述方法包括将预处理后的地面地震记录和VSP地震记录分别划分为多个时窗，并对每个时窗的地震记录执行以下步骤：

a、在地面地震记录的双谱的复赛谱域中提取第一混合相位子波；

b、在VSP地震记录的双谱的复赛谱域中提取第二混合相位子波；

c、从第一混合相位子波向第二混合相位子波进行匹配滤波，计算得到子波替换反褶积算子；

d、将所述子波替换反褶积算子与所述地面地震记录进行子波反褶积，得到提高分辨率后的地面地震记录。

2.根据权利要求1所述的提高地震记录分辨率的反褶积方法，其特征在于，所述步骤a包括：

地面地震记录基于褶积原理表示为：

x(t)＝w(t)*r(t)+n(t)   (1)

在等式(1)中，x(t)为地面地震记录，w(t)为地震子波，r(t)为地层脉冲响应，n(t)为噪声信号；

假设地层脉冲响应r(t)是平稳随机独立非高斯过程，噪声n(t)为高斯白噪声，则根据高阶累积量的性质，地层脉冲响应和噪声的三阶累积量分别为：

${\mathrm{ACCUM}}_{3r}({\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2)={\mathrm{\&gamma;}}_{3r}{\mathrm{\&delta;}}_{{\mathrm{\&tau;}}_1,{\mathrm{\&tau;}}_2}---\left(\text{2}\right)$

ACCUM_3n(τ₁,τ₂)＝0   (3)

在等式(2)和(3)中，γ_3r为地层脉冲响应的斜度，γ_3r为常数，为多维脉冲函数，ACCUM_3r(τ₁,τ₂)为地层脉冲响应的三阶累积量，ACCUM_3n(τ₁,τ₂)为噪声的三阶累积量，τ₁、τ₂为三阶累加量中时间延迟，而地面地震记录的三阶累积量ACCUM_3x(τ₁,τ₂)为：

ACCUM_3x(τ₁,τ₂)＝γ_3r∑w(t)w(t+τ₁)w(t+τ₂)   (4)

对等式(4)做二维傅里叶变换可得地面地震记录的双谱：

$B_x({\mathrm{\&omega;}}_1,{\mathrm{\&omega;}}_2)={\mathrm{\&gamma;}}_{3r}W\left(e^{{\mathrm{i\&omega;}}_1}\right)W\left(e^{{\mathrm{i\&omega;}}_2}\right)W\left(e^{-i({\mathrm{\&omega;}}_1+{\mathrm{\&omega;}}_2)}\right)---\left(5\right)$

在等式(5)中，ω₁和ω₂表示频率，B_x(ω₁,ω₂)为地面地震记录的双谱，为地震子波频率为ω₁的傅里叶谱，为地震子波频率为ω₂的傅里叶谱，为地震子波频率为ω₁+ω₂的傅里叶谱；

令z＝e^iω，并定义双谱域中复赛谱表达式为：

C_x(z₁,z₂)＝lnB_x(z₁,z₂)   (6)

可利用Z变换将混合相位地震子波分解为最小相位分量和最大相位分量：

W(z)＝W_min(z^-1)W_max(z)   (7)

在等式(7)中，W_min(z^-1)为地震子波最小相位多项式，W_max(z)为地震子波最大相位多项式，其中：

$\left\{\begin{array}{c}{W}_{\min }\left(z^{-1}\right)=\underset{i=1}{\overset{p}{\mathrm{\&Pi;}}}(1-a_i{z}^{-1}),\left|a_i\right|<1\\ W_{\max }\left(z\right)=\underset{i=1}{\overset{q}{\mathrm{\&Pi;}}}(1-b_{i}z),\left|b_i\right|<1\end{array}\right.---\left(\text{8}\right)$

在等式(8)中，a_i和b_i表示多项式系数，p和q表示多项式次数；

将等式(5)、(7)和(8)代入等式(6)得：

$\begin{array}{c}{C}_x(z_1,z_2)=\ln B_x(z_1,z_2)\\ =\ln \mathrm{\&gamma;}3r+\ln W\left(z_1\right)+\ln W\left(z_2\right)+\ln W(z_1^{-1},z_2^{-1})\\ =\ln \mathrm{\&gamma;}3r+\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_1^{-1})+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_i)\\ +\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_2^{-1})+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_2)\\ +\ln \underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}(1-a_i{z}_1{z}_2)+\ln \underset{i=1}{\overset{q}{\mathrm{\&Sigma;}}}(1-b_i{z}_1^{-1}{z}_2^{-1})\end{array}---\left(9\right)$

对等式(9)作逆Z变换获得时间域复赛谱：

令等式(10)中对等式(9)取关于z₁或者z₂的偏微分，并把等式(10)代入，以获得三阶累积量和复赛谱的关系：

$\begin{array}{c}\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}\left\{A^{\left(i\right)}\right[{\mathrm{ACCUM}}_x(m-i,n)-{\mathrm{ACCUM}}_x(m+i,n+i)\left]\right\}\\ +\underset{j=1}{\overset{q}{\mathrm{\&Sigma;}}}\left\{B^{\left(j\right)}\right[{\mathrm{ACCUM}}_x(m-j,n-j)-{\mathrm{ACCUM}}_x(m+j,n)\left]\right\}\\ =-{\mathrm{mACCUM}}_x(m,n)\end{array}---\left(11\right)$

对于等式(11)，令ω＝max[p,q]取整数，选择m＝-ω,...,0,...ω，n＝-ω/2,...,0,...,ω/2构建以下线性方程：

Da＝y   (12)

在等式(12)中，D是由ACCUM_3x不同时间延迟三阶累积量差构成的矩阵，a＝[A⁽¹⁾,...,A^(p),B⁽¹⁾,...,B^(q)]^T是未知向量，y是由ACCUM_3x不同时间延迟三阶累积量构成的列向量，则等式(12)的最小平方解为：

a＝[D^TD]^-1D^Ty   (13)

求出未知向量a并利用等式(7)得到所述第一混合相位子波。

3.根据权利要求2所述的提高地震记录分辨率的反褶积方法，其特征在于，所述步骤b对时窗内的VSP地震记录采用与所述步骤a相同的操作得到所述第二混合相位子波。

4.根据权利要求1所述的提高地震记录分辨率的反褶积方法，其特征在于，所述步骤c包括根据维纳滤波原理，得到等式(14)并通过求解等式(14)得到所述子波替换反褶积算子：

$\left[\begin{array}{cccc}{\mathrm{\&phi;}}_{11}\left(0\right)& {\mathrm{\&phi;}}_{11}\left(1\right)& ...& {\mathrm{\&phi;}}_{11}\left(k\right)\\ {\mathrm{\&phi;}}_{11}\left(1\right)& {\mathrm{\&phi;}}_{11}\left(0\right)& ...& {\mathrm{\&phi;}}_{11}(k-1)\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\&phi;}}_{11}\left(k\right)& {\mathrm{\&phi;}}_{11}(k-1)& ...& {\mathrm{\&phi;}}_{11}\left(0\right)\end{array}\right]\left[\begin{array}{c}f\left(0\right)\\ f\left(1\right)\\ .\\ .\\ .\\ f\left(k\right)\end{array}\right]=\left[\begin{array}{c}{\mathrm{\&phi;}}_{12}\left(0\right)\\ {\mathrm{\&phi;}}_{12}\left(1\right)\\ .\\ .\\ .\\ {\mathrm{\&phi;}}_{12}\left(k\right)\end{array}\right]---\left(14\right)$

在等式(14)中，φ₁₁(0),φ₁₁(1),...,φ₁₁(k-1),φ₁₁(k)为所述第一混合相位子波的自相关序列，φ₁₂(0),φ₁₂(1),...,φ₁₂(k-1),φ₁₂(k)为所述第一混合相位子波与所述第二混合相位子波的互相关序列，f(0),f(1),...,f(k)为所述子波替换反褶积算子，k为子波长度。