一种基于动态匹配的抗假频地震数据插值方法、电子设备及存储介质

1.本发明属于地震勘探技术领域，具体涉及一种基于动态匹配的抗假频地震数据插值方法、电子设备及存储介质。


背景技术：

2.目前，随着国家经济建设的日益增长，对能源的需求量也是越来越大，特别是石油、天然气等资源的开采都离不开勘探物理学。而在野外勘探油气资源工作中，地震勘探一直发挥着极其重要的作用。但是，随着野外的环境和地质条件正在变得越来越复杂，导致地震采集工作越来越困难。由于自然环境和经济成本的约束，采集到的数据往往会出现不规则或规则缺失的现象，得到的频谱图中出现不同程度的空间假频现象，一般不规则缺失的数据转换到频域引起的假频信息转化为低幅值不相干的随机噪声，可通过设置阈值和稀疏迭代的方法进行消除。但是，规则缺失数据在频域中产生的假频是与真实频谱相近的，很难辨别出现真实频谱信息，一般的插值算法很难压制这类假频的干扰。这种假频现象会严重降低地震剖面图的分辨率，同时影响着地震后续处理工作，所以如何有效的进行抗假频插值尤其重要。
3.传统的插值算法如基于滤波器的方法，此类方法的原理是根据线性同向轴在f-k域中具有可预测这一性质，从低频信号中取出具有高频成分的滤波算子来实现地震数据的重构，比较常用的方法有spitz方法。但此类方法都是基于线性同相轴假设，对于弯曲同相轴其抗假频效果欠佳，在实际的复杂数据插值中更是无法起到抗假频的效果。对于弯曲同相轴地震数据目前大都是基于curvelet变换作为稀疏算子，并结合凸集投影算法来进行插值重构，但此方法在规则采样数据中其抗假频效果并不好，不能有效区分真实频谱与假频信息，同样也无法在实际数据中应用。
4.近年来，动态规划相关的算法在语音识别、模块匹配等领域广泛应用。借助于动态规划，该模型在数据的匹配上具有显著优势。而在地震数据插值技术中，现有的算法缺乏对数据抗假频插值的能力，能够适用于实际数据抗假频插值算法的研究是有必要的。


技术实现要素：

5.本发明要解决的问题是现有插值算法进行抗假频插值不准确的问题，提供一种基于动态匹配的抗假频地震数据插值方法、电子设备及存储介质。
6.为实现上述目的，本发明通过以下技术方案实现：
7.一种基于动态匹配的抗假频地震数据插值方法，包括如下步骤：
8.s1、获取规则采样的二维地震数据；
9.s2、利用动态时间规整方法对步骤s1得到的地震数据中相邻两列地震数据进行匹配，得到相邻两列地震数据的匹配结果数据；
10.s3、根据步骤s2得到的相邻两列地震数据的匹配结果数据计算局部斜率；
11.s4、根据步骤s3得到的局部斜率在局部方向进行方向线性插值；
12.s5、重复步骤s1-s4，遍历所有相邻两列地震数据进行方向线性插值，得到重构的地震数据。
13.进一步的，步骤s2的具体实现方法包括如下步骤：
14.s2.1、设置步骤s1得到的地震数据中相邻两列地震数据为两个一维地震道的时间序列s＝[s1,s2,...sn]∈rn，t＝[t1,t2,...,tn]∈rn，其中n为地震数据的采样点的个数；
[0015]
s2.2、设置p为一维地震道的时间序列s和一维地震道的时间序列t之间的匹配点序列，匹配公式为：
[0016]
p＝{p(1)＝(a1,b1),p(2)＝(a2,b2),...,p(k)＝(ak,bk),...,p(n)＝(an,bn)}
[0017]
其中，(ak,bk)表示s序列中的点和t序列中的点相匹配，ak,bk为1,...n中的任意一个数，k＝1,...,n；
[0018]
s2.3、利用动态时间规整方法对s和t进行匹配求解，输出两个序列的匹配结果，假设d为点和点间的距离，动态时间规整方法的模型表达式d(s,t,p)为：
[0019][0020]
s2.4、对于所有的1≤i≤n和1≤j≤n构建值的dtw表格，其中i是行位置，j是列位置，表示(s1,...,si)和(t1,...,tj)在dtw约束下的匹配值；
[0021]
s2.5、设置动态时间规整方法的模型约束条件，然后在约束条件下进行最小化求解，得到最优匹配值
[0022]
s2.6、根据步骤s2.5得到的最优匹配值在步骤s2.4构建的dtw表格中进行反向追踪，通过回溯处理得到最优匹配结果p
*
，计算公式如下：
[0023][0024]
进一步的，步骤s2.5的动态时间规整方法的模型约束条件为：
[0025]
s2.5.1、边界约束：p(1)＝(1,1)，p(n)＝(n,n)，则有s序列和t序列的a1＝b1＝1，an＝bn＝n；
[0026]
s2.5.2、单调性约束：设置p(k)＝(ak,bk)，p(k+1)＝(a
k+1
,b
k+1
)，则有
[0027][0028]
s2.5.3、步长约束：设置p(k)＝(ak,bk),p(k+1)＝(a
k+1
,b
k+1
),则有：
[0029][0030]
表示在匹配过程中，每次匹配一个时间步。
[0031]
进一步的，步骤s3的计算公式为：
[0032][0033]
其中，q(k)为s序列中的点和t序列中的点间的斜率，δx为地震数据的道间距。
[0034]
进一步的，步骤s4的根据步骤s3得到的局部斜率在在两列地震数据之间局部方向进行沿着斜率方向的插值，插值公式为：
[0035]
o(i)＝(s
(i+q(i)
·
δx)/2
+t
(i-q(i)
·
δx)/2
)/2
[0036]
其中，o(i)为插值后的i点的地震数据，s
(i+q(i)
·
δx)/2
为s序列中的第(i+q(i)
·
δx)/2个点，t
(i-q(i)
·
δx)/2
为t序列中的第(i-q(i)
·
δx)/2个点。
[0037]
电子设备，包括存储器和处理器，存储器存储有计算机程序，所述的处理器执行所述计算机程序时实现所述的一种基于动态匹配的抗假频地震数据插值方法的步骤。
[0038]
计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现所述的一种基于动态匹配的抗假频地震数据插值方法。
[0039]
本发明的有益效果：
[0040]
本发明所述的一种基于动态匹配的抗假频地震数据插值方法，将动态规划算法应用于地震数据插值中，将采集到的数据，相邻两列之间进行匹配，利用匹配信息得到局部斜率方向，然后在局部进行方向线性插值。有效的利用了原始数据的信息和结构特征，进而提高了插值的准确性和地震数据的重构效果。
[0041]
本发明所述的一种基于动态匹配的抗假频地震数据插值方法，在对地震数据进行插值时，借助于所计算出的方向信息，能够有效的控制插值的方向，从而达到很好的抗假频效果。
[0042]
本发明所述的一种基于动态匹配的抗假频地震数据插值方法，在模型求解过程中，由于只需要进行动态匹配然后利用方向信息进行线性插值，整个过程计算较为简单，相比传统的方法，在计算时间上有所下降，在保证抗假频插值效果的同时还提升了插值效率。
附图说明
[0043]
图1为本发明所述的一种基于动态匹配的抗假频地震数据插值方法的流程图；
[0044]
图2为具体实施方式一中的模拟地震数据图；
[0045]
图3为具体实施方式一中的模拟地震数据的频谱图；
[0046]
图4为具体实施方式一中的采样率为50％的模拟地震数据图；
[0047]
图5为具体实施方式一中的采样率为50％的模拟地震数据的频谱图；
[0048]
图6为具体实施方式一中的利用spitz方法插值后的地震数据图；
[0049]
图7为具体实施方式一中的利用spitz方法插值后的地震数据的频谱图；
[0050]
图8为具体实施方式一中的利用curvelet变换方法插值后的地震数据图；
[0051]
图9为具体实施方式一中的利用curvelet变换方法插值后的地震数据的频谱图；
[0052]
图10为具体实施方式一中的利用本发明方法插值后的地震数据图；
[0053]
图11为具体实施方式一中的利用本发明方法插值后的地震数据的频谱图；
[0054]
图12为具体实施方式一中的利用本发明方法插值后的地震数据的斜率图；
[0055]
图13为具体实施方式二中的实际采样的地震数据图；
[0056]
图14为具体实施方式二中的实际采样的地震数据的频谱图；
[0057]
图15为具体实施方式二中的利用spitz方法插值后的地震数据图；
[0058]
图16为具体实施方式二中的利用spitz方法插值后的地震数据的频谱图；
[0059]
图17为具体实施方式二中的利用curvelet变换方法插值后的地震数据图；
[0060]
图18为具体实施方式二中的利用curvelet变换方法插值后的地震数据的频谱图；
[0061]
图19为具体实施方式二中的利用本发明方法插值后的地震数据图；
[0062]
图20为具体实施方式二中的利用本发明方法插值后的地震数据的频谱图；
[0063]
图21为具体实施方式二中的利用本发明方法插值后的地震数据的斜率图。
具体实施方式
[0064]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及具体实施方式，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施方式仅用以解释本发明，并不用于限定本发明，即所描述的具体实施方式仅仅是本发明一部分实施方式，而不是全部的具体实施方式。通常在此处附图中描述和展示的本发明具体实施方式的组件可以以各种不同的配置来布置和设计，本发明还可以具有其他实施方式。
[0065]
因此，以下对在附图中提供的本发明的具体实施方式的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定具体实施方式。基于本发明的具体实施方式，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他具体实施方式，都属于本发明保护的范围。
[0066]
为能进一步了解本发明的发明内容、特点及功效，兹例举以下具体实施方式，并配合附图1-附图21详细说明如下：
[0067]
具体实施方式一：
[0068]
一种基于动态匹配的抗假频地震数据插值方法，包括如下步骤：
[0069]
s1、获取规则采样的二维地震数据；
[0070]
s2、利用动态时间规整方法对步骤s1得到的地震数据中相邻两列地震数据进行匹配，得到相邻两列地震数据的匹配结果数据；
[0071]
进一步的，步骤s2的具体实现方法包括如下步骤：
[0072]
s2.1、设置步骤s1得到的地震数据中相邻两列地震数据为两个一维地震道的时间序列s＝[s1,s2,...sn]∈rn，t＝[t1,t2,...,tn]∈rn，其中n为地震数据的采样点的个数；
[0073]
s2.2、设置p为一维地震道的时间序列s和一维地震道的时间序列t之间的匹配点序列，匹配公式为：
[0074]
p＝{p(1)＝(a1,b1),p(2)＝(a2,b2),...,p(k)＝(ak,bk),...,p(n)＝(an,bn)}
[0075]
其中，(ak,bk)表示s序列中的点和t序列中的点相匹配，ak,bk为1,...n中的任意一个数，k＝1,...,n；
[0076]
s2.3、利用动态时间规整方法对s和t进行匹配求解，输出两个序列的匹配结果，假设d为点和点间的距离，动态时间规整方法的模型表达式d(s,t,p)为：
[0077][0078]
s2.4、对于所有的1≤i≤n和1≤j≤n构建值的dtw表格，其中i是行位置，j是列位置，表示(s1,...,si)和(t1,...,tj)在dtw约束下的匹配值；
[0079]
s2.5、设置动态时间规整方法的模型约束条件，然后在约束条件下进行最小化求
解，得到最优匹配值
[0080]
进一步的，步骤s2.5的动态时间规整方法的模型约束条件为：
[0081]
s2.5.1、边界约束：p(1)＝(1,1)，p(n)＝(n,n)，则有s序列和t序列的a1＝b1＝1，an＝bn＝n；
[0082]
s2.5.2、单调性约束：设置p(k)＝(ak,bk)，p(k+1)＝(a
k+1
,b
k+1
)，则有
[0083][0084]
s2.5.3、步长约束：设置p(k)＝(ak,bk),p(k+1)＝(a
k+1
,b
k+1
),则有：
[0085][0086]
表示在匹配过程中，每次匹配一个时间步；
[0087]
s2.6、根据步骤s2.5得到的最优匹配值在步骤s2.4构建的dtw表格中进行反向追踪，通过回溯处理得到最优匹配结果p
*
，计算公式如下：
[0088][0089]
s3、根据步骤s2得到的相邻两列地震数据的匹配结果数据计算局部斜率；步骤s3的计算公式为：
[0090][0091]
其中，q(k)为s序列中的点和t序列中的点间的斜率，δx为地震数据的道间距；
[0092]
s4、根据步骤s3得到的局部斜率在局部方向进行方向线性插值；
[0093]
进一步的，步骤s4的根据步骤s3得到的局部斜率在在两列地震数据之间局部方向进行沿着斜率方向的插值，插值公式为：
[0094]
o(i)＝(s
(i+q(i)
·
δx)/2
+t
(i-q(i)
·
δx)/2
)/2
[0095]
其中，o(i)为插值后的i点的地震数据，s
(i+q(i)
·
δx)/2
为s序列中的第(i+q(i)
·
δx)/2个点，t
(i-q(i)
·
δx)/2
为t序列中的第(i-q(i)
·
δx)/2个点。
[0096]
s5、重复步骤s1-s4，遍历所有相邻两列地震数据进行方向线性插值，得到重构的地震数据。
[0097]
分别使用spitz算法和curvelet变换的插值方法与本实施方式的插值方法对50％规则采样的模拟地震数据进行插值操作，分别得到重建地震数据。
[0098]
附图2-附图3是一个大小为256*202的模拟地震数据(道间距＝10m,时间采样点t＝0.004s)及其频谱图，
[0099]
附图4-附图5是对其进行50％均匀采样后的地震数据及其频谱图，可以看到该数据具有严重的假频干扰。附图6-附图7、附图8-附图9、附图10-附图11分别是利用spitz方法、curvelet变换方法和本实施方式方法插值后的数据及其频谱图，从频谱图上可以看到本实施方式方法的效果更好，附图12是利用本实施方式方法计算出的该地震数据的斜率图。
[0100]
根据信噪比量化恢复质量，信噪比计算公式如下：
[0101][0102]
其中，x为完整地震数据，xn为重建地震数据。三种方法插值后的信噪比以及运行时间如表1所示。
[0103]
表1信噪比对比
[0104][0105]
从表1中可以看到，本实施方式方法在信噪比和运行时间上都优于另外两种方法。
[0106]
具体实施方式二：
[0107]
一种基于动态匹配的抗假频地震数据插值方法，包括如下步骤：
[0108]
s1、获取规则采样的二维地震数据；
[0109]
s2、利用动态时间规整方法对步骤s1得到的地震数据中相邻两列地震数据进行匹配，得到相邻两列地震数据的匹配结果数据；
[0110]
进一步的，步骤s2的具体实现方法包括如下步骤：
[0111]
s2.1、设置步骤s1得到的地震数据中相邻两列地震数据为两个一维地震道的时间序列s＝[s1,s2,...sn]∈rn，t＝[t1,t2,...,tn]∈rn，其中n为地震数据的采样点的个数；
[0112]
s2.2、设置p为一维地震道的时间序列s和一维地震道的时间序列t之间的匹配点序列，匹配公式为：
[0113]
p＝{p(1)＝(a1,b1),p(2)＝(a2,b2),...,p(k)＝(ak,bk),...,p(n)＝(an,bn)}
[0114]
其中，(ak,bk)表示s序列中的点和t序列中的点相匹配，ak,bk为1,...n中的任意一个数，k＝1,...,n；
[0115]
s2.3、利用动态时间规整方法对s和t进行匹配求解，输出两个序列的匹配结果，假设d为点和点间的距离，动态时间规整方法的模型表达式d(s,t,p)为：
[0116][0117]
s2.4、对于所有的1≤i≤n和1≤j≤n构建值的dtw表格，其中i是行位置，j是列位置，表示(s1,...,si)和(t1,...,tj)在dtw约束下的匹配值；
[0118]
s2.5、设置动态时间规整方法的模型约束条件，然后在约束条件下进行最小化求解，得到最优匹配值
[0119]
进一步的，步骤s2.5的动态时间规整方法的模型约束条件为：
[0120]
s2.5.1、边界约束：p(1)＝(1,1)，p(n)＝(n,n)，则有s序列和t序列的a1＝b1＝1，an＝bn＝n；
[0121]
s2.5.2、单调性约束：设置p(k)＝(ak,bk)，p(k+1)＝(a
k+1
,b
k+1
)，则有
[0122][0123]
s2.5.3、步长约束：设置p(k)＝(ak,bk),p(k+1)＝(a
k+1
,b
k+1
),则有：
[0124][0125]
表示在匹配过程中，每次匹配一个时间步；
[0126]
s2.6、根据步骤s2.5得到的最优匹配值在步骤s2.4构建的dtw表格中进行反向追踪，通过回溯处理得到最优匹配结果p
*
，计算公式如下：
[0127][0128]
s3、根据步骤s2得到的相邻两列地震数据的匹配结果数据计算局部斜率；步骤s3的计算公式为：
[0129][0130]
其中，q(k)为s序列中的点和t序列中的点间的斜率，δx为地震数据的道间距。
[0131]
s4、根据步骤s3得到的局部斜率在局部方向进行方向线性插值；
[0132]
进一步的，步骤s4的根据步骤s3得到的局部斜率在在两列地震数据之间局部方向进行沿着斜率方向的插值，插值公式为：
[0133]
o(i)＝(s
(i+q(i)
·
δx)/2
+t
(i-q(i)
·
δx)/2
)/2
[0134]
其中，o(i)为插值后的i点的地震数据，s
(i+q(i)
·
δx)/2
为s序列中的第(i+q(i)
·
δx)/2个点，t
(i-q(i)
·
δx)/2
为t序列中的第(i-q(i)
·
δx)/2个点。
[0135]
s5、重复步骤s1-s4，遍历所有相邻两列地震数据进行方向线性插值，得到重构的地震数据。
[0136]
分别使用spitz算法和curvelet变换的插值方法与本实施方式的插值方法对50％规则采样的实际地震数据进行插值操作，分别得到重建地震数据。
[0137]
附图13-附图14是一个大小为1251*386(50％均匀采样)的实际地震数据(道间距＝40m,时间采样点t＝0.004s)及其频谱图，可以看到该数据的假频现象更加严重。
[0138]
附图15-附图16、附图17-附图18、附图19-附图20分别是利用spitz方法、curvelet变换方法和本实施方式方法插值后的数据及其频谱图，从频谱图上可以看到spitz和curvelet变换的方法基本不能够去除假频的干扰，插值后的数据仍然具体严重的假频干扰，从附图19-附图20可以看到本实施方式方法对实际地震数据也具有良好的抗假频效果。并且三种方法的运算时间分别为16.73s、184.95s和5.29s，本实施方式方法更快。附图21是利用本实施方式方法计算出的该地震数据的斜率图。
[0139]
具体实施方式三：
[0140]
电子设备，包括存储器和处理器，存储器存储有计算机程序，所述的处理器执行所述计算机程序时实现具体实施方式一或二所述的一种基于动态匹配的抗假频地震数据插值方法的步骤
[0141]
本发明的计算机装置可以是包括有处理器以及存储器等装置，例如包含中央处理器的单片机等。并且，处理器用于执行存储器中存储的计算机程序时实现上述的基于creo
软件的可修改由关系驱动的推荐数据的推荐方法的步骤。
[0142]
所称处理器可以是中央处理单元(central processing unit，cpu)，还可以是其他通用处理器、数字信号处理器(digital signal processor，dsp)、专用集成电路(application specific integrated circuit，asic)、现成可编程门阵列(field-programmable gate array，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等。
[0143]
所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序(比如声音播放功能、图像播放功能等)等；存储数据区可存储根据手机的使用所创建的数据(比如音频数据、电话本等)等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(smart media card,smc)，安全数字(secure digital,sd)卡，闪存卡(flash card)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。
[0144]
具体实施方式四：
[0145]
计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现具体实施方式一或二所述的一种基于动态匹配的抗假频地震数据插值方法。
[0146]
本发明的计算机可读存储介质可以是被计算机装置的处理器所读取的任何形式的存储介质，包括但不限于非易失性存储器、易失性存储器、铁电存储器等，计算机可读存储介质上存储有计算机程序，当计算机装置的处理器读取并执行存储器中所存储的计算机程序时，可以实现上述的基于creo软件的可修改由关系驱动的建模数据的建模方法的步骤。
[0147]
所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、u盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(rom，read-only memory)、随机存取存储器(ram，random access memory)、电载波信号、电信信号以及软件分发介质等。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。
[0148]
需要说明的是，术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
[0149]
虽然在上文中已经参考具体实施方式对本技术进行了描述，然而在不脱离本技术的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，本技术所披露的具体实施方式中的各项特征均可通过任意方式相互结合起来使用，在本说明书中未对这些组合的情况进行穷举性的描述仅仅是出于省略篇幅和节
约资源的考虑。因此，本技术并不局限于文中公开的特定具体实施方式，而是包括落入权利要求的范围内的所有技术方案。