CN103675912B - 一种基于vsp时变偏振分析的矢量波场合成及分离方法 - Google Patents

Abstract

本发明提供了一种基于VSP时变偏振分析的矢量波场合成及分离方法，属于石油地球物理勘探领域。所述方法包括以下步骤：(1)读入水平H_p分量和垂直z分量这两个记录；(2)在频率域进行时窗自动选取得到偏振分析时窗；(3)在步骤(2)得到的偏振分析时窗内利用协方差矩阵法判断所述两个记录的波场类型，如果判断结果是单一波场，则采用能量准则法进行偏振分析，得到分离及合成后的P波和SV波的波场，如果判断结果是混合波场，则采用Cho和Spencer法进行偏振分析，得到分离及合成后的P波和SV波的波场。利用本发明方法可以使P波及SV波波场较好的分离，并且能够从标量波场还原矢量波场。

Description

一种基于VSP时变偏振分析的矢量波场合成及分离方法

技术领域

本发明属于石油地球物理勘探领域，具体涉及一种基于VSP时变偏振分析的矢量波场合成及分离方法。

背景技术

偏振分析又称极化分析，是通过研究三分量地震记录的波的偏振特性来分离纵、横波和提高信噪比的一种重要处理技术。三分量VSP记录所提供的重要特征是地震波的空间偏振特性，而偏振分析主要是通过利用各种类型的波的偏振方向来确定其线性偏振信号方位，从而进行波场分离，之后由水平方向和垂直方向的标量记录合成不同波的矢量波场。只有在矢量波场中，才能从各种波的空间偏振分量上得到波真实的矢量振幅和波场，才能为资料处理和解释提供保证。

目前，在VSP偏振分析中常用的方法是用能量准则解析法来确定各种波的偏振主方向，然后进行偏振投影。这种时不变的偏振分析方法在VSP三分量资料处理中被广泛的应用，并取得了较好的效果。

现阶段的VSP偏振分析方法主要适用的方法有：能量准则解析法、协方差矩阵分析法，Cho和Spencer估算混合波场中各种波极化方向偏振分析法，以及根据模型利用射线和波动方程等求取偏振方向的方法。时变偏振分析方法的分析研究仍处于理论实验的阶段；在实际资料中的应用较少，且效果并不是很好。

发明内容

本发明的目的在于解决上述现有技术中存在的难题，提供一种基于VSP时变偏振分析的矢量波场合成及分离方法，将能量准则、协方差、Cho和Spencer三种方法相融合，在时窗内根据波场类型的不同，选取最优的偏振分析方法，从而进行VSP复杂波场分离，实现矢量波场合成及分离。

本发明是通过以下技术方案实现的：

一种基于VSP时变偏振分析的矢量波场合成及分离方法，所述方法包括以下步骤：

(1)读入水平H_p分量和垂直z分量这两个记录；

(2)在频率域进行时窗自动选取得到偏振分析时窗；

(3)在步骤(2)得到的偏振分析时窗内利用协方差矩阵法判断所述两个记录的波场类型，如果判断结果是单一波场，则采用能量准则法进行偏振分析，得到分离及合成后的P波和SV波的波场，如果判断结果是混合波场，则采用Cho和Spencer法进行偏振分析，得到分离及合成后的P波和SV波的波场。

所述步骤(2)具体如下：选取高斯函数作为时窗函数，给定N进行时窗的自动选取得到偏振分析时窗，然后在此偏振分析时窗内分别对步骤(1)读入的两个记录进行时频域短时离散傅里叶变换，将两个记录从时域转化为频域；

所述高斯函数为：

$g^0\left(t\right)=2^{1/4}{e}^{2\mathrm{\π}{t}^2}$

所述时频域短时离散傅里叶变换采用短时离散傅里叶变换公式：

$\mathrm{STFT}(m,{\mathrm{\ω}}_k)=\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}x\left(n\right)g*(n-\mathrm{mN})e^{-j\frac{\mathrm{\π}}{M}\mathrm{nk}}$

式中：频率k为整数，k＝0至M，x(n)为给定的水平H_p分量或垂直z分量，n＝0至L-1，L为数据的长度，N为时间轴上的时窗函数移动的步长，M是将频域的一个周期2π分成了M个点，m是任意正整数，是时窗函数移动步长的系数，*为复共轭。

与现有技术相比，本发明的有益效果是：

(1)现有的时变偏振分析方法的时窗选取方法为：传统简单的三点法或过零点法自动选取时窗；这种方法不能很好的根据波场的实际特点和周期划分时窗的，也不能适应各种方法对时窗中波场类型的要求。然而本发明方法选取在频率域根据波场的频率特性选择时窗的方法，在一定程度上克服了现有技术的缺点；

(2)利用本发明方法可以使P波及SV波波场较好的分离，并且能够从标量波场还原矢量波场。

附图说明

图1是本发明地质模型与观测系统示意图。

图2是在图1所示的观测系统下的反射波及转换波的射线路径图。

图3-1是理论合成的下行P波和上行SV波的模拟记录。

图3-2是理论合成的下行SV波和上行P波的模拟记录。

图4-1是上行波场理论合成记录中的水平H_p分量。

图4-2是上行波场理论合成记录中的垂直z分量。

图5-1是下行波场理论合成记录中的水平H_p分量。

图5-2是下行波场理论合成记录中的垂直z分量。

图6-1是全波场波场分离的记录中的下行P波和上行SV波。

图6-2是全波场波场分离的记录中的下行SV波和上行P波。

图7-1是上行波场波场分离的记录中的上行SV波。

图7-2是上行波场波场分离的记录中的上行P波。

图8-1是下行波场波场分离的记录中的下行P波。

图8-2是下行波场波场分离的记录中的下行SV波。

图9是本发明方法的步骤框图。

具体实施方式

下面结合附图对本发明作进一步详细描述：由于能量准则法和协方差矩阵法对单一波场的计算精度较高，而Cho和Spencer方法恰恰适合在有两种不同偏振属性的波存在的混合波场中使用的优点和适用范围，因此本发明将能量准则法(请参考朱光明，垂直地震剖面法.北京[M]：石油工业出版社，1988)、协方差矩阵法(请参考朱衍镛.二分量记录的空间方向滤波[J].石油地球物理勘探，1995，30(S2).116-125)及Cho和Spencer法(请参考ChoWH，Spencer T W.Estimation of polarization and slowness in mixedwavefields.Geophysics，1992，57(6)：805-814)相结合形成一种新的混合波场时变偏振分析方法。

本发明方法首先在频率域选择适当的时窗，通过协方差矩阵法判断时窗内的波场类型，即单一波场或者混合波场，然后根据不同的波场类型，选取最优的偏振分析方法，从而进行VSP复杂波场分离，实现矢量波场合成及分离。

如图9所示，本发明方法的具体实现的步骤为：

(1)读入水平H_p分量和垂直z分量这两个记录；

(2)频率域的时窗自动选取：

时变偏振分析的基本单元为时窗选取，合理和精确的时窗是偏振分析的前提和保证，只有准确的偏振角才能进行彻底的波场分离和偏振合成。因此选择在频率域进行时窗自动选取，即根据波形在频率域的特性自动判断一个波长作为一个时窗(通过选择为高斯函数作为时窗函数，给定N(时间轴上的窗函数移动的步长)进行时窗的自动选取，后进行短时离散傅里叶变换)；

对两个记录进行时域和频率域(即时频域)短时离散傅里叶变换，窗数选择为高斯函数，即Gabor基函数：

$g^0\left(t\right)=2^{1/4}{e}^{2\mathrm{\π}{t}^2}---\left(1.1\right)$

进行短时离散傅里叶变换采用短时离散傅里叶变换公式：

$\mathrm{STFT}(m,{\mathrm{\ω}}_k)=\underset{n=-\∞}{\overset{\∞}{\mathrm{\Σ}}}x\left(n\right)g*(n-\mathrm{mN})e^{-j\frac{\mathrm{\π}}{M}\mathrm{nk}}---\left(1.2\right)$

式中：频率k为整数，k＝0至M，x(n)为给定的信号，即水平H_p分量和垂直z分量这两个记录，n＝0，1，...，L-1，L为数据(序列)的长度，N为时间轴上的窗函数移动的步长，M为频率域的采样点数，即将频域的一个周期2π分成了M个点，m是任意正整数，是窗函数移动步长的系数，*为复共轭。进行该变换的目的是将数据从时域转化为频域。

步骤(2)通过选取高斯函数(1.1)作为偏振分析的时窗，输入的地震数据在此时窗内进行短时离散傅里叶变化，将数据从时域转化为频域从而进行步骤(3)的偏振分析，之后移动时窗，进行下个时窗内的偏振分析。

(3)在时窗内利用协方差矩阵判断其波场的类型(单一波场或混合波场)；

协方差矩阵法：是一种解析信号的分析方法；采用矩阵的特征参数来反映质点轨迹图的特点，从而求解质点运动轨迹的偏振参数。对于VSP三分量记录而言，可以选取三分量记录中N个样点作为一个时窗T(t₁，t₂)，而且在这个T(t₁，t₂)时窗内，每个样点i的三个分量值x_j，y_i，z_i的均值或数学期望分别为：

$m_x=\frac{1}{N}\underset{i=N_1}{\overset{N_2}{\mathrm{\Σ}}}{x}_i,m_y=\frac{1}{N}\underset{i=N_1}{\overset{N_2}{\mathrm{\Σ}}}{y}_i,m_z=\frac{1}{N}\underset{i=N_1}{\overset{N_2}{\mathrm{\Σ}}}{z}_i---\left(2.1\right)$

其中，i＝N₁…N₂；N＝N₂-N₁+1；(N₁-N₂)Δt＝t₂-t₁；Δt是时间采样率。

由此，定义协方差矩阵就可以写作：

$M_C=\left[\begin{array}{ccc}\mathrm{\Σ}{A}^2& \mathrm{\ΣAB}& \mathrm{\ΣAC}\\ \mathrm{\ΣBA}& \mathrm{\Σ}{B}^2& \mathrm{\ΣBC}\\ \mathrm{\ΣCA}& \mathrm{\ΣCB}& C^2\end{array}\right]---\left(2.2\right)$

式中∑为A＝(x_i-m_x)；B＝(y_i-m_y)；C＝(z_i-m_z)。

通过对协方差矩阵求解计算特征值和对应的特征向量。通常情况下，可以得到λ₁、λ₂、λ₃三个特征值，它们关系为(λ₁≥λ₂≥λ₃)，并且可以求出相应的特征向量。如果时窗内为单一波场，则只有一个非零的特征值；如果时窗内记录为混合波场，则有两个非零特征值。

(4)对于单一波场，选取能量准则法进行偏振分析；

如果时窗内只有一个非零的特征值，可以判断此时窗内为单一波场，采用能量准则法进行偏振分析。

能量准则法：VSP三分量记录的能量是随着方位角而变化，能使能量达到最大的方位角就是偏振角。所以，设能量式可以用下式表示：

E(θ)＝∑_i(x_i cosθ+y_i sinθ)² (3.1)

i为从时窗起点到时窗终点所编的序号，使能量取极大的必要条件是：

$\frac{\∂E\left(\mathrm{\θ}\right)}{\∂\mathrm{\θ}}=0---\left(3.2\right)$

从而可求得计算偏振角的公式：

$\tan 2{\mathrm{\θ}}_{\max }=\frac{2{\mathrm{\Σ}}_i{x}_i{y}_i}{{\mathrm{\Σ}}_i(y_i^2-x_i^2)}---\left(3.3\right)$

即可求得偏振角θ＝θ_max。

根据偏振投影公式(3.4)式：

$\left\{\begin{array}{c}P=z\cos \mathrm{\θ}+H_p\sin \mathrm{\θ}\\ \mathrm{SV}=-z\sin \mathrm{\θ}+H_p\cos \mathrm{\θ}\end{array}\right.---\left(3.4\right)$

进行波场分离及合成，得到最终的P波和SV波场。

(5)对于混合波场，选取Cho和Spencer法进行偏振分析；

如果时窗内有两个非零的特征值，则判断此时窗内为混合波场，采用Cho和Spencer法进行偏振分析。

Cho和Spencer法是在选择了一个包含多个道集的空间滑动窗口的情况下。假设在相邻道集波场的振幅、相位变化不大的情况下，根据纵波(即P波)和转换波(SV是转化波的一种，文中的SV波是转化的SV波)在传播中相对于不同偏振极化矢量和慢度得到的垂向和水平向分量，在最小方差约束下从多分量输入数据中获取纵波和转换波的偏振方向。

多分量检波器接收波场信息时，水平和垂向检波器记录到的波场分别是纵波波场和转换波波场在水平和垂向上的投影组合。检波器记录可以表示为波场偏振极化方向与纵波和转换波组合矩阵的乘积，即第n道检波器的记录U_i ⁿ(t)时域可以表示为

$\left\{\begin{array}{c}{U}_3^{\left(n\right)}\left(t\right)=a_3^P{w}^P\left(t\right)+a_3^{\mathrm{SV}}{w}^{\mathrm{SV}}\left(t\right)\\ U_1^{\left(n\right)}\left(t\right)=a_1^P{w}^P\left(t\right)+a_1^{\mathrm{SV}}{w}^{\mathrm{SV}}\left(t\right)\end{array}\right.---\left(4.1\right)$

其中a_i ^P是P波在i方向分量上的偏振极化方向，w^p(t)是P的波场，U_i ⁽ⁿ⁾(t)代表第n道检波器在t时刻的记录。下角标i＝1时代表水平方向记录，i＝3时代表垂向记录，上角标P代表P波、SV代表转换的SV波、(n)代表第n道记录。

第n+1道波场和记录频率域可以表示为：

其中φ_P＝2πfq^PΔ和φ_SV＝2πfq^SVΔ分别是P波和SV波相邻检波器之间的相位差。

从式(4.1)和式(4.2)可得第n道与第n+1道之间在频率域的关系；

$\left\{\begin{array}{c}{U}_3^{(n+1)}\left(f\right)\\ U_1^{(n+1)}\left(f\right)\end{array}\right\}=\mathrm{A\Λ}{A}^{-1}\left\{\begin{array}{c}{U}_3^{\left(n\right)}\left(f\right)\\ U_1^{\left(n\right)}\left(f\right)\end{array}\right\}---\left(4.3\right)$

其中： $A=\left[\begin{array}{cc}{a}_3^P& a_3^{\mathrm{SV}}\\ a_1^P& a_1^{\mathrm{SV}}\end{array}\right],\mathrm{\Λ}=\left[\begin{array}{cc}{g}_P{e}^{j2\mathrm{\π}{f_q}^P\mathrm{\Δ}}& 0\\ 0& g_S{e}^{j2\mathrm{\π}{f_q}^{\mathrm{SV}}\mathrm{\Δ}}\end{array}\right].$

因为P波和SV波的偏振方向不在同一方向上，那么矩阵A是可逆的，即A^-1是矩阵A的逆矩阵。在二维复向量空间中，(4.3)式可以写成一个线性变换，可以表示为：

Uⁿ⁺¹＝YUⁿ (4.4)

其中：Y＝AΛA^-1

如果选取的同一个窗口，相邻道集波场的两种波的特性如偏振极化方向、慢度和振幅变化率等变化差异不大时，通过最小方差约束条件(4.5)可以从式(4.4)中求出传递矩阵，即

||Uⁿ⁺¹-YUⁿ||＝min (4.5)

通过垂向和水平方向波场记录可以计算出Y后，对其做特征值和特征向量分解，可以得出波场分解的最后表达式：

$\left\{\begin{array}{c}{W}^P\left(f\right)\\ W^{\mathrm{SV}}\left(f\right)\end{array}\right\}=A^{-1}\left\{\begin{array}{c}{U}_3^{\left(n\right)}\left(f\right)\\ U_1^{\left(n\right)}\left(f\right)\end{array}\right\}---\left(4.6\right)$

利用式(4.6)最终可以得到分离及合成后的两个矢量波场，即P波和SV波的波场。

波场分离和合成是一起的过程，原始的记录为H_P和z两个记录，其中都含有两种波场(P波和SV波)，先将原始的两个记录中的P波和SV波分离出，最终将两个记录中的P波和SV波合成只含有P波和只含有SV波的两个矢量波场。

具体实施时，采用C语言编程实现了本发明的方法。选取模型参数为：震源位于地面，偏移距为300m。井位于震源的左侧，检波器放置在井中，道距5m，接收点深度范围为200m～1500m，共261道接收。地层选择七层的平层地质模型，地质模型和观测系统示意图如图1所示，此观测系统下的反射波及转换波的射线路径图如图2所示。通过地质模型合成的水平H_P分量和垂直z分量的VSP记录。为了便于与偏振分析方法波场分离合成的结果进行对比，同时生成了全波场进行分离后的模拟记录(图3-1中包含下行P波和上行SV波，图3-2中包含下行SV波和上行P波)。

图4-1和4-2是上行波场进行分离后的模拟记录，图5-1和5-2是下行波场分离后的模拟记录。应用本方法分别对理论的全波场记录，上行波记录和下行波记录做了测试。图6-1及6-2是全波场记录的测试结果，从分析效果看，本发明方法把六种波，按它们的偏振象限，较好的分离在两个剖面上。图7-1及7-2是用只有上行波场的记录分离的结果，在除去了下行波的干扰后，在多种波混合的区域分离效果有了一定改善，分离效果较好。图8-1及8-2是只有下行波场的记录分离的结果，由于下行波场比较单一，只有P波和SV波，没有上行波和转换波，所以两种波分离的很彻底。测试结果证明本发明方法对于波场分离以及矢量波场还原均有较好的效果。

上述技术方案只是本发明的一种实施方式，对于本领域内的技术人员而言，在本发明公开了应用方法和原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施方式所描述的方法，因此前面描述的方式只是优选的，而并不具有限制性的意义。

Claims (1)

1.一种基于VSP时变偏振分析的矢量波场合成及分离方法，其特征在于：所述方法包括以下步骤：

(1)读入水平H_p分量和垂直z分量这两个记录；

(2)在频率域进行时窗自动选取得到偏振分析时窗；

(3)在步骤(2)得到的偏振分析时窗内利用协方差矩阵法判断所述两个记录的波场类型，如果判断结果是单一波场，则采用能量准则法进行偏振分析,得到分离及合成后的P波和SV波的波场，如果判断结果是混合波场，则采用Cho和Spencer法进行偏振分析，得到分离及合成后的P波和SV波的波场；

其中，所述步骤(2)具体如下：选取高斯函数作为时窗函数，给定N进行时窗的自动选取得到偏振分析时窗，然后在此偏振分析时窗内分别对步骤(1)读入的两个记录进行时频域短时离散傅里叶变换，将两个记录从时域转化为频域；

所述高斯函数为：

$g^0\left(t\right){=2}^{1/4}{e}^{{-\mathrm{\πt}}^2}$

所述时频域短时离散傅里叶变换采用短时离散傅里叶变换公式：

$STFT(m,{\mathrm{\ω}}_k)=\underset{n=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\mathrm{\Σ}}}x\left(n\right)g*(n-mN)e^{-j\frac{\mathrm{\π}}{M}nk}$

式中：频率k为整数，k＝0至M，x(n)为给定的水平H_p分量或垂直z分量，n＝0至L-1，L为数据的长度,N为时间轴上的时窗函数移动的步长，M是将频域的一个周期2π分成了M个点，m是任意正整数，是时窗函数移动步长的系数，*为复共轭。