CN102830430B - 一种层位速度建模方法 - Google Patents

Abstract

本发明公开了一种层位速度建模方法，具体包括步骤：导入样品点的坐标及属性值；计算实验变差函数；选取理论变差函数模型并进行拟合，得到三个轴上的变差函数；构造权重矩阵，将三个方向上的变差函数统一套合得到各向同性变差函数；求解克里金方程组，进行速度插值生成速度模型。本发明的层位速度建模方法采用Nelder-Mead算法进行拟合，可以处理大数量级的数据；该方案考虑了不同尺度上的变异，采用混合模型代替原有的单一模型，用几种单一的变差函数混合的方式来描述变异性，大大提高拟合精度；统一套合方法结构简单，可以解决三维空间中不同方向变差函数类型不同时的套合，适用于多种复杂变差函数组合，具有普适性。

Description

一种层位速度建模方法

技术领域

本发明属于地质勘探技术领域，具体涉及其中的层位速度建模方法。

背景技术

在地质勘探过程中，地震速度是进行时深转换、岩性预测、构造建模等构造解释工作的不可缺少参数，精确的建立地震速度模型是地质勘探过程中的重要工作之一。三维地质建模以反映地质体及其特征的数据为基础，利用计算机技术模拟地质体的表面形态和内部属性特征，并以图形图像的方式再现地质体的三维构造，最后进行地质体的空间形态和内部属性分析。三维地质建模按技术主要分为两类：一类是确定性建模，一类是随机建模。随机建模是指以已知的信息为基础，以随机函数为理论，应用随机模拟方法，产生可选的、等可能的储层模型的方法。确定性建模是指从具有确定性资料的控制点出发，推测出点间确定的、唯一的储层参数，其建模的方法主要有储层地震学方法、储层沉积学方法及地质统计学克里金插值方法。

克里金方法在国外的固体矿产资源估算中有广泛的应用，它是一种综合考虑变量结构性和随机性后，对未采样点属性进行最优无偏估计的估值方法。利用克里金方法进行插值，必须完成以下两个步骤：（1）根据数据的空间相关性，计算变差函数并进行拟合，得到插值模型；（2）选择一种克里金方法（例如普通克里金，指示克里金），进行插值。

在应用克里金插值进行层位速度建模的过程中，预测未知点速度值的方法已经很成熟，而获取空间变异性（即求取变差函数）的方法还存在很多问题，获取空间变异性有三个重要的内容，选取理论变差函数模型、拟合算法、套合成统一模型，现有的研究方法都存在一些不足。在三维克里金中，变差函数可能会出现不同尺度上变异性，这种情况下的变异性由混合模型来描述，现有的拟合理论变差函数的方法为最小二乘法，用最小二乘法拟合这种混合模型的变差函数会很复杂，此外，变差函数也会出现不同方向上的变异性，即在不同方向上，变差函数的模型都可能会不相同，这种情况下，变差函数的组合情况有很多种，目前还没有一个较好的方法可以解决所有组合情况下的三维变差函数套合。

发明内容

本发明的目的是为了解决在层位速度建模中现有的克里金方法存在的上述问题，提出了一种层位速度建模方法。

本发明的技术方案是：一种层位速度建模方法，包括如下步骤：

步骤1.导入样品点的坐标及属性值；

步骤2.分别计算X轴、Y轴、Z轴上的实验变差函数γ*_x(h)、γ*_y(h)、γ*_z(h)；

步骤3.利用步骤2得到的三个方向上的实验变差函数，选取理论变差函数模型并进行拟合，得到X轴、Y轴、Z轴上的变差函数γ_x(h)、γ_y(h)、γ_z(h)；

步骤4.利用三个轴向的比例函数构造权重矩阵，将步骤3得到的三个方向上的变差函数统一套合得到最终的各向同性变差函数γ(h)；

步骤5.利用步骤4得到的变差函数γ(h)，求解克里金方程组，对未知点进行速度插值，最终生成速度模型。

进一步的，步骤3所述的拟合具体采用Nelder-Mead算法。

本发明的有益效果：本发明的层位速度建模方法采用Nelder-Mead算法进行拟合，可以处理大数量级的数据；该方案考虑了不同尺度上的变异，采用混合模型代替原有的单一模型，用几种单一的变差函数混合的方式来描述变异性，大大提高拟合精度；统一套合方法结构简单，可以解决三维空间中不同方向变差函数类型不同时的套合，适用于多种复杂变差函数组合，具有普适性，并在工程实现中应用，具有可行性。本发明的方法一方面可以改善变差函数拟合效果，另一方面可以处理三维空间中不同类型的各向异性套合，能处理几何各向异性、带状各向异性，以及多种不同类型变差函数模型的组合。

附图说明

图1为本发明的层位速度建模方法流程示意图。

图2为计算实验变差函数所需参数的描述示意图。

图3为一幅实验变差函数示意图。

图4为高斯模型变差函数示意图。

图5为指数模型变差函数示意图。

图6为高斯-指数混合模型变差函数示意图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步的说明。

本发明的层位速度建模方法流程示意图如图1所示，具体包括如下步骤：

步骤1.导入样品点的坐标及层位速度值。

三维层位建模的数据源是一系列已知层位速度值的空间坐标点，这些点一般来源于打井得到的采样点。建模首要的任务便是收集采样点数据，由于克里金方法要求采样点的层位速度值必须服从正态分布，因此需要对收集整理后的实验数据进行统计分析，分析主要包括计算算术平均数，方差，相关系数等，统计分析用以判断数据是否满足正态分布，并可以确定哪些点是极端值，加以处理，对采样点数据进行处理后，便得到样品点。

步骤2.分别计算X轴、Y轴、Z轴上的实验变差函数γ_x(h)、γ_y(h)、γ_z(h)，具体过程如下：

在三维空间中，实验变差函数的计算公式如下：

${\mathrm{\γ}}^{*}\left(\stackrel{\→}{h}\right)=\frac{1}{2N}\underset{i,j,k=1}{\overset{N}{\mathrm{\Σ}}}{[Z(x_i,y_j,z_k)-Z({x_i}^{,},{y_j}^{,},{z_k}^{,})]}^2---\left(1\right)$

其中， 表示距离为的实验变差函数值，N为空间网格上相距为的点对数；Z(x_i,y_j,z_k)表示在点(x_i,y_j,z_k)上的Z值；Z(x_i’,y_j’,z_k’)表示在点(x_i’y_j’,z_k’)处的Z值。

在计算实验变差函数前，需要知道以下参数：滞后距个数Num：由用户给定，一般设为20；最大搜索距离max_d：由用户给定，最大值不能超过列线长度的一半；基本滞后距h：最大搜索距离除以滞后距个数；距离容差tolerance_h：由用户根据实际情况给定；角度容差tolerance_angle：由用户根据实际情况给定。上述参数的描述如图2所示。

对于某一特定方向上的实验变差函数，计算步骤如下：

a)计算所有已知点对之间的距离、角度以及属性值之差的平方。

b)确定滞后距：lag＝h*K（K＝0,1,2,3…）。

c)每个滞后距处，寻找落在其tolerance_h和tolerance_angle范围内的点对，根据式（1）计算实验变差函数值，图3给出了一幅实验变差函数图。

通过上述过程，可以得到X轴、Y轴、Z轴上的实验变差函数γ*_x(h)、γ*_y(h)、γ*_z(h)。

步骤3.利用步骤2得到的三个方向上的实验变差函数，选取理论变差函数模型并进行拟合，得到X轴、Y轴、Z轴上的变差函数γ_x(h)、γ_y(h)、γ_z(h)。

现有的理论变差函数模型主要有线性模型，球状模型，指数模型，高斯模型，五球形模型等等，但是在实际中，某个方向上的变差函数会出现不同尺度上的变异，这种情况下，选取单一的模型可能无法很好的拟合实验变差函数，本发明的方案考虑了混合模型来拟合实验变差函数。

下面给出高斯模型，指数模型，高斯-指数混合模型的表达式。

高斯模型：

$\mathrm{\γ}\left(h\right)=C_0+C(1-e^{-\frac{{\left(3h\right)}^2}{a^2}})---\left(2\right)$

指数模型：

$\mathrm{\γ}\left(h\right)=C_0+C(1-e^{-\frac{3h}{a}})---\left(3\right)$

高斯-指数混合模型：

$\mathrm{\γ}\left(h\right)=C_0+C_1(1-e^{-\frac{3h}{a_1}})+C_2(1-e^{-\frac{{\left(3h\right)}^2}{{a_2}^2}})---\left(4\right)$

其中，C₀表示块金常数，它反映了层位速度的连续性很差，即使在很短的距离内，速度的变化也可能很大；a、a₁、a₂表示变程，变差函数通过变程来反映速度的影响范围，变差函数随着距离的增大而增大，但当距离大于变程值后，变差函数不再增大，稳定在γ(∞)附近，γ(∞)称为基台值，C₀+C表示基台值，C、C₁、C₂表示拱高。

在本发明的方法中，考虑到不同尺度上的变异，采用混合模型代替原有的单一模型，用几种单一的变差函数混合的方式来描述变异性，因而可以提高拟合精度，在此具体选择高斯-指数混合模型。

在拟合实验变差函数时，这里采用Nelder-Mead（单纯性）算法。以未知估计参数的初始值为基础，Nelder-Mead算法根据迭代的控制变量，在每次迭代中修正待估参数，计算出目标函数的极小值并进行比较，进入下一次迭代，算法的收敛是由迭代次数或者待估参数精度控制。

对于式（4）所示的混合模型，计算变差函数残差目标函数表达式如下：

$f(a_1,a_2,c_0,c_1,c_2)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{w}_i{[\mathrm{\γ}\left(h_i\right)-{\mathrm{\γ}}^{*}\left(h_i\right)]}^2$

$=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{w}_i{[c_0+c_1(1-e^{-\frac{3h}{a_1}})+c_2(1-e^{-\frac{{\left(3h\right)}^2}{{a_2}^2}})-{\mathrm{\γ}}^{*}\left(h_i\right)]}^2---\left(5\right)$

其中，w_i表示权重值，Nelder-Mead迭代算法的控制参量取值如表1所示。

表1

标量参数	性质	取值
			ρ	反射系数	1
χ	延伸系数	2
			γ	缩小系数	1/2
σ	缩减系数	1/2

这里先定义一个变量x，x代表混合模型的待估参数：x=(a₁,a₂,c₀,c₁,c₂)。Nelder-Mead算法迭代的控制变量由表1给出，下面介绍迭代过程：

首先给出初始n+1（n指参数个数，参数为a₁,a₂,c₀,c₁,c₂，n=5）个x值，即给出n+1套a₁,a₂,c₀,c₁,c₂的值，给定的方法可以由用户预先设定。

按照下列步骤进行迭代：

第一步：排序，按照f(x₁)≤f(x₂)≤…≤f(x_n+1)进行排序，这里的f(x)由式（5）求得。

第二步：反射，计算反射点x_r： $x_r=\stackrel{\&OverBar;}{x}+\mathrm{\&rho;}(\stackrel{\&OverBar;}{x}-x_{n+1})=(1+\mathrm{\&rho;})\stackrel{\&OverBar;}{x}-{\mathrm{\&rho;x}}_{n+1},$ 其中计算f_r＝f(x_r)，若f₁≤f_r<f_n，用反射点x_r代替x_n+1，回到第一步；若f_r<f₁，跳到第三步；若f_r≥f_n，跳到第四步。

第三步：延伸，计算点x_e： $x_e=\stackrel{\&OverBar;}{x}+\mathrm{\&chi;}(x_r-\stackrel{\&OverBar;}{x})=(1+\mathrm{\&rho;\&chi;})\stackrel{\&OverBar;}{x}-\mathrm{\&rho;\&chi;}{x}_{n+1}$

若f_e<f_r，用点x_e代替点x_n+1，回到第一步；反之，用反射点x_r代替x_n+1，回到第一步。

第四步：缩小，计算点x_c：若f_n≤f_r<f_n+1，设点若f_c≤f_r，用x_c代替x_n+1，回到第一步；反之，跳到第五步；若f_r≥f_n+1，设点若f_cc<f_n+1，用x_cc代替x_n+1，回到第一步；反之，跳到第五步。

第五步：缩减，对于除去最优点外的所有点，用下面的点代替：

x_i＝x₁+σ(x_i-x₁)foralli∈{2,…,n+1},回到第一步。

为了验证在实际应用中，混合模型比单一模型能更好的描述空间变异性，这里进行了变差函数拟合仿真，分别选用高斯-指数混合模型、高斯模型、指数模型作为理论变差函数模型，用Nelder-Mead（单纯性）算法分别拟合，并比较拟合效果。

实例的原始数据来源于工程地质体打井数据，取样数据共112个，井数据的格式如下：

HorizonXYDepthTimeVelocity

其中，X和Y表示井的坐标；Horizon表示当前的数据对应的层位；Depth表示当前层位的深度值；Time表示当前层位的层位时间；Velocity表示层速度。在仿真过程中，只用到X，Y，Velocity三个值。

变差函数拟合的仿真结果如图4、5、6所示，其中，图4为高斯模型变差函数示意图，图5为指数模型变差函数示意图，图6为高斯-指数混合模型变差函数示意图。比较三幅图中的拟合结果，可以发现，高斯-指数混合模型拟合之后的变差函数曲线更接近原始的实验变差函数，除了距离较大的几个点，拟合后的曲线几乎可以和实验变差函数重合。

步骤4.利用三个轴向的比例函数构造权重矩阵，将步骤3得到的三个方向上的变差函数统一套合得到最终的各向同性变差函数γ(h)；

在工程实现中，会遇到不同方向变差函数模型不同的情况，在此提出一种近似的统一套合方案，得到一个各向同性的模型，来解决遇到的问题。下面介绍这种套合形式的具体过程：

已知三个轴向方向（X，Y，Z）上的变差函数为γ_x(h)、γ_y(h)、γ_z(h)，定义比例函数如下：

$W_i\left(h\right)=\frac{{\mathrm{\&gamma;}}_i\left(h_i\right)}{\underset{i=x,y,z}{\mathrm{\&Sigma;}}{\mathrm{\&gamma;}}_i\left(h_i\right)},i=x,y,z---\left(6\right)$

比例函数W_i(h)表示某个方向上的变异性所占的比重，其中，h_i表示矢量距离h在X，Y，Z方向上的距离分量。

根据比例函数构造权重系数矩阵，如式（7）所示：

$D=\left[\begin{array}{ccc}{W}_x\left(h\right)& 0& 0\\ 0& W_y\left(h\right)& 0\\ 0& 0& W_z\left(h\right)\end{array}\right]---\left(7\right)$

统一套合后将得到各向同性变差函数γ(h)，此时三个坐标轴的各项同性变差函数只和距离矢量h的模||h||有关，形式如下：

$\mathrm{\&gamma;}\left(h\right)=D\&CenterDot;\left[\begin{array}{c}{\mathrm{\&gamma;}}_x\left(\right|\left|h\right|\left|\right)\\ {\mathrm{\&gamma;}}_y\left(\right|\left|h\right|\left|\right)\\ {\mathrm{\&gamma;}}_z\left(\right|\left|h\right|\left|\right)\end{array}\right]---\left(8\right)$

其中，参数h表示空间矢量距离，||h||表示标量距离值。

由式（8）可以看出：只要轴向上的变差函数为非负，则套合后的变差函数模型也为非负，满足理论变差函数模型的要求，套合后的模型结构随不同方向动态变化，其比例由各个轴向的距离和变差函数确定。

这里，利用各轴向上比例函数可以很好的将各轴向上的变异性所占的变异比重表示出来，为空间统一套合提供思路，通过权重系数矩阵将各轴向上相互不同的变差函数模型转变为最终的各向同性的变差函数，实现了套合。

步骤5.利用步骤4得到的变差函数γ(h)，求解克里金方程组，对未知点进行速度插值，最终生成速度模型。

在得到套合后的变差函数统一模型之后，就可以预测待估点的速度值。克里金预测方法包括普通克里金、简单克里金、泛克里金、指示克里金等，这里采用的是普通克里金，普通克里金方程组如（9）式所示：

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

通过求解方程组，在得到一系列权值λ_i后，通过式（10）得到待估点的速度值：

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

其中，Z(x₀)表示在x₀点处的层位速度值（未知），Z(x_i)表示在x_i点处的层位速度值（已知），λ_i表示权重。

最后通过插值得到各个未知点的速度，然后建立三维速度模型。

本方案解决了获取层位速度空间变异性时的不精确问题，同时为实现变异性的空间统一套合提供了一种可行的方法。本发明的层位速度建模方法采用Nelder-Mead算法进行拟合，可以处理大数量级的数据；该方案考虑了不同尺度上的变异，采用混合模型代替原有的单一模型，用几种单一的变差函数混合的方式来描述变异性，提高了拟合精度；统一套合方法结构简单，解决了三维空间中不同方向变差函数类型不同时的套合，适用于多种复杂变差函数组合，具有普适性，并在工程实现中应用，具有可行性。

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

Claims (1)

1.一种层位速度建模方法，包括如下步骤：

步骤1.导入样品点的坐标及属性值；

步骤2.分别计算X轴、Y轴、Z轴上的实验变差函数γ*_x(h)、γ*_y(h)、γ*_z(h)；

步骤3.利用步骤2得到的三个方向上的实验变差函数，选取理论变差函数模型并进行拟合，得到X轴、Y轴、Z轴上的变差函数γ_x(h)、γ_y(h)、γ_z(h)；

所述的理论变差函数模型具体为：高斯-指数混合模型，所述的拟合具体采用Nelder-Mead算法；

所述的高斯-指数混合模型具体如下：

$G\left(h\right)=C_0+C_1(1-e^{-\frac{3h}{a_1}})+C_2(1-e^{-\frac{{\left(3h\right)}^2}{{a_2}^2}})$

其中，C₀表示块金常数；C₁、C₂表示拱高；a₁、a₂表示变程；

步骤4.利用三个轴向的比例函数构造权重矩阵，将步骤3得到的三个方向上的变差函数统一套合得到最终的各向同性变差函数γ(h)；

所述的得到最终的各向同性变差函数γ(h)的具体过程如下：

定义比例函数如下：

$W_i\left(h\right)=\frac{{\mathrm{\&gamma;}}_i\left(h_i\right)}{\underset{i=x,y,z}{\&Sigma;}{\mathrm{\&gamma;}}_i\left(h_i\right)},i=x,y,z$

比例函数W_i(h)表示某个方向上的变异性所占的比重，其中，h_i表示矢量距离h在X，Y，Z方向上的距离分量；

根据比例函数构造权重系数矩阵：

$D=\left[\begin{array}{ccc}{W}_x\left(h\right)& 0& 0\\ 0& W_y\left(h\right)& 0\\ 0& 0& W_z\left(h\right)\end{array}\right]$

统一套合后得到的各向同性变差函数γ(h)形式如下：

$\mathrm{\&gamma;}\left(h\right)=D\&CenterDot;\left[\begin{array}{c}{\mathrm{\&gamma;}}_x\left(\right|\left|h\right|\left|\right)\\ {\mathrm{\&gamma;}}_y\left(\right|\left|h\right|\left|\right)\\ {\mathrm{\&gamma;}}_z\left(\right|\left|h\right|\left|\right)\end{array}\right]$

其中，||h||表示标量距离值；

步骤5.利用步骤4得到的变差函数γ(h)，求解克里金方程组，对未知点进行速度插值，最终生成速度模型。