基于空间克里金插值的高精度多波匹配方法
【技术领域】
[0001] 本发明属于地质解释领域中多波匹配问题，特别涉及一种基于空间克里金插值的 高精度多波匹配方法。
【背景技术】
[0002] 多波地震勘探是进行岩性油气藏和隐蔽油气藏勘探的一种非常有潜力的手段，但 是，由于诸多原因，多波多分量理论研究和油气田实际勘探地质需求的结合问题、复杂条件 下的转换波地震资料处理问题和多波综合解释、全波属性的地质应用等问题一直没有取得 显著进展，并且已经成为制约多波地震勘探技术进一步发展的"瓶颈"。而解决这些问题 的基础是做好多波多分量资料处理，提供高质量的各向同性和各向异性处理成果。其中， 多波传播机理的基础研究、多波资料中的纵横波匹配方法研究是目前多波地震资料后续处 理的重点和难点，是多波精确成像和叠前纵横波联合反演及岩性识别、储层预测和含气性 识别的重要基础，是体现多波多分量地震勘探技术实际勘探开发应用价值的关键。因此， 基于多波传播机理，研究纵横波高精度的匹配新方法，有利于充分利用多波多分量地震资 料、准确认识多波地质响应特征，突出多波多分量地震资料解决地质问题的能力，具有重 大的意义。
[0003] 目前多波匹配有基于反射特征的匹配方法和基于多波层位的匹配方法，前者通过 PP波和PS波地震数据的波形和波组特征进行对比生成γ。值，然后基于该γ。值实现两者 的时间域匹配。后者首先分别基于ΡΡ波和PS波地震数据追踪解释出对应的层位，然后通 过标定对应的层位产生时移体，最后将时移体应用于PS地震数据，实现PS与ΡΡ地震数据 的匹配。
[0004] 随着勘探目标要求的提高，多波匹配技术的研究越来越受到人们的重视，纵横波 匹配技术已经成为地球物理学的研究热点。JamesE.G(1996)详细介绍了纵横波速度比 (γ。）的求取方法，并用最大相关法求取γ<：、平均γ<：、层间Yc等，使用VSP资料验证从PP 波和PS波剖面中求取γ。，且γ。可以用短波长振幅反演。Wai-KinChan等（1997)利用常 数值多次试算法，在时间对数域实现纵、横波匹配。由于浅层和深层的γc值不同，因此 该匹配方法只能适应于特定的目的层。RichardV.D和JamesG等（2001)利用最大相似 性原理，扫描PP波和PS波的γ。谱，拾取平均γ。值，将二者在时间域进行匹配。VanDok 和Kristiansen(2003)采用人工手动拾取强同相轴和其他明显的地质特征，实现了ΡΡ波和 PS波在时间域的粗略匹配。MichaelV.D等（2003)从ΡΡ波和PS波中求得纵横波速度比 和泊松比，基于建立的横波模型，在深度域实现纵横波匹配。就纵横波匹配技术软件现状而 言，Transform软件的匹配功能相对完善，提供的基于反射特征的匹配方法和基于多波层位 的匹配方法都能实现PP数据和PS数据的匹配。纵观国内外纵横波匹配技术的发展，几乎 所有的方法均是直接以PP波和PS波同一地层反射同相轴相似性最大的原则求取纵横波的 速度比，直接对PP波和PS波在时间域进行匹配，其匹配程度是比较粗略的。
【发明内容】
[0005] 本发明的目的在于克服现有技术，提供一种通过空间克里金插值算法对初始的 ga_a数据体进行插值，采用重采样算法对PS数据进行初始匹配；然后通过控制层位点以 及划定时窗的方式完成PS数据和PP数据的精细匹配，完成最终多波匹配，提高了多波匹配 的精度以及匹配的效率的基于空间克里金插值的高精度多波匹配方法。
[0006] 本发明的目的是通过以下技术方案来实现的：基于空间克里金插值的高精度多波 匹配方法，包括初始匹配和精细匹配两个步骤，通过基于空间克里金插值的算法完成对原 始PS数据进行初始匹配，包括以下步骤：
[0007]S1、对最初的gamma数据体进行三角剖分，构建三维Delaunay三角网；
[0008] S2、对于构建的Delaunay三角网，通过空间索引的方法对待插值点进行空间定 位；
[0009] S3、对每一个待插值点进行克里金插值运算，计算出每个点的8曰_3值，根据数据 体内所有点的ga_a值，采用重采样算法对PS数据进行初始匹配；
[0010] 所述的精细匹配包括以下步骤：
[0011] S4、对初始匹配后的PS数据和原始PP数据，通过基于层位约束的精细匹配算法， 通过控制层位点以及划定时窗的方式完成精细匹配。
[0012] 进一步地，所述的步骤S1中，以人机交互设置的数据点和自动追踪算法计算出来 的数据点作为最初的gamma数据体；其实现方法具体包括以下子步骤：
[0013] SI1、根据给定的点集的区域构建一个初始三角形，使该三角形能够包含所有的数 据点，作为初始Delaunay三角网；
[0014]S12、从点集中选择一个点进行插入操作，在已有的Delaunay三角网中找到包含 该点的三角形，将该点与包含它的三角形的顶点连接形成三个新的三角形；
[0015] S13、根据Deluanay三角网的外接空圆特性，利用Lawson提出的L0P算法对当前 Delaunay三角网进行优化，使之满足Delaunay特性；
[0016] S14、重复步骤S12和S13,当点集中所有点都包含在Delaunay三角网中时结束；
[0017] S15、将与初始三角形连接的所有三角形删除。
[0018] 进一步地，所述的步骤S2包括以下子步骤：
[0019] S21、对于当前四面体，计算它的四个顶点的最大和最小Z值，进而计算出Z方向的 索引号范围；
[0020] S22、按照Z值从小到大的顺序，将Z固定，依次对该四面体进行切割，计算该Z平 面与四面体的交点；
[0021] S23、若求得交点个数为1，则该点定位完成，若交点个数为3或4,则按照Y值大小 将交点排序，交点构成三角形或者四边形；
[0022] S24、计算出Y方向的索引号范围，按照Y值从小到大的顺序，将Y值固定，依次对 三角形或四边形进行切割，计算该平行于X轴的直线与三角形或四边形的交点Xi，X2,则落 在Xi，X2范围内的点即为下一个待处理的四面体；
[0023]S25、对于下一个四面体重复步骤S21~S24,直到遍历完所有四面体。
[0024] 进一步地，所述的步骤S3包括以下子步骤：
[0025]S31、进行克里金插值运算：
[0027] 其中区域化变量为Z(x)，待插值点为x。，样本点记为Xi(i= 1，2,......，η)，η为 样本点的个数，在点^处的属性值记为Ζ(χJ，则点χ。处的属性值Ζ(χ。）是由各个样本点属 性值加权求和得到的，其中，λ 待定权系数；
[0028]S32、在无偏估计和区域化变量满足二阶平稳的条件下，为了使得到的估计方差值 达到最小，得到加权系数的方程组，即克里金方程组：
[0030] 在实际插值过程中，为了对方程组求解得到权重系数，将方程组转换为矩阵的形 式进行表达，之后利用矩阵运算进行求解，便于计算，矩阵形式如式（3)所示：
[0031] [Κ]·[λ] = [Μ] (3)
[0032] 在式（3)中，Κ矩阵的表达式如下：
[0034] 其中，γ(νη，νη)表示第η个已知点与第η个已知点之间的变差函数值，在Κ矩阵 中，前η行η列的值是两两已知点之间的变差函数值；
[0035] Μ矩阵的表达式如下：
[0037] γ(νη，ν)表示第η个已知点与未知点之间的变差函数值，Μ矩阵的前η行表示待 插值点与每个已知点之间的变差函数值；
[0038]λ矩阵的表达式如下：
[0040]λ矩阵即权重系数矩阵，由λ= [Κ] + [Μ]计算得到，[Κ]+表示的是Κ矩阵的逆矩 阵；
[0041]S33、将λ矩阵带入公式（1)中，计算点χ。处的属性值Z(x。)，即为该点的gamma 值；
[0042]S34、重复步骤S31~S33,计算所有待插值点的gamma值；
[0043]S35、根据所有待插值点的gamma值，利用重采样算法重新计算该道数据上每个点 的数据值，完成初始匹配。
[0044] 进一步地，所述的步骤S4具体包括以下子步骤：
[0045]S41、记录每道PS数据的层位点坐标；
[0046]S42、抽取一道PP数据和对应的一道初始匹配后的PS数据构成一个N行2列的数 据面，其中N为每道PS数据上采样点的个数；
[0047]S43、设定时窗参数：凡、乂、队，其中，Nw为波形对比窗口的大小，N"为波形对比窗口 的移动量，队为波形对比窗口的搜索半径，Nw，N"，Ns均为整数；以Nw为窗口大小，以N"为窗 口移动量，从上到下，滑动窗口，以窗口中心点作为种子点，计算每个种子点移动量；具体计 算方法为：
[0048]设is= (k-1) *Nm+Nw/2,ie=iS+Nw，定义集合f= …,义丨丨，k为波形对比窗口 个数；Z={-Λζ，-Μ+1，··.，0,…，-1，Λζ}，
[0049] 定义长度为（Nw+1)的向量f，f的第i个元素/(/) =ΖΧ(十
[0050] 定义2NS+1个长度为（Nw+1)的向量gl，gl的第i个元素 g,(i) =D(i^+/+z-1, /2), /e?,Ι
[0051] 最优移动量为Sa,j2)，其中jJPj2代表第jJPj2道；
[0052] 则S(j\，j2)的计算方法为：
[0054] 对于第k个波形对比窗口，计算拉平种子点移动量，定义长度为J的向量 m用来存储每道PS数据拉平种子点移动量，所述的拉平种子点移动量为：m(j2)= sa，j2)-sa，j);
[0055]S44、增加对层位点