CN105242306A - 基于空间克里金插值的高精度多波匹配方法 - Google Patents

Abstract

本发明公开了一种基于空间克里金插值的高精度多波匹配方法，包括初始匹配和精细匹配两个步骤，初始匹配包括以下步骤：S1、对最初的gamma数据体进行三角剖分，构建三维Delaunay三角网；S2、对于构建的Delaunay三角网，通过空间索引的方法对待插值点进行空间定位；S3、对每一个待插值点进行克里金插值运算，计算出每个点的gamma值，采用重采样算法对PS数据进行初始匹配；精细匹配包括以下步骤：S4、对初始匹配后的PS数据和原始PP数据，通过控制层位点以及划定时窗的方式完成精细匹配。本发明通过空间克里金插值算法对初始的gamma数据体进行插值，采用重采样算法对PS数据进行初始匹配；然后进行精细匹配，提高了多波匹配的精度以及匹配的效率。

Description

基于空间克里金插值的高精度多波匹配方法

技术领域

本发明属于地质解释领域中多波匹配问题，特别涉及一种基于空间克里金插值的高精度多波匹配方法。

背景技术

多波地震勘探是进行岩性油气藏和隐蔽油气藏勘探的一种非常有潜力的手段，但是，由于诸多原因，多波多分量理论研究和油气田实际勘探地质需求的结合问题、复杂条件下的转换波地震资料处理问题和多波综合解释、全波属性的地质应用等问题一直没有取得显著进展，并且已经成为制约多波地震勘探技术进一步发展的“瓶颈”。而解决这些问题的基础是做好多波多分量资料处理，提供高质量的各向同性和各向异性处理成果。其中,多波传播机理的基础研究、多波资料中的纵横波匹配方法研究是目前多波地震资料后续处理的重点和难点，是多波精确成像和叠前纵横波联合反演及岩性识别、储层预测和含气性识别的重要基础,是体现多波多分量地震勘探技术实际勘探开发应用价值的关键。因此，基于多波传播机理,研究纵横波高精度的匹配新方法，有利于充分利用多波多分量地震资料、准确认识多波地质响应特征,突出多波多分量地震资料解决地质问题的能力，具有重大的意义。

目前多波匹配有基于反射特征的匹配方法和基于多波层位的匹配方法，前者通过PP波和PS波地震数据的波形和波组特征进行对比生成γ₀值，然后基于该γ₀值实现两者的时间域匹配。后者首先分别基于PP波和PS波地震数据追踪解释出对应的层位，然后通过标定对应的层位产生时移体，最后将时移体应用于PS地震数据，实现PS与PP地震数据的匹配。

随着勘探目标要求的提高，多波匹配技术的研究越来越受到人们的重视，纵横波匹配技术已经成为地球物理学的研究热点。JamesE.G(1996)详细介绍了纵横波速度比(γ₀)的求取方法，并用最大相关法求取γ₀、平均γ₀、层间γ₀等，使用VSP资料验证从PP波和PS波剖面中求取γ₀，且γ₀可以用短波长振幅反演。Wai-KinChan等(1997)利用常数γ₀值多次试算法，在时间对数域实现纵、横波匹配。由于浅层和深层的γ₀值不同，因此该匹配方法只能适应于特定的目的层。RichardV.D和JamesG等(2001)利用最大相似性原理，扫描PP波和PS波的γ₀谱，拾取平均γ₀值，将二者在时间域进行匹配。VanDok和Kristiansen(2003)采用人工手动拾取强同相轴和其他明显的地质特征，实现了PP波和PS波在时间域的粗略匹配。MichaelV.D等(2003)从PP波和PS波中求得纵横波速度比和泊松比，基于建立的横波模型，在深度域实现纵横波匹配。就纵横波匹配技术软件现状而言，Transform软件的匹配功能相对完善，提供的基于反射特征的匹配方法和基于多波层位的匹配方法都能实现PP数据和PS数据的匹配。纵观国内外纵横波匹配技术的发展，几乎所有的方法均是直接以PP波和PS波同一地层反射同相轴相似性最大的原则求取纵横波的速度比，直接对PP波和PS波在时间域进行匹配，其匹配程度是比较粗略的。

发明内容

本发明的目的在于克服现有技术，提供一种通过空间克里金插值算法对初始的gamma数据体进行插值，采用重采样算法对PS数据进行初始匹配；然后通过控制层位点以及划定时窗的方式完成PS数据和PP数据的精细匹配，完成最终多波匹配，提高了多波匹配的精度以及匹配的效率的基于空间克里金插值的高精度多波匹配方法。

本发明的目的是通过以下技术方案来实现的：基于空间克里金插值的高精度多波匹配方法，包括初始匹配和精细匹配两个步骤，通过基于空间克里金插值的算法完成对原始PS数据进行初始匹配，包括以下步骤：

S1、对最初的gamma数据体进行三角剖分，构建三维Delaunay三角网；

S2、对于构建的Delaunay三角网，通过空间索引的方法对待插值点进行空间定位；

S3、对每一个待插值点进行克里金插值运算，计算出每个点的gamma值，根据数据体内所有点的gamma值，采用重采样算法对PS数据进行初始匹配；

所述的精细匹配包括以下步骤：

S4、对初始匹配后的PS数据和原始PP数据，通过基于层位约束的精细匹配算法，通过控制层位点以及划定时窗的方式完成精细匹配。

进一步地，所述的步骤S1中，以人机交互设置的数据点和自动追踪算法计算出来的数据点作为最初的gamma数据体；其实现方法具体包括以下子步骤：

S11、根据给定的点集的区域构建一个初始三角形，使该三角形能够包含所有的数据点，作为初始Delaunay三角网；

S12、从点集中选择一个点进行插入操作，在已有的Delaunay三角网中找到包含该点的三角形，将该点与包含它的三角形的顶点连接形成三个新的三角形；

S13、根据Deluanay三角网的外接空圆特性，利用Lawson提出的LOP算法对当前Delaunay三角网进行优化，使之满足Delaunay特性；

S14、重复步骤S12和S13，当点集中所有点都包含在Delaunay三角网中时结束；

S15、将与初始三角形连接的所有三角形删除。

进一步地，所述的步骤S2包括以下子步骤：

S21、对于当前四面体，计算它的四个顶点的最大和最小Z值，进而计算出Z方向的索引号范围；

S22、按照Z值从小到大的顺序，将Z固定，依次对该四面体进行切割，计算该Z平面与四面体的交点；

S23、若求得交点个数为1，则该点定位完成，若交点个数为3或4，则按照Y值大小将交点排序，交点构成三角形或者四边形；

S24、计算出Y方向的索引号范围，按照Y值从小到大的顺序，将Y值固定，依次对三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点X₁，X₂，则落在X₁，X₂范围内的点即为下一个待处理的四面体；

S25、对于下一个四面体重复步骤S21～S24，直到遍历完所有四面体。

进一步地，所述的步骤S3包括以下子步骤：

S31、进行克里金插值运算：

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_{i}Z\left(x_i\right)---\left(1\right)$

其中区域化变量为Z(x)，待插值点为x₀，样本点记为x_i(i＝1,2,......,n)，n为样本点的个数，在点x_i处的属性值记为Z(x_i)，则点x₀处的属性值Z(x₀)是由各个样本点属性值加权求和得到的，其中，λ_i为待定权系数；

S32、在无偏估计和区域化变量满足二阶平稳的条件下，为了使得到的估计方差值达到最小，得到加权系数的方程组，即克里金方程组：

$\left\{\begin{array}{c}\underset{j=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_{j}C(x_i,x_j)-\mathrm{\μ}=\stackrel{\‾}{C}(x_i,V),(i=1,2...,n)\\ \underset{i=1}{\overset{n}{\Σ}}{\mathrm{\λ}}_i=1\end{array}\right.---\left(2\right)$

在实际插值过程中，为了对方程组求解得到权重系数，将方程组转换为矩阵的形式进行表达，之后利用矩阵运算进行求解，便于计算，矩阵形式如式(3)所示：

[K]·[λ]＝[M](3)

在式(3)中，K矩阵的表达式如下：

$\[K\]=\left[\begin{array}{ccccc}\mathrm{\γ}(v_1,v_1)& \mathrm{\γ}(v_1,v_2)& ...& \mathrm{\γ}(v_1,v_n)& 1\\ \mathrm{\γ}(v_2,v_1)& \mathrm{\γ}(v_2,v_2)& ...& \mathrm{\γ}(v_2,v_n)& 1\\ .& .& & .& .\\ .& .& & .& .\\ .& .& ...& .& .\\ \mathrm{\γ}(v_n,v_1)& \mathrm{\γ}(v_n,v_2)& ...& \mathrm{\γ}(v_n,v_n)& 1\\ 1& 1& ...& 1& 0\end{array}\right]---\left(4\right)$

其中，γ(v_n,v_n)表示第n个已知点与第n个已知点之间的变差函数值，在K矩阵中，前n行n列的值是两两已知点之间的变差函数值；

M矩阵的表达式如下：

$\[M\]=\left[\begin{array}{c}\mathrm{\γ}(v_1,v)\\ \mathrm{\γ}(v_2,v)\\ .\\ .\\ .\\ \mathrm{\γ}(v_n,v)\\ 1\end{array}\right]---\left(5\right)$

γ(v_n,v)表示第n个已知点与未知点之间的变差函数值，M矩阵的前n行表示待插值点与每个已知点之间的变差函数值；

λ矩阵的表达式如下：

$\[\mathrm{\λ}\]=\left[\begin{array}{c}{\mathrm{\λ}}_1\\ {\mathrm{\λ}}_2\\ .\\ .\\ .\\ {\mathrm{\λ}}_n\\ \mathrm{\μ}\end{array}\right]---\left(6\right)$

λ矩阵即权重系数矩阵，由λ＝[K]⁺[M]计算得到，[K]⁺表示的是K矩阵的逆矩阵；

S33、将λ矩阵带入公式(1)中，计算点x₀处的属性值Z(x₀)，即为该点的gamma值；

S34、重复步骤S31～S33，计算所有待插值点的gamma值；

S35、根据所有待插值点的gamma值，利用重采样算法重新计算该道数据上每个点的数据值，完成初始匹配。

进一步地，所述的步骤S4具体包括以下子步骤：

S41、记录每道PS数据的层位点坐标；

S42、抽取一道PP数据和对应的一道初始匹配后的PS数据构成一个N行2列的数据面，其中N为每道PS数据上采样点的个数；

S43、设定时窗参数：N_w、N_m、N_s，其中，N_w为波形对比窗口的大小，N_m为波形对比窗口的移动量，N_s为波形对比窗口的搜索半径，N_w，N_m，N_s均为整数；以N_w为窗口大小，以N_m为窗口移动量，从上到下，滑动窗口，以窗口中心点作为种子点，计算每个种子点移动量；具体计算方法为：

设i_s＝(k-1)*N_m+N_w/2，i_e＝i_s+N_w，定义集合k为波形对比窗口个数； $\hat{L}=\{-N_s,-N_s+1,...,0,...,N_s-1,N_s\},$

定义长度为(N_w+1)的向量f，f的第i个元素

定义2N_s+1个长度为(N_w+1)的向量g_l，g_l的第i个元素

最优移动量为S(j₁,j₂)，其中j₁和j₂代表第j₁和j₂道；

则S(j₁,j₂)的计算方法为：

$S(j_1,j_2)=\arg \underset{l\∈\hat{L}}{max}\frac{f\·g_l}{\sqrt{f\·f}\sqrt{g_l\·g_l}}$

对于第k个波形对比窗口，计算拉平种子点移动量，定义长度为J的向量m用来存储每道PS数据拉平种子点移动量，所述的拉平种子点移动量为：m(j₂)＝S(j₁,j₂)-S(j₁,j₁)；

S44、增加对层位点的控制，对m(j₂)进行处理，此时每个窗口中心点的坐标为i_s＝(k-1)*N_m+N_w/2，如果该坐标是层位点的坐标，则设置m(j₂)＝0，否则不处理；

S45、重复步骤S41～S44，求出整个PS数据面上的种子点的初始移动量，这种方法下，每道PS数据上种子点的位置都一样，然后对每interval道数据同一位置的种子点进行相加求取平均值作为第interval/2道数据种子点的移动量，其中interval为道间隔；求出每第interval/2+i*interval道PS数据上种子点的移动量，其中0<＝i<＝(M-1)/interval，M为该PS数据面上PS数据的道数；

S46、依次对所求出的每interval/2+i*interval道PS数据上同一位置的种子点移动量进行横向平滑处理，然后线性插值出该PS数据面上其他道种子点的最终移动量；

S47、根据步骤S46得到的种子点的最终移动量，采用线性插值算法，分别计算出每道PS数据上每个采样点的移动量；

S48、每个采样点的新的采样坐标为X^new(i,j)＝i+m_i(j)；设置原始采样坐标，具体地，X(i,j)＝i,i＝1,2,...,I；j＝1,2；对于任意一个定值以X的第j列X_j为自变量、以D的第j列D_j为函数值，构建三次样条差值函数，则得到精细匹配后的PS数据。

本发明的有益效果是：本发明通过初始匹配和精细匹配两个步骤实现多波精细匹配，初始匹配中，通过空间克里金插值算法对初始的gamma数据体进行插值，根据插值过后的gamma数据体，采用重采样算法对ps数据进行初始匹配；然后利用初始匹配后的PS数据，采用精细匹配算法，通过控制层位点以及划定时窗的方式完成PS数据和PP数据的精细匹配；从而完成最终多波匹配，提高了多波匹配的精度以及匹配的效率。

附图说明

图1为本发明的构建三维Delaunay三角网的流程图；

图2为本发明的构建的Delaunay三角网示意图；

图3为本发明的克里金插值运算的流程图；

图4为本发明的空间离散数据点示意图；

图5为本发明的两道数据进行匹配时划的定时窗窗口示意图；

图6为本发明的原始地震数据；

图7为本发明进行初始匹配后的地震数据；

图8为本发明进行精细匹配后的地震数据。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

基于空间克里金插值的高精度多波匹配方法，包括初始匹配和精细匹配两个步骤，通过基于空间克里金插值的算法完成对原始PS数据进行初始匹配，包括以下步骤：

S1、对最初的gamma数据体进行三角剖分，构建三维Delaunay三角网；

S2、对于构建的Delaunay三角网，通过空间索引的方法对待插值点进行空间定位，定位待插值点所在的四面体，实现邻域点的快速搜索；

S3、对每一个待插值点进行克里金插值运算，计算出每个点的gamma值，根据数据体内所有点的gamma值，采用重采样算法对PS数据进行初始匹配；

所述的精细匹配包括以下步骤：

S4、对初始匹配后的PS数据和原始PP数据，通过基于层位约束的精细匹配算法，通过控制层位点以及划定时窗的方式完成精细匹配。

进一步地，本发明所述的步骤S1中，以人机交互设置的数据点和自动追踪算法计算出来的数据点作为最初的gamma数据体；如图1所示，其实现方法具体包括以下子步骤：

S11、根据给定的点集的区域构建一个初始三角形，使该三角形能够包含所有的数据点，作为初始Delaunay三角网；

S12、从点集中选择一个点进行插入操作，在已有的Delaunay三角网中找到包含该点的三角形，将该点与包含它的三角形的顶点连接形成三个新的三角形；

S13、根据Deluanay三角网的外接空圆特性，利用Lawson提出的LOP算法对当前Delaunay三角网进行优化，使之满足Delaunay特性；

S14、重复步骤S12和S13，当点集中所有点都包含在Delaunay三角网中时结束；

S15、将与初始三角形连接的所有三角形删除，得到如图2所示的三角网。

进一步地，所述的步骤S2包括以下子步骤：

S21、对于当前四面体，计算它的四个顶点的最大和最小Z值，进而计算出Z方向的索引号范围；

S22、按照Z值从小到大的顺序，将Z固定，依次对该四面体进行切割，计算该Z平面与四面体的交点；

S23、若求得交点个数为1，则该点定位完成，若交点个数为3或4，则按照Y值大小将交点排序，交点构成三角形或者四边形；

S24、计算出Y方向的索引号范围，按照Y值从小到大的顺序，将Y值固定，依次对三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点X₁，X₂，则落在X₁，X₂范围内的点即为下一个待处理的四面体；

S25、对于下一个四面体重复步骤S21～S24，直到遍历完所有四面体。

进一步地，如图3所示，所述的步骤S3包括以下子步骤：

S31、进行克里金插值运算：

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)---\left(1\right)$

其中区域化变量为Z(x)，待插值点为x₀，样本点记为x_i(i＝1,2,......,n)，n为样本点的个数，在点x_i处的属性值记为Z(x_i)，则点x₀处的属性值Z(x₀)是由各个样本点属性值加权求和得到的，其中，λ_i为待定权系数；

S32、在无偏估计和区域化变量满足二阶平稳的条件下，为了使得到的估计方差值达到最小，得到加权系数的方程组，即克里金方程组：

$\left\{\begin{array}{c}\underset{j=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{j}C(x_i,x_j)-\mathrm{\&mu;}=\stackrel{\&OverBar;}{C}(x_i,V),(i=1,2...,n)\\ \underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_i=1\end{array}\right.---\left(2\right)$

在实际插值过程中，为了对方程组求解得到权重系数，将方程组转换为矩阵的形式进行表达，之后利用矩阵运算进行求解，便于计算，矩阵形式如式(3)所示：

[K]·[λ]＝[M](3)

在式(3)中，K矩阵的表达式如下：

$\&lsqb;K\&rsqb;=\left[\begin{array}{ccccc}\mathrm{\&gamma;}(v_1,v_1)& \mathrm{\&gamma;}(v_1,v_2)& ...& \mathrm{\&gamma;}(v_1,v_n)& 1\\ \mathrm{\&gamma;}(v_2,v_1)& \mathrm{\&gamma;}(v_2,v_2)& ...& \mathrm{\&gamma;}(v_2,v_n)& 1\\ .& .& & .& .\\ .& .& & .& .\\ .& .& ...& .& .\\ \mathrm{\&gamma;}(v_n,v_1)& \mathrm{\&gamma;}(v_n,v_2)& ...& \mathrm{\&gamma;}(v_n,v_n)& 1\\ 1& 1& ...& 1& 0\end{array}\right]---\left(4\right)$

其中，γ(v_n,v_n)表示第n个已知点与第n个已知点之间的变差函数值，在K矩阵中，前n行n列的值是两两已知点之间的变差函数值；

M矩阵的表达式如下：

$\&lsqb;M\&rsqb;=\left[\begin{array}{c}\mathrm{\&gamma;}(v_1,v)\\ \mathrm{\&gamma;}(v_2,v)\\ .\\ .\\ .\\ \mathrm{\&gamma;}(v_n,v)\\ 1\end{array}\right]---\left(5\right)$

γ(v_n,v)表示第n个已知点与未知点之间的变差函数值，M矩阵的前n行表示待插值点与每个已知点之间的变差函数值；

λ矩阵的表达式如下：

$\&lsqb;\mathrm{\&lambda;}\&rsqb;=\left[\begin{array}{c}{\mathrm{\&lambda;}}_1\\ {\mathrm{\&lambda;}}_2\\ .\\ .\\ .\\ {\mathrm{\&lambda;}}_n\\ \mathrm{\&mu;}\end{array}\right]---\left(6\right)$

λ矩阵即权重系数矩阵，由λ＝[K]⁺[M]计算得到，[K]⁺表示的是K矩阵的逆矩阵；

S33、将λ矩阵带入公式(1)中，计算点x₀处的属性值Z(x₀)，即为该点的gamma值；

S34、重复步骤S31～S33，计算所有待插值点的gamma值；

S35、根据所有待插值点的gamma值，利用重采样算法重新计算该道数据上每个点的数据值，完成初始匹配。

进一步地，所述的步骤S4具体包括以下子步骤：

S41、记录每道PS数据的层位点坐标；

S42、抽取一道PP数据和对应的一道初始匹配后的PS数据构成一个N行2列的数据面，其中N为每道PS数据上采样点的个数；

S43、设定时窗参数：N_w、N_m、N_s，如图5所示，其中，N_w为波形对比窗口的大小，N_m为波形对比窗口的移动量，N_s为波形对比窗口的搜索半径，N_w，N_m，N_s均为整数；以N_w为窗口大小，以N_m为窗口移动量，从上到下，滑动窗口，以窗口中心点作为种子点，计算每个种子点移动量；具体计算方法为：

设i_s＝(k-1)*N_m+N_w/2，i_e＝i_s+N_w，定义集合k为波形对比窗口个数； $\hat{L}=\{-N_s,-N_s+1,...,0,...,N_s-1,N_s\},$

定义长度为(N_w+1)的向量f，f的第i个元素

定义2N_s+1个长度为(N_w+1)的向量g_l，g_l的第i个元素

最优移动量为S(j₁,j₂)，其中j₁和j₂代表第j₁和j₂道；

则S(j₁,j₂)的计算方法为：

$S(j_1,j_2)=\arg \underset{l\&Element;\hat{L}}{max}\frac{f\&CenterDot;g_l}{\sqrt{f\&CenterDot;f}\sqrt{g_l\&CenterDot;g_l}}$

对于第k个波形对比窗口，计算拉平种子点移动量，定义长度为J的向量m用来存储每道PS数据拉平种子点移动量，所述的拉平种子点移动量为：m(j₂)＝S(j₁,j₂)-S(j₁,j₁)；

S44、增加对层位点的控制，对m(j₂)进行处理，此时每个窗口中心点的坐标为i_s＝(k-1)*N_m+N_w/2，如果该坐标是层位点的坐标，则设置m(j₂)＝0，否则不处理；

S45、重复步骤S41～S44，求出整个PS数据面上的种子点的初始移动量，这种方法下，每道PS数据上种子点的位置都一样，然后对每interval道数据同一位置的种子点进行相加求取平均值作为第interval/2道数据种子点的移动量，其中interval为道间隔；求出每第interval/2+i*interval道PS数据上种子点的移动量，其中0<＝i<＝(M-1)/interval，M为该PS数据面上PS数据的道数；

S46、依次对所求出的每interval/2+i*interval道PS数据上同一位置的种子点移动量进行横向平滑处理，然后线性插值出该PS数据面上其他道种子点的最终移动量；

S47、根据步骤S46得到的种子点的最终移动量，采用线性插值算法，分别计算出每道PS数据上每个采样点的移动量；

S48、每个采样点的新的采样坐标为X^new(i,j)＝i+m_i(j)；设置原始采样坐标，具体地，X(i,j)＝i,i＝1,2,...,I；j＝1,2；对于任意一个定值以X的第j列X_j为自变量、以D的第j列D_j为函数值，构建三次样条差值函数，则得到精细匹配后的PS数据。

图6为原始的地震数据，其中，左侧为PP数据，右侧为PS数据。图7为采用本发明的方法进行初始匹配后的地震数据，图8为采用本发明的方法进行精细匹配后的地震数据。从图中可以看出，采用本发明的方法可以达到良好的匹配效果。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.基于空间克里金插值的高精度多波匹配方法，其特征在于，包括初始匹配和精细匹配两个步骤，通过基于空间克里金插值的算法完成对原始PS数据进行初始匹配，包括以下步骤：

S1、对最初的gamma数据体进行三角剖分，构建三维Delaunay三角网；

S2、对于构建的Delaunay三角网，通过空间索引的方法对待插值点进行空间定位；

S3、对每一个待插值点进行克里金插值运算，计算出每个点的gamma值，根据数据体内所有点的gamma值，采用重采样算法对PS数据进行初始匹配；

所述的精细匹配包括以下步骤：

S4、对初始匹配后的PS数据和原始PP数据，通过基于层位约束的精细匹配算法，通过控制层位点以及划定时窗的方式完成精细匹配。

2.根据权利要求1所述的基于空间克里金插值的高精度多波匹配方法，其特征在于，所述的步骤S1中，以人机交互设置的数据点和自动追踪算法计算出来的数据点作为最初的gamma数据体；其实现方法具体包括以下子步骤：

S11、根据给定的点集的区域构建一个初始三角形，使该三角形能够包含所有的数据点，作为初始Delaunay三角网；

S12、从点集中选择一个点进行插入操作，在已有的Delaunay三角网中找到包含该点的三角形，将该点与包含它的三角形的顶点连接形成三个新的三角形；

S13、根据Deluanay三角网的外接空圆特性，利用Lawson提出的LOP算法对当前Delaunay三角网进行优化，使之满足Delaunay特性；

S14、重复步骤S12和S13，当点集中所有点都包含在Delaunay三角网中时结束；

S15、将与初始三角形连接的所有三角形删除。

3.根据权利要求1所述的基于空间克里金插值的高精度多波匹配方法，其特征在于，所述的步骤S2包括以下子步骤：

S21、对于当前四面体，计算它的四个顶点的最大和最小Z值，进而计算出Z方向的索引号范围；

S22、按照Z值从小到大的顺序，将Z固定，依次对该四面体进行切割，计算该Z平面与四面体的交点；

S23、若求得交点个数为1，则该点定位完成，若交点个数为3或4，则按照Y值大小将交点排序，交点构成三角形或者四边形；

S24、计算出Y方向的索引号范围，按照Y值从小到大的顺序，将Y值固定，依次对三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点X₁，X₂，则落在X₁，X₂范围内的点即为下一个待处理的四面体；

S25、对于下一个四面体重复步骤S21～S24，直到遍历完所有四面体。

4.根据权利要求1所述的基于空间克里金插值的高精度多波匹配方法，其特征在于，所述的步骤S3包括以下子步骤：

S31、进行克里金插值运算：

$Z\left(x_0\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{i}Z\left(x_i\right)---\left(1\right)$

其中区域化变量为Z(x)，待插值点为x₀，样本点记为x_i(i＝1,2,......,n)，n为样本点的个数，在点x_i处的属性值记为Z(x_i)，则点x₀处的属性值Z(x₀)是由各个样本点属性值加权求和得到的，其中，λ_i为待定权系数；

S32、在无偏估计和区域化变量满足二阶平稳的条件下，为了使得到的估计方差值达到最小，得到加权系数的方程组，即克里金方程组：

$\left\{\begin{array}{c}\underset{j=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_{j}C(x_i,x_j)-\mathrm{\&mu;}=\stackrel{\&OverBar;}{C}(x_i,V),(i=1,2...,n)\\ \underset{i=1}{\overset{n}{\&Sigma;}}{\mathrm{\&lambda;}}_i=1\end{array}\right.---\left(2\right)$

在实际插值过程中，为了对方程组求解得到权重系数，将方程组转换为矩阵的形式进行表达，之后利用矩阵运算进行求解，便于计算，矩阵形式如式(3)所示：

[K]·[λ]＝[M](3)

在式(3)中，K矩阵的表达式如下：

$\&lsqb;K\&rsqb;=\left[\begin{array}{ccccc}\mathrm{\&gamma;}(v_1,v_1)& \mathrm{\&gamma;}(v_1,v_2)& ...& \mathrm{\&gamma;}(v_1,v_n)& 1\\ \mathrm{\&gamma;}(v_2,v_1)& \mathrm{\&gamma;}(v_2,v_2)& ...& \mathrm{\&gamma;}(v_2,v_n)& 1\\ .& .& & .& .\\ .& .& & .& .\\ .& .& ...& .& .\\ \mathrm{\&gamma;}(v_n,v_1)& \mathrm{\&gamma;}(v_n,v_2)& ...& \mathrm{\&gamma;}(v_n,v_n)& 1\\ 1& 1& ...& 1& 0\end{array}\right]---\left(4\right)$

其中，γ(v_n,v_n)表示第n个已知点与第n个已知点之间的变差函数值，在K矩阵中，前n行n列的值是两两已知点之间的变差函数值；

M矩阵的表达式如下：

$\&lsqb;M\&rsqb;=\left[\begin{array}{c}\mathrm{\&gamma;}(v_1,v)\\ \mathrm{\&gamma;}(v_2,v)\\ .\\ .\\ .\\ \mathrm{\&gamma;}(v_n,v)\\ 1\end{array}\right]---\left(5\right)$

γ(v_n,v)表示第n个已知点与未知点之间的变差函数值，M矩阵的前n行表示待插值点与每个已知点之间的变差函数值；

λ矩阵的表达式如下：

$\&lsqb;\mathrm{\&lambda;}\&rsqb;=\left[\begin{array}{c}{\mathrm{\&lambda;}}_1\\ {\mathrm{\&lambda;}}_2\\ .\\ .\\ .\\ {\mathrm{\&lambda;}}_n\\ \mathrm{\&mu;}\end{array}\right]---6)$

λ矩阵即权重系数矩阵，由λ＝[K]⁺[M]计算得到，[K]⁺表示的是K矩阵的逆矩阵；

S33、将λ矩阵带入公式(1)中，计算点x₀处的属性值Z(x₀)，即为该点的gamma值；

S34、重复步骤S31～S33，计算所有待插值点的gamma值；

S35、根据所有待插值点的gamma值，利用重采样算法重新计算该道数据上每个点的数据值，完成初始匹配。

5.根据权利要求1所述的基于空间克里金插值的高精度多波匹配方法，其特征在于，所述的步骤S4具体包括以下子步骤：

S41、记录每道PS数据的层位点坐标；

S42、抽取一道PP数据和对应的一道初始匹配后的PS数据构成一个N行2列的数据面，其中N为每道PS数据上采样点的个数；

S43、设定时窗参数：N_w、N_m、N_s，其中，N_w为波形对比窗口的大小，N_m为波形对比窗口的移动量，N_s为波形对比窗口的搜索半径，N_w，N_m，N_s均为整数；以N_w为窗口大小，以N_m为窗口移动量，从上到下，滑动窗口，以窗口中心点作为种子点，计算每个种子点移动量；具体计算方法为：

设i_s＝(k-1)*N_m+N_w/2，i_e＝i_s+N_w，定义集合k为波形对比窗口个数；

$\hat{L}=\{-N_s,-N_s+1,...,0,...,N_s-1,N_s\},$

定义长度为(N_w+1)的向量f，f的第i个元素

定义2N_s+1个长度为(N_w+1)的向量g_l，g_l的第i个元素

最优移动量为S(j₁,j₂)，其中j₁和j₂代表第j₁和j₂道；

则S(j₁,j₂)的计算方法为：

$S(j_1,j_2)=\arg \underset{l\&Element;\hat{L}}{max}\frac{f\&CenterDot;g_l}{\sqrt{f\&CenterDot;f}\sqrt{g_l\&CenterDot;g_l}}$

对于第k个波形对比窗口，计算拉平种子点移动量，定义长度为J的向量m用来存储每道PS数据拉平种子点移动量，所述的拉平种子点移动量为：m(j₂)＝S(j₁,j₂)-S(j₁,j₁)；

S44、增加对层位点的控制，对m(j₂)进行处理，此时每个窗口中心点的坐标为i_s＝(k-1)*N_m+N_w/2，如果该坐标是层位点的坐标，则设置m(j₂)＝0，否则不处理；

S45、重复步骤S41～S44，求出整个PS数据面上的种子点的初始移动量，这种方法下，每道PS数据上种子点的位置都一样，然后对每interval道数据同一位置的种子点进行相加求取平均值作为第interval/2道数据种子点的移动量，其中interval为道间隔；求出每第interval/2+i*interval道PS数据上种子点的移动量，其中0<＝i<＝(M-1)/interval，M为该PS数据面上PS数据的道数；

S46、依次对所求出的每interval/2+i*interval道PS数据上同一位置的种子点移动量进行横向平滑处理，然后线性插值出该PS数据面上其他道种子点的最终移动量；

S47、根据步骤S46得到的种子点的最终移动量，采用线性插值算法，分别计算出每道PS数据上每个采样点的移动量；

S48、每个采样点的新的采样坐标为X^new(i,j)＝i+m_i(j)；设置原始采样坐标，具体地，X(i,j)＝i,i＝1,2,...,I；j＝1,2；对于任意一个定值以X的第j列X_j为自变量、以D的第j列D_j为函数值，构建三次样条差值函数则得到精细匹配后的PS数据。