技术领域
[0001] 本
发明涉及地球物理勘探领域，更特别地，涉及石油地球物理勘探领域的一种
地震波分解方法。
背景技术
[0002]
地震勘探技术是利用弹性振动研究
地层中不同
岩石的弹性差异来认识地层的分布情况和地质构造状况的一种地球物理勘探方法，主要应用于石油、
天然气勘探开发、矿产资源普查、工程地质勘查等领域，也是石油勘探中寻找储油圈闭最普遍最有效的一种途径。由地震所提供的地下信息不仅可以有效地划分圈闭范围，而且还可以研究判断地层的岩石性质和
流体性质，甚至还可以扩展到油藏特征描述等方面。目前，存在多种对地震
信号进行处理的方法，其中地震信号多子波分解方法可以将采集到的地震信号(
地震道数据)分解为各种不同的子波，这些子波具有不同主频、不同子波宽度以及不同时间
位置。将这些子波
叠加在一起即可精确恢复原有的地震信号以用于地震信号分析。
[0003] 近年来匹配追踪技术应用到地震
信号处理中，将一个地震波分解成一系列子波，这些子波在子波时频分解中被称为时频单位。在匹配追踪过程中，为了合理地表征连续小波，需要利用以下五个参数：振幅A、
时移u、尺度(时间轴长度)σ、平均
频率ωm和
相位φ。常规的匹配追踪
算法是一个很麻烦的
迭代过程，需要从大量丰富的小波字典中选取出合适参数的子波。
[0004] 对于连续小波而言，重要参数之一是尺度σ，它控制了频域带宽和时域子
波长度。但是基于窗式傅式变换或者Gabor变换时进行谱分解却忽略了尺度这个可调参数，甚至于在一些匹配追踪技术应用中也被忽略了，所以现有的子波长度是不变的。比如基于Liu&Marfurt思想的则将小波宽度当作时中心频率的函数，在这种情况下实际就是常规
小波变换。因此现有的采用匹配追踪技术的地震波分解方法存在问题：一、由于大量翻阅小波字典，导致计算速度很慢；二、尺度σ为固定值，子波为固定长度。
发明内容
[0005] 为解决上述问题，本发明提供一种地震波分解方法，该方法包括：
[0006] 1)选择子
波形式，该子波形式为由公式(1)定义的Gabor小波形式：
[0007] 公式(1)，
[0008] 其中gr(t)是Gabor小波函数、t是时间变量、 是高斯窗函数、u是时移、σ是时间轴长度、ω是中心频率、φ是
相移、i表示虚数单位；
[0009] 2)获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式(2)所定义的形式：
[0010] 公式(2)，
[0011] 其中f(t)是地震信号的函数，n为分解次数，an是与第n次分解的小波 对应的(N+1) (1)振幅值，N为正整数，R f为第N次分解后的残差， 且R f＝f，标记
j为1至N的正整数；
[0012] 根据所述地震信号的复数属性来计算第j次分解时的参数集γj＝{uj，σj，ωj，φj}；
[0013] 在计算出所述参数集γj＝{uj，σj，ωj，φj}的
基础上计算第j次分解的小波的振幅值aj；
[0014] 3)将所述计算出的γj＝{uj，σj，ωj，φj}代入到公式(1)以得到在第j次分解出的小波 ，将计算出的振幅值aj与该小波 相乘以得到第j次分解出的子波。
[0015] 根据本发明的可替换方法，其中公式(1)可以由以下公式来替代：
[0016] 公式(6)，
[0017] 其中m(t)是Morlet小波函数、u是时移、σ是时间轴长度、ωm是平均
角频率、φ是相移。
[0018] 本发明提供的地震波分解方法，对传统的匹配追踪方法进行了改进，其尺度参数σ是可自适应调节的，通过σ的变化可以控制子波的长度。在对分解后得到的多子波进行分析中，还可以通过限制σ的大小来抑制脉冲或者正弦信号等
干扰信号。
附图说明
[0019] 图1是根据本发明的一个实施方式的一种地震波分解的方法的
流程图；
[0020] 图2是根据本发明的另一个实施方式的一种地震波分解的方法的流程图。
具体实施方式
[0021] 下面结合附图对本发明提供的地震波分解方法做进一步描述。
[0022] 如图1所示，根据本发明的一个实施方式，提供一种地震波分解方法，包括：
[0023] 1)选择子波形式，该子波形式为由公式(1)定义的Gabor小波形式：
[0024] 公式(1)，
[0025] 其中gr(t)是Gabor小波函数，w(·)表示高斯窗函数、u是时移、σ是时间轴长度、ω是中心频率、φ是相移、i表示虚数单位；选择子波形式可以是从本领域技术人员公知的Gabor小波字典、离散Dirac函数和基于复指数的离散傅里叶函数中选择匹配地震信号的成一系列小波或小波
原子，优选地本发明提供的地震波分解方法选择由公式(1)定义的Gabor小波形式作为子波形式。Gabor小波是一种时域的复函数，波形类似于复正弦信号与高斯信号包络相乘的结果。w(t)可以例如为
[0026] 2)获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式(2)所定义的形式：
[0027] 公式(2)，
[0028] 其中f(t)是地震信号的函数，n为分解次数，an是与第n次分解的小波 对应的(N+1) (1)振幅值，N为正整数，R f为第N次分解后的残差， 且R f＝f，标记
(j+1)
j为1-N的正整数；对于每一次分解，都有 其中 与R f
正交。
[0029] 根据所述地震信号的复数属性来计算第j次分解时的参数集γj＝{uj，σj，ωj，φj}，其中计算参数集γj＝{uj，σj，ωj，φj}包括：对初始信号进行Hilbert变换以构建复数信号，将该复数信号的最大振幅包络处所对应的时间作为时移uj，并将该最大振幅包络处对应的时刻处的瞬时频率和瞬时相位分别作为中心频率ωj和相移 根据公式(3)估算参数σj：
[0030] 公式(3)，
[0031] 其中D是小波集合， 表示R(j)f和 的内积， 为第j次分解的小波 的标准化；其中，如果j＝1，则所述初始信号为所获取的地震信号；如果j＞1，(j)
则所述初始信号为在第j-1次分解后的残差R f。
[0032] 具体地，在根据Hilbert变换得到参数uj，ωj，φj后，给定参数σj的范围，在该给定范围中选取多个值作为σj分别与参数uj，ωj，φj一起代入到公式(1)中，由此得到一组小波集合，设该小波集合为D。也就是说小波集合D由相同值的参数uj，ωj，φj和一组不同值的参数σj确定的一组小波。然后根据公式(3)在小波集合D中寻找满足公式(3)的小波 ，该小波 中的σj即为估算出的σj。这样就避免了传统方法中需要翻阅小波字典。所述给定范围可以是例如1-20。在给定范围中多个值的选取可以根据估算
精度的需求而有所不同，例如，在给定范围中选择多个值，相邻两个值之间的间隔可以是0.5-2，优选为1。
[0033] 在计算出参数集γj＝{uj，σj，ωj，φj}的基础上计算第j次分解时的小波的振幅值aj，计算振幅值aj包括根据公式(4)来估算aj：
[0034] 公式(4)，
[0035] 其中
[0036] 且公式(4)满足条件 (公式5)；
[0037] 其中Γ表示复数域。
[0038] 这里，将计算出的参数集γj＝{uj，σj，ωj，φj}代入到公式(1)得到 ，然后根据以上公式(4)可以估算出aj。
[0039] 3)将所述计算出的γj＝{uj，σj，ωj，φj}代入到公式(1)以得到在第j次分解出的小波 ，将计算出的振幅值aj与该小波 相乘以得到第j次分解出的子波；
[0040] 其中所述分解完成条件可以为所述残差R(N+1)f的
能量小于所述地震信号平均能量的2％，优选为小于1％。地震信号的能量的计算方法可以采用本领域技术人员公知的方法，这里不再赘述。
[0041] 对于当前分解而言， ，其中 与R(j+1)f是正交函数，满足对于正交投影(公式3)，选取使
最大，即‖R(j+1)f‖2最小时的 ，因此，公式(4)满足残差最小的条件：
[0042] 公式(5)。
[0043] 本发明提供的地震波分解方法是多次迭代分解的过程。对于第一次迭代分解，有即 此时将f(t)作为初始信号并根据上述方法确定γ1＝{u1，σ1，ω1，φ1}和a1，由此得到第一次迭代分解时分解出的子波 。接(2)
下来，由于 则有 在第二次迭代分解，将R f作为
初始信号，根据上述方法确定γ2＝{u2，σ2，ω2，φ2}和a2，由此得到第二次迭代分解时分(3)
解出的子波 在形式上有 在第三次迭代分解，将R
f作为初始信号，并根据上述方法确定γ3＝{u3，σ3，ω3，φ3}和a3。依次类推，在经过N(N+1)
次迭代分解后就能得到N个子波和残差R f，即如公式(2)描述的形式。当残差满足一定的条件后，不再对残差进行分解，该地震波分解结束。残差满足的条件根据处理地震信号的类型、环境、分解精度的不同而有所不同，例如待分解的地震信号的类型、特性以及所需的分解精度等这些因素发生变化，残差满足的条件也可以发生变化。残差满足的条件可以是如上所述残差的能量小于所述地震信号平均能量的2％，优选地，可以是1％。在分解结束(N+1)
后，残差R f可以被当作噪音数据对其进行处理，例如将其过滤掉。
[0044] 可替换地，可以采用Morlet小波作为时频分解的子波单位，因为Morlet小波能够较好地表征弹性波在孔隙介质中传播所产生的能量衰减和速度频散。因此可以选择Morlet小波形式作为子波形式，该Morlet小波由公式(6)定义：
[0045] 公式(6)，
[0046] 其中m(t)是Morlet小波函数、t是时间变量、u是时移、σ是时间轴长度、ωm是平均角频率、φ是相移、i表示虚数单位。采用Morlet小波对地震信号进行分解的方法与采用Gabor小波的相同，根据上述的方法计算参数集γ＝{u，σ，ωm，φ}和an，即计算参数集γ＝{u，σ，ωm，φ}的方法与计算参数集γ＝{u，σ，ω，φ}的方法相同，具体为：
[0047] 1)选择子波形式，该子波形式为由(6)定义的Morlet小波形式：
[0048] 公式(6)，
[0049] 其中m(t)是Morlet小波函数、t是时间变量、u是时移、σ是时间轴长度、ωm是平均角频率、φ是相移、i表示虚数单位；
[0050] 2)获取地震信号，并对获取的地震信号进行多次分解，直至满足分解完成条件，分解后的地震信号的形式为公式(7)所定义的形式：
[0051] 公式(7)，
[0052] 其中f(t)是地震信号的函数，n为分解次数，an是与第n次分解的小波mn对应的(N+1) (0)振幅值，N为正整数，R f为第N次分解后的残差， 且R f＝f，标记
j为1至N的正整数；
[0053] 根据所述地震信号的复数属性来计算第j次分解时的参数集 ；其中计算第j次分解时的参数集 包括：对初始信号进行Hilbert变换以
构建复数信号，将该复数信号的最大振幅包络处所对应的时间作为时移uj，并将该最大振幅包络处对应的时刻处的瞬时频率和瞬时相位分别作为平均角频率 和相移 根据公式(8)估算参数σj：
[0054] 公式(8)，(j) (j)
[0055] 其中D是小波集合，〈R f，mj〉表示R f和mj的内积， 为第j次分解的小波mj的标准化；
[0056] 其中，如果j＝1，则所述初始信号为所获取的地震信号；如果j＞1，则所述初始(j)信号为在第j-1次分解后的残差R f。
[0057] 与Gabor小波的情况类似，在根据Hilbert变换得到参数uj， φj后，给定参数σj的范围，在该给定范围中选取多个值作为σj分别与参数uj， φj一起代入到公式(6)中，由此得到一组小波集合，设该小波集合为D。也就是说小波集合D由相同值的参数uj， φj和一组不同值的参数σj，确定的一组小波。然后根据公式(8)在小波集合D中寻找满足公式(8)的小波mj，该小波mj中的σj即为估算出的σj。所述给定范围可以例如是1-20。在给定范围中多个值的选取可以根据估算精度的需求而有所不同，例如，相邻两个值之间的间隔可以是0.5-2，优选为1。
[0058] 在计算出所述参数集 的基础上计算第j次分解时的小波的振幅值aj；计算第j次分解的小波的振幅值aj包括根据公式(9)来计算该振幅值aj：
[0059] 公式(9)，
[0060] 其中R(i)f＝aimi+R(i+1)f，
[0061] 且公式(9)满足条件 (公式10)；
[0062] 这里，将计算出的参数集 代入到公式(6)得到mj，然后根据以上公式(9)可以估算出σj。
[0063] 3)将所述计算出的 代入到公式(6)以得到在第j次分解出的小波mj，将计算出的振幅值aj与该小波mj相乘以得到第j次分解出的子波。
(j)
[0064] 与采用Gabor小波的情形类似，如果所述残差R f满足一定条件，则分解完成。残差满足的条件可以是残差的能量小于所述地震信号平均能量的2％，优选为1％。
[0065] 另外，采用Morlet小波将获得的原始地震信号(地震道数据)分解为多个子波，这样能够很明显地看出能量在时频域的分布。
[0066] 在该方法中除了可输出子波重构的数据外还可以输出子波的时
频谱，进而进行时频分析：
[0067] 将地震信号f(t)分解为一系列小波 n＝0，1，2，…，N-1，这样能够很明显地看出能量在时频域的分布。在时频域信号f(t)的能量
密度表达式为：
[0068] 公式(11)
[0069] 其中 是某一
选定子波 的Wigner分布。同样地，我们将时频域振幅谱定义为：
[0070] 公式(12)
[0071] 对于公式(6)所给出的Morlet小波m(t)，其时频域的Wigner分布可定义为：
[0072] 公式(13)
[0073] 其 中 为 m(t) 的 复 共 轭。 再 次 利 用 积 分 变 换 式可得：
[0074] 公式(14)
[0075] 因此其时频域振幅谱解析式可写为：
[0076]
[0077] 公式(15)
[0078] 其中 为第n个小波的平均频率， 为归一化因数。
[0079] 为了提高对地震信号进行分解时的分解速度，一个方面是要提高模值 或的计算速度，因为在上述计算参数σj和aj的公式(3)(公式(8))和公式(4)(公式(9))2 2
中都需要对 或 (‖mj‖ 或‖mj‖)进行计算。如果能提高 或 (‖mj‖
或‖mj‖)的计算速度，那么对地震波分解的速度会有很大提高。因为地震信号是实值函数，因此公式(1)可以改写成或替换成公式(16)：
[0080] 公式(16)
[0081] 其中，gr(t)是Gabor小波函数，w(·)是高斯窗函数、u是时移、σ是时间轴长度、ω是中心频率、φ是相移；
[0082] 同样地公式(6)可以改写成或替换成公式(17)：
[0083] 公式(17)
[0084] 其中m(t)是Morlet小波函数、u是时移、σ是时间轴长度、ωm是平均角频率、φ是相移。
[0085] 通过下面的公式(18)和公式(19)
[0086] 公式(18)
[0087] 公式(19)
[0088] 可以计算：
[0089] 公式(20)。
[0090] 根据公式(18)和公式(19)计算：
[0091]
[0092] 由此可以得出
[0093] 公式(21)
[0094] 以及类似地
[0095] 公式(22)
[0096] 在Gabor小波的情况中，将上述计算得到的参数集γj＝{uj，σj，ωj，φj}代入到公式(1)和公式(21)；将根据公式(1)得到的 (计算的振幅能量)和根据公式(21)得到的 利用公式(3)以重新计算σj；将根据公式(1)得到的 和根据公式(21)得到的利用公式(4)以重新计算aj，这一步就是匹配；用重新计算得到的σj和aj分别替换之前估算的σj和aj。根据Hilbert变换得到的参数uj，ωj，φj以及这里重新计算得到的参数σj和aj即为进行匹配后的参数，利用参数uj，ωj，φj以及这里重新计算得到的参数σj和aj进行更精确的追踪方法以重新计算γj＝{uj，σj，ωj，φj}和aj。所述追踪方法可以是本领域技术人员公知的
匹配追踪算法中的追踪方法。例如，如果对γj＝{uj，σj，ωj，φj}中的uj进行追踪，则可以先保持其它三个参数σj，ωj，φj的值不变，在给定范围(根据半时长表达式 计算宽度，其中lt为半时长，ε为最小振幅值，单位为dB，可以设定ε＝-120dB，该宽度即为给定范围)内进行扫描，扫描的间隔例如可以是1，在扫描范围中计算地震信号与gabor小波的内积，寻找其能量与小波模平方 的比值最大的点，此时对应的uj就是最优的uj。追踪其它三个参数σj，ωj，φj的方法与追踪uj的类似。在进行追踪得到γj＝{uj，σj，ωj，φj}后，可以根据该γj＝{uj，σj，ωj，φj}和公式(4)得到追踪后的aj。最后利用追踪后的γj＝{uj，σj，ωj，φj}和aj可以得到分解出的子波。
[0097] 与Gabor子波的情况类似，在Morlet小波的情况中，将上述计算得到的参数集代入到公式(6)和公式(22)；将根据公式(6)得到的mj(计算的振幅能量)和根据公式(22)得到的‖mj‖利用公式(8)以重新计算σj；将根据公式(6)得到的mj和根据公式(22)得到的‖mj‖2利用公式(9)以重新计算aj，这一步就是匹配；用重新计算得到的σj和aj分别替换之前估算的σj和aj。根据Hilbert变换得到的参数uj，φj以及这里重新计算得到的参数σj和aj即为进行匹配后的参数，再利用根据Hilbert变换得到的参数uj， φj以及这里重新计算得到的参数σj和aj进行更精确的追踪方法以重新计算 和aj。这里的追踪方法可以是本领域技术人员公知的匹配
追踪算法中的追踪方法。例如，如果对 中的uj进行追踪，则可以先保持
其它三个参数σj， φj的值不变，在给定范围(例如，根据半时长表达式
计算宽度，其中lt为半时长，ε为最小振幅值，单位为dB，可以设定ε＝-120dB，该宽度即为给定范围)内进行扫描，扫描的间隔例如可以是1，在扫描范围中计算地震信号与
2
morlet小波mj的内积，寻找其能量与小波模平方‖mj‖ 的比值最大的点，此时对应的uj就是最优的uj。追踪其它三个参数σj， φj的方法与追踪uj的类似，也是先保持参数集中的三个参数的值固定，对剩下的一个参数在给定范围内进行扫描，找
到满足公式(8)的最优的该参数值。在进行追踪方法得到 后，可以根据
该 和公式(9)得到追踪后的aj。最后利用追踪后的 和
aj可以得到分解出的子波。
[0098] 在追踪过程中也应用了公式(21)或(22)进行模值的计算。
[0099] 这样分别采用公式(21)和(22)在匹配和追踪过程中都可以提高公式(3)和公式(4)以及公式(8)和公式(9)的计算速度。
[0100] 在传统方法中由于没有小波模的解析表达式，需要大量查阅小波字典，通过多次的匹配来完成，所以运算效率较低，而且尺度参数σ为一固定值使得子波长度是固定的。使用本发明提供的地震波分解方法，应用小波模及时频谱的解析表达式进行计算，不需要大量查阅小波字典，从而提高分解速度；而且尺度参数σ可以调节，从而子波的长度是可以变化的；不仅能减少分解后的残差，还能丰富小波字典。
[0101] 在对地震信号进行分解之后，可以对分解出的各个子波进行处理，例如滤波、去噪等。对这些子波进行叠加就能重构地震信号。但是由于地震波分解时的频率间隔不可能无限的小，这就会使得一些信号缺失，由此使得重构的地震信号反射同相轴不连续。为此，本发明提供的一种地震波分解的方法还包括对多个子波进行重构，在地震波重构时进行反褶积处理以保证同相轴的连续性。所采用的反褶积处理可以是本领域技术人员公知的反褶积方法，优选为多道预测反褶积方法。
[0102] 另外，当尺度参数σ的值极小时，对应的子波为脉冲信号；当σ值很大时，对应的子波为正弦信号。在重构时可以通过限定尺度参数σ的范围来剔除尖脉冲、单频谐振等噪声。例如可以将σ＜0.4以及σ＞10所对应的子波过滤掉，利用σ滤波，可以使重构后的地震信号的时频谱更为清晰，能更有效地进行分析。
[0103] 本发明提供的方法可以通过
硬件、
软件或者硬件与软件的结合来实现，这些硬件为所属领域技术人员公知的器件，例如
检波器、
滤波器、处理器等；所述软件也为所属领域技术人员公知的载有采用计算机语言编程的程序的介质，例如
软盘、
硬盘、光盘以及移动硬盘等；采用的计算机语言为所属领域公知的语言。
