CN104932016A - 一种Viterbi-BMC算法的自动速度分析方法 - Google Patents

Abstract

本发明公开了一种Viterbi-BMC算法的自动速度分析方法，在一个确定的t₀下，采用短排列CMP道集，进行常规速度分析，得到一个初始的动校正速度V_nmo；将V_nmo和t₀代入非双曲时差公式，对大排列地震数据进行扫描，得到一个初始的η值代入非双曲时差公式，对大排列地震数据进行扫描，重新确定动校正速度V_nmo；反复迭代V_nmo和η两到三次，获得稳定的V_nmo和η；对下一个t₀进行扫描，执行前面迭代过程，直到全部的t₀被扫描完成，输出速度V_nmo谱和η谱。本发明的有益效果是能够处理大排列地震数据，处理结果准确。

Description

一种Viterbi-BMC算法的自动速度分析方法

技术领域

本发明属于地震测量技术领域，涉及一种Viterbi-BMC算法的自动速度分析方法。

背景技术

在勘探地震学中，速度信息是一个非常重要的参数，不仅在地震，测井，岩石物性之间起到了桥梁的作用，而且是整个地震数据处理的基石(NMO，DMO，偏移，时深转换等)。速度分析具有如下的三个作用：a.主要提供叠加速度，从而提高叠加剖面的信噪比和分辨率；b.为偏移和时深转换，VVAZ等提供速度场信息；c.为岩性的划分提高可靠的速度信息。

自Dix(1955)提出基于各项同性介质常规反射P波的双曲时差公式后，该理论一直被地震勘探界广泛的应用。但是随着地震勘探技术的发展，地震数据由短排列向大排列发展，介质的模型由各向同性向各向异性发展，常规速度分析的精度已经不能满足工业界的需求。基于这个问题，我们对于三维的大排列地震数据，采用基于VTI介质的Akhaliaf非双曲速度分析计算，形成Vnmo谱和η谱，然后采用viterbi自动拾取优化算法，拾取v-t曲线，采用BMS动校正法，进行动校正处理，提高对浅层远偏移距的地震数据的信息保留，提高速度分析及动校正的精度和效率。目前常规速度分析技术有：

1.非双曲速度分析；非双曲线速度分析是建立在各向异性地质模型上的一种地震波速度反演方法，在常规速度分析的基础上采用了非双曲时差公式进行迭代反演方法。在传统的速度分析技术中是将地层假设为各向同性的，地震波在水平层状介质中传播时，其反射波的走时曲线是成双曲线方程的。实际地层大多是各向异性的，如果运用常规速度分析对地震资料处理，所获得的速度信息会产生极大的误差，使得动校正无法将远道拉平，从而得不到可靠与满意的效果。根据各向异性介质模型，建立非双曲线速度分析方法，不仅更符合实际地质情况，而且也能得到更精确的速度，从而得到更高精度的地震资料。

2.Viterbi自动拾取；Viterbi算法是基于随机过程的Morkov链理论的，寻找最短路径的一种最优化搜索算法，它被通讯行业里被广泛的应用于编译卷积码。该算法的原理是：如果点和间的最短路径通过中间点，那么，和间的这个路径段也是该路径段间的最短路径。这个算法包括两个步骤：向前最短路径累积(积分)计算和向后递减跟踪。本技术在速度自动拾取的应用中,我们称为“向前最大速度谱能量团累积计算和向后递减跟踪”。

3.BMS动校正法；自从Buchholtz首次揭示常规NMO校正的拉伸效应以及Dunkin等对其进行研究以来,致力于消除拉伸效应的研究一直十分活跃。Rupert(1969)等提出了称作BMS的校正方法,该方法根据同相轴的分布交错划分一系列数据块,对不同的块实行静态平移,然后对平移后的块再进行合并的动校正技术。

发明内容

本发明的目的在于提供一种Viterbi-BMC算法的自动速度分析方法，解决了目前常规速度分析对于大排列地震数据失效的问题。

本发明所采用的技术方案是按照以下步骤进行：

步骤1：输入地震数据，对于P波地震数据采用非双曲速度分析的双谱算法：

Step one：在一个确定的t₀下，采用短排列CMP道集，进行常规速度分析，得到一个初始的动校正速度V_nmo；

Step two：将Step one中得到初始动校正速度V_nmo和t₀代入非双曲时差公式，对大排列地震数据进行扫描，得到一个初始的η值；

Step three：用初始η值代入非双曲时差公式，对大排列地震数据进行扫描，重新确定动校正速度V_nmo；

Step four：重复执行Step two和Step three反复迭代_Vnmo和η两到三次，获得稳定的V_nmo和η；

Step five：对下一个t₀进行扫描，执行前面四步，直到全部的t₀被扫描完成，输出速度V_nmo谱和η谱。

步骤2：对步骤1中得到速度V_nmo谱和η谱进行自动拾取得到V_nmo-t₀曲线和η-t₀曲线。

步骤3：对步骤2中的v-t曲线和η-t曲线进行线性插值得到新的v-t0曲线和η-t0曲线；

步骤4：基于非双曲时差的BMS动校正法：采用BMS动校正方法进行插值搬移；

先假设在CMP道集中任意反射波脉冲的波形是完全相同的，首先将CMP道集中零偏移距的地震道分成若干时窗长度为BL的有重叠单元，假设每一段的长度都可以容下一个完整反射地震波脉冲，沿着非双曲时差曲线投影到其他的地震道上形成等价单元，然后进行非双曲搬移，对于其它单元重复操作，校平整个道集。

进一步，所述步骤2中，自动拾取的步骤为：

Step one：速度谱和η谱的圆滑处理；

设共中心点道集(CMP)有N道地震记录，该CMP道集数据表达为u(i,j)，根据速度谱计算原理和Taner的相似系数准则可获得速度谱S(t_0,V_nmo)和η谱S(t_0,η)表达如下式：

$S\left(t_0,V_{nmo}\right)=\frac{\underset{j=0}{\overset{M}{\Σ}}\[\underset{i=1}{\overset{N}{\Σ}}u{\left(i,j+r_i\right)}^2\]}{\underset{j=0}{\overset{M}{\Σ}}\underset{i=1}{\overset{N}{\Σ}}u{\left(i,j+r_i\right)}^2}---\left(2\right)$

式中i＝(1,2,3…N)—地震道道号

j＝(0,1,2,…,M)—地震道内采样点序号

—地震波偏移距延迟时间的序号

t_i—偏移距引起的延迟时间

△t—时间采样间隔

S(t_0,V_nmo)—相似系数

t_0,—零偏移距的双程旅行时

V_nmo—动校正速度

现在为了方便表达，设μ(x_i,y_j)＝S(t_0,V_nmo)，式中x_i＝t₀，i＝(0,1,2…m)，y_j＝V_nmo(j＝0,1,…,n)，x_i为零偏移距双程旅行时序列，y_j为动校正速度序列，引入另一个量w(i,j)来表示圆滑后的速度谱其关系如下：

$w\left(i,j\right)=\frac{\underset{jw=w_1}{\overset{w_2}{\Σ}}\underset{iw=l_1}{\overset{l_2}{\Σ}}\mathrm{\μ}\left(x_{iw},y_{jw}\right)}{\left(w_2-w_1\right)\×\left(l_2-l_1\right)}---\left(3\right)$

式中w₁＝i-l_w/2，w₂＝i+l_w/2，l₁＝j-l_w/2，l₂＝j+l_w/2，l_w为滑动窗口的长度；

Step two：使用Viterbi算法进行速度自动拾取；

由Viterbi算法可知，V＝V₁ ^T＝(V₁,V₂,…V_T)为一个所观测的动校正速度的序列，则其联合概率分布总可表示为：

$P\left({V_1}^T\right)=P\left(V_1\right)\underset{t=2}{\overset{T}{\Π}}P\left(V_t\right|{V_1}^{t-1})---\left(4\right)$

式中，P(V₁)为初始动校正速度观测序列，为观测序列条件分布概率，t＝1,2…,T，上式表明：在t时刻，观测序列V_T的条件分布概率依赖于V_T-1前的所有的值，P(V₁ ^t-1)为其先验概率，为其后验概率；

假设在t时刻的动校正速度状态(V_nmo)_t是{1,…,M}内一个有限的数，该状态的向量表示为：

V_nmo＝{(V_nmo)₁}^T＝((V_nmo)₁,(V_nmo)₂,…(V_nmo)_T)       (5)

根据随机过程一阶Markov链理论和Bayes规则，前面的观测动校正速度序列V₁ ^T和状态序列之间的概率关系为：

$P\left(V_t|{\left(V_{nmo}\right)}_1^T,V_1^{t-1}\right)=P{\left(V_t|\left(V_{nmo}\right)\right)}_t---\left(6\right)$

$P\left(V_t|{\left(V_{nmo}\right)}_1^{t+1},{V_1}^t\right)=P\left(V_{t+1}\right|{\left(V_{\mathrm{nmo}}\right)}_t)---\left(7\right)$

式(6)中P(V_t|(V_nmo)_t)为发射概率，即两种序列联合条件分布概率，P(V₁ ^t-1)为先验概率，联合条件分布概率P(V_t|(V_nmo)_t)为后验概率。式(7)中，P(V_t+1|(V_nmo)_t)为传递概率，即状态序列分布概率，P(V₁ ^t)为先验概率，P(V_t+1|(V_nmo)_t)为后验概率；

由式(6)和式(7)可知，当要预测动校正速度序列，动校正速度状态值是过去所有相关的动校正速度状态值和观测序列值的函数，根据独立概率分布的假设可知，观测序列和状态变量的联合概率分布的关系表示为：

$P\left({V_1}^t,{\left(V_{nmo}\right)}_1^T\right)=P\left({\left(V_{nmo}\right)}_1\right)\underset{t=1}{\overset{T-1}{\Π}}P\left({\left(V_{nmo}\right)}_{t+1}|{\left(V_{nmo}\right)}_t\right)\underset{t=1}{\overset{T}{\Π}}P\left(V_t|{\left(V_{nmo}\right)}_t\right)-\left(8\right)$

通过最大概率状态这可以利用下面的最大算法得到：

$q_1^{T*}=\arg \max {|}_{q1}^{T}P<\left({\left(V_{nmo}\right)}_1^T\right|V_1^T)=\arg \max {|}_{q_1}^{T}P\left({\left(V_{nmo}\right)}_1^T,V_1^T\right)---\left(9\right)$

式中为条件分布概率的最大值，

Viterbi算法求解最大值首先定义

$M\left(iu,t\right)=P\left(V_t|{\left(V_{nmo}\right)}_t\right)\underset{jv}{\max }\left\{P\left({\left(V_{nmo}\right)}_t|{\left(V_{nmo}\right)}_{t-1}\right)M\left(jv,t-1\right)\right\}---\left(10\right)$

iu＝(V_nmo)_t，jv＝(V_nmo)_t-1，t＝1,2,3,…,T，P(V_t|(V_nmo)_t＝iv)表示在时刻t的观测动校正速度观测序列和状态函数的联合概率分布，M(iu,t)表示传递序列变量携带了V_t-1之前所有信息；

M(iu,t)的初始值为M(iu,t)＝P(V₁|(V_nmo)₁)P((V_nmo)₁)，因此获得最后的动校正速度序列

${\max }_{q_1}^{T}P\left(V_1^T,{\left(V_{nmo}\right)}_1^T\right)={\max }_{iu}V\left(iu,T\right)---\left(11\right)$

将预处理的速度谱w(i,j)为观测的动校正速度观测序列V₁ ^T，将动校正速度状态函数序列用于记录路径的传递变量记为L(i,j)＝i_k|j-1；通过自动拾取动校正速度V_nmo和η，得到v-t曲线和η-t曲线。

进一步，所述步骤4中：

Step one：对其输入v-t曲线和η-t曲线，按下式计算走时；

$t\left(k\right)=\sqrt{t_{p0}^2\left(k\right)+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}\left(k\right)x^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+\left(1+2{\mathrm{\&eta;}}_{eff}\right)x^2\&rsqb;}};$

Step two：对于零偏移距每个t0开分析时窗形成搬移块，找到控制点，根据Akhaliaf的非双曲各向异性时距方程计算轨迹窗在窗口里计算时差进行整体的搬移；

Step three：由于地震信号是离散的，如果t(k)在两个离散点之间，那么我们采用一维lanczons插值算法进行8点插值：

$A\left(t\left(k\right)\right)=\underset{n=k-3}{\overset{k+4}{\&Sigma;}}A\left(n\right)L\left(t\left(k\right)-k\right)---\left(20\right)$

式中L表示如下计算：

$L\left(x\right)=\left\{\begin{array}{cc}1& \left(x=0\right)\\ \frac{\sin \left(\mathrm{\&pi;}x\right)}{\mathrm{\&pi;}x}\frac{\sin \left(\mathrm{\&pi;}x/N\right)}{\mathrm{\&pi;}x/N}& \left(x\&NotEqual;0^\left|x\right|<N\right)\\ 0& \left(\left|x\right|\&GreaterEqual;N\right)\end{array}\right.---\left(21\right)$

A(t(k))——是动校正后t(k)的振幅值；

$t\left(k\right)=\sqrt{t_{p0}^2\left(k\right)+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}\left(k\right)x^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+\left(1+2{\mathrm{\&eta;}}_{eff}\right)x^2\&rsqb;}}.---\left(22\right)$

本发明的有益效果是能够处理大排列地震数据，处理结果准确。

附图说明

图1是本发明方法步骤流程示意图；

图2是基于Akhalif时距曲线非双曲速度分析实现的技术流程图；

图3是Viterbi自动拾取的计算流程图；

图4是viterbi自动拾取算法实现流程图；

图5是BMS计算示意图；

图6是BMS动校正技术实现流程图；

图7是原始CMP数据；

图8是速度谱和η谱；

图9是BMS动校正结果和某商业软件动校正结果对比图。

具体实施方式

下面结合具体实施方式对本发明进行详细说明。

本发明是一套基于各向异性介质的理论的叠前处理算法方案，对大排列的CMP(CDP)地震数据进行非双曲速度分析形成Vnmo谱和η谱，然后采用viterbi自动拾取优化算法，拾取v-t曲线，采用BMS动校正法，进行动校正处理。

如图1所示，为本发明方法的步骤流程。图2是基于Akhalif时距曲线非双曲速度分析实现的技术流程图。

双谱法非双曲速度分析

1、原理

利用Akhalifah非双曲时差分析公式作为各向异性非双曲速度的基础，如下式：

$t^2\left(x\right)=t_{p0}^2+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}{x}^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+\left(1+2{\mathrm{\&eta;}}_{eff}\right)x^2\&rsqb;}---\left(1\right)$

式中：t(x)—偏移距x时的旅行时

t_p0—纵波双程旅行时

V_nmo—动校正速度

η_eff—等效η各向异性参数

x—偏移距。

该双曲时差分析公式由两部分组成：第一部分是由前两项组成的双曲时差，第二部分是由第三项构成的各向异性的非双曲时差。通过理论分析可知，该式对大偏移距(大排列)的地震数据具有良好处理效果。

非双曲时距曲线公式主要有两个参数：一个是(动校正速度)，另一个是(等效各向异性参数η)。这两个参数对于非双曲旅行时的影响在偏移距上的分布是不一样的，动校正速度对于非双曲旅行时的影响是全局性的，对于整个排列起决定性的作用，而等效各向异性参数只在大偏移距的情况下才会对非双曲旅行时有较明显的影响，所以我们只有在对(X/D>1)远道数据进行分析的情况下才能得到比较精确的值。再结合常规速度分析相似系数的判断准则，就形成了基于各向异性非双曲速度分析。

步骤1：如图2所示，输入地震数据(道集格式为CMP道集或CDP道集)。对于P波地震数据采用如下非双曲速度分析的双谱算法，实现基于各向异性非双曲速度分析，减少了运算量：

Step one：在一个确定的t₀下，采用短排列CMP道集，进行常规速度分析，得到一个初始的动校正速度V_nmo；

Step two：将Step one中得到初始动校正速度和t₀代入Alkhalifah非双曲时差公式，对大排列地震数据进行扫描，得到一个初始的η值；

Step three：我们将Step two中用初始η值在代入非双曲时差公式，对大排列地震数据进行扫描，重新确定动校正速度V_nmo；

Step four：重复执行Step two和Step three反复迭代V_nmo和η两到三次，我们将获得稳定的V_nmo和η；

Step five：对下一个t₀进行扫描，执行前面四步，直到全部的t₀被扫描完成，输出速度V_nmo谱和η谱。

步骤2：对步骤1中得到速度V_nmo谱和η谱进行自动拾取得到V_nmo-t₀曲线和η-t₀曲线。

拾取有效动校正速度是速度分析中一个很重要的工作环节。拾取有效动校正速度的方式一般可以分为人机交互拾取，人工拾取和计算机自动拾取。一般而言，人机交互拾取和人工拾取具有较高的灵活性和容错性，因为其人眼对聚焦性一般不是很高的速度谱具有较高的拾取能力，但是前两者同时也具有较大的局限性：效率低下，人工成本较高且精度不高等。本发明就是将自动优化拾取技术应用到对动校正速度的拾取，不但提高了拾取的精度，而且还能大大提高速度分析的效率。

在本发明中，我们将采用Viterbi算法进行优化拾取。如图3所示为Viterbi自动拾取的计算流程图，Viterbi算法是一种按照约束条件自动拾取目标函数最优解的算法。该算法最早被意大利电讯工程师Viterbi用于卷积解码问题而得名。

自动拾取问题同其他的优化问题一样，包括两个步骤：1.定义优化目标函数；2.优化搜索(优化拾取)。本发明所使用Viterbi方法就是使用在第二步骤中进行优化搜索最优解。

由于Viterbi算法在理论上要求其目标函数在数据源上连续光滑，所以速度谱连续光滑性直接影响影响Viterbi算法的搜索的准确性。如果我们所得到速度谱和η谱的聚焦性较差，自动拾取将无法得到准确的结果，所以要对速度谱和η谱进行圆滑处理。

综上所述，基于Viterbi算法的速度谱自动拾取的两个步骤总结为：1.对速度谱和η谱进行圆滑处理；2.对圆滑处理后速度谱和η谱进行Viterbi优化拾取。

Step one：速度谱和η谱的圆滑处理。

设共中心点道集(CMP)有N道地震记录,该CMP道集数据表达为u(i,j)，根据速度谱计算原理和Taner的相似系数准则可获得速度谱S(t_0,V_nmo)和η谱S(t_0,η)表达如下式：

$S\left(t_0,V_{nmo}\right)=\frac{\underset{j=0}{\overset{M}{\&Sigma;}}{\&lsqb;\underset{i=1}{\overset{N}{\&Sigma;}}u\left(i,j+r_i\right)\&rsqb;}^2}{N\underset{j=0}{\overset{M}{\&Sigma;}}\underset{i=1}{\overset{N}{\&Sigma;}}{u}^2\left(i,j+r_i\right)}---\left(2\right)$

式中i＝(1,2,3…N)—地震道道号

j＝(0,1,2,…,M)—地震道内采样点序号

—地震波偏移距延迟时间的序号

t_i—偏移距引起的延迟时间

△t—时间采样间隔

S(t_0,V_nmo)—相似系数

t_0,—零偏移距的双程旅行时

V_nmo—动校正速度

现在为了方便表达，设μ(x_i,y_j)＝S(t_0,V_nmo)式中x_i＝t₀，i＝(0,1,2…m)，y_j＝V_nmo(j＝0,1,…,n)，x_i为零偏移距双程旅行时序列，y_j为动校正速度序列。现在我们在引入另一个量w(i,j)来表示圆滑后的速度谱其关系如下，在这里我们采用空间均值滤波方法：

$w\left(i,j\right)=\frac{\underset{jw=w_1}{\overset{w_2}{\&Sigma;}}\underset{iw=l_1}{\overset{l_2}{\&Sigma;}}\mathrm{\&mu;}\left(x_{iw},y_{jw}\right)}{\left(w_2-w_1\right)\&times;\left(l_2-l_1\right)}---\left(3\right)$

式中w₁＝i-l_w/2，w₂＝i+l_w/2，l₁＝j-l_w/2，l₂＝j+l_w/2，l_w为滑动窗口的长度。

Step two：使用Viterbi算法进行速度自动拾取。

由于Viterbi算法可知，V＝V₁ ^T＝(V₁,V₂,…V_T)为一个所观测的动校正速度的序列，则其联合概率分布总可表示为：

$P\left({V_1}^T\right)=P\left(V_1\right)\underset{t=2}{\overset{T}{\&Pi;}}P\left(V_t\right|{V_1}^{t-1})---\left(4\right)$

式中，P(V₁)为初始动校正速度观测序列，为观测序列条件分布概率，t＝1,2…,T。式(4)表明：在t时刻，观测序列V_T的条件分布概率依赖于V_T-1前的所有的值，P(V₁ ^t-1)为其先验概率，为其后验概率。

假设在t时刻的动校正速度状态(V_nmo)_t是{1,…,M}内一个有限的数，该状态的向量表示为

V_nmo＝{(V_nmo)₁}^T＝((V_nmo)₁,(V_nmo)₂,…(V_nmo)_T)      (5)

根据随机过程一阶Markov链理论和Bayes规则，前面的观测动校正速度序列V₁ ^T和状态序列之间的概率关系为

$P\left(V_t|{\left(V_{nmo}\right)}_1^T,V_1^{t-1}\right)=P\left(V_t|{\left(V_{nmo}\right)}_t\right)---\left(6\right)$

$P\left(V_t|{\left(V_{nmo}\right)}_1^{t+1},V_1^t\right)=P{\left(V_{t+1}|\left(V_{nmo}\right)\right)}_t---\left(7\right)$

式(6)中P(V_t|(V_nmo)_t)为发射概率，即两种序列联合条件分布概率，P(V₁ ^t-1)为先验概率，联合条件分布概率P(V_t|(V_nmo)_t)为后验概率。式(7)中，P(V_t+1|(V_nmo)_t)为传递概率，即状态序列分布概率，P(V₁ ^t)为先验概率，P(V_t+1|(V_nmo)_t)为后验概率。

由于式(6)和式(7)可知，当要预测动校正速度序列，动校正速度状态值是过去所有相关的动校正速度状态值和观测序列值的函数。根据独立概率分布的假设可知，观测序列和状态变量的联合概率分布的关系还可以表示为：

$P\left({V_1}^t,{\left(V_{nmo}\right)}_1^T\right)=P\left({\left(V_{nmo}\right)}_1\right)\underset{t=1}{\overset{T-1}{\&Pi;}}P\left({\left(V_{nmo}\right)}_{t+1}|{\left(V_{nmo}\right)}_t\right)\underset{t=1}{\overset{T}{\&Pi;}}P\left(V_t|{\left(V_{nmo}\right)}_t\right)---\left(8\right)$

通过最大概率状态这可以利用下面的最大算法可以得到

$q_1^{T*}=\arg \max {|}_{q1}^{T}P<\left({\left(V_{nmo}\right)}_1^T|V_1^T\right)=\arg \max {|}_{q_1}^{T}P\left({\left(V_{nmo}\right)}_1^T,V_1^T\right)---\left(9\right)$

式中为条件分布概率的最大值。

Viterbi算法可以求解最大值首先定义

$M\left(iu,t\right)=P\left(V_t|{\left(V_{nmo}\right)}_t\right)\underset{jv}{\max }\left\{P\left({\left(V_{nmo}\right)}_t|{\left(V_{nmo}\right)}_{t-1}\right)M\left(jv,t-1\right)\right\}---\left(10\right)$

式(10)，iu＝(V_nmo)_t，jv＝(V_nmo)_t-1，t＝1,2,3,…,T。P(V_t|(V_nmo)_t＝iv)表示在时刻t的观测动校正速度观测序列和状态函数的联合概率分布，或称为传递函数，主要起到递归计算作用，M(iu,t)表示传递序列变量携带了V_t-1之前所有信息。

M(iu,t)的初始值为M(iu,t)＝P(V₁|(V_nmo)₁)P((V_nmo)₁)，因此获得最后的动校正速度序列

${\max }_{q_1}^{T}P({V_1}^T,{\left(V_{nmo}\right)}_1^T)={\max }_{iu}V(iu,T)---\left(11\right)$

如果上述递归的最大增量ju^*(iu,t)得以确定，那么最佳的动校正速度

由于上面的算法我们可以把将我们预处理的速度谱w(i,j)为观测的动校正速度观测序列V₁ ^T，而结合上述的Viterbi算法，可以将动校正速度状态函数序列用于记录路径的传递变量记为L(i,j)＝i_k|j-1；我们就可以自动拾取动校正速度V_nmo和η，得到v-t曲线和η-t曲线。

步骤3：对步骤2中的v-t曲线和η-t曲线进行线性插值得到新的v-t0曲线和η-t0曲线。由于速度谱零偏移距双程旅行时间t0的间隔与地震资料时间采样间隔不同，所以要进行重新采样，才能进行动校正。

设离散v-t曲线序列为v(t₀)(t＝0,dt₀,2dt₀,…,kdt₀,…,ndt₀)，dt₀为速度分析时设定的双程旅行时的扫描采样间隔，为了方便表达，我们将其记成v(kdt₀)＝v(k)；设现行插值后离散V-t₀曲线序列为

V(t)(t＝0,dt,2dt,…,idt,…,ndt)，假设任意确定的t值在kdt₀和(k+1)dt₀之间。线性插值表达式如下：

$V\left(t\right)=\frac{v(k+1)-v\left(k\right)}{{\mathrm{dt}}_0}(t-{\mathrm{kdt}}_0)+v\left(k\right)---\left(12\right)$

由此处理可以得到重新采样的v-t曲线，我们也可以将v-t曲线换成η-t曲线得到重新采样的η-t曲线。

步骤4：基于非双曲时差的BMS动校正法。

采用BMS动校正方法进行插值搬移。(流程图如下图4)

在BMS法中，先假设在CMP道集中任意反射波脉冲的波形是完全相同的。BMS流程(如图5所示)首先将CMP道集中零偏移距的地震道分成若干时窗长度为BL的有重叠单元，例如图5中AB，CD段。我们假设每一段的长度都可以容下一个完整反射地震波脉冲。AB单元，沿着非双曲时差曲线投影到其他的地震道上形成等价单元，然后进行非双曲搬移，对于其它单元重复操作，校平整个道集。图6为BMS动校正技术实现流程图。

实现步骤：

Step one：对其输入v-t曲线和η-t曲线，按下式计算走时；

$t\left(k\right)=\sqrt{t_{p0}^2\left(k\right)+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}\left(k\right)x^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+(1+2{\mathrm{\&eta;}}_{eff})x^2\&rsqb;}};$

Step two：对于零偏移距每个t0开分析时窗形成“搬移块”，找到控制点，根据Akhaliaf的非双曲各向异性时距方程计算轨迹窗，按照图5中所示的搬移方式在窗口里计算时差进行整体的搬移；

Step three：由于地震信号是离散的，如果t(k)在两个离散点之间，那么我们采用一维lanczons插值算法进行8点插值(见下式(20))。

$A\left(t\left(k\right)\right)=\underset{n=k-3}{\overset{k+4}{\&Sigma;}}A\left(n\right)L\left(t\left(k\right)-k\right)---\left(20\right)$

式中L表示如下计算：

$L\left(x\right)=\left\{\begin{array}{cc}1& (x=0)\\ \frac{sin\left(\mathrm{\&pi;}x\right)}{\mathrm{\&pi;}x}\frac{sin(\mathrm{\&pi;}x/N)}{\mathrm{\&pi;}x/N}& (x\&NotEqual;0^|x|<N)\\ 0& \left(\right|x|\&GreaterEqual;N)\end{array}\right.---\left(21\right)$

其中：A(t(k))——是动校正后t(k)的振幅值；

$t\left(k\right)=\sqrt{t_{p0}^2\left(k\right)+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}\left(k\right)x^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+(1+2{\mathrm{\&eta;}}_{eff})x^2\&rsqb;}}---\left(22\right)$

本发明的优点：

(a)以非双曲线速度分析为基础，同时利用viterbi优化算法，对速度谱和η进行高精度的自动拾取，极大的提高非双曲速度分析效率；

(b)采用BMS法联合Literilysinc非线性插值，有效的压制了拉伸畸变，提高了动校正的精度，提高了叠加剖面的质量。

本发明的实例

我们应用该方案对下面的一个大偏移距CMP道集进行处理(数据参数如下表1)。原始数据如下：

表1：原始超CMP道集参数表

参数名称	参数大小
		道间距	25m
覆盖次数	581
		采样率	4ms
最大跑间距	8km

原始CMP数据如图7所示(横坐标为偏移距，纵坐标为零偏移距旅行时)，速度谱和η谱如图8所示，该道集的速度分析结果可以得到我们拾取速度曲线和η曲线与速度谱和η相似系数最大处较为吻合。

图9：左边为BMS动校正结果、右边为某商业软件动校正结果。

由上图可以看出BMS的动校正方法对浅层远偏移距的同相轴的拉平效果较好，较好的压制拉伸畸变。

本发明对于三维的大排列地震数据，采用基于VTI介质的Akhaliaf非双曲速度分析计算，形成Vnmo谱和η谱，然后采用viterbi自动拾取优化算法，拾取v-t曲线，采用BMS动校正法，进行动校正处理，提高对浅层远偏移距的地震数据的信息保留，提高速度分析及动校正的精度和效率，提高叠加剖面的精度。

以上所述仅是对本发明的较佳实施方式而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施方式所做的任何简单修改，等同变化与修饰，均属于本发明技术方案的范围内。

Claims (3)

1.一种Viterbi-BMC算法的自动速度分析方法，其特征在于按照以下步骤进行：

步骤1：输入地震数据，对于P波地震数据采用非双曲速度分析的双谱算法：

Step one：在一个确定的t₀下，采用短排列CMP道集，进行常规速度分析，得到一个初始的动校正速度v_nmo；

Step two：将Step one中得到初始动校正速度v_nmo和t₀代入非双曲时差公式，对大排列地震数据进行扫描，得到一个初始的η值；

Step three：用初始η值代入非双曲时差公式，对大排列地震数据进行扫描，重新确定动校正速度v_nmo；

Step four：重复执行Step two和Step three反复迭代v_nmo和η两到三次，获得稳定的v_nmo和η；

Step five：对下一个t₀进行扫描，执行前面四步，直到全部的t₀被扫描完成，输出速度v_nmo谱和η谱。

步骤2：对步骤1中得到速度v_nmo谱和η谱进行自动拾取得到V_nmo-t₀曲线和η-t₀曲线；

步骤3：对步骤2中的v-t曲线和η-t曲线进行线性插值得到新的v-t0曲线和η-t0曲线；

步骤4：基于非双曲时差的BMS动校正法：采用BMS动校正方法进行插值搬移；

先假设在CMP道集中任意反射波脉冲的波形是完全相同的，首先将CMP道集中零偏移距的地震道分成若干时窗长度为BL的有重叠单元，假设每一段的长度都可以容下一个完整反射地震波脉冲，沿着非双曲时差曲线投影到其他的地震道上形成等价单元，然后进行非双曲搬移，对于其它单元重复操作，校平整个道集。

2.按照权利要求1所述一种Viterbi-BMC算法的自动速度分析方法，其特征在于：所述步骤2中，自动拾取的步骤为：

Step one：速度谱和η谱的圆滑处理；

设共中心点道集(CMP)有N道地震记录，该CMP道集数据表达为u(i,j)，根据速度谱计算原理和Taner的相似系数准则可获得速度谱S(t_0,V_nmo)和η谱S(t_0,η)表达如下式：

$S\left(t_{0,}{V}_{nmo}\right)=\frac{\underset{j=0}{\overset{M}{\&Sigma;}}{\&lsqb;\underset{i=1}{\overset{N}{\&Sigma;}}u\left(i,j+r_i\right)\&rsqb;}^2}{N\underset{j=0}{\overset{M}{\&Sigma;}}\underset{i=1}{\overset{N}{\&Sigma;}}{u}^2\left(i,j+r_i\right)}---\left(2\right)$

式中i＝(1,2,3…N)—地震道道号

j＝(0,1,2,…,M)—地震道内采样点序号

—地震波偏移距延迟时间的序号

t_i—偏移距引起的延迟时间

△t—时间采样间隔

S(t_0,V_nmo)—相似系数

t_0,—零偏移距的双程旅行时

V_nmo—动校正速度

现在为了方便表达，设μ(x_i,y_j)＝S(t_0,V_nmo)，式中x_i＝t₀，i＝(0,1,2…m)，y_j＝V_nmo(j＝0,1,…,n)，x_i为零偏移距双程旅行时序列，y_j为动校正速度序列，引入另一个量w(i,j)来表示圆滑后的速度谱其关系如下：

$w\left(i,j\right)=\frac{\underset{jw=w_1}{\overset{w_2}{\&Sigma;}}\underset{iw=l_1}{\overset{l_2}{\&Sigma;}}\mathrm{\&mu;}\left(x_{iw},y_{jw}\right)}{\left(w_2-w_1\right)\&times;\left(l_2-l_1\right)}---\left(3\right)$

式中w₁＝i-l_w/2，w₂＝i+l_w/2，l₁＝j-l_w/2，l₂＝j+l_w/2，l_w为滑动窗口的长度；

Step two：使用Viterbi算法进行速度自动拾取；

由Viterbi算法可知，V＝V₁ ^T＝(V₁,V₂,…V_T)为一个所观测的动校正速度的序列，则其联合概率分布总可表示为：

$P\left(V_1^T\right)=P\left(V_1\right)\underset{t=2}{\overset{T}{\&Pi;}}P\left(V_t\right|V_1^{t-1})---\left(4\right)$

式中，P(V₁)为初始动校正速度观测序列，为观测序列条件分布概率，t＝1,2…,T，上式表明：在t时刻，观测序列V_T的条件分布概率依赖于V_T-1前的所有的值，P(V₁ ^t-1)为其先验概率，为其后验概率；

假设在t时刻的动校正速度状态(V_nmo)_t是{1,…,M}内一个有限的数，该状态的向量表示为

V_nmo＝{(V_nmo)₁}^T _＝((V_nmo)₁,(V_nmo)₂,…(V_nmo)_T)       (5)

根据随机过程一阶Markov链理论和Bayes规则，前面的观测动校正速度序列V₁ ^T和状态序列之间的概率关系为

$P\left(V_t\right|{\left(V_{nmo}\right)}_1^T,V_1^{t-1})=P\left(V_t\right|{\left(V_{nmo}\right)}_t)---\left(6\right)$

$P\left(V_t\right|{\left(V_{nmo}\right)}_1^{t+1},V_1^t)=P\left(V_{t+1}\right|{\left(V_{nmo}\right)}_t)---\left(7\right)$

式(6)中P(V_t|(V_nmo)_t)为发射概率，即两种序列联合条件分布概率，P(V₁ ^t-1)为先验概率，联合条件分布概率P(V_t|(V_nmo)_t)为后验概率，式(7)中，P(V_t+1|(V_nmo)_t)为传递概率，即状态序列分布概率，P(V₁ ^t)为先验概率，P(V_t+1|(V_nmo)_t)为后验概率；

由于式(6)和式(7)可知，当要预测动校正速度序列，动校正速度状态值是过去所有相关的动校正速度状态值和观测序列值的函数，根据独立概率分布的假设可知，观测序列和状态变量的联合概率分布的关系表示为：

$P(V_1^t,{\left(V_{nmo}\right)}_1^T)=P\left({\left(V_{nmo}\right)}_1\right)\underset{t=1}{\overset{T-1}{\&Pi;}}P\left({\left(V_{nmo}\right)}_{t+1}|{\left(V_{nmo}\right)}_t\right)\underset{t=1}{\overset{T}{\&Pi;}}P\left(V_t|{\left(V_{nmo}\right)}_t\right)---\left(8\right)$

通过最大概率状态这可以利用下面的最大算法得到：

$q_1^{T*}=\mathrm{argmax}{|}_{q1}^{T}P\left({\left(V_{nmo}\right)}_1^T|V_1^T\right)=\mathrm{argmax}{|}_{q_1}^{T}P\left({\left(V_{nmo}\right)}_1^T,V_1^T\right)---\left(9\right)$

式中为条件分布概率的最大值，

Viterbi算法求解最大值首先定义

$M(iu,t)=P\left(V_t|{\left(V_{nmo}\right)}_t\right)\underset{jv}{max}\left\{P\right({\left(V_{nmo}\right)}_t|{\left(V_{nmo}\right)}_{t-1}\left)M(jv,t-1)\right\}---\left(10\right)$

iu＝(V_nmo)_t，jv＝(V_nmo)_t-1，t＝1,2,3,…,T，P(V_t|(V_nmo)_t＝iv)表示在时刻t的观测动校正速度观测序列和状态函数的联合概率分布，M(iu,t)表示传递序列变量携带了V_t-1之前所有信息；

M(iu,t)的初始值为M(iu,t)＝P(V₁|(V_nmo)₁)P((V_nmo)₁)，因此获得最后的动校正速度序列

${\max }_{q_1}^{T}P\left(V_1^T,{\left(V_{nmo}\right)}_1^T\right)={\max }_{iu}V(iu,T)---\left(11\right)$

将预处理的速度谱w(i,j)为观测的动校正速度观测序列V₁ ^T，将动校正速度状态函数序列用于记录路径的传递变量记为L(i,j)＝i_k|_j-1；通过自动拾取动校正速度V_nmo和η，得到v-t曲线和η-t曲线。

3.按照权利要求1所述一种Viterbi-BMC算法的自动速度分析方法，其特征在于：所述步骤4中：

Step one：对其输入v-t曲线和η-t曲线，按下式计算走时；

$t\left(k\right)=\sqrt{t_{p0}^2\left(k\right)+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}\left(k\right)x^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+(1+2{\mathrm{\&eta;}}_{eff})x^2\&rsqb;}};$

Step two：对于零偏移距每个t0开分析时窗形成搬移块，找到控制点，根据Akhaliaf的非双曲各向异性时距方程计算轨迹窗在窗口里计算时差进行整体的搬移；

Step three：由于地震信号是离散的，如果t(k)在两个离散点之间，那么我们采用一维lanczons插值算法进行8点插值：

$A\left(t\left(k\right)\right)=\underset{n=k-3}{\overset{k+4}{\&Sigma;}}A\left(n\right)L(t\left(k\right)-k)---\left(20\right)$

式中L表示如下计算：

A(t(k))——是动校正后t(k)的振幅值；

$t\left(k\right)=\sqrt{t_{p0}^2\left(k\right)+\frac{x^2}{V_{nmo}^2\left(0\right)}-\frac{2{\mathrm{\&eta;}}_{eff}\left(k\right)x^4}{V_{nmo}^2\left(0\right)\&lsqb;t_{p0}^2{V}_{nmo}^2\left(0\right)+(1+2{\mathrm{\&eta;}}_{eff})x^2\&rsqb;}}.---\left(22\right)$