CN105277974A - 一种地层数据插值方法 - Google Patents

Abstract

本发明公开了一种地层数据插值方法，该方法包括：拟合点获取步骤，在未知点区域的四周获取预设数量的地层数据已知的已知点，作为拟合点；变异函数确定步骤，根据拟合点和预设变异模型，确定变异函数；加权系数确定步骤，根据已知点和变异函数，确定加权系数；地层数据插值步骤，利用加权系数对已知点的地层数据加权求和，确定得到未知点的地层数据。本发明提高了最终得到的未知点的地层数据的准确性，有效消除了“野值”现象。

Description

一种地层数据插值方法

技术领域

本发明涉及石油勘探技术领域，具体地说，涉及一种地层数据插值方法。

背景技术

数据插值是现代数据估计与分析中的一个重要手段，它是根据已知点数值估计未知点数值的方法称为数据插值，其目的是把无原始值或者原始值是异常值(也叫野值)的情况下，通过插值估计方法插入或者替换原始值，从而获得符合要求的数据体并进行后续处理和分析。

对层位速度进行编辑分析可以获得更平滑的速度模型，为接下来的偏移得到好的成像剖面提供好的前提条件。但是现有的线性插值方法因为应用数据的局限性不能得到精度较高的估计值。如利用线性或者双线性插值方法只能用到两个方向的值来对未知数据进行估计，从而得到新的层位速度。

然而在使用线性插值方法对层位进行估计的时候，仍然会出现类似“野值”的现象，这对层位速度的插值结果影响很大，严重不利于进一步分析应用。

基于上述情况，亟需一种能够对诸如层位速度等地层数据进行准确插值的方法。

发明内容

为解决上述问题，本发明提供了一种地层数据插值方法，所述方法包括：

拟合点获取步骤，在未知点区域的四周获取预设数量的地层数据已知的已知点，作为拟合点；

变异函数确定步骤，根据所述拟合点和预设变异模型，确定变异函数；

加权系数确定步骤，根据所述已知点和变异函数，确定加权系数；

地层数据插值步骤，利用所述加权系数对已知点的地层数据加权求和，确定得到未知点的地层数据。

根据本发明的一个实施例，所述地层数据包括层位速度。

根据本发明的一个实施例，所述变异函数确定步骤包括：

根据所述拟合点构建实验变异函数；

根据所述实验变异函数和预设变异模型，确定变异函数。

根据本发明的一个实施例，根据如下公式构建实验变异函数：

2γ₀(h)＝E[f(x)-f(x+h)]²

其中，γ₀(h)表示实验变异函数，f(x)表示点x的地层数据，f(x+h)表示沿x轴方向距离为h的点x+h的地层数据。

根据本发明的一个实施例，所述预设变异模型包括球状模型，所述球状模型表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(\frac{3}{2}\frac{h}{a}-\frac{1}{2}\frac{h^3}{a^3})& 0<h\&le;a\\ C_0+c& h>a\end{array}\right.$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，a表示变程。

根据本发明的一个实施例，所述预设变异模型包括高斯模型，所述高斯模型表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(1-e^{-\frac{h^2}{a^2}})& h>0\end{array}\right.$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，表示变程。

根据本发明的一个实施例，所述预设变异模型包括指数模型，所述指数模型表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+C\&CenterDot;(1-e^{-\frac{h}{a}})& h>0\end{array}\right.$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，3a表示变程。

根据本发明的一个实施例，所述加权系数确定步骤包括：

1)根据如下公式计算待构建矩阵K和矩阵M的元素C_ij：

C_ij＝C(|S_i-S_j|)＝C₀+c-γ(|S_i-S_j|)

其中，|S_i-S_j|表示点S_i与点S_j之间的距离，C₀表示块金常数，c表示拱高；

2)根据C_ij构建矩阵K和矩阵M；

3)根据如下公式计算矩阵R：

R＝K^-1·M

4)从矩阵R中提取对应于各个已知点的加权系数。

根据本发明的一个实施例，在所述地层数据插值步骤中，根据如下公式计算所述未知点的地层数据：

F＝λ^T·D

其中，F表示未知点的地层数据矩阵，λ^T表示加权系数矩阵λ的转置，D表示已知点的地层数据矩阵。

本发明提供的地层数据插值方法通过选取未知点区域四周的已知点来作为拟合点以用于确定变异函数。相较于现有的地层数据插值方法，本发明能够有效利用未知点周围的全部信息，从而克服了现有地层数据插值方法的弊端，提高了最终得到的未知点的地层数据的准确性，有效消除了“野值”现象。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要的附图做简单的介绍：

图1是根据本发明一个实施例的地层数据插值方法的流程图；

图2是根据本发明一个实施例的待分析区域层位速度的示意图；

图3是根据本发明一个实施例的拟合点的分布示意图；

图4是根据本发明一个实施例的利用拟合点构建的实验变异函数曲线；

图5是根据本发明一个实施例的插值后的未知点区域的层位速度的示意图；

图6是根据本发明一个实施例的插值后的待分析区域的层位速度示意图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

同时，在以下说明中，出于解释的目的而阐述了许多具体细节，以提供对本发明实施例的彻底理解。然而，对本领域的技术人员来说显而易见的是，本发明可以不用这里的具体细节或者所描述的特定方式来实施。

另外，在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行，并且，虽然在流程图中示出了逻辑顺序，但是在某些情况下，可以以不同于此处的顺序执行所示出或描述的步骤。

众所周知，许多实际问题都用函数y＝f(x)来表示某种内在规律的数量关系，其中相当一部分函数是通过试验或观测得到的。虽然f(x)在[a,b]上是存在的，有的还是连续的。但在实际应用中，例如对于层位速度，只能给出[a,b]上的一系列点x_i的函数值y_i，即：

y_i＝f(x_i)i＝0,1,......,n(1)

但是为了研究函数以及层位速度的变化规律，往往需要求出不在以上点的函数值。因此就需要一个既能反映函数f(x)的特性，又便于计算的简单函数P(x)。用P(x)近似f(x)，通常选一类简单的函数作为P(x)，P(x)具有如下特性：

P(x_i)＝f(x_i)i＝1,2,......,n(2)

这个函数P(x)就是希望得到的插值估计函数，通过函数P(x)便可以计算出未知点的层位速度，这也就是层位速度的插值。

现有的层位速度插值方法主要为线性插值法，而线性插值法一般指的是双线性插值(也叫双线性内插)。这种插值方法是有两个变量的插值函数的拓展，其核心思想是在两个方向上分别进行一次线性插值。

如期望得到函数f在点P＝(x,y)的值，假设已知函数f在Q₁₁＝(x₁,y₁)、Q₁₂＝(x₁,y₂)、Q₂₁＝(x₂,y₁)、Q₂₂＝(x₂,y₂)四个点的值。那么首先在x方向进行线性插值，然后在y方向进行线性插值，最终便得到双线性插值的结果。

针对于线性插值，其误差R_T可以表示为：

R_T＝f(x)-p(x)(3)

其中，p(x)可以根据如下公式计算得到：

$p\left(x\right)=f\left(x_2\right)+\frac{f\left(x_1\right)-f\left(x_2\right)}{x_1-x_2}(x-x_2)---\left(4\right)$

根据罗尔定理，可以知道，如果f有二阶连续导数，那么误差R_T的范围可以根据如下公式计算得到：

$\left|R_T\right|\&le;\frac{{(x_2-x_1)}^2}{8}\underset{x_1\&le;x{\&le;x}_2}{\max }\left|f^{\&prime;\&prime;}\left(x\right)\right|---\left(5\right)$

从公式(5)中可以看出，函数上两点之间的近似随着所近似的函数的二阶导数的增大而逐渐变差，即函数的曲率越大，简单线性插值近似的误差也就越大。所以在进行层位速度插值的过程中，需要得到近似误差较小的结果，显然线性插值算法无法满足

本发明通过对线性插值方法的分析，发现由于线性插值方法只能用到估计点附近一个或几个方向上的已知信息，而不能利用周围的全部信息。这种方向性的限制是线性插值算法的主要弊端，同时也造成了估计值的误差较大，导致出现“野值”现象。

针对线性插值方法的弊端，本发明提供了一种能够有效克服该弊端的层位速度插值方法，图1示出了本实施例中该方法的流程图。同时，为了更清楚的阐述本方法的目的、原理以及优点，以下结合图2所示出的待分析区域的示意图，以地层数据中的层位速度为例来进行说明。需要说明的是，以层位速度为例仅仅是为了更加清楚地阐述本发明的目的、原理以及优点，其并不是作为对本发明的限定，在本发明的其他实施例中，地层数据还可以为其他合理的数据。

如图2所示，本实施例所提供的待分析区域中，在线号区间为[1951，2000]、CDP号区间为[1101，1150]的区域中所包含的2500个点是层位速度未知的未知点，而待分析区域中的其他点则为层位速度已知的已知点，而本方法的目的即是通过已知点来确定未知点的层位速度。

如图1所示，本方法在步骤S101中从未知点区域的四周获取预设数量的已知点，并将获取到的已知点作为拟合点。其中，图3示出了本实施例中获取的拟合点的分布示意图。从图3中可以看出，本实施例从未知点的四周共获取了1700个已知点来作为拟合点。

在步骤S102中，根据步骤S101获取的拟合点构建实验变异函数。设点x的层位速度为f(x)，沿x轴方向距离为h的点x+h的层位速度为f(x+h)，则拟合点沿x轴方向可以被相同的距离h分割成许多对点，从而得到相应的一组差值f(x)-f(x+h)，而此组差值的平方期望即为实验变异函数γ₀(h)。其中实验变异函数γ₀(h)可以根据如下公式计算得到：

2γ₀(h)＝E[f(x)-f(x+h)]²(6)

公式(6)还可以表示为：

${\mathrm{\&gamma;}}_0\left(h\right)=\frac{\underset{i=1}{\overset{N\left(h\right)}{\mathrm{\&Sigma;}}}{[f\left(x_i\right)-f(x_i+h)]}^2}{2\&CenterDot;N\left(h\right)}---\left(7\right)$

其中，N(h)为距离等于h的点对的总数，x_i和x_i+h构成第i对点对。以距离h为横坐标，γ₀(h)为纵坐标作图，即可绘制出实验变异函数曲线。本实施例中，利用拟合点所得到的实验变异函数曲线如图4所示。

地质统计学中的Kriging方法可以用于空间数据的最优线性无偏内插估计，该方法在地质、土壤、气象、生态、地球化学等众多领域都有着广泛的应用。在进行Kriging空间插值时，关键是对区域化变量的空间相关性进行分析，即需要对选用的变异函数模型进行比较取舍，以及对模型参数进行优化估计。通常是采用一定的理论模型对实验变异函数进行拟合，从而求得最佳的参数。

所以再次如图1所示，本实施例在步骤S103中利用预设变异模型对步骤S102得到的实验变异函数进行拟合，得到最佳的模型参数，进而得到变异模型。其中，本实施例中，预设变异模型选用球状模型，该模型可以用如下公式表示：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(\frac{3}{2}\frac{h}{a}-\frac{1}{2}\frac{h^3}{a^3})& 0<h\&le;a\\ C_0+c& h>a\end{array}\right.---\left(8\right)$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，a表示变程。

为了求解最佳参数，即块金常数C₀、拱高c、变程a的最佳值，就需要采用一定的拟合标准。而考虑到滞后距离分组后，各组h_j的点对数N(h_j)(j＝1,2,3,...,m，其中j为组号，即实验变异函数曲线上的点号，m为组数)不一，对建立拟合标准的影响程度也不同。所以各组点对数对拟合标准的影响程度可以采用权重系数来进行衡量，其中，本实施例中，权重系数可以采用如下公式计算得到：

其拟合标准为：

$C_{\mathrm{\&omega;}}=\min \left\{\underset{j=1}{\overset{m}{\mathrm{\&Sigma;}}}\mathrm{\&omega;}\left(i\right){[r\left(h_j\right)-r_0\left(h_j\right)]}^2\right\}---\left(10\right)$

其中，r(h_j)表示变异曲线上点h_j的取值，C_ω即是模型优化问题的目标函数。由于球状模型为非连续的的非连续可导，求解参数时需要采用直接搜索法，具体地，可以采用遗传算法与模式搜索法相结合的方法来求解模型参数。

需要说明的是，在本发明的其他实施例中，预设变异模型还可以采用其他合理的变异模型，例如高斯模型或指数模型，本发明不限于此。其中，高斯模型可以用如下公式表示：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(1-e^{-\frac{h^2}{a^2}})& h>0\end{array}\right.---\left(11\right)$

与球状模型不同的是，在高斯模型下，表示变程。

而指数模型可以表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(1-e^{-\frac{h}{a}})& h>0\end{array}\right.---\left(12\right)$

与球状模型不同的是，在高斯模型下，3a表示变程。

此外，在本发明的其他实施例中，也可以采用其他合理方式来求解各个变异模型的最佳参数，本发明不限于此。

得到变异函数后，如图1所示，本实施例在步骤S104中根据已知点和变异函数确定加权系数。具体地，首先根据如下公式计算各个点之间的C_ij：

C_ij＝C(|S_i-S_j|)＝C₀+c-γ(|S_i-S_j|)(13)

其中，|S_i-S_j|表示点S_i与点S_j之间的距离，C₀表示块金常数，c表示拱高。

随后根据C_ij构建矩阵K和矩阵M。本实施例中，假设已知点的数量为5000个，所以构建得到的矩阵K可以表示为：

$K=\left[\begin{array}{cc}\left[C_{\mathrm{kl}}\right]& 1\\ 1& 0\end{array}\right],k=\mathrm{1,2},...,5000,l=\mathrm{1,2},...,5000---\left(14\right)$

其中，[C_kl]表示由各个已知点之间的C_kl构成的矩阵，本实施例中它的维度为5000×5000，那么构建得到的矩阵K的维度即为5001×5001。

未知点的数量为2500个，那么构建得到的矩阵M可以表示为：

$M=\left[\begin{array}{c}\left[C_{\mathrm{mn}}\right]\\ 1\end{array}\right],m=\mathrm{1,2},...,5000,n=\mathrm{1,2},...,2500---\left(15\right)$

其中，[C_mn]表示由各个未知点与各个已知点之间的C_mn构成的矩阵，本实施例中，它的维度为5000×2500，那么构建得到的矩阵M的维度则为5001×2500。

那么矩阵R则可以根据如下公式计算得到：

R＝K^-1·M(16)

其中， $R=\left[\begin{array}{c}\mathrm{\&lambda;}\\ -\mathrm{\&mu;}\end{array}\right].$ 那么根据公式(16)计算得到的矩阵R即可以提取得到加权系数矩阵λ，从而得到各个已知点对未知点的加权系数。本实施例中，加权系数矩阵λ的维度为5000×2500

最后在步骤S105中，利用加权系数对已知点的层位速度进行加权求和，得到未知点的层位速度。其可以用如下公式表示：

F＝λ^T·D(17)

其中，F表示未知点的层位速度矩阵，λ^T表示加权系数矩阵λ的逆矩阵，D表示已知点的层位速度矩阵，其维度为5000×1。根据公式(17)即可得到维度为2500×1的矩阵，其中各个元素即为各个未知点的层位速度。

至此完成层位速度的插值。

图5示出了利用上述方法所得到的未知点区域的层位速度的示意图，而图6则示出了整个待分析区域的层位速度示意图。从图5和图6中均可以看出，本方法所得到的层位速度插值结果中不存在“野值”现象。

从以上描述中可以看出，本发明提供的层位速度插值方法通过选取未知点区域四周的已知点来作为拟合点以用于确定变异函数。相较于现有的层位速度插值方法，本发明能够有效利用未知点周围的全部信息，从而克服了现有层位速度插值方法的弊端，提高了最终得到的未知点的层位速度的准确性，有效消除了“野值”现象。

应该理解的是，本发明所公开的实施例不限于这里所公开的特定结构、处理步骤或材料，而应当延伸到相关领域的普通技术人员所理解的这些特征的等同替代。还应当理解的是，在此使用的术语仅用于描述特定实施例的目的，而并不意味着限制。

说明书中提到的“一个实施例”或“实施例”意指结合实施例描述的特定特征、结构或特性包括在本发明的至少一个实施例中。因此，说明书通篇各个地方出现的短语“一个实施例”或“实施例”并不一定均指同一个实施例。

虽然上述示例用于说明本发明在一个或多个应用中的原理，但对于本领域的技术人员来说，在不背离本发明的原理和思想的情况下，明显可以在形式上、用法及实施的细节上作各种修改而不用付出创造性劳动。因此，本发明由所附的权利要求书来限定。

Claims (9)

1.一种地层数据插值方法，其特征在于，所述方法包括：

拟合点获取步骤，在未知点区域的四周获取预设数量的地层数据已知的已知点，作为拟合点；

变异函数确定步骤，根据所述拟合点和预设变异模型，确定变异函数；

加权系数确定步骤，根据所述已知点和变异函数，确定加权系数；

地层数据插值步骤，利用所述加权系数对已知点的地层数据加权求和，确定得到未知点的地层数据。

2.如权利要求1所述的方法，其特征在于，所述地层数据包括层位速度。

3.如权利要求1或2所述的方法，其特征在于，所述变异函数确定步骤包括：

根据所述拟合点构建实验变异函数；

根据所述实验变异函数和预设变异模型，确定变异函数。

4.如权利要求3所述的方法，其特征在于，根据如下公式构建实验变异函数：

2γ₀(h)＝E[f(x)-f(x+h)]²

其中，γ₀(h)表示实验变异函数，f(x)表示点x的地层数据，f(x+h)表示沿x轴方向距离为h的点x+h的地层数据。

5.如权利要求4所述的方法，其特征在于，所述预设变异模型包括球状模型，所述球状模型表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(\frac{3}{2}\frac{h}{a}-\frac{1}{2}\frac{h^3}{a^3})& 0<h\&le;a\\ C_0+c& h>a\end{array}\right.$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，a表示变程。

6.如权利要求4所述的方法，其特征在于，所述预设变异模型包括高斯模型，所述高斯模型表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+c\&CenterDot;(1-e^{-\frac{h^2}{a^2}})& h>0\end{array}\right.$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，表示变程。

7.如权利要求4所述的方法，其特征在于，所述预设变异模型包括指数模型，所述指数模型表示为：

$\mathrm{\&gamma;}\left(h\right)=\left\{\begin{array}{cc}0& h=0\\ C_0+C\&CenterDot;(1-e^{-\frac{h}{a}})& h>0\end{array}\right.$

其中，h表示距离，γ(h)表示变异函数，C₀表示块金常数，c表示拱高，3a表示变程。

8.如权利要求1～7中任一项所述的方法，其特征在于，所述加权系数确定步骤包括：

1)根据如下公式计算待构建矩阵K和矩阵M的元素C_ij：

C_ij＝C(|S_i-S_j|)＝C₀+c-γ(|S_i-S_j|)

其中，|S_i-S_j|表示点S_i与点S_j之间的距离，C₀表示块金常数，c表示拱高；

2)根据C_ij构建矩阵K和矩阵M；

3)根据如下公式计算矩阵R：

R＝K^-1·M

4)从矩阵R中提取对应于各个已知点的加权系数。

9.如权利要求8所述的方法，其特征在于，在所述地层数据插值步骤中，根据如下公式计算所述未知点的地层数据：

F＝λ^T·D

其中，F表示未知点的地层数据矩阵，λ^T表示加权系数矩阵λ的转置，D表示已知点的地层数据矩阵。