CN104297782A - 一种采用频变的复数相关系数估算瞬时慢度的方法 - Google Patents

Abstract

本发明是一种采用频变的复数相关系数估算瞬时慢度的方法，首先通过频谱分析来确定地震剖面的主频，并在时间方向把地震道投影到以该主频为频率的富里叶变换的基函数上，然后利用投影后的地震剖面按空间道序再按时间采样点顺序逐点估算每个采样点的瞬时慢度。估算瞬时慢度采用的是慢度扫描方式，首先对每一扫描慢度利用当前采样点对应的当前处理道及其相邻道的当前处理时窗内的数据计算出总的复数相关系数，然后利用每一扫描慢度的总的复数相关系数的实部和虚部估算出当前采样点的瞬时慢度。本发明抗噪能力更强，可以利用复数相关系数的实部和虚部的性质更精确地估算瞬时慢度。

Description

一种采用频变的复数相关系数估算瞬时慢度的方法

技术领域

本发明涉及地球物理勘探技术，属于地震资料处理中属性提取范畴，是一种采用频变的复数相关系数估算瞬时慢度的方法。

背景技术

在地震资料处理中，瞬时慢度是地震数据一个非常重要的属性。瞬时慢度可广泛应用于地震数据插值或外推、噪声压制、VSP数据上下行波分离、波束偏移、速度层析反演和地震解释中断层预测等。

估算瞬时慢度的方法一般分为三大类。第一类是基于相关性分析的慢度扫描方法；第二类基于局部Radon变换的慢度方法；第三类是基于平面波微分方程的慢度估算方法。

基于相关性分析的慢度扫描方法以当前地震道的当前采样点为中心点，按扫描慢度确定信号在其相邻地震道上的轨迹，然后以信号在地震道上的轨迹点为中心取时窗，计算所有时窗内数据之间的互相关系数并求和，这样就得到了以扫描慢度为变量的互相关函数，最后求互相关函数的最大值对应的扫描慢度就可以得到当前地震道的当前采样点处的瞬时慢度。这种方法的优点是算法简单，缺点是得到的慢度精度低。

发明内容

本发明目的在于提供一种可更准确估算瞬时慢度的采用频变的复数相关系数估算瞬时慢度的方法。

本发明的具体实施步骤如下：

1）采集得到二维地震剖面；

2）对地震剖面进行频谱分析并确定地震剖面的主频f_c；

对地震剖面进行F-K谱分析并确定地震剖面的最小和最大视速度v_min和v_max；

3）在时间方向，把地震道投影到以地震剖面的主频f_c为频率的富里叶变换的基函数上，得到投影后的地震剖面A（x，t，f_c）：

$A(x,t,f_c)=a(x,t)e^{-2\mathrm{\πi}{f}_{c}t}$

式中，a（x，t）是二维地震剖面，x和t分别为地震道的空间和时间坐标，f_c为地震剖面的主频；A（x，t，f_c）为投影后的地震剖面；

4）对投影后的地震剖面A（x，t，f_c）先按空间道序再按时间采样点顺序逐点按步骤5）-18）估算瞬时慢度；

5）根据步骤2)确定的地震剖面的主频f_c，用下式计算以采样点为中心的沿时间方向的处理时窗的长度：

$w_t=\frac{1}{f_c\mathrm{\Δt}}$

式中，f_c为地震剖面的主频，Δt为地震剖面的时间采样率，w_t为处理时窗的长度；

6）根据步骤2)确定的地震剖面的最小和最大视速度v_min和v_max，计算出慢度扫描的最少值和最大值；

p_min=1/v_max

p_max=1/v_min

式中，p_min和p_max表示慢度扫描的最少值和最大值，v_min和v_max是由步骤2）确定的地震剖面的最小和最大视速度；

7）根据慢度扫描个数，计算慢度扫描的增量：

dp=（p_max-p_min）/（n_p-1）

式中，n_p为用户给出的慢度扫描个数；dp为慢度扫描的增量；p_min和p_max表示慢度扫描的最少值和最大值；

8）以当前处理地震道为中心道，确定参与瞬时慢度估算的相邻地震道数w_x；

所述的相邻地震道数w_x少于9大于3。

9）从慢度扫描的最少值开始到慢度扫描的最大值为止，按慢度扫描增量计算扫描慢度p_i：

p_i=（i-1）*dp（i=1,2,…,n_p）

其中，i表示扫描慢度的编号，p_i是第i个扫描慢度，dp表示由步骤7）计算的慢度增量；

10）以当前处理地震道和处理时窗的中心点为基点，对每一扫描慢度按照下式计算其对应的相邻地震道上的处理时窗的中心点时间τ_i,j和时窗对应的信号的延迟时Δτ_i,j；

τ_i,j=τ₀+（x_j-x₀）*p_i  （-w_x/2≤j≤w_x/2,i=1,2…..n_p）

Δτ_i,j=（x_j-x₀）*p_i  （-w_x/2≤j≤w_x/2，i=1,2…..n_p）

式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0编号表示当前处理道；p_i表示第i个扫描慢度。x₀表示当前处理地震道的空间作标，x_j表示第j个相邻地震道的空间作标；τ₀表示当前处理地震道的处理时窗的中心点时间；τ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δτ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；w_x表示由步骤8）定义的相邻地震道数；n_p表示由步骤7）定义的慢度扫描个数；

11）对每一扫描慢度根据步骤10）计算的当前道及其相邻道的处理时窗的中心点时间以及步骤5）定义的时窗长度，提取当前道及其相邻道的处理时窗内的采样点的数值

${\tilde{B}}_{j,n}\left(p_i\right)=A(x_j,m_{i,j}+n,f_c)$

m_i,j=τ_i,j/Δt

1≤i≤n_p，-w_x/2≤j≤w_x/2，-w_t/2≤n≤w_t/2

式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度。τ_i,j表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δt为地震剖面的时间采样率；m_i,j表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点对应的采样点号；A(x_j,m_i,j+n,fc)表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗内的编号为n的采样点的数值；表示当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；w_x表示步骤8）定义的相邻地震道数；n_p表示步骤7）定义的慢度扫描个数。w_t表示步骤5）定义的时窗长度；f_c表示由步骤2）确定的地震剖面的主频；

12）对每一扫描慢度，根据10）计算的相邻道处理时窗对应的信号延迟时，对提取的处理时窗内的采样点的数值进行相位校正；

$B_{j,n}\left(p_i\right)={\tilde{B}}_{j,n}\left(p_i\right)e^{-2\mathrm{\πi}{f}_c{\mathrm{\Δ\τ}}_{i,j}}$

式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；表示当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；Δτ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；f_c表示由2）确定的地震剖面的主频；

13）对每一扫描慢度，计算相位校正后的每个相邻地震道的处理时窗内的信号的能量E_j(p_i)：

$E_j\left(p_i\right)=\underset{n=-w_t/2}{\overset{w_t/2}{\mathrm{\Σ}}}{B}_{j,n}\left(p_i\right){\stackrel{\‾}{B}}_{j,n}\left(p_i\right)$

式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；为B_j,n(p_i)的共轭；E_j(p_i)表示当用慢度p_i进行扫描时第j个相邻地震道处理时窗内信号的能量；w_t表示由步骤5）定义的时窗长度；

14）对每一扫描慢度,计算参与慢度估算的各地震道处理时窗内的信号的复数互相关系数C_j,k(p_i)；

$C_{j,k}\left(p_i\right)=\underset{n=-w_t/2}{\overset{w_t/2}{\mathrm{\Σ}}}{B}_{j,n}\left(p_i\right){\stackrel{\‾}{B}}_{k,n}\left(p_i\right)/\sqrt{E_j\left(p_i\right)E_k\left(p_i\right)}$

-w_x/2≤j≤w_x/2，k>j

式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；B_k,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第k个相邻地震道的处理时窗内提取的编号为n的采样点的数值；表示B_k,n(p_i)的共轭。E_j(p_i)表示当用慢度p_i进行扫描时第j个相邻地震道处理时窗内的信号的能量；E_k(p_i)表示当用慢度p_i进行扫描时第k个相邻地震道处理时窗内信号的能量；

C_j,k(p_i)表示当用慢度p_i进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；w_t表示由步骤5）定义的时窗长度；w_x表示由步骤8）定义的相邻地震道数；

15）对每一扫描慢度,把其对应的各地震道处理时窗内信号的复数互相关系数进行求和，得到该慢度的总复数相关系数C(p_i)：

$C\left(p_i\right)=\underset{j=-w_x/2}{\overset{w_x/2}{\mathrm{\Σ}}}\underset{k>j}{\overset{w_x/2}{\mathrm{\Σ}}}{C}_{j,k}\left(p_i\right)$

式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；p_i表示第i个扫描慢度；C_j,k(p_i)表示当用慢度p_i进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；w_x表示由步骤8）定义的相邻地震道数；

16）对每一扫描慢度,计算该慢度对应的总复数相关系数的实部和虚部；

RC(p_i)=real(C(p_i))

IC(p_i)=image(C(p_i))

式中，i表示扫描慢度的编号；p_i表示第i个扫描慢度；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；RC(p_i)表示C(p_i)的实部；IC(p_i)表示C(p_i)的虚部；

17）根据步骤16）的计算结果，在步骤6）慢度扫描的最小值和最大值之间寻找使复数相关系数虚部为零的绝对值最小的扫描慢度q₁：

IC(q₁)=0

式中，IC(p_i)表示C(p_i)的虚部。C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；q₁表示使复数相关系数虚部为零的绝对值最小的扫描慢度；

18）在慢度0和使复数相关系数虚部为零的绝对值最小的零点之间，寻找使复数相关系数实部为最大的扫描慢度q₂，这个扫描慢度即为当前处理点的瞬时慢度；

RC(q₂)=maxRC(p_i)

式中，RC(p_i)表示C(p_i)的虚部；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和，q₂表示使复数相关系数实部为最大的扫描慢度。

本发明一方面由于是在投影后的地震数据上进行瞬时慢度估算，因此抗噪能力更强，另一方面由于采用了复数相关系数，因此可以利用复数相关系数的实部和虚部的性质更精确地估算瞬时慢度。

附图说明

图1地震剖面（左上）、投影后的地震剖面的实部（左下）、虚部（右下）和模（左上）。在该图例中，投影所采用的地震剖面的主频为25Hz。

图2当前处理地震道的空间坐标x₀，相邻地震道的空间坐标x_j、相邻地震道数w_x、沿时间方向的处理时窗的长度w_t、当前处理时窗的中心点τ₀，当用慢度p_i进行扫描时第j个相邻地震道处理时窗的中心点时间τ_i,j和时窗对应的信号的延迟时Δτ_i,j示意图。

图3沿某一地震剖面某一时间的各采样点的复数互相关函数的实部（上）和虚部（下）

图4合成地震数据（左）和由本发明方法估算的瞬时慢度（右）

具体实施方案

以下结合附图详细说明是本发明具体实施方案。

1）采集得到二维地震剖面；

2）对地震剖面进行频谱分析并确定地震剖面的主频f_c；

对地震剖面进行F-K谱分析并确定地震剖面的最小和最大视速度v_min和v_max；

3）在时间方向，把地震道投影到以地震剖面的主频f_c为频率的富里叶变换的基函数上，得到投影后的地震剖面A（x，t，f_c）：

$A(x,t,f_c)=a(x,t)e^{-2\mathrm{\πi}{f}_{c}t}$

式中，a（x，t）是二维地震剖面，x和t分别为地震道的空间和时间坐标，f_c为地震剖面的主频；A（x，t，f_c）为投影后的地震剖面；

图1展示的是a（x，t）、投影后的地震剖面A（x，t，f_c）的实部real（A（x，t，f_c））、虚部image（A（x，t，f_c））和模||A（x，t，f_c）||。在该图例中投影所使用的地震剖面的主频f_c为25Hz。

4）对投影后的地震剖面A（x，t，f_c）先按空间道序再按时间采样点顺序逐点按步骤5）-18）估算瞬时慢度；

5）根据步骤2)地震剖面的主频f_c，用下式计算以采样点为中心的沿时间方向的处理时窗的长度：

$w_t=\frac{1}{f_c\mathrm{\Δt}}$

式中，f_c为地震剖面的主频，Δt为地震剖面的时间采样率，w_t为处理时窗的长度；

图2中两个矩形方框示意的是沿时间方向的处理时窗。

6）根据步骤2)确定的地震剖面的最小和最大视速度v_min和v_max，计算出慢度扫描的最少值和最大值；

p_min=1/v_max

p_max=1/v_min

式中，p_min和p_max表示慢度扫描的最少值和最大值，v_min和v_max是由步骤2）确定的地震剖面的最小和最大视速度；

7）根据慢度扫描个数，计算慢度扫描的增量：

dp=（p_max-p_min）/（n_p-1）

式中，n_p为用户给出的慢度扫描个数；dp为慢度扫描的增量；p_min和p_max表示慢度扫描的最少值和最大值；

8）以当前处理地震道为中心道，确定参与瞬时慢度估算的相邻地震道数w_x；

所述的相邻地震道数w_x少于9大于3。

图2中x₀和x_j示意的分别是当前地震道和相邻地震道的空间坐标。

9）从慢度扫描的最少值开始到慢度扫描的最大值为止，按慢度扫描增量计算扫描慢度p_i：

p_i=（i-1）*dp（i=1,2,…,n_p）

其中，i表示扫描慢度的编号，p_i是第i个扫描慢度，dp表示由步骤7）计算的慢度增量；

10）以当前处理地震道和处理时窗的中心点为基点，对每一扫描慢度按照下式计算其对应的相邻地震道上的处理时窗的中心点时间τ_i,j和时窗对应的信号的延迟时Δτ_i,j；

τ_i,j=τ₀+（x_j-x₀）*p_i  （-w_x/2≤j≤w_x/2,i=1,2…..n_p）

Δτ_i,j=（x_j-x₀）*p_i  （-w_x/2≤j≤w_x/2，i=1,2…..n_p）

式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0编号表示当前处理道；p_i表示第i个扫描慢度。x₀表示当前处理地震道的空间作标，x_j表示第j个相邻地震道的空间作标；τ₀表示当前处理地震道的处理时窗的中心点时间；τ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δτ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；w_x表示由步骤8）定义的相邻地震道数；n_p表示由步骤7）定义的慢度扫描个数；

图2中两个实心圆点示意的分别是当前处理地震道的处理时窗的中心点时间τ₀和当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点时间τ_i,j。两个两个实心圆点之间的时差为第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时。

11）对每一扫描慢度根据步骤10）计算的当前道及其相邻道的处理时窗的中心点时间以及步骤5）定义的时窗长度，提取当前道及其相邻道的处理时窗内的采样点的数值

${\tilde{B}}_{j,n}\left(p_i\right)=A(x_j,m_{i,j}+n,f_c)$

m_i,j=τ_i,j/Δt

1≤i≤n_p，  -w_x/2≤j≤w_x/2，-w_t/2≤n≤w_t/2

式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度。τ_i,j表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δt为地震剖面的时间采样率；m_i,j表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点对应的采样点号；A(x_j,m_i,j+n,fc)表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗内的编号为n的采样点的数值；表示当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；w_x表示步骤8）定义的相邻地震道数；n_p表示步骤7）定义的慢度扫描个数。w_t表示步骤5）定义的时窗长度；f_c表示由步骤2）确定的地震剖面的主频；

12）对每一扫描慢度，根据10）计算的相邻道处理时窗对应的信号延迟时，对提取的处理时窗内的采样点的数值进行相位校正；

$B_{j,n}\left(p_i\right)={\tilde{B}}_{j,n}\left(p_i\right)e^{-2\mathrm{\πi}{f}_c{\mathrm{\Δ\τ}}_{i,j}}$

式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；表示当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；Δτ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；f_c表示由2）确定的地震剖面的主频；

13）对每一扫描慢度，计算相位校正后的每个相邻地震道的处理时窗内的信号的能量E_j(p_i)：

$E_j\left(p_i\right)=\underset{n=-w_t/2}{\overset{w_t/2}{\mathrm{\Σ}}}{B}_{j,n}\left(p_i\right){\stackrel{\‾}{B}}_{j,n}\left(p_i\right)$

式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；为B_j,n(p_i)的共轭；E_j(p_i)表示当用慢度p_i进行扫描时第j个相邻地震道处理时窗内信号的能量；w_t表示由步骤5）定义的时窗长度；

14）对每一扫描慢度,计算参与慢度估算的各地震道处理时窗内的信号的复数互相关系数C_j,k(p_i)；

$C_{j,k}\left(p_i\right)=\underset{n=-w_t/2}{\overset{w_t/2}{\mathrm{\Σ}}}{B}_{j,n}\left(p_i\right){\stackrel{\‾}{B}}_{k,n}\left(p_i\right)/\sqrt{E_j\left(p_i\right)E_k\left(p_i\right)}$

-w_x/2≤j≤w_x/2，k>j

式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；B_k,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第k个相邻地震道的处理时窗内提取的编号为n的采样点的数值；表示B_k,n(p_i)的共轭。E_j(p_i)表示当用慢度p_i进行扫描时第j个相邻地震道处理时窗内的信号的能量；E_k(p_i)表示当用慢度p_i进行扫描时第k个相邻地震道处理时窗内信号的能量；

C_j,k(p_i)表示当用慢度p_i进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；w_t表示由步骤5）定义的时窗长度；w_x表示由步骤8）定义的相邻地震道数；

15）对每一扫描慢度,把其对应的各地震道处理时窗内信号的复数互相关系数进行求和，得到该慢度的总复数相关系数C(p_i)：

$C\left(p_i\right)=\underset{j=-w_x/2}{\overset{w_x/2}{\mathrm{\Σ}}}\underset{k>j}{\overset{w_x/2}{\mathrm{\Σ}}}{C}_{j,k}\left(p_i\right)$

式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；p_i表示第i个扫描慢度；C_j,k(p_i)表示当用慢度p_i进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；w_x表示由步骤8）定义的相邻地震道数；

16）对每一扫描慢度,计算该慢度对应的总复数相关系数的实部和虚部；

RC(p_i)=real(C(p_i))

IC(p_i)=image(C(p_i))

式中，i表示扫描慢度的编号；p_i表示第i个扫描慢度；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；RC(p_i)表示C(p_i)的实部；IC(p_i)表示C(p_i)的虚部；

图3展示的是在慢度扫描的最少值到慢度扫描的最大值之间各扫描慢度对应的总复数相关系数的实部RC(p_i)和虚部IC(p_i)。

17）根据步骤16）的计算结果，在步骤6）慢度扫描的最小值和最大值之间寻找使复数相关系数虚部为零的绝对值最小的扫描慢度q₁：

IC(q₁)=0

式中，IC(p_i)表示C(p_i)的虚部。C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；q₁表示使复数相关系数虚部为零的绝对值最小的扫描慢度；

18）在慢度0和使复数相关系数虚部为零的绝对值最小的零点之间，寻找使复数相关系数实部为最大的扫描慢度q₂，这个扫描慢度即为当前处理点的瞬时慢度；

RC(q₂)=maxRC(p_i)

式中，RC(p_i)表示C(p_i)的虚部；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和。q₂表示使复数相关系数实部为最大的扫描慢度。

图4展示的是合成地震记录和按步骤5）-18）估算的瞬时慢度剖面。

Claims (2)

1.一种采用频变的复数相关系数估算瞬时慢度的方法，特点是具体实施步骤如下：

1）采集得到二维地震剖面；

2）对地震剖面进行频谱分析并确定地震剖面的主频f_c；对地震剖面进行F-K谱分析并确定地震剖面的最小和最大视速度v_min和v_max；

3）在时间方向，把地震道投影到以地震剖面的主频f_c为频率的富里叶变换的基函数上，得到投影后的地震剖面A(x,t,f_c)：

$A(x,t,f_c)=a(x,t)e^{-2\mathrm{\πi}{f}_{c}t}$

式中，a（x,t）是二维地震剖面，x和t分别为地震道的空间和时间坐标，f_c为地震剖面的主频；A(x,t,f_c)为投影后的地震剖面；

4）对投影后的地震剖面A(x,t,f_c)先按空间道序再按时间采样点顺序逐点按步骤5）-18）估算瞬时慢度；

5）根据步骤2)确定的地震剖面的主频f_c，用下式计算以采样点为中心的沿时间方向的处理时窗的长度：

$w_t=\frac{1}{f_c\mathrm{\Δt}}$

式中，f_c为地震剖面的主频，Δt为地震剖面的时间采样率，w_t为处理时窗的长度；

6）根据步骤2)确定的地震剖面的最小和最大视速度v_min和v_max，计算出慢度扫描的最少值和最大值；

p_min=1/v_max

p_max=1/v_min

式中，p_min和p_max表示慢度扫描的最少值和最大值，v_min和v_max是由步骤2）确定的地震剖面的最小和最大视速度；

7）根据慢度扫描个数，计算慢度扫描的增量：

dp=（p_max-p_min）/（n_p-1）

式中，n_p为用户给出的慢度扫描个数；dp为慢度扫描的增量；p_min和p_max表示慢度扫描的最少值和最大值；

8）以当前处理地震道为中心道，确定参与瞬时慢度估算的相邻地震道数w_x；

9）从慢度扫描的最少值开始到慢度扫描的最大值为止，按慢度扫描增量计算扫描慢度p_i：

p_i=（i-1）*dp（i=1,2,…,n_p）

其中，i表示扫描慢度的编号，p_i是第i个扫描慢度，dp表示由步骤7）计算的慢度增量；

10）以当前处理地震道和处理时窗的中心点为基点，对每一扫描慢度按照下式计算其对应的相邻地震道上的处理时窗的中心点时间τ_i,j和时窗对应的信号的延迟时Δτ_i,j；

τ_i,j=τ₀+（x_j-x₀）*p_i  （-w_x/2≤j≤w_x/2,i=1,2…..n_p）

Δτ_i,j=（x_j-x₀）*p_i（-w_x/2≤j≤w_x/2，i=1,2…..n_p）

式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0编号表示当前处理道；p_i表示第i个扫描慢度。x₀表示当前处理地震道的空间作标，x_j表示第j个相邻地震道的空间作标；τ₀表示当前处理地震道的处理时窗的中心点时间；τ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δτ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；w_x表示由步骤8）定义的相邻地震道数；n_p表示由步骤7）定义的慢度扫描个数；

11）对每一扫描慢度根据步骤10）计算的当前道及其相邻道的处理时窗的中心点时间以及步骤5）定义的时窗长度，提取当前道及其相邻道的处理时窗内的采样点的数值

${\tilde{B}}_{j,n}\left(p_i\right)=A(x_j,m_{i,j}+n,f_c)$

m_i,j=τ_i,j/Δt

1≤i≤n_p，-w_x/2≤j≤w_x/2，-w_t/2≤n≤w_t/2

式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；τ_i,j表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δt为地震剖面的时间采样率；m_i,j表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗的中心点对应的采样点号；A(x_j,m_i,j+n,fc)表示当用慢度p_i进行扫描时第j个相邻地震道的处理时窗内的编号为n的采样点的数值；表示当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；w_x表示步骤8）定义的相邻地震道数；n_p表示步骤7）定义的慢度扫描个数。w_t表示步骤5）定义的时窗长度；f_c表示由步骤2）确定的地震剖面的主频；

12）对每一扫描慢度，根据10）计算的相邻道处理时窗对应的信号延迟时，对提取的处理时窗内的采样点的数值进行相位校正；

$B_{j,n}\left(p_i\right)={\tilde{B}}_{j,n}\left(p_i\right)e^{-2\mathrm{\πi}{f}_c{\mathrm{\Δ\τ}}_{i,j}}$

式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；表示当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；Δτ_i,j表示当用慢度p_i进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；f_c表示由2）确定的地震剖面的主频；

13）对每一扫描慢度，计算相位校正后的每个相邻地震道的处理时窗内的信号的能量E_j(p_i)：

$E_j\left(p_i\right)=\underset{n=-w_t/2}{\overset{w_t/2}{\mathrm{\Σ}}}{B}_{j,n}\left(p_i\right){\stackrel{\‾}{B}}_{j,n}\left(p_i\right)$

式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；为B_j,n(p_i)的共轭；E_j(p_i)表示当用慢度p_i进行扫描时第j个相邻地震道处理时窗内信号的能量；w_t表示由步骤5）定义的时窗长度；

14）对每一扫描慢度,计算参与慢度估算的各地震道处理时窗内的信号的复数互相关系数C_j,k(p_i)；

$C_{j,k}\left(p_i\right)=\underset{n=-w_t/2}{\overset{w_t/2}{\mathrm{\Σ}}}{B}_{j,n}\left(p_i\right){\stackrel{\‾}{B}}_{k,n}\left(p_i\right)/\sqrt{E_j\left(p_i\right)E_k\left(p_i\right)}$

-w_x/2≤j≤w_x/2，k>j

式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；p_i表示第i个扫描慢度；B_j,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；B_k,n(p_i)表示经过相位校正后的当用慢度p_i进行扫描时从第k个相邻地震道的处理时窗内提取的编号为n的采样点的数值；表示B_k,n(p_i)的共轭；E_j(p_i)表示当用慢度p_i进行扫描时第j个相邻地震道处理时窗内的信号的能量；E_k(p_i)表示当用慢度p_i进行扫描时第k个相邻地震道处理时窗内信号的能量；C_j,k(p_i)表示当用慢度p_i进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；w_t表示由步骤5）定义的时窗长度；w_x表示由步骤8）定义的相邻地震道数；

15）对每一扫描慢度,把其对应的各地震道处理时窗内信号的复数互相关系数进行求和，得到该慢度的总复数相关系数C(p_i)：

$C\left(p_i\right)=\underset{j=-w_x/2}{\overset{w_x/2}{\mathrm{\Σ}}}\underset{k>j}{\overset{w_x/2}{\mathrm{\Σ}}}{C}_{j,k}\left(p_i\right)$

式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；p_i表示第i个扫描慢度；C_j,k(p_i)表示当用慢度p_i进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；w_x表示由步骤8）定义的相邻地震道数；

16）对每一扫描慢度,计算该慢度对应的总复数相关系数的实部和虚部；

RC(p_i)=real(C(p_i))

IC(p_i)=image(C(p_i))

式中，i表示扫描慢度的编号；p_i表示第i个扫描慢度；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；RC(p_i)表示C(p_i)的实部；IC(p_i)表示C(p_i)的虚部；

17）根据步骤16）的计算结果，在步骤6）慢度扫描的最小值和最大值之间寻找使复数相关系数虚部为零的绝对值最小的扫描慢度q₁：

IC(q₁)=0

式中，IC(p_i)表示C(p_i)的虚部。C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；q₁表示使复数相关系数虚部为零的绝对值最小的扫描慢度；

18）在慢度0和使复数相关系数虚部为零的绝对值最小的零点之间，寻找使复数相关系数实部为最大的扫描慢度q₂，得到当前处理点的瞬时慢度；

RC(q₂)=maxRC(p_i)

式中，RC(p_i)表示C(p_i)的虚部；C(p_i)表示当用慢度p_i进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和，q₂表示使复数相关系数实部为最大的扫描慢度。

2.根据权利要求1所述的方法，特点是步骤8）所述的相邻地震道数w_x少于9大于3。