[发明专利]一种采用频变的复数相关系数估算瞬时慢度的方法

摘要

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

主权项

一种采用频变的复数相关系数估算瞬时慢度的方法，特点是具体实施步骤如下：1)采集得到二维地震剖面；2)对地震剖面进行频谱分析并确定地震剖面的主频fc；对地震剖面进行F‑K谱分析并确定地震剖面的最小和最大视速度vmin和vmax；3)在时间方向，把地震道投影到以地震剖面的主频fc为频率的富里叶变换的基函数上，得到投影后的地震剖面A(x,t,fc)：A(x,t,fc)=a(x,t)e-2πifct]]>式中，a(x,t)是二维地震剖面，x和t分别为地震道的空间和时间坐标，fc为地震剖面的主频；A(x,t,fc)为投影后的地震剖面；4)对投影后的地震剖面A(x,t,fc)先按空间道序再按时间采样点顺序逐点按步骤5)‑18)估算瞬时慢度；5)根据步骤2)确定的地震剖面的主频fc，用下式计算以采样点为中心的沿时间方向的处理时窗的长度：wt=1fcΔt]]>式中，fc为地震剖面的主频，Δt为地震剖面的时间采样率，wt为处理时窗的长度；6)根据步骤2)确定的地震剖面的最小和最大视速度vmin和vmax，计算出慢度扫描的最少值和最大值；pmin＝1/vmaxpmax＝1/vmin式中，pmin和pmax表示慢度扫描的最少值和最大值，vmin和vmax是由步骤2)确定的地震剖面的最小和最大视速度；7)根据慢度扫描个数，计算慢度扫描的增量：dp＝(pmax‑pmin)/(np‑1)式中，np为用户给出的慢度扫描个数；dp为慢度扫描的增量；pmin和pmax表示慢度扫描的最少值和最大值；8)以当前处理地震道为中心道，确定参与瞬时慢度估算的相邻地震道数wx；9)从慢度扫描的最少值开始到慢度扫描的最大值为止，按慢度扫描增量计算扫描慢度pi：pi＝(i‑1)*dp i＝1,2,…,np其中，i表示扫描慢度的编号，pi是第i个扫描慢度，dp表示由步骤7)计算的慢度增量；10)以当前处理地震道和处理时窗的中心点为基点，对每一扫描慢度按照下式计算其对应的相邻地震道上的处理时窗的中心点时间τi,j和时窗对应的信号的延迟时Δτi,j；τi,j＝τ0+(xj‑x0)*pi‑wx/2≤j≤wx/2,i＝1,2…..npΔτi,j＝(xj‑x0)*pi‑wx/2≤j≤wx/2，i＝1,2…..np式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0编号表示当前处理道；pi表示第i个扫描慢度；x0表示当前处理地震道的空间作标，xj表示第j个相邻地震道的空间作标；τ0表示当前处理地震道的处理时窗的中心点时间；τi,j表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δτi,j表示当用慢度pi进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；wx表示由步骤8)定义的相邻地震道数；np表示由步骤7)定义的慢度扫描个数；11)对每一扫描慢度根据步骤10)计算的当前道及其相邻道的处理时窗的中心点时间以及步骤5)定义的时窗长度，提取当前道及其相邻道的处理时窗内的采样点的数值B~j,n(pi)=A(xj,mi,j+n,fc)]]> mi,j＝τi,j/Δt 1≤i≤np，‑wx/2≤j≤wx/2，‑wt/2≤n≤wt/2式中，i表示扫描慢度的编号；j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；pi表示第i个扫描慢度；τi,j表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗的中心点时间；Δt为地震剖面的时间采样率；mi,j表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗的中心点对应的采样点号；A(xj,mi,j+n,fc)表示当用慢度pi进行扫描时第j个相邻地震道的处理时窗内的编号为n的采样点的数值；表示当用慢度pi进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；wx表示步骤8)定义的相邻地震道数；np表示步骤7)定义的慢度扫描个数；wt表示步骤5)定义的时窗长度；fc表示由步骤2)确定的地震剖面的主频；12)对每一扫描慢度，根据10)计算的相邻道处理时窗对应的信号延迟时，对提取的处理时窗内的采样点的数值进行相位校正；Bj,n(pi)=B~j,n(pi)e-2πifcΔτi,j]]>式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；pi表示第i个扫描慢度；表示当用慢度pi进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；Bj,n(pi)表示经过相位校正后的当用慢度pi进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；Δτi,j表示当用慢度pi进行扫描时第j个相邻地震道与当前处理道的处理时窗中心点之间的时间延迟时；fc表示由2)确定的地震剖面的主频；13)对每一扫描慢度，计算相位校正后的每个相邻地震道的处理时窗内的信号的能量Ej(pi)：Ej(pi)=Σn=-wt/2wt/2Bj,n(pi)B‾j,n(pi)]]>式中，i表示扫描慢度的编号，j表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；pi表示第i个扫描慢度；Bj,n(pi)表示经过相位校正后的当用慢度pi进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；为Bj,n(pi)的共轭；Ej(pi)表示当用慢度pi进行扫描时第j个相邻地震道处理时窗内信号的能量；wt表示由步骤5)定义的时窗长度；14)对每一扫描慢度,计算参与慢度估算的各地震道处理时窗内的信号的复数互相关系数Cj,k(pi)；Cj,k(pi)=Σn=-wt/2wt/2Bj,n(pi)B‾k,n(pi)/Ej(pi)Ek(pi)]]> ‑wx/2≤j≤wx/2，k>j式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；n表示处理时窗内的采样点的编号，0编号表示处理时窗的中心采样点；pi表示第i个扫描慢度；Bj,n(pi)表示经过相位校正后的当用慢度pi进行扫描时从第j个相邻地震道的处理时窗内提取的编号为n的采样点的数值；Bk,n(pi)表示经过相位校正后的当用慢度pi进行扫描时从第k个相邻地震道的处理时窗内提取的编号为n的采样点的数值；表示Bk,n(pi)的共轭；Ej(pi)表示当用慢度pi进行扫描时第j个相邻地震道处理时窗内的信号的能量；Ek(pi)表示当用慢度pi进行扫描时第k个相邻地震道处理时窗内信号的能量；Cj,k(pi)表示当用慢度pi进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；wt表示由步骤5)定义的时窗长度；wx表示由步骤8)定义的相邻地震道数；15)对每一扫描慢度,把其对应的各地震道处理时窗内信号的复数互相关系数进行求和，得到该慢度的总复数相关系数C(pi)：C(pi)=Σj=-wx/2wx/2Σk>jwx/2Cj,k(pi)]]>式中，i表示扫描慢度的编号；j和k表示相邻地震道的编号，0序号表示当前处理道；pi表示第i个扫描慢度；Cj,k(pi)表示当用慢度pi进行扫描时第j个和第k个相邻地震道的处理时窗内信号的复相关系数；C(pi)表示当用慢度pi进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；wx表示由步骤8)定义的相邻地震道数；16)对每一扫描慢度,计算该慢度对应的总复数相关系数的实部和虚部；RC(pi)＝real(C(pi))IC(pi)＝image(C(pi))式中，i表示扫描慢度的编号；pi表示第i个扫描慢度；C(pi)表示当用慢度pi进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；RC(pi)表示C(pi)的实部；IC(pi)表示C(pi)的虚部；17)根据步骤16)的计算结果，在步骤6)慢度扫描的最小值和最大值之间寻找使复数相关系数虚部为零的绝对值最小的扫描慢度q1：IC(q1)＝0式中，IC(pi)表示C(pi)的虚部；C(pi)表示当用慢度pi进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和；q1表示使复数相关系数虚部为零的绝对值最小的扫描慢度；18)在慢度0和使复数相关系数虚部为零的绝对值最小的零点之间，寻找使复数相关系数实部为最大的扫描慢度q2，得到当前处理点的瞬时慢度；RC(q2)＝maxRC(pi)式中，RC(pi)表示C(pi)的虚部；C(pi)表示当用慢度pi进行扫描时所有相邻道的处理时窗内的信号的复数相关系数的和，q2表示使复数相关系数实部为最大的扫描慢度。