一种基于改进的多道匹配追踪算法的地震时频分析方法与流程

1.本发明涉及油气田地震勘探领域，具体涉及一种基于改进的多道匹配追踪算法的地震时频分析方法。


背景技术：

2.地震勘探以其覆盖范围广、探测深度大等优势，一直是油气资源勘探与开发的重要指导信息。通过炸药或可控震源在地表激发地震波，然后再经地层界面反射并由检波器接收，经过一定技术手段处理后，地震资料能够对地下空间实现有效成像。尤其是三维地震勘探引入后，大幅度提高了地震资料对地下空间的成像能力，为了解地下构造特征提供了重要数据。地震波在地下空间传播过程中，其波形频率、能量、相位等特征是随时间不断变化的。当遇到不同弹性或物性特征的地层，地震反射波会在多个方面表现出属性变化，例如瞬时频率等。为此，基于地震资料开展时频分析，获取的地震时频谱属性能够指导油气资源的勘探开发，其在开展地下储层，尤其是气藏识别中具有重要指导意义。
3.地震资料时频解释技术利用地震信号丰富的频率信息，减少了常规解释的不确定性，在薄储层及小尺度地质目标的识别中起到了重要作用。地震波在地层介质中传播，散射和吸收衰减等均与频率有关，信号的功率谱密度是时变的，严格意义上属于非平稳信号。时频分析(time frequency analysis)是这类非平稳信号分析的重要手段，在地震勘探数据分析中已有重要应用，如地震波能量衰减补偿、时变滤波、层序检测、提高分辨率、地震频谱分解、地震旋回分析、瞬时属性提取等。谱分解技术可对地震数据进行全频扫描，得到一定频率的描述目的层地质特征的数据体，已被广泛应用于岩性识别及烃类检测。时频分析是谱分解的基础，由于时间分辨率和频率分辨率的相互制约，常规时频分析技术如短时傅里叶变换、小波变换等不能同时获得高精度的时域和频域谱信息；s变换虽然结合了短时傅里叶变换与小波变换的优点，但其在频率间断处刻画的频率信息并不准确；而匹配追踪算法是一种具有自适应性的地震信号分解技术，可自适应搜索最佳匹配原子，并获得高精度时频谱。该算法的原理如下：首先建立冗余的子波库(原子核)，然后将信号自适应性地分解为一些原子的核，使得残余能量快速衰减，得到最佳匹配。该算法能够根据信号特点自适应地选择匹配原子，具有自适应性，可取得较高的时频分析精度。
4.然而，传统的基于匹配追踪的地震时频方法亦存在一定问题亟待解决。比如，其对于地震数据开展逐道处理，由于资料噪声和算法稳定性的影响，各道获取的时频谱可能存在差异性较大的问题，降低了结果横向连续性。此外，由于传统匹配追踪在计算地震瞬时属性方面仍然存在稳定性较差的问题，从而影响地震时频谱精度。为此，需要对传统匹配追踪算法进行优化改善，以提高地震时频分析结果的可靠性。


技术实现要素：

5.针对传统地震时频分析方法所存在的问题，本发明的目的是提供一种基于改进的多道匹配追踪算法的地震时频分析方法。
6.为了实现以上目的，本发明采用以下技术方案：
7.一种基于改进的多道匹配追踪算法的地震时频分析方法，该地震时频分析方法包括以下步骤：
8.获得叠后地震资料；
9.根据所述叠后地震资料提取当前地震道及其邻道信息，以计算对应的加权平均地震道数据；
10.对所述加权平均地震道数据采用优化算法提取其对应的时频原子的瞬时时频属性；
11.针对所述加权平均地震道数据开展迭代求解以获取其对应的时频属性，迭代直至残差小于给定阈值；
12.针对所述叠后地震资料逐道计算其加权平均地震道数据所对应的地震时频属性，迭代直至所有地震道遍历完毕；
13.输出所述叠后地震资料对应的地震时频信息并开展地震时频特征分析。
14.本发明在匹配追踪理论的基础上，针对传统单道时频分析存在横向连续性不佳的问题，引入多道匹配追踪算法以提高地震时频分析结果的横向连续性和稳定性，为高精度地震时频分析奠定基础。此外，对于传统基于匹配追踪的地震时频分析中瞬时属性(如瞬时频率、瞬时相位等)计算精度不高，稳定性较差的问题，本发明采用了一种改进的地震时频原子瞬时属性计算方法，大幅度提高了瞬时频率、瞬时相位等属性的计算精度和稳定性，抑制了频率负数的出现，使计算结果更具有地球物理意义。本发明联合了多道匹配追踪和瞬时属性优化计算算法以改进传统的基于匹配追踪的地震时频分析方法，大幅度提高了地震时频分析的精度，亦提高了地震时频谱的横向连续性和稳定性，对油气藏勘探和开发具有重要指导意义。
15.根据本发明的地震时频分析方法，优选地，所述叠后地震资料经过地震资料处理环节，包括：静校正、去噪、振幅补偿、速度分析、动校正和偏移。
16.根据本发明的地震时频分析方法，优选地，所述根据所述叠后地震资料提取当前地震道及其邻道信息，以计算对应的加权平均地震道数据的步骤中，对于任意一道地震数据，记为s0(t)，记其邻近地震道分别为[s1(t),s2(t),...,s
n-1
(t)]；所述加权平均地震道数据为通过如下公式计算：
[0017][0018]
其中，ai为加权系数，i＝0,1,...,n-1；0＜ai＜1，且所有ai的和为1。
[0019]
根据本发明的地震时频分析方法，优选地，所述对所述加权平均地震道数据采用优化算法提取其对应的时频原子的瞬时时频属性的步骤中，加入一定的阻尼因子以降低地震噪声影响。
[0020]
根据本发明的地震时频分析方法，优选地，所述对所述加权平均地震道数据采用优化算法提取其对应的时频原子的瞬时时频属性的步骤中，瞬时频率的计算公式如下：
[0021]
insfreq(t)＝(d λe)-1b[0022]
[0023]
其中，insfreq(t)表示瞬时频率，λ为阻尼因子，0＜λ＜1，通常远小余1接近0；e表示单位矩阵；r(t)为复数地震道的虚部，即地震道的希尔伯特变换结果。
[0024]
根据本发明的地震时频分析方法，优选地，所述针对所述加权平均地震道数据开展迭代求解以获取其对应的地震时频属性，迭代直至残差小于给定阈值的步骤具体包括：
[0025]
针对所述加权平均地震道数据，通过迭代求解的方法在每次迭代中计算地震道残余波形在特定时间位置处所对应的地震原子；并在此基础上，获取该原子的时频谱信息；每次迭代确定一个原子，并将该原子的时频谱与之前所有原子的时频谱相加，然后从当前地震残余波形中减去当前地震原子，并继续迭代直至地震道残余能量小于给定的阈值；此时，该道地震数据的时频谱即为其分解的所有原子的时频谱的总和。
[0026]
根据本发明的地震时频分析方法，优选地，所述针对所述叠后地震资料逐道计算其加权平均地震道数据所对应的地震时频属性，迭代直至所有地震道遍历完毕的步骤具体包括：
[0027]
对每一道地震数据，首先计算所述加权平均地震道数据；在此基础上，结合所述对所述加权平均地震道数据采用优化算法提取其对应的时频原子的瞬时时频属性的步骤，以及所述针对当前加权平均地震道数据开展迭代求解以获取其对应的时频属性，迭代直至残差小于给定阈值的步骤，迭代求解当前地震道的时频谱信息，并针对各道地震数据开展逐道迭代以求取整个地震剖面或地震体的时频属性特征。
[0028]
根据本发明的地震时频分析方法，优选地，所述输出所述叠后地震资料对应的地震时频信息并开展地震时频特征分析的步骤包括：
[0029]
将每道地震数据所计算的地震时频属性按照一定规律(即是指地震道的相对空间位置)组合后即可获取该叠后地震数据所对应的时频属性体。
[0030]
根据本发明的地震时频分析方法，优选地，所述时频属性体比对应的叠后地震数据高一个维度。
[0031]
根据本发明的地震时频分析方法，优选地，对于地震时频属性的分析，对比同一属性在不同地震主频情况下的变换特征，通过对比分析以确定地下某一空间位置上所具有的油气勘探潜力。
[0032]
本发明的有益效果包括：
[0033]
1)本发明采用多道匹配追踪技术，不仅优化了地震时频属性横向连续性和稳定性，亦提高了地震时频分析的计算效率；
[0034]
2)本发明采用改进的地震时频原子瞬时属性计算方法，提高了地震资料瞬时相位、瞬时频率等参数的计算精度，有利于压制地震时频属性异常值，提高结果稳定性；
[0035]
3)本发明通过联合多道匹配追踪技术和改进的瞬时属性计算方法，与传统地震时频分析技术相比，具有运行效率高、时频属性横向稳定性好、结果精度高等优点；
[0036]
4)本发明算法简单易行，运行效率高，易于形成软件功能模块推广应用。
附图说明
[0037]
图1是本发明实施例提供的一种基于改进的多道匹配追踪算法的地震时频分析方法流程图。
[0038]
图2a是本发明实施例1输入的两道实际地震数据之一。
[0039]
图2b是本发明实施例1输入的两道实际地震数据之二。
[0040]
图3a是本发明实施例1以图2a中地震道为输入所获取的瞬时频率属性。
[0041]
图3b是本发明实施例1以图2a中地震道为输入所获取的瞬时频率属性。
[0042]
虚线为传统方法计算的地震道瞬时频率，实线为改进算法计算的地震道瞬时频率。
[0043]
图4a-图4d是本发明实施例1以图2a和图2b中地震道为输入所获取的地震时频属性；
[0044]
图4a-图4b和图4c-图4d分别对应于图2a和图2b的地震道数据；
[0045]
其中图4a、图4c为传统方法计算的地震时频属性，图4b、图4d为本发明算法计算的地震时频属性。
[0046]
图5是本发明实施例2输入的叠后二维地震剖面。
[0047]
图6a-图6d是本发明实施例2以图5中叠后二维地震剖面为输入所获取的地震分频能量剖面；
[0048]
图6a-图6b对应主频为10hz，c-d)对应主频为20hz，其中图6a、图6c为本发明算法所获取的地震分频能量剖面，图6b、图6d为传统方法所获取的地震分频能量剖面。
[0049]
图7是本发明实施例2以图5中叠后二维地震剖面的某道数据为输入所获取的剩余能量百分比变化趋势；白框为当前算法在每次迭代中的能量变化规律；黑框为传统方法在迭代过程中的能量变化规律。
[0050]
图8是本发明实施例3输入的中国西部某工区叠后三维地震数据体。
[0051]
图9a-图9d是本发明实施例3以图8中叠后三维地震数据体为输入所获取的地震分频能量体；
[0052]
图9a-图9d中地震分频能量体所对应的地震主频分别为10hz、20hz、30hz以及40hz。
具体实施方式
[0053]
为了更清楚地说明本发明，下面结合优选实施例对本发明做进一步的说明。本领域技术人员应当理解，下面所具体描述的内容是说明性的而非限制性的，不应以此限制本发明的保护范围。
[0054]
本发明是在研究了以下问题的基础之上提出的：
[0055]
(1)传统匹配追踪方法忽视了地震资料的横向相关性特征，其所计算的地震时频属性通常横向稳定性较差，影响地震时频分析精度；本发明通过引入多道匹配追踪方法，在计算当前地震道时频谱过程中引入其地震邻域信息，可提高地震时频属性计算的横向稳定性；
[0056]
(2)传统方法未考虑地震相位解缠过程中相位的单调递增性，其所计算的地震瞬时频率存在负值频率现象，不符合实际地震勘探的物理学特征；因此本发明在地震相位计算过程中通过开展相位校正，确保解缠过程中相位的单调递增性，可消除地震瞬时频率负值出现的可能性；
[0057]
(3)在实际地震数据中，由于噪声等因素的影响存在相位急剧变换的特征，从而引起地震瞬时相位属性的突变，进而产生频率异常值；本发明通过在瞬时频率计算过程中加
入一定的阻尼因子，提高瞬时属性计算结果的稳定性，从而达到压制异常值，提高地震时频属性提取精度和稳定性的目的。
[0058]
本发明具体采取以下工作步骤来实现上述技术方案：
[0059]
获取叠后地震资料；
[0060]
根据所述叠后地震资料提取当前地震道及其邻道信息以计算对应的加权平均地震道数据；
[0061]
对上述获取的加权平均地震道数据采用优化算法提取其对应的时频原子的瞬时时频属性；
[0062]
针对当前加权平均地震道数据开展迭代求解以获取其对应的时频属性，迭代直至残差小于给定阈值；
[0063]
针对当前地震资料逐道计算其平均地震道数据所对应的地震时频属性，迭代直至所有地震道遍历完毕；
[0064]
输出该叠后地震资料对应的地震时频信息并开展地震时频特征分析。
[0065]
如图1所示，本发明实施例提供的基于改进的多道匹配追踪算法的地震时频分析方法具体包括：
[0066]
101、采集叠后地震资料：
[0067]
采集经过一系列地震资料处理环节的叠后地震资料，这些环节包括但不限于：静校正、去噪、振幅补偿、速度分析、动校正、偏移等，最终处理的叠后地震资料等同于自激自收条件下的采集数据，可记为s(t,x)，其中t表示地震波双程旅行时，x为地震道号。
[0068]
102、根据所述叠后地震资料提取当前地震道及其邻道信息，以计算对应的加权平均地震道数据：
[0069]
传统地震时频分析技术基于单道地震数据逐道计算各个地震道时频属性，由于地震噪声及算法稳定性的影响，其所计算的地震时频属性通常会出现横向不连续特征，影响高精度地震时频分析效果。对此，提出了一种基于多道匹配追踪的地震时频分析方法，其通过在当前地震道时频属性计算过程中引入邻域地震道信息，从而增强地震时频属性的横向稳定性和精度。
[0070]
对于任意一道地震数据，可记为s0(t)，记其邻近地震道分别为[s1(t),s2(t),...,s
n-1
(t)]。由于s0(t)与[s1(t),s2(t),...,s
n-1
(t)]的数据记录距离较近，其与邻道地震波形之间存在较强的相关性。为此可以通过加权平均融合当道数据与邻道信息以提高地震时频分析结果的横向稳定性。设加权平均数据为其可通过如下公式计算：
[0071][0072]
其中，ai为加权系数(i＝0,1,...,n-1)。对于a，可通过调整其大小以控制任意邻道对当前道时频属性的贡献权重；权系数值的大小一般与其和邻道距离有关，距离越小，权系数越大；反之亦然。
[0073]
103、对上述获取的加权平均地震道数据采用优化算法提取其对应的时频原子的瞬时时频属性：
[0074]
地震时频属性计算存在一个相位解缠的过程。实际资料中，由于地震噪声等影响，相位解缠所获取的地震瞬时相位不能保证满足单调递增的规律，计算获得的地震瞬时频率
存在负值，不合符实际地震勘探的物理学特征。此外，由于地震瞬时相位在某些时间位置存在相位突变，因此导致出现瞬时频率异常值，其远超出实际地震资料频率的变化范围，因此有必要改进当前地震瞬时属性的计算方法，以提高地震时频分析的精度和稳定性。
[0075]
设地震道所对应的复数地震道hs(t)为：
[0076][0077]
其中r(t)为复数地震道的虚部，即地震道的希尔伯特变换结果；j表示虚数单位。在此基础上，地震道的瞬时相位可表示为：
[0078][0079]
其中，unwrap(.)表示相位解缠。在获取瞬时相位insphase(t)的基础上，开展相位校正，满足瞬时相位在时间维度上的单调递增性：
[0080]
insphase(ti)≤insphase(t
i 1
)
ꢀꢀꢀ
(4)
[0081]
对于瞬时频率insfreq(t)，其与瞬时相位的关系可写为:
[0082][0083]
其中，insphase'(t)表示瞬时相位的导数，可通过如下公式计算：
[0084][0085]
故瞬时频率insfreq(t)可直接写为：
[0086][0087]
在此基础上，设：
[0088][0089]
则瞬时频率的计算表达式可进一步写为：
[0090]
insfreq(t)＝d-1bꢀꢀꢀ
(9)
[0091]
传统方法通常基于上式计算瞬时频率；当地震数据存在随机噪声等因素干扰的情况下，局部会出现相位值突变，从而导致瞬时频率计算的异常值。对此，通过在瞬时频率计算公式中加入一定的阻尼因子以降低地震噪声影响，从而提高地震时频分析的精度：
[0092]
insfreq(t)＝(d λe)-1bꢀꢀꢀ
(10)
[0093]
其中，λ为阻尼因子，e表示单位矩阵。通过调节阻尼因子的大小，可以有效提高瞬时频率提取的稳定性。在此基础上，进一步计算地震原子的尺度因子、振幅等参数，从而确定当前地震道在该位置处所对应的地震原子。
[0094]
104、针对当前加权平均地震道数据开展迭代求解以获取其对应的时频属性，迭代直至残差小于给定阈值：
[0095]
针对上述地震道，通过迭代求解的方法在每次迭代中计算地震道残余波形在特定时间位置处所对应的地震原子；并在此基础上，获取该原子的时频谱信息。每次迭代确定一
个原子，并将该原子的时频谱与之前所有原子的时频谱相加，然后从当前地震残余波形中减去当前地震原子，并继续迭代直至地震道残余能量小于给定的阈值。此时，该道地震数据的时频谱即为其分解的所有原子的时频谱的总和。
[0096]
105、针对当前地震资料逐道计算其平均地震道数据所对应的地震时频属性，迭代直至所有地震道遍历完毕：
[0097]
对每一道地震数据，首先通过结合其邻道数据计算加权平均地震道数据以提高时频谱的横向稳定性。在此基础上，结合103-104步骤迭代求解当前地震道的时频谱信息，并针对各道地震数据开展逐道迭代以求取整个地震剖面或地震体的时频属性特征。
[0098]
106、输出该叠后地震资料对应的地震时频信息并开展地震时频特征分析：
[0099]
将每道地震数据所计算的时频谱按照一定规律(即是指地震道的相对空间位置)组合后即可获取该地震数据所对应的时频属性体，其通常比对应的地震数据高一个维度：一道地震数据的时频谱为一个二维剖面，二维地震剖面对应的时频属性一般为一个三维数据体，而三维地震数据体所对应的时频属性则通常为一个四维数据，即为多个三维数据体。对于地震时频属性的分析，通常可以对比同一属性在不同地震主频情况下的变换特征，通过对比分析以确定地下某一空间位置上所具有的油气勘探潜力。例如，当地下存在含气储层时，高频地震波能量在此位置通常衰减较快，而低频能量则相对则衰减较慢，对比分析地震高低频能量衰减差异，可以精确确定地下气藏的位置和品质。
[0100]
实施例1
[0101]
本实施例输入如图2a和图2b所示的两道实际地震数据。
[0102]
图3a和图3b是以图2a和图2b中地震道为输入所获取的瞬时频率属性：虚线为传统方法计算的地震道瞬时频率，实线为改进算法计算的地震道瞬时频率。由图3a和图3b可知，由于传统方法不考虑相位解缠过程中相位的单调递增性，其所获取的地震瞬时频率在某些时间位置上出现频率负值；此外，由于个别位置处相位的急剧变换，传统方法获取的地震瞬时频率脱离实际参数范围出现异常值，不符合地震波传播的物理学特征。在改进方法中，通过保持解缠过程中相位的单调递增性，同时在计算瞬时频率过程中加入一定的阻尼因子(阻尼因子根据测试和经验确定，经过测试后再一个合适的范围内根据经验确定；本实施例中阻尼因子取值为0.03)以杜绝异常值的出现，其所获取的瞬时频率属性与传统方法相比，既不存在负值频率现象亦极大地抑制了异常值，故当前方案更有利用保证地震时频属性提取的精度和稳定性。
[0103]
图4a-图4d是以图2a和图2b中地震道为输入所获取的地震时频属性；图4a-图4b和图4c-图4d分别对应于图2a和图2b的地震道数据；其中图4a、图4c为传统方法计算的地震时频属性，图4b、图4d为本发明算法计算的地震时频属性。如图中箭头所示，在本发明实施例算法所获取的地震时频谱(图4b、图4d)中，频率维度的异常值被明显地压制，同时地震时频谱在水平方向上稳定展布且垂向变化亦稳定变化，其精度和稳定性均优于传统算法结果(图4a、图4c)。
[0104]
实施例2
[0105]
本实施例输入如图5所示的叠后二维地震剖面。
[0106]
图6a-图6d是本发明实施例2以图5中叠后二维地震剖面为输入所获取的地震分频能量剖面：图6a-图6b对应主频为10hz，图6c-图6d对应主频为20hz，其中图6a、图6c为本发
明算法所获取的地震分频能量剖面，图6b、图6d为传统方法所获取的地震分频能量剖面。如图中椭圆和箭头所示，由于采用了多道匹配追踪技术，其所获取的地震分频能量剖面(图6a、图6c)的横向稳定性明显优于传统算法结果(图6b、图6d)。此外，由于当前算法能够压制地震频率异常值，图(图6a、图6c)中的地震分频能量值较为稳定，在一定区间内呈规律稳定变化；而传统方法获取的地震分频能量在水平和垂向上都呈现出不稳定变化特征，体现了基于当前算法开展高精度地震时频分析的优势。
[0107]
图7是以图5中叠后二维地震剖面的某道数据为输入所获取的剩余能量百分比变化趋势：白框为当前算法在每次迭代中的能量变化规律；黑框为传统方法在迭代过程中的能量变化规律。如图所示，与传统方法相比，当前算法收敛速度较快，能够提高地震时频分析的效率。
[0108]
实施例3
[0109]
本实施例输入如图8所示的中国西部某工区叠后三维地震数据体。
[0110]
图9a-图9d是本发明实施例3以图8中叠后三维地震数据体为输入所获取的地震分频能量体；图9a-图9d中地震分频能量体所对应的地震主频分别为10hz、20hz、30hz以及40hz。如图所示，当前算法所获取的多个地震分频能量体的能量展布特征与地震资料变化趋势基本一致，横向稳定性好，无杂乱特征；其获取的地震分频能量体能够精确表征地震信号在不同地震主频上的能量变化趋势，对于油气资源勘探具有重要指导意义。
[0111]
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定，对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动，这里无法对所有的实施方式予以穷举，凡是属于本发明的技术方案所引伸出的显而易见的变化或变动仍处于本发明的保护范围之列。