CN102062873A - 一种纵横波匹配方法 - Google Patents

Abstract

本发明提供的纵横波匹配方法包括对纵波地震数据和转换横波地震数据进行匹配，所述匹配为反射时间匹配、反射能量匹配和波形匹配中的至少一者，其中反射时间匹配包括使用时间校正系数α对转换横波剖面的反射时间进行校正，以使得校正后的转换横波剖面的反射时间能够与相应的纵波剖面的反射时间匹配；反射能量匹配包括使用能量校正系数β对转换横波的反射能量进行校正，以使得校正后的转换横波的反射能量能够与相应的纵波的反射能量匹配；波形匹配包括使用计算出的最优维纳滤波器对转换横波的波形进行校正，以使得校正后的转换横波的波形能够与纵波的波形匹配。本发明提供的纵横波匹配方法精度高且效率高。

Description

一种纵横波匹配方法

技术领域

本发明涉及地震勘探领域，尤其涉及多分量地震勘探中的地震数据处理。

背景技术

在多波多分量地震勘探中，通常以纵波激发，采用多分量检波器在每个接收点上记录反射的地震波信号，进而对接收到的地震波数据进行进一步处理和分析。多波多分量地震勘探之所以能够提供比纯纵波或纯横波勘探更丰富、更准确的地质信息，是由于同时获得了反射纵波和转换横波地震数据。由于纵波反射之后产生的纵波和转换横波在反射时间、反射能量和波形上均存在明显差异，因此需要在同一数据范围内对纵波和转换横波在反射时间、反射能量和波形方面进行匹配处理，才能利用它们进行后续的地质解释和油气预测等。

现有的纵横波匹配方法主要有两种：一种是依据纵波层位和横波层位利用手工拉伸或压缩来实现纵横波匹配；另一种是手工拾取Gamma谱(纵横波速度比值)，逐点完成纵横波匹配。这两种匹配方法都是通过人为手工实现，主观因素影响极大、精度低、效率低，且只能对大套的标志性地层进行匹配处理，无法对薄层、复杂小断层等进行匹配，并且也无法实现海量地震数据的精确匹配。

发明内容

本发明的目的是针对现有技术中的纵横波匹配方法均通过手工实现，精度低、效率低等问题，提供一种高精度、高效率的纵横波匹配方法。

本发明提供的纵横波匹配方法包括对纵波地震数据和转换横波地震数据进行匹配，所述匹配为反射时间匹配、反射能量匹配和波形匹配中的至少一者。

其中所述反射时间匹配包括以下步骤：确定纵波传播速度V_P与转换横波传播速度V_S的比γ₀，以得到时间校正系数α，其中

以及使用所述时间校正系数α对转换横波剖面的反射时间进行校正，以使得校正后的转换横波剖面的反射时间能够与相应的纵波剖面的反射时间匹配。

所述反射能量匹配包括以下步骤：计算纵波反射系数R_PP和转换横波反射系数R_PS，以得到能量校正系数β，其中

以及使用所述能量校正系数β对转换横波的反射能量进行校正，以使得校正后的转换横波的反射能量能够与相应的纵波的反射能量匹配。

所述波形匹配包括以下步骤：计算最优维纳滤波器，该最优维纳滤波器为将转换横波地震数据作为输入并使得该滤波器的实际输出与所述纵波地震数据之间的均方误差最小的维纳滤波器；以及使用计算出的最优维纳滤波器对转换横波的波形进行校正，以使得校正后的转换横波的波形能够与纵波的波形匹配。

根据本发明提供的纵横波匹配方法，用于反射时间匹配的时间校正系数α是通过确定精确的纵横波传播速度之比而计算得出的，用于反射能量匹配的能量校正系数β是利用实际的纵横波地震数据通过公式计算纵横波反射系数而获得的，并且用于校正转换横波波形、使之与纵波波形相匹配的最优维纳滤波器也是通过求解维纳-霍夫方程实现的，本发明提供的纵横波匹配方法均是基于实际的地震数据通过计算实现的，未受人为主观判断的影响，与现有技术中通过主观判断采用人为手工实现的匹配方法相比，精确性和效率都得到了很大的提高。

附图说明

图1是本发明提供的纵横波匹配方法中纵波入射之后的反射及透射示意图；以及

图2是本发明提供的纵横波匹配方法中采用的维纳滤波的基本原理图。

具体实施方式

下面结合附图对本发明提供的纵横波匹配方法做进一步的详细描述。

本发明提供的纵横波匹配方法包括对纵波地震数据和转换横波地震数据进行匹配，所述匹配为反射时间匹配、反射能量匹配和波形匹配中的至少一者，其中所述反射时间匹配包括以下步骤：确定纵波传播速度V_P与转换横波传播速度V_S的比γ₀，以得到时间校正系数α，其中

以及使用所述时间校正系数α对转换横波剖面的反射时间进行校正，以使得校正后的转换横波剖面的反射时间能够与相应的纵波剖面的反射时间匹配；所述反射能量匹配包括以下步骤：计算纵波反射系数R_PP和转换横波反射系数R_PS，以得到能量校正系数β，其中

以及使用所述能量校正系数β对转换横波的反射能量进行校正，以使得校正后的转换横波的反射能量能够与相应的纵波的反射能量匹配；所述波形匹配包括以下步骤：计算最优维纳滤波器，该最优维纳滤波器为将转换横波地震数据作为输入并使得该滤波器的实际输出与所述纵波地震数据之间的均方误差最小的维纳滤波器；以及使用计算出的最优维纳滤波器对转换横波的波形进行校正，以使得校正后的转换横波的波形能够与纵波的波形匹配。为了实现更精确的纵横波匹配，优选情况下，所述匹配为反射时间匹配、反射能量匹配和波形匹配。

下面分别从反射时间匹配、反射能量匹配和波形匹配三个方面对本发明提供的纵横波匹配方法进行详细描述。

1.反射时间匹配

纵波在地质界面入射之后的反射波包括反射纵波和转换横波，为本领域技术人员所公知。由于反射得到的转换横波的传播速度低于反射纵波的传播速度，导致同一反射层在转换横波剖面和纵波剖面上的到达时间不一致。为了利用纵波和反射横波进行联合反演以及地质解释、油气预测等，需要对纵波和转换横波进行反射时间匹配，使得反映同一地质层位的纵波和转换横波在时间轴上处于相同的时间位置。在本领域中，对纵波和转换横波的时间匹配指的是以纵波剖面的到达时为标准，校正转换横波剖面时间，使校正后的转换横波剖面时间与纵波剖面时间在同一时间轴上对齐，为本领域技术人员所公知。

根据本发明提供的方法，假设纵波垂直入射，距地震剖面为零偏移，纵波的传播时间为：

$t_{\mathrm{PP}}=\frac{2H}{V_P}---\left(1\right)$

在式(1)中，t_PP为纵波的传播时间，H为纵波传播垂直距离，V_P为纵波的传播速度。

转换横波的传播时间为：

$t_{\mathrm{PS}}=\frac{H}{V_P}+\frac{H}{V_S}---\left(2\right)$

在式(2)中，t_PS为转换横波的传播时间，V_S为转换横波的传播速度。

由式(1)和式(2)可得：

$\frac{t_{\mathrm{PP}}}{t_{\mathrm{PS}}}=\frac{2}{1+\frac{V_P}{V_S}}---\left(3\right)$

令时间校正系数

由式(3)可知，只要得到准确的纵横波传播速度比

由于

就能得到时间校正系数α。从而可以通过该系数α将转换横波剖面的反射时间校正到纵波剖面的反射时间，实现对纵波传播时间t_PP和转换横波传播时间t_PS的匹配。通常在离散的地震数据序列中，可以通过模拟退火反演法来获取精确的纵横波传播速度比

模拟退火反演法是本领域中求取最优解的常规非线性反演方法，为本领域技术人员所公知，例如在《地球物理反演的理论与方法》(1996，杨文采著)专著中有详细描述。在本发明中，通过将采用人机交互解释拾取的粗略的纵横波速度比作为该算法迭代的起点(初始模型值x₀)，以迭代的方式反演可以获得目标模型值x_N，该值即为所要求取的精确的垂直纵横波速度比γ₀。采用模拟退火反演法确定γ₀的方法原理如下：

首先，利用纵波数据与匹配后转换横波数据互相关能量最大准则，确定目标函数E(x_N)：

$E\left(x_N\right)=\left|\frac{\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{y}_n{y}_n^{\′}}{\sqrt{\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{y}_n^2\underset{n=0}{\overset{N}{\mathrm{\Σ}}}{y}_n^{\′2}}}\right|$

上式中，y_n为匹配时窗内的纵波数据，N为时窗长度，y′_n为给定初始γ₀后进行时间和能量匹配后的转换波数据。反演的具体步骤如下：

①给定模型的每一参数(匹配时窗内每一样点处γ₀值)变化范围，在这个范围内随机选择一个初始模型x₀(即给匹配时窗内的每一样点赋一初始值)，并计算相应的目标函数值E(x₀)；

②对当前模型x₀进行扰动产生一个新模型x，计算相应的目标函数值E(x)，得到新模型目标函数与当前模型目标函数之差ΔE＝E(x)-E(x₀)；

③若ΔE＜0，则新模型x被接受；若ΔE＞0，则新模型x按概率P＝exp(-ΔE/T)进行接受，其中T为温度。当模型被接受时，置x₀＝x，E(x₀)＝E(x)

④在温度T下，重复一定次数(次数根据模型被接受的实际情况而设定，当模型易被接受时，次数可减少，反之则增加)的扰动和接受过程，即重复步骤②和③；

⑤缓慢地降低温度T；

⑥重复步骤②-⑤，直至收敛条件满足为止(即目标函数的值小于或等于事先给定的精度(例如1×10^-6)为止)，则当收敛条件满足时得到的目标模型值x_N即是所要求取的精确的纵横波传播速度比

在以上模拟退火反演γ₀过程中涉及三个方面的技术问题：

①模型扰动

模拟退火中新模型的产生是对当前模型进行扰动得到的，可采用依赖于温度的似Cauchy分布产生新模型，即

$\left.\begin{array}{c}{x}_i^{\′}=x_i+y_i(B_i-A_i)\\ y_i=\mathrm{T\; sgn}(u-0.5)[{(1+1/T)}^{|2u-1|}-1]\end{array}\right\}$

式中，x′_i为扰动后的模型，x_i为当前模型中的第i个变量；方程组中第二个方程中的y_i为第一个方程的中间系数表达式，u为[0，1]均匀分布的随机数；[A_i，B_i]为x_i的取值范围，且要求扰动后的x′_i∈[A_i，B_i]；sgn(·)为符号函数。

②接受概率

非常快速模拟退火算法用广义Gibbs分布产生接受概率P，即

P＝[1-(1-h)ΔE/T]^1/(1-h)

在上式中，ΔE为新模型目标函数与当前模型目标函数之差；T为温度；h为实数。

③降温方式

Ingber(1989)给出了非常快速模拟退火方法的降温方式：

T(K)＝T₀exp(-CK^1/N)

式中，T₀为初始温度；K为迭代次数；C为给定常数；N为待反演的参数个数。它可以改写为

$T\left(K\right)=T_0{\mathrm{\α}}^{K^{1/N}}$

通常选择0.7≤α≤1，在实际应用中，常采用0.5或1代替式子中的1/N。

由此，将采用模拟退火反演法获得的精确的纵横波速度比γ₀代入上式

即可得到时间校正系数

从而能够使用所述时间校正系数α对转换横波剖面的反射时间进行校正，以使得校正后的转换横波剖面的反射时间与相应的纵波剖面的反射时间匹配。

2.反射能量匹配

当纵波在地质界面入射之后，会同时发生反射、透射和转换。其中，反射波包括反射纵波和转换横波，透射波包括透射纵波和转换横波。

在实际多波多分量地震勘探中，由于纵波和转换横波的传播速度、震动方式等的差异，导致接收到的转换横波分量的能量相对纵波分量的能量而言较弱，且后续对多波地震数据的处理也将进一步增强纵波和转换横波能量之间的差异。因此，需要在联合应用纵横波资料之前对转换横波和纵波的能量进行匹配。在本领域中，对纵波和转换横波的能量匹配指的是以纵波的反射能量为标准，校正转换横波的反射能量，使校正后的转换横波的反射能量与纵波的反射能量相匹配，为本领域技术人员所公知。

图1为纵波入射之后的反射及透射示意图。如图1所示，P为入射纵波，PP为反射纵波，PS为反射转换横波，PP’为透射纵波，PS′为透射转换横波；其中，θ₁为纵波入射角，θ₂为纵波透射角，则纵波平均反射角θ＝(θ₁+θ₂)/2；φ₁为转换横波反射角，φ₂为转换横波透射角，则转换横波近似反射角φ＝sin^-1(γsinθ)；V_S1为转换横波在介质1中的传播速度，V_S2为转换横波在介质2中的传播速度，则转换横波平均速度V_S＝(V_S1+V_S2)/2，转换横波速度差ΔV_S＝V_S2-V_S1；V_P1为纵波在介质1中的传播速度，V_P2为纵波在介质2中的传播速度，则纵波平均速度V_P＝(V_P1+V_P2)/2，纵波速度差ΔV_P＝V_P2-V_P1，转换横波与纵波的速度比γ＝V_S/V_P；ρ₁为介质1的密度，ρ₂为介质2的密度，则介质密度差Δρ＝ρ₂-ρ₁，地层介质平均密度ρ＝(ρ₁+ρ₂)/2。

在多波多分量地震勘探中，对于任意地层界面，当上下地层的纵横波传播速度与地层介质平均密度一定时，由Zoeppritz方程式(4)可以推导出式(5)中的纵波反射系数和转换横波反射系数：

$\left[\begin{array}{c}{R}_{\mathrm{pp}}\\ R_{\mathrm{ps}}\\ T_{\mathrm{pp}}\\ T_{\mathrm{ps}}\end{array}\right]={\left[\begin{array}{cccc}-\sin {\mathrm{\θ}}_1& -\cos {\mathrm{\φ}}_1& \sin {\mathrm{\θ}}_2& \cos {\mathrm{\φ}}_2\\ \cos {\mathrm{\θ}}_1& -\sin {\mathrm{\φ}}_1& \cos {\mathrm{\θ}}_2& -\sin {\mathrm{\φ}}_2\\ \sin {2\mathrm{\θ}}_1& \frac{V_{P1}}{V_{S1}}\cos {2\mathrm{\φ}}_1& \frac{{\mathrm{\ρ}}_2{V}_{S2}^2{V}_{P1}}{{\mathrm{\ρ}}_1{V}_{S1}^2{V}_{P2}}\cos {2\mathrm{\φ}}_1& \frac{{\mathrm{\ρ}}_2{V}_{S2}{V}_{P1}}{{\mathrm{\ρ}}_1{V}_{S1}^2}\cos {2\mathrm{\φ}}_2\\ -\cos {2\mathrm{\φ}}_1& \frac{V_{S1}}{V_{P1}}\sin {2\mathrm{\φ}}_1& \frac{{\mathrm{\ρ}}_2{V}_{P2}}{{\mathrm{\ρ}}_1{V}_{P1}}\cos 2{\mathrm{\φ}}_2& -\frac{{\mathrm{\ρ}}_2{V}_{S2}}{{\mathrm{\ρ}}_1{V}_{P1}}\sin {2\mathrm{\φ}}_2\end{array}\right]}^{-1}\left[\begin{array}{c}\sin {\mathrm{\θ}}_1\\ \cos {\mathrm{\θ}}_1\\ \sin 2{\mathrm{\θ}}_1\\ \cos {2\mathrm{\φ}}_1\end{array}\right]---\left(4\right)$

$\left\{\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)\≈\frac{1+{\tan }^2\mathrm{\θ}}{2}(\frac{\mathrm{\Δ}{V}_P}{V_P}+\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}})-\frac{8{\mathrm{\γ}}^2{\sin }^2\mathrm{\θ}}{2}(\frac{\mathrm{\Δ}{V}_S}{V_S}\\ +\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}})+(2{\mathrm{\γ}}^2{\sin }^2\mathrm{\θ}-\frac{{\tan }^2\mathrm{\θ}}{2})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ R_{\mathrm{PS}}(\mathrm{\θ},\mathrm{\φ})\≈\frac{\tan \mathrm{\φ}}{2\mathrm{\γ}}(4{\sin }^2\mathrm{\φ}-4\mathrm{\γ}\cos \mathrm{\θ}\cos \mathrm{\φ})(\frac{\mathrm{\Δ}{V}_S}{V_S}\\ +\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}})-\frac{\tan \mathrm{\φ}}{2\mathrm{\γ}}(1+2{\sin }^2\mathrm{\φ}-2\mathrm{\γ}\cos \mathrm{\θ}\cos \mathrm{\φ})\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}\right.---\left(5\right)$

其中，R_PP为纵波反射系数，R_PS为转换横波反射系数，T_PP为纵波透射系数，T_PS为转换横波反射系数。

由于地层界面的纵波反射系数R_PP与转换横波反射系数R_PS均为确定值，反射系数的大小仅与上下地层弹性参数及纵波入射角有关，且由于地震道是地层反射系数与子波的褶积，由此纵波反射系数R_PP与转换横波反射系数R_PS的比例关系等于纵波反射能量与转换横波反射能量的比例关系。令P_PP为纵波反射能量，P_PS为转换横波反射能量，则能量校正系数

由此，只要计算出纵波反射系数R_PP和转换横波反射系数R_PS，就能确定纵波和转换横波反射能量之间的对应关系，得到能量校正系数β，从而可以通过该系数β对转换横波的反射能量进行校正，以使得校正后的转换横波的反射能量与相应的纵波的反射能量匹配。

可以利用已知的地层弹性参数，通过波场正演模拟来确定纵波和转换横波反射能量之间的对应关系。

首先，利用已知测井数据建立地层模型；根据实际资料中的主频及子波特征，采用数值模拟方法，分别计算纵横波叠前道集；计算纵波CMP道集与转换波CCP道集的总能量，其计算方法如下：

利用实际测井数据和模拟数据计算纵波反射系数R_PP和转换横波反射系数R_PS：

$R_{\mathrm{PP}}(t_p,H_p)=\frac{\underset{0}{\overset{t_{p_{\max }}}{\∫}}\left|S_{\mathrm{mode}\mathrm{lPP}}({t_p}^{\′},H_p){{\mathrm{dt}}_p}^{\′}\right|}{\underset{0}{\overset{t_{p_{\max }}}{\∫}}\left|S_{\mathrm{dataPP}}({t_p}^{\′},H_p){{\mathrm{dt}}_p}^{\′}\right|}{S}_{\mathrm{dataPP}}(t_p,H_p)$

$R_{\mathrm{PS}}(t_s,H_s)=\frac{\underset{0}{\overset{t_{s_{\max }}}{\∫}}\left|S_{\mathrm{mode}\mathrm{lPS}}({t_s}^{\′},H_s){{\mathrm{dt}}_s}^{\′}\right|}{\underset{0}{\overset{t_{s_{\max }}}{\∫}}\left|S_{\mathrm{dataPS}}({t_s}^{\′},H_s){{\mathrm{dt}}_s}^{\′}\right|}{S}_{\mathrm{dataPS}}(t_s,H_s)---\left(6\right)$

在式(6)中，t_p、t_s分别是纵波和转换横波的双程旅行时，t_p′、t_s′分别是对时间的积分变量，t_pmax、t_smax分别是纵波和转换横波的匹配时窗内的最大时间长度，H_p、H_s分别是纵波和转换横波的偏移距，S_modelPP、S_modelPS分别是利用测井资料合成的纵波和转换横波叠前道集理论数据、S_dataPP、S_dataPS分别是纵波和转换横波叠前道集实际数据。

根据采用公式(6)计算得到的纵波反射系数R_PP和转换横波反射系数R_PS，能够确定能量校正系数β，利用该系数β对CCP道集各样点进行标量乘，则可以得到能够与相应的CMP道集能量相匹配的CCP道集能量，从而实现转换波与纵波的反射能量匹配。

3.波形匹配

纵波与转换横波地震数据经传播时间匹配处理后，将引起转换横波地震数据的波形变化。这种波形变化对多分量地震数据的影响主要表现为转换横波地震数据的振幅谱异常。假设地层反射系数的频谱满足“白噪”假设条件，则纵波与转换横波经传播时间匹配处理后，转换横波地震数据的振幅谱异常可用转换横波记录子波变化来表示。由此，为了使得纵波与转换横波地震数据波形一致，以利用它们进行后续的地质解释和油气预测等，就需要对转换横波地震数据子波及相位进行整形处理。在本领域中，对纵波和转换横波的波形匹配指的是以纵波的地震数据波形为标准，校正转换横波的地震数据波形，使校正后的转换横波的地震数据波形与纵波的地震数据波形一致，为本领域技术人员所公知。

本发明提供的方法中，波形匹配是将纵波地震数据作为希望输出，将转换横波地震数据作为输入，通过褶积模型与最小能量准则，采用维纳(Wiener)最佳滤波器设计方法，找出最优滤波器，使得实际输出与希望输出之间的均方误差为最小，由此实现转换横波的地震数据波形与纵波的波形的匹配。

滤波是指从连续的(或离散的)输入数据中滤除噪声和干扰以提取有用信息的过程，其相应的装置称为滤波器。最佳滤波器是指能够根据某一最佳准则进行滤波的滤波器。假定线性滤波器的输入为有用信号和噪声之和，两者均为广义平稳过程且知它们的二阶统计特性，根据最小均方误差准则(滤波器的输出信号与需要信号之差的均方值最小)，求得最佳线性滤波器的参数的滤波器称为维纳滤波器，为本领域技术人员所公知。

设维纳滤波器的输入为含噪声的随机信号，期望输出与实际输出之间的差值为误差，对该误差求均方，即为均方误差。因此均方误差越小，噪声滤除效果就越好。为了使均方误差最小，关键在于求脉冲响应。如果能够满足维纳-霍夫(wiener-Hoof)方程，就可使维纳滤波器达到最佳。根据维纳-霍夫方程，最佳维纳滤波器的脉冲响应完全由输入自相关函数以及输入与期望输出的互相关函数所决定。

因此，本发明提供的波形匹配处理包括两个部分，即波形匹配与地震数据的相位匹配处理。

详细的实现原理如下：

在一定条件下，地震数据记录时间序列的总体平均可以用时间平均来代替。如果平稳过程(其统计特性不随时间而变的随机过程)对任一函数的有限时间平均或有限样本(以概率为1)收敛于总体平均，则称这个过程为各态历经过程或称为遍历过程。根据维纳-辛钦定理(各态历经定理)，对于平稳的各态历经过程，其自相关函数为其功率谱密度的傅里叶反变换。运用这个定理可以得到，在一定条件下，处在统计平衡的时间序列的时间平均等于相平均。有了这个前提，就可以从时间序列的过去数据推知未来和预测未来。如前所述，滤波问题就是尽可能地恢复一个被噪声干扰了的信息流问题。实质上，就是预测一个被噪声或其他因素干扰的时间序列的问题，而维纳滤波所解决的就是在最小均方误差准则下的滤波问题。

设计维纳滤波器的过程就是寻求最佳滤波器脉冲响应的明确表达式，其实质就是求解维纳-霍夫方程。通过根据概率论中的协方差函数描述的滤波器中的不同信号的统计特性，在滤波器的随机输入(包括信号与噪声)条件下，可以通过维纳-霍夫方程求解。这个方程是积分方程，其解是维纳滤波器的最优脉冲响应。

求解维纳-霍夫方程，即求解最佳维纳滤波器的过程是现有数字信号处理技术中的惯用技术，为本领域技术人员所公知。例如在《现代信号处理》(2003，张贤达著)专著中有详细描述。根据本发明提供的方法，求解维纳-霍夫方程以找到最佳维纳滤波器的过程如下：

①最佳维纳滤波求解关系

输入信号：x(t)＝(x(0)，x(1)，…，x(n-1))

滤波因子：h(t)＝(h(0)，h(1)，…，h(m))

实际输出：y(t)＝h(t)*x(t)＝∑_τh(τ)x(t-τ)

其中，τ为延迟时。

          y(t)＝(y(0)，y(1)，…y(M))，M＝m+n

期望输出：d(t)＝(d(0)，d(1)，…，d(M))

输入误差：e(t)＝d(t)-y(t)

误差能量：Q＝∑_te²(t)＝∑t(y(t)-d(t))²

为了求取适合的滤波因子h(t)，令误差平方和(总误差能量)Q为最小，滤波因子的求解即成为误差能量Q对滤波因子h(t)求极值的问题。

②滤波方程的推导

对滤波因子h(t)求偏导数，化简，并令其为零：

$\frac{\∂Q}{\∂h\left(s\right)}=\frac{\∂}{\∂h\left(s\right)}{\mathrm{\Σ}}_t[{\mathrm{\Σ}}_{t}h\left(\mathrm{\τ}\right)x(t-\mathrm{\τ})-d\left(t\right){]}^2$

$={\mathrm{\Σ}}_t\frac{\∂}{\∂h\left(s\right)}{[{\mathrm{\Σ}}_{\mathrm{\τ}}h\left(\mathrm{\τ}\right)x(t-\mathrm{\τ})-d\left(t\right)]}^2$

由此得

∑_th(t)∑_tx(t-τ)x(t-s)＝∑_td(t)x(t-s)            (7)

其中，S为滤波因子样点序列号(s＝0，1，…，m)

令

r_xxx(t-τ)＝∑_tx(t-τ)x(t-s)                      (8)

r_dxx(t-τ)＝∑_tx(t-τ)x(t-s)                      (9)

由自相关和互相关函数的定义可知，式(8)中r_xx是关于x(t)的自相关，式(9)中r_dx则是关于d(t)与x(t)的互相关。由式(8)和式(9)，方程组(7)可写成

${\mathrm{\Σ}}_{\mathrm{\τ}=0}^m{r}_{\mathrm{xx}}(\mathrm{\τ}-s)h\left(\mathrm{\τ}\right)=r_{\mathrm{dx}}\left(s\right),$ (s＝0，1，…，m)                        (10)

由于自相关函数r_xx为一偶函数，则式(10)可以写成矩阵形式

得到滤波因子h(t)之后，输入信号x(t)褶积，就可得到最佳输出。这一过程如图2所示。

式(11)中左端自相关矩阵为托不里兹(Tpeplitz)矩阵，它是信号处理中的一个常用矩阵，为本领域技术人员所公知。该托不里兹矩阵的各元素为实型值，对角元素相等，其它元素对称于对角线。该方程可用莱文森递归算法快速求解。

③最小能量误差的计算

将求取误差能量的方程Q＝∑_te²(t)＝∑_t(y(t)-d(t))²表示为

$Q_{\min }={\mathrm{\Σ}}_t{d}_t^2-2\left({\mathrm{\Σ}}_t{\mathrm{\Σ}}_{\mathrm{\τ}}{d}_t{h}_{\mathrm{\τ}}{x}_{t-\mathrm{\τ}}\right)+{\mathrm{\Σ}}_t{\left({\mathrm{\Σ}}_{\mathrm{\τ}}{h}_{\mathrm{\τ}}{x}_{t-\mathrm{\τ}}\right)}^2---\left(12\right)$

Q_min代表最小能量误差。

将关系式

∑_td_tx_t-τ＝r_dx                    (13)

∑_tx_tx_t-τ＝r_xx                    (14)

带入方程(12)中，得到

$Q_{\min }={\mathrm{\Σ}}_t{d}_t^2-2{\mathrm{\Σ}}_{\mathrm{\τ}}{h}_{\mathrm{\τ}}{r}_{\mathrm{dx}}+{\mathrm{\Σ}}_{\mathrm{\τ}}{h}_{\mathrm{\τ}}{\mathrm{\Σ}}_i{h}_i{\mathrm{\Σ}}_t{x}_{t-\mathrm{\τ}}{x}_{t-i}---\left(15\right)$

$Q_{\min }={\mathrm{\Σ}}_t{d}_t^2-2{\mathrm{\Σ}}_{\mathrm{\τ}}{h}_{\mathrm{\τ}}{r}_{\mathrm{dx}}+{\mathrm{\Σ}}_{\mathrm{\τ}}{h}_{\mathrm{\τ}}{\mathrm{\Σ}}_i{h}_i{\mathrm{\Σ}}_t{h}_i{x}_{\mathrm{\τ}-i}---\left(16\right)$

最后利用方程(10)，得到

$Q_{\min }={\mathrm{\Σ}}_t{d}_t^2-{\mathrm{\Σ}}_{\mathrm{\τ}}{h}_{\mathrm{\τ}}{r}_{\mathrm{dx}}---\left(17\right)$

④莱文森递归算法

为了确定维纳滤波系数，需要求解正规方程(11)，采用莱文森递归法可以快速求解该方程中的托不里兹矩阵。莱文森递归法是先做一个两点滤波器，由此推导出一个三点滤波器，如此继续下去，直至推导出n点滤波器，莱文森递归法是信号处理中求解托不里兹方程组的常用算法，为本领域技术人员所公知。下面分析期望输出是零之后脉冲d_t＝(1，0，0，…)的情况下利用莱文森递归算法求解托不里兹方程的过程。

假定输入信号序列x_t＝(x₀，x₁，x₂，...)，方程(11)可以写成下面形式，即

等式两边同除以h(0)，得

式中，a_i＝h(i)/h(0)，i＝1，2，…，n-1；

v＝x₀/h(0)。

根据方程(19)求解未知量v和未知滤波器系数(a₁，a₂，...，a_n-1)。由于期望输出为零滞后脉冲，因此滤波器(1，a₁，a₂，...，a_n-1)表示了脉冲反褶积过程。采用莱文森递归算法求解方程(19)，先从两项滤波器(1，a₁)开始，接着求解滤波器(1，a₁，a₂)，依此类推。当n＝2时，方程(19)为

$\left[\begin{array}{cc}{r}_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)\\ r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1\end{array}\right]=\left[\begin{array}{c}v\\ 0\end{array}\right]---\left(20\right)$

写出联立方程

r_xx(0)+r_xx(1)a₁＝v                        (21)

和

r_xx(1)+r_xx(0)a₁＝0                        (22)

第一次迭代求得滤波器系数a₁和未知变量v，即

a₁＝-r_xx(1)/r_xx(0)                        (23)

和

v＝r_xx(0)+r_xx(1)a₁                    (24)

现在将方程(20)写成三项滤波器(1，a’₁，a’₂)，即

$\left[\begin{array}{ccc}{r}_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(2\right)\\ r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)\\ r_{\mathrm{xx}}\left(2\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1^{\′}\\ a_2^{\′}\end{array}\right]=\left[\begin{array}{c}{v}^{\′}\\ 0\\ 0\end{array}\right]---\left(25\right)$

这一步骤将用到前面各步骤的结果(方程(21)和(22))求解滤波系数(a’₁，a’₂)。

首先，用上面建立联立方程组的方法，重写方程(20)为：

$\left[\begin{array}{ccc}{r}_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(2\right)\\ r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)\\ r_{\mathrm{xx}}\left(2\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ 0\end{array}\right]=\left[\begin{array}{c}v\\ 0\\ e\end{array}\right]---\left(26\right)$

式中

e＝r_xx(2)+r_xx(1)a₁                    (27)

交换1、3行顺序，将方程变形为：

$\left[\begin{array}{ccc}{r}_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(2\right)\\ r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)\\ r_{\mathrm{xx}}\left(2\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}0\\ a_1\\ 1\end{array}\right]=\left[\begin{array}{c}e\\ 0\\ v\end{array}\right]---\left(28\right)$

方程(28)两边都乘以变量c，然后从方程(26)中减去上述结果得

$\left[\begin{array}{ccc}{r}_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(2\right)\\ r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)& r_{\mathrm{xx}}\left(1\right)\\ r_{\mathrm{xx}}\left(2\right)& r_{\mathrm{xx}}\left(1\right)& r_{\mathrm{xx}}\left(0\right)\end{array}\right]\left[\begin{array}{c}1\\ a_1-c^{a_1}\\ -c\end{array}\right]=\left[\begin{array}{c}v-\mathrm{ce}\\ 0\\ e-\mathrm{vc}\end{array}\right]---\left(29\right)$

最后，比较方程(25)和(29)，可发现

$\left[\begin{array}{c}1\\ a_1^{\′}\\ a_2^{\′}\end{array}\right]=\left[\begin{array}{c}1\\ a_1-c^{a_1}\\ -c\end{array}\right]---\left(30\right)$

和

$\left[\begin{array}{c}{v}^{\′}\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}v-\mathrm{ce}\\ 0\\ e-\mathrm{vc}\end{array}\right]---\left(31\right)$

求解c和v′(变形后的未知变量)得

$v^{\′}=v[1-{\left(\frac{e}{v}\right)}^2]---\left(32\right)$

和

C＝e/v                            (33)

因此新的滤波系数(a′₁，a′₂)为方程(30)

$a_1^{\′}=a_1-\frac{c}{v}{a}_1---\left(34\right)$

和

$a_2^{\′}=-\frac{c}{v}---\left(35\right)$

重复递推模式，确定下一组滤波器系数：

a.用输入序列的自相关滞后和当前的滤波器系数(方程(28))计算v和e；

b.计算下一组滤波器系数(方程(32))；

c.计算新的v值(方程(34))。

这一递推过程得到期望长度为n的维纳滤波器(1，a₁，a₂，...，a_n-1)，该维纳滤波器即为使得实际输出与希望输出之间的均方误差为最小的最佳滤波器。由此，通过将纵波地震记录作为希望输出，将转换横波地震记录作为输入，采用得到的最佳滤波器，可以实现纵波与转换横波地震数据的波形匹配。

Claims (5)

1.一种纵横波匹配方法，该方法包括对纵波地震数据和转换横波地震数据进行匹配，所述匹配为反射时间匹配、反射能量匹配和波形匹配中的至少一者，其中：

所述反射时间匹配包括以下步骤：

确定纵波传播速度V_P与转换横波传播速度V_S的比

以得到时间校正系数α，其中

以及

使用所述时间校正系数α对转换横波剖面的反射时间进行校正，以使得校正后的转换横波剖面的反射时间能够与相应的纵波剖面的反射时间匹配；

所述反射能量匹配包括以下步骤：

计算纵波反射系数R_PP和转换横波反射系数R_PS，以得到能量校正系数β，其中

以及

使用所述能量校正系数β对转换横波的反射能量进行校正，以使得校正后的转换横波的反射能量能够与相应的纵波的反射能量匹配；

所述波形匹配包括以下步骤：

计算最优维纳滤波器，该最优维纳滤波器为将转换横波地震数据作为输入并使得该滤波器的实际输出与所述纵波地震数据之间的均方误差最小的维纳滤波器；以及

使用计算出的最优维纳滤波器对转换横波的波形进行校正，以使得校正后的转换横波的波形能够与纵波的波形匹配。

2.根据权利要求1所述的方法，其中，所述匹配为反射时间匹配、反射能量匹配和波形匹配。

3.根据权利要求1所述的方法，其中纵波传播速度V_P与转换横波传播速度V_S的比

采用模拟退火反演法进行确定。

4.根据权利要求1所述的方法，其中纵波反射系数R_PP和转换横波反射系数R_PS通过下式计算获得：

$R_{\mathrm{PP}}(t_p,H_p)=\frac{\underset{0}{\overset{t_{p_{\max }}}{\∫}}\left|S_{\mathrm{mode}\mathrm{lPP}}({t_p}^{\′},H_p){{\mathrm{dt}}_p}^{\′}\right|}{\underset{0}{\overset{t_{p_{\max }}}{\∫}}\left|S_{\mathrm{dataPP}}({t_p}^{\′},H_p){{\mathrm{dt}}_p}^{\′}\right|}{S}_{\mathrm{dataPP}}(t_p,H_p)$

$R_{\mathrm{PS}}(t_s,H_s)=\frac{\underset{0}{\overset{t_{s_{\max }}}{\∫}}\left|S_{\mathrm{mode}\mathrm{lPS}}({t_s}^{\′},H_s){{\mathrm{dt}}_s}^{\′}\right|}{\underset{0}{\overset{t_{s_{\max }}}{\∫}}\left|S_{\mathrm{dataPS}}({t_s}^{\′},H_s){{\mathrm{dt}}_s}^{\′}\right|}{S}_{\mathrm{dataPS}}(t_s,H_s)$

其中t_p、t_s分别是纵波和转换横波的双程旅行时，t_p′、t_s′分别是对时间的积分变量，t_pmax、t_smax分别是纵波和转换横波的匹配时窗内的最大时间长度，H_p、H_s分别是纵波和转换横波的偏移距，S_modelPP、S_modelPS分别是纵波和转换横波叠前道集理论数据、S_dataPP、S_dataPS分别是纵波和转换横波叠前道集实际数据。

5.根据权利要求1所述的方法，其中所述最优维纳滤波器的计算是通过求解维纳-霍夫方程实现的。

Images