CN108897039B - 基于层位控制的dtw多波匹配方法 - Google Patents

Abstract

本发明公开一种基于层位控制的DTW多波匹配方法，针对传统的DTW算法由于多波数据的复杂性且当有一些信息的缺失时效果往往不理想的问题，本发明的方法在DTW算法的基础上加入了层位控制点，为算法匹配时的最佳匹配路径提供参照点，使得通过HCDTW算法最后能够得到较为理想的匹配的结果，并采用自适应的网格，以及二阶累加函数使得求出来的压缩后的PS波未见畸变，且与PP波匹配的匹配精度有较大提高，使最终得到的匹配路径更加平滑和符合真实的情况。

Description

基于层位控制的DTW多波匹配方法

技术领域

本发明属于多波勘探领域，特别涉及一种多波地震资料匹配技术。

背景技术

多波勘探技术最早起源于上世纪的30年代的前苏联，他们使用了多波勘探中的横波分量来进行油气勘探。但之后，由于理论基础的匮乏，以及采集设备，解释方法等的落后，在很长的一段时间里，多波勘探技术无法在理论上有着更进一步的前进。直到20世纪90年代，在采集技术不断提升的情况下，更多的研究人员投入到多波勘探的技术的研究中，也使得该技术不断进步，完善，也逐渐在实际中得到应用。

多波地震资料能否进行高精度的匹配是整个多波勘探技术的核心所在，随着对地震勘探以及反演等的准确度要求的提高，对于多波地震资料匹配的精确度也就越来越高。但还是受限制于数据采集的局限性，以及地质的复杂性，使得高精度的多波地震资料匹配成为了该技术的难点之一。怎样在上文所述的情况下，使得横纵波对应的层位能够相互匹配，并且能够获得隐藏的层位信息，是其技术的核心要点。伴随着而来的是，多波地震资料匹配相关的理论不断涌现。

1996年James E.G提出多波勘探技术的核心所在就是能不能将横波和纵波进行联合解释。1996年James E.G详细介绍了纵横波速度比的求取方法，并用最大相关法求取γ、平均γ、层间γ等，且γ可以用于短波长振幅反演。2001年，Richard等，使用了最大相关系数法将横波和纵波进行匹配。他们使用的方法主要是通过求取平均速度比，然后将横波和纵波在时间域进行匹配。2002年，Garotta等通过Gaiser提出的最大相关系数法来对纵横波进行匹配，并联合二者AVO属性来进行建模。最大相关系数法匹配时关注的是局部信息，因此只是局部的最优不能达到整体的最优。

2003年，Sergey F为了实现纵横波自动匹配，认为剖面到剖面的时间映射弯曲函数可以通过波时间函数除去振幅增益补偿函数来得到，利用最小二乘法优化全局最优解，并估算逐点的互相关函数、推导纵横波速度比。2003年Backus等为了减少人工干预，他们采用多次快速迭代的方法来实现层位的匹配，与之前不同的是他们采用的是最小二乘法来进行横纵波的匹配。2003年Kristiansen等通过人工拾取层位的方式，实现了多波匹配。他们将地震专家的经验成功带入到匹配流程中。而由于多波数据的复杂性以及一些信息的缺失，解释人员的参与对于匹配的效果有着极大的提升。最小二乘法与传统的最大相关系数法相比，虽然可以从整体上去优化匹配，得到一个全局的较优解，但是计算复杂度较高。

2004年Michael Nickle提出了一种匹配算法，为了估算PS波和PP波的时差，该算法采用了多次迭代来计算PS波和PP波属性，然后进行低通滤波；迭代足够次数可以求得高精度的纵横波速度比，进而在时间域高精度地匹配PS波和PP波。2008年，Yuan将模拟退火算法与最大相关系数法结合，在满足局部最优的情况下使得整体的匹配有着较好的效果。此方法无论是理论数据和实际应用下都取得了相当不错的结果。模拟退火算法优点是局部搜索能力强，运行时间较短，缺点是全局搜索能力差，容易受参数的影响。

1978年，Sakoe和Chiba在语音识别的背景下开发了一种不同的动态规划解决方案，并且今天被广泛地称为动态时间规整(DTW)。1983年，Anderson和Gaby提出了DTW对地球物理问题的几种应用，他们称这种算法为“动态波形匹配”。2003年，Keysers和Unger已经证明DTW对于估计多维图像中的位移最直接的扩展没有计算上可行的解决方案。具体而言，计算精确解所需的时间随着图像尺寸呈指数增长。因此，许多作者提出了解决这个棘手问题的实际解决方案。2010年Pishchulin提供了一个最近的总结。2013年，Hale提出了动态图像规整算法(DIW)，该算法将信号匹配中常用的动态时间规整算法(DTW)应用到地震数据处理中。这无疑为后来的研究者打开了新的思路。2014年，Hale基于Mottle开发的近似解决方案提出了平滑的动态时间规整算法，与原算法相比，新的算法对原数据进行采样，排除掉一些匹配系数较低的点，使得最后的结果更加平滑。在通过已知的，大的且快速变化的移位相关的图像对的测试中，Hale证明了DIW可以估计已知移位的准确性。在进一步的测试中，Hale表明DIW可以比基于局部互相关的方法更精确，特别是当时间或空间变化很快时。诸如Hall和Hale提出的用于估计时间推移地震图像中的偏移的互相关方法仅在偏移变化较慢的情况下才是准确的。传统的DTW算法由于多波数据的复杂性且当有一些信息的缺失时效果往往不理想。

发明内容

为解决上述技术问题，本发明提出一种基于层位控制的DTW多波匹配方法，在原有的DTW算法基础上加入了层位控制点，结合了地质专家解释，最后通过自适应的网格平滑匹配结果，使得匹配更加准确。

本发明采用的技术方案为：一种基于层位控制的DTW多波匹配方法，包括层位控制点约束。

进一步地，具体包括以下步骤：

S1、对PP波和PS波分别进行预处理；

S2、获取待匹配的PP波和PS波层位控制点，并根据层位控制点对待匹配的PP波和PS波进行分段；

S3、采用改进的距离累加公式计算每段匹配的累加距离矩阵；

S4、根据累加距离矩阵得到每个分段的匹配路径，通过将每个分段的匹配路径整合得到全局匹配路径；

S5、自适应选取种子点。

更进一步地，包括以下分步骤：

S21、在PP波上确定大层位位置，设置大网格道间隔，根据大层位位置和大网格间隔确定大网格的位置；

S22、在PS波上确定与PP波对应的大层位位置，并采用与步骤S21中相同的大网格道间隔来确定相应的大网格位置；

S23、根据大网格点位置，选取PP波上的大网格点作为基准，在对应的PS波上选择相应的大网格作为计算中心，取其上下各W个点，并选取相应的一段长度t与对应的PP波上一段长度t来计算这2W+1个点的最大相关系数并选择其最大的点，作为与PP波上的大网格点相匹配的层位点；

S24、使用大网格上的点将各层位插满，并用平滑函数对其做平滑处理，得到每一道数据的层位控制点。

进一步地，长度t＝2W+1。

进一步地，步骤S3改进的距离累加公式具体为：

for a＞＝1,b＜＝1,

i∈[2,N-1],j∈[2,M-1]

其中，d[i,j]表示累加距离矩阵，a和b表示权重因子，N表示匹配时一道PS波的点数，M表示匹配时一道PP波的点数。

进一步地，步骤S5具体为：

S51、记录下预处理波形中振幅极值点；

S52、筛选出间隔大于设置间隔阈值的极值点作为路径种子点；

S53、将全局匹配路径对应的值赋值给这些路径种子点；

S54、通过插值算法将路径补全，得到最终的匹配路径；

S55、根据最终的匹配路径对PS波进行压缩，得到压缩后的PS波，记作PS’；

S56、对PS’进行整数时间的重采样，得到PS”作为最终的结果。

本发明的有益效果：本发明的方法在DTW算法的基础上加入了层位控制点，为算法匹配时的最佳匹配路径提供参照点，使得通过HCDTW算法最后能够得到较为理想的匹配的结果，最后对匹配路径进行了自适应网格的筛选，使得得到的匹配路径更加平滑和符合真实的情况。

附图说明

图1为本发明的方案流程图；

图2为本发明实施例提供的本发明方法与传统DTW算法所得匹配路径的对比图；

图3为本发明实施例提供的PP波实际数据；

图4为本发明实施例提供的PS波实际数据；

图5为本发明实施例提供的实际数据下本发明方法与传统DTW算法路径比较；

图6为本发明实施例提供的未经自适应路径筛选的DTW匹配结果；

图7为本发明实施例提供的未经自适应路径筛选的本发明方法的匹配结果；

图8为本发明实施例提供的经自适应网格筛选的DTW匹配结果；

图9为本发明实施例提供的层位匹配结果；

图10为本发明实施例提供的经自适应网格筛选的本发明方法的匹配结果。

具体实施方式

为便于本领域技术人员理解本发明的技术内容，下面结合附图对本发明内容进一步阐释。

如图1所示为本发明的方案流程图，本发明的技术方案为：一种基于层位控制的DTW多波匹配方法，加入了层位控制，使得匹配更加准确。DTW为现有技术，本发明中不做详细阐述，本发明的方法具体包括以下步骤：

S1、对PP波和PS波分别进行预处理；

S2、获取待匹配的PP波和PS波层位控制点，并根据层位控制点对待匹配的PP波和PS波进行分段；

S3、采用改进的距离累加公式计算每段匹配的累加距离矩阵；

S4、根据累加距离矩阵得到每个分段的匹配路径，通过将每个分段的匹配路径整合得到全局匹配路径；

S5、自适应选取种子点。

步骤S1具体为：

由于地震数据的采集不可避免的导致信息缺失以及PP波和PS波传播机理的不同，所以在对其进行匹配前需要对数据进行一些预处理，以减少其它因素对于匹配的影响。本发明采用D.Hale中的预处理公式对PP波和PS波分别进行预处理：

其中，为原始数据的希尔伯特变换，f(t)为原始的波形数据，a(t)为求得的预处理后的波形。

步骤S2具体为：

在原始的PP波数据上的大层位上设置大网格标签，然后让地质专家对这些大网格标签来进行匹配，即是在对应的PS波数据上寻求相应的匹配点，构成层位控制点对。可以根据最大相关系数值(R.V.Dok,J.Gaiser.Stratigraphic description of the MorrowFormation using mode-converted shear waves:Interpretation tools andtechniques for three land surveys[J].Leading Edge,2001,20(9).)的大小来确定与PP波层位点所对应匹配的PS波上的数据点。之后对得到的层位点进行插值处理，并用平滑函数将其平滑，使得每一道地震数据上都有着匹配精准的层位控制点。这样，即保证了得到的层位控制点有着手工解释的准确度，也同时能够大大的缩短解释的时间，也使得各道之间的层位控制点呈现出道与道间的平滑性，使得分道匹配后的得到的压缩后的PS波不会出现跳变。本发明采用如下层位控制点获得的流程(参考：徐天吉.纵横波高精度匹配方法研究与应用[D].成都理工大学,2012.)：

S21、在PP波上确定大层位位置，设置大网格道间隔，根据大层位位置和大网格间隔确定大网格的位置；

S22、在PS波上确定与PP波对应的大层位位置，并采用与步骤S21中相同的大网格道间隔来确定相应的大网格位置；

S23、根据大网格点位置，选取PP波上的大网格点作为基准，在对应的PS波上选择相应的大网格作为计算中心，取其上下各W个点，并选取相应的一段长度t与对应的PP波上一段长度t来计算这个2W+1个点的最大相关系数并选择其最大的点，作为与PP波上的大网格点相匹配的层位点。

S24、使用大网格上的点将各层位插满，并用平滑函数对其做平滑处理，得到每一道数据的层位控制点。

步骤S3具体为：

现有的二阶对称递归求和函数公式：

将式(2)的递推公式中第五项进行改进：

增加了两个权重因子a和b，当a＝1,b＝1时，与原来的二阶求和公式(2)相同。当a＞1时，相当于加大了d[i-2,j-1]到d[i,j]的权值，由于这个方向的变化是PP波压缩两个位移量，PS波压缩一个位移量，而这种情况应该尽量避免。由于路径是求最小值因此系数a可以取2，相当于加大了压缩PP波这一方向的权值，回溯时路径就不会往这一方向走。a值可以根据实际数据的情况加以抉择。b＜1时，可以看出相当于降低了d[i-1,j-2]到d[i,j]的值，而这个变化方向则是想要匹配进行的方向，b的取值一般为PP波的长度与待匹配PS波的长度的比值，这样使得在匹配是充分考虑到匹配序列长度不一的问题，以及对匹配的方向上根据数据特性进行加权处理。

由公式(3)可以看出，当前位置的累加量并非仅仅由上一个位置的点定，还是由上一个位置点的再上一个位置决定的。即最小可以是两个位置点变化一个单位的移动量，可以认为最小的位移量为半个单位，这在精度上有较大的提升。在通过改进后的公式对多波进行匹配后，也不会出现PP波被压缩，PS波不被压缩的违背理论的情况。并且当前点累加量的计算与之前的五个点都有关系，而不再是使用一阶累加公式时的三个点。虽然这会略微增加计算的时间开销，但是这会使得得到当前点的累加量时更为可靠，使得最后的结果更加稳定，符合原始数据的特征。

步骤S4具体为:

假设获得的层位控制点对为H₁,H₂,…H_l,H_l+1,…,H_n，其中H_l表示第l个层位控制点在路径中的编号，令两个层位点H_l,H_l+1之间的最小代价路径为W_l，则带层位约束最小代价路径D为：

其中W_begin为起始点[0,0]到第一个层位控制点H₁之间的路径，W_l为两个层位点H_l,H_l+1之间的最小代价路径，W_end为最后一个层位控制点H_n到终点[N-1,M-1]之间的路径，w_i表示路径中第i个点的匹配关系，D为最终的路径。如图2所示，DTW算法得到的匹配路径，以及本发明方法层位点控制下的HCDTW算法所求得的匹配路径，图2中contol表示层位控制点。

步骤S5具体为：

HCDTW算法采用D.Hale(D.Hale.Dynamic warping of seismic images:Geophysics,2013,78,no.2,S105–S115.)中一种较为自适应的方式来对种子点进行筛选。在波形预处理时，通过公式(1)可以得到PS波处理之后的波形a(t)。处理之后的波形a(t)均为正振幅，可以选择其中波形振幅的极值点作为种子点，对这些点设置间隔阈值Δt，使它们不至于过于密集而相互影响，即满足：

将满足条件的a(t₀)路径种子点，这样得到的路径种子点可以最大程度涵盖原有路径的特征，之后，通过插值的方式得出最终的路径。路径优化筛选的流程如下：

1.记录下预处理波形a(t)中振幅极值点。

2.筛除间隔小于Δt的极值点。

3.选择筛选出来的点作为路径种子点，并将匹配路径对应的值赋给这些点。

通过插值算法将路径完善，得到最后结果。

以下通过具体仿真数据对本发明的内容做进一步阐述：

图3、图4为某地区的多波数据，作为对比，图3为PP波数据，图4为PS波数据。两组数据均截取相同的剖面进行说明对比，PP波数据和PS波数据均为第400inline号剖面，并截取1道至100道共100道数据，且PP波长度为441个采样点，采样率为2ms，时间为880ms，PS波长度为831个采样点，采样率为2ms，转换为时间即为1660ms。其中用黑色线标识出由专家解释的层位线，PP波层位分别记为ppHor1，ppHor2，ppHor3，与之对应的PS波层位记为psHor1，psHor2，psHor3。本实施例对PP波和PS波先进行振幅匹配，即将PS波数据按照其与PP波形比例进行压缩。本实施例选取第50道数据作为对比，在其上通过真实层位导出的数据作为层位控制点，之后分别对算法进行DTW匹配，和HCDTW匹配，并对匹配路径是否进行自适应路径筛选进行对比。

图5为第50道数据在DTW算法和HCDTW算法下的求得的匹配路径，其中3号实线的路径为HCDTW算法求得的结果，4号实线路径为DTW算法求得的结果。以层位控制点(Horizencontrol)为基准，可以看到DTW算法求得路径与层位控制点相比，离得较远，并且整体呈现上凸的趋势，这与理论分析的结果并不相符，因为随着地层深度的增加，纵波和横波的速度比是呈增大趋势，在图中大概可以用匹配路径的切线的斜率变化去描述速度比的变大或是变小。故使用DTW算法得到匹配路径不太符合横纵波地质特性，而HCDTW算法不仅符合层位控制点的趋势，也符合实际的地质特性。1号和2号虚线分别为未使用自适应网格的HCDTW算法(HCDTW Without Gird)和DTW算法(DTW Without Gird)得到的曲线，可以看见实线的使用自适应网格进行筛选后的路径较虚线路径更为平滑。自适应大网格间隔为25个位移量，在进行路径筛选优化后，得到了基于自适应大网格的匹配路径。由图5可见，经自适应网格得到的路径较每点进行匹配后得到路径更为平滑，两条路径在前面具有较大的相似性，符合预期。在最后一个控制点后，HCDTW算法与大网格优化后的路径比差距较大，这是由于随着地层深度的增加，在PS波和PP波数据的末端，采集到的数据已经不再具有高分辨率，再对每一个进行匹配势必会得到较大的误差，故采用基于自适应大网格的路径优化去克服这个问题。图5中的PS Sample Number表示PS波样本数，PP Sample Number表示PP波样本数。

图6、图7分别是未经自适应路径筛选的DTW算法和HCDTW算法的匹配结果。从图中可以看出，无论是图6还是图7中，右边PS波波形已经明显畸变，实验结果也说明对所有数据点进行DTW算法匹配会使得波形得到畸变，这与之前D.Hale(S.Compton,D.Hale.Estimating VP/VS ratios using smooth dynamic image warping[J].Geophysics,2014,79(6):1639-1643.)的论证结果一致，也说明对匹配后路径进行筛选的重要性。且在图6中层位错动较大，由于算法采用的是回溯的方式去寻求最佳路径，即是从下至上来进行路径的选择，所以在图6中层位底层层位错动较小，上层由于误差累计较大。图7中层位对齐，但由于未使用路径筛选，使得层位出现畸变的情况。图6、图7中的横坐标Trace number表示地震记录道数目；

图8为经过自适应网格筛选后DTW算法得到的匹配结果，从图中可以看出波形未见畸变，但层位与图6出现类似的不匹配的问题，说明在DTW算法中，仅仅对路径进行优化筛选并不能得到准确的匹配结果，还需要加入层位进行控制和约束。

图9为DTW算法层位匹配后的结果，其中图左边为原始的1至50道PP波数据，右边为经压缩并重采样后PS波数据，其排列方式为从第51道至100道排列。图10为本发明HCDTW算法匹配的结果。图9右边的压缩后PS波与压缩前的PS波波形图相比波形未见畸变，重采样保持了原有数据的特征。但与左边的PP波相比，矩形实线方框区域的层位之间的局域并未进行较好的匹配，存在同向轴错动的问题，与之进行对比的本发明的HCDTW算法得到的匹配结果，在相应的矩形实线方框区域内匹配效果较好，将强反射层位都进行了对齐。并且压缩后的PS波未见畸变，从整体上来看达到匹配的要求。本发明提出的HCDTW算法由于使用了自适应的网格，以及二阶累加函数使得求出来的压缩后的PS波未见畸变，且与PP波匹配的匹配精度有较大提高。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。

Claims (4)

1.一种基于层位控制的DTW多波匹配方法，其特征在于，包括层位控制点约束；具体包括以下步骤：

S1、对PP波和PS波分别进行预处理；

S2、获取待匹配的PP波和PS波层位控制点，并根据层位控制点对待匹配的PP波和PS波进行分段；

S3、采用改进的距离累加公式计算每段匹配的累加距离矩阵；步骤S3改进的距离累加公式具体为：

for a＞＝1,b＜＝1,

i∈[2,N-1],j∈[2,M-1]

其中，d[i,j]表示累加距离，a和b表示权重因子，N表示匹配时一道PS波的点数，M表示匹配时一道PP波的点数；

S4、根据累加距离矩阵得到每个分段的匹配路径，通过将每个分段的匹配路径整合得到全局匹配路径；

S5、自适应选取种子点。

2.根据权利要求1所述的一种基于层位控制的DTW多波匹配方法，其特征在于，包括以下分步骤：

S21、在PP波上确定大层位位置，设置大网格道间隔，根据大层位位置和大网格间隔确定大网格的位置；

S22、在PS波上确定与PP波对应的大层位位置，并采用与步骤S21中相同的大网格道间隔来确定相应的大网格位置；

S23、根据大网格点位置，选取PP波上的大网格点作为基准，在对应的PS波上选择相应的大网格作为计算中心，取其上下各W个点，并选取相应的一段长度t与对应的PP波上一段长度t来计算这2W+1个点的最大相关系数并选择其最大的点，作为与PP波上的大网格点相匹配的层位点；

S24、使用大网格上的点将各层位插满，并用平滑函数对其做平滑处理，得到每一道数据的层位控制点。

3.根据权利要求2所述的一种基于层位控制的DTW多波匹配方法，其特征在于，长度t＝2W+1。

4.根据权利要求3所述的一种基于层位控制的DTW多波匹配方法，其特征在于，步骤S5具体为：

S51、记录下预处理波形中振幅极值点；

S52、筛选出间隔大于设置间隔阈值的极值点作为路径种子点；

S53、将全局匹配路径对应的值赋值给这些路径种子点；

S54、通过插值算法将路径补全，得到最终的匹配路径；

S55、根据最终的匹配路径对PS波进行压缩，得到压缩后的PS波，记作PS’；

S56、对PS’进行整数时间的重采样，得到PS”作为最终的结果。