一种基于相控非参数各向异性变差函数构建方法与流程

本发明涉及地质统计学领域，特别涉及一种变差函数构建方法。

背景技术：

随机建模的主要工具是地质统计学，地质统计学是利用变差函数去统计采样数据间的空间关系规律，然后利用这个统计规律去模拟未采样位置处随机变量的值的一种统计方法。地质统计学在研究区域化变量的空间分布结构特征规律性的基础上，选择各种合适的克立格法，以达到更精确的估计或对区域化变量进行条件模拟的目的。它一般包括三大基本组成部分：空间函数的相关性分析、克立格估计和条件模拟。空间函数相关性分析是对变差函数和协方差函数的分析，包括它们的定义和估计方差，它是克立格估计和条件模拟的基础。

自地质统计学问世以来，对地球科学中有关资源信息的空间变异性用区域化变量理论进行表征，不仅对地质预测和评估理论有很大发展，也对许多非地质科学如水土资源、水文气象、生物环境和工程技术产生了广泛影响，目前已经形成对空间信息分布特征或模拟其离散性和波动性的研究，均可用地质统计学理论进行，而地质统计学也已经成为评估各种区域性自然现象、自然资源及再现其波动性过程的新的工程科学。

变差函数作为地质统计学的一个基本工具，是随机建模的关键因素，其主要作用是刻画空间变量的空间关系结构，在建模中可以帮助人们分析储层特征、获得特征。变差函数可以充分描述储层参数空间相关性，在不同的沉积条件下，储层的空间分布特征各不相同，尤其是非均质性较强的地区，更需要用变差函数来准确估计储层的各项异性特征。因此变差函数的构建是随机模拟和随机反演的核心技术。

变差函数

地质统计学中用区域化变量来表示所研究的对象。区域化变量是以空间点X的三维直角坐标(x_u,x_v,x_w)为自变量的随机场Z(x_u,x_v,x_w)＝Z(x)，它反映了地质变量特有的双重特征，即随机特征和空间结构特征。

在点x处的变差函数定义为x处相距为h的区域化随机变量之差的方差之半，记为γ(x,h)，即

r(x,h)＝0.5*D²[Z(x)-Z(x+h)]

＝0.5*E[Z(x)-Z(x+h)]²-0.5*{E[Z(x)-Z(x+h)]}² (1-1)

变差函数有三个基本参数：变程a、块金值c₀及基台值c。

①变程

变程a反映了区域化变量的影响范围，在变程范围之内具有相关性，且随着距离的增大相关性逐渐减弱，超出该范围则不具有相关性。

②块金值

块金值c₀＝γ(0)，理论上变差函数在原点处的值应为零，但是由于测量的误差和微观变异性，使得在很小的距离之内随机区域化变量之间出现了较大的变化。

③基台值

基台值反映了区域化变量在研究区域范围内的变异强度。

变差函数能很好地描述地质变量的特征，如变程反映了变量的影响范围，基台值反映了变量的变异强度，变差函数在零点附近起始段的斜率反映了在较小的距离内，变量的变化是剧烈还是平缓，不同方向上的变差图可反映变量的各向异性等等^[1]。

随机模拟

变差函数的构建是进行随机模拟的前提。随机建模是指根据收集到的信息，运用随机函数理论，采用随机模拟方法，从而产生可选的、等概率的、高精度的储层模型的方法。这种方法承认控制点以外的储层参数具有一定的不确定性，即具有一定的随机性。其流程如图1所示。

技术实现要素：

本发明为解决上述技术问题，提出了一种基于相控非参数各向异性变差函数构建方法，结合工区不同的沉积相相位，在不同沉积相分别构建非参数各向异性变差函数，拟合参数时提出一种基于蚁群算法的方式，可以较好的估计出变差函数对应参数。为进一步的模拟工作提供良好的技术支持。

本发明采用的技术方案是：一种基于相控非参数各向异性变差函数构建方法，包括：

S1、工区沉积相划分，采用指示模拟对工区沉积相进行划分；

S2、构建不同沉积相的变差函数，具体包括以下分步骤：

S21、初始化：导入已经反演好的波阻抗数据和测井数据；

S22、相控建模：利用PETREL软件针测井数据进行相控建模，得到每个点对应的标签phase_label；

S23、对角度进行等分，设置角度区间数量为M等分，记等分M＝degree_num，间隔记为θ，θ＝360/M，则第n个角度区间对应斜率为[arctan((n-1)·θ),arctan(n·θ))；

S24、随机选取点对，将不同标签下的点进行分类，任意两点随机组合成点对，将每组点对都进行标号；

S25、计算所选取点对的斜率，根据其取值将该点对放入对应的角度区间；

S26、对各的沉积相分别重复步骤S33至步骤S35，直到各沉积相的各个角度区间都有足够的采样点；

S27、根据不同沉积相下的采样数据分别计算变差函数；

S3、拟合变差函数参数求解，包括：

S31、对候选解{x₁,x₂,…}的每一变量x_i用字长为N的二进制码串{b_Nb_N-1…b₁b₀}进行编码，根据下面的公式进行解码：

其中，b∈{0,1},j＝1,2…N，b_N-1为最高位，b₀为最低位，变量x_i的左边界为实数值x_imin，右边界为实数值x_imax，z表示二进制码串对应的十进制整数值的左边界为实数值；

S32、将待拟合参数转化为有向图的形式；

S33、采用指数模型进行参数拟合。

进一步地，所述步骤S1之前还包括：指示变换，具体为：根据门限值，把连续确定的原始数据离散成布尔量0或1。

进一步地，所述步骤S32具体为：

定义有向图G＝(C,L)，其中顶点集C为

其中，v_s为起始点，顶点和分别表示二进制码串中位b_j取值为0和1的状态，j＝1,2…N，c₀(v_s),均表示顶点集C中的元素，e＝1,2,…,2N。

进一步地，所述步骤S33具体包括以下分步骤：

S331、nc＝0，各τ_ij和Δτ_ij的初始化，将m个蚂蚁置于n个顶点上；

其中，nc为迭代步数或搜索次数，τ_ij为t时刻在i,j连线上残留的信息量，为蚂蚁k在边弧(i,j)上留下的单位长度轨迹信息素数量；

S332、将各蚂蚁的初始出发点置于当前解集中，对每个蚂蚁k(k＝1,…,m)，按概率移至下一顶点j；将顶点j置于当前解集；

S333、计算各蚂蚁的目标函数值Z_k，记录使得目标函数值Z_k最小的解；k＝1,…m；

S334、设Path^*(t)为第t搜索周期内的最佳路径，该最佳路径对应的目标函数值为f^*(t)，边弧(i,j)中顶点i对应候选解的第Κ位，则蚁群搜索的信息素按照下式更新

其中，f^*(t+1)表示第t+1搜索周期内的最佳路径对应的目标函数值，τ_ij(t,k)表示t时刻的信息素，τ_ij(t+1,k)表示t+1时刻的信息素，L为候选解的二进制编码的编码长度，为正整数，k表示第k个蚂蚁；

S335、对各弧边(i,j)，置Δτ_ij＝0；nc＝nc+1；

S336、若nc<预定的迭代次数，转步骤S332。

本发明的有益效果：本发明基于地层不同沉积相进行不同的各向异性变差函数拟合，实现了对指数型变差函数的构建和求解，通过合理的相建模，具有以下优点：

(1)本发明所构建的变差函数是基于不同沉积相的各向异性变差函数，与传统的变差函数构建方式相比，本发明拟合出的变差函数更接近于真实的地质情况，有效避免了因变差函数构建不合理而在随后的随机模拟工作中所产生的误差；

(2)本发明建立了一种基于反演数据的点对随机选取方法，通过在各个角度区间点对的随机选取作为变差函数构建时的参数输入，能更好的反应变差函数的各向异性特征；

(3)本发明提出了一种基于蚁群算法的变差函数参数拟合方法，通过将问题的转化，能够提高变差函数参数拟合的准确性，为以后的工作提供了技术保障。

附图说明

图1为随机模拟的流程图。

图2为本发明提供的变差函数的构建流程图。

具体实施方式

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

如图1所示为本发明的方案流程图，本发明的技术方案为：一种基于相控非参数各向异性变差函数构建方法，包括：

S1、工区沉积相划分，采用指示模拟对工区沉积相进行划分；指示模拟是划分工区沉积相的有效方法。在进行模拟之前，本申请先要进行指示变换，即根据一系列门限值，把连续确定的原始数据离散成一个布尔量0或1的过程。设{Z(x_α),α＝1,…,N}是所研究的储集层上的一组观测数据，N表示进行了N次观测，得到N个不同值，给定一个门限值Z，可把所有的观测数据转换成指示值I(x_α；Z)；

则以x₀为中心的待估区的区域化变量I^*(x,Z)的概率估计为

式中的权系数λ_α可由指示Kriging方程组求出，最终可以求得待估区域的估计值Z^*(X)；

式中的权系数λ_b可由指示Kriging方程组求出，Z(X)表示区域变化量，是一个随机变量。

在PETREL软件中可以选择利用指示模拟的方法进行沉积相的划分。在导入工区测井数据后，PETREL软件会根据测井数据进行聚类，一般情况下，三种不同的沉积相就能满足需要，因此，手动调节将其分成三类，之后选择指示模拟方式对其进行建模，就可以得到工区中每个点所对应的相标签，将每个点的标签都保存在phase_label的矩阵中，标签值取为1，2，3。

S2、构建不同沉积相的变差函数，传统的变差函数仅仅利用测井数据进行构建，不能很好的拟合出变差函数曲线。本发明采用已经反演好的数据，利用较为丰富的数据，提出一种点对随机选择的方法，针对不同的沉积相标签，分别进行变差函数的构建，如图2所示，具体包括以下分步骤：

S21、初始化：导入已经反演好的波阻抗数据和测井数据；

S22、相控建模：利用PETREL软件针测井数据进行相控建模，得到每个点对应的标签phase_label；

S23、对角度进行等分，设置角度区间数量为M等分，M＝degree_num，间隔记为θ，θ＝360/M，则第n个角度区间对应斜率为[arctan((n-1)·θ),arctan(n·θ))；

S24、随机选取点对，将不同标签下的点进行分类，任意两点随机组合成点对，将每组点对都进行标号；本申请利用rsenne Twister随机数发生器为每组点对产生一个随机数，rsenne Twister随机数发生器递推公式为

其中，初始值为(x₀,x₁,...,x_n-1)，是一个n维行向量，每一个x＝x^(w)x^(w-1)...x⁽⁰⁾是w位的二进制表示，左边为高位；表示x_s的前w-r位(u为upper)，表示x_s+1的后r位(l为lower)，是由x_s的前w-r位和x_s+1的后r位组成的新的w位数据；1≤m<n，r为字符分割点，0≤r≤w-1；“|”为或运算符，“⊕”为异或运算符，“>>”为右移，“<<”为左移；A是MT算法的回旋变换矩阵，大小为w×w维；“:＝”是一种表示方式，意思是定义等于。

按照随机数大小对点对进行排序，选择N组点对作为变差函数计算输入。

S25、计算所选取点对的斜率，根据其取值将该点对放入对应的角度区间；

S26、对各的沉积相分别重复步骤S33至步骤S35，直到各沉积相的各个角度区间采样点个数大于或等于N；按照划分的角度区间大小进行判断，角度区间较大的话就多一些，较小的话就可以少点，本申请中设定采样点个数N为20个左右。

S27、根据不同沉积相下的采样数据分别计算变差函数；

S3、拟合变差函数参数求解，虽然目前已提出了多种变异函数的拟合方法，但这些方法用起来有很多不便和缺点。如加权多项式拟合法很大程度上依赖地质人员的经验，其权系数的确定也较困难；线性规划法与加权多项式回归方法相比，仅在求取多项式系数方法上有所改进，没有考虑加权问题；目标规划法则又太复杂；加权线性规划法考虑到了不同滞后距h下所得实验变异函数值进行加权，还保证了拟合的成功，同时进行了人工干预，但该方法是在加权多项式拟合法和线性规划法上提出的，目标函数太复杂，很难求解。本发明根据这些拟合方法的特点，提出了一种运用蚁群算法来拟合变差函数参数的方法。包括：

S31、对候选解{x₁,x₂,…}的每一变量x_i用字长为N的二进制码串{b_Nb_N-1…b₁b₀}进行编码，根据下面的公式进行解码：

其中，b∈{0,1},j＝1,2…N，b_N-1为最高位，b₀为最低位，变量x_i的左边界为实数值x_imin，右边界为实数值x_imax，z表示二进制码串对应的十进制整数值的左边界为实数值；

S32、将待拟合参数转化为有向图的形式；定义有向图G＝(C,L)，其中顶点集C为

v_s为起始点，顶点和分别表示二进制码串中位b_j取值为0和1的状态，j＝0，1，2，…，N，c₁是实际的顶点，括号里的表示它的生成方式，上标0,1表示它的两种状态，下标表示从N开始，则有向弧L为

即，对于j＝2,3,…N，在所有顶点和处，分别有且只有指向的两条有向弧；

S33、采用指数模型进行参数拟合，其公式为

待拟合参数为a和c。

由试验变差函数可得a∈(40,100),c∈(8000,16000)，通过步骤S32分别对其进行编码并转化为有向图G的形式。则对于问题

其最小值问题可以转化为以参数a和c为图的蚁群算法最小路径求解问题，其中γ*(h_i)为对应试验变差函数滞后距为h_i的函数值。其求解步骤如下：

S331、nc＝0，各τ_ij和Δτ_ij的初始化，将m个蚂蚁置于n个顶点上；

其中，nc为迭代步数或搜索次数，τ_ij为t时刻在i,j连线上残留的信息量，为蚂蚁k在弧边(i,j)上留下的单位长度轨迹信息素数量；

S332、将各蚂蚁的初始出发点置于当前解集中，对每个蚂蚁k，k＝1,…,m，按概率移至下一顶点j；将顶点j置于当前解集；

S333、计算各蚂蚁的目标函数值Z_k，记录使得目标函数值Z_k最小的解；

S334、设Path^*(t)为第t搜索周期内的最佳路径，该最佳路径对应的目标函数值为f^*(t)，弧边(i,j)中顶点i对应候选解的第Κ位，则蚁群搜索的信息素按照下式更新

其中，f^*(t+1)表示第t+1搜索周期内的最佳路径对应的目标函数值，τ_ij(t,k)表示t时刻的信息素，τ_ij(t+1,k)表示t+1时刻的信息素，L为候选解的二进制编码的编码长度，为正整数，k表示第k个蚂蚁；

S335、对各弧边(i,j)，置Δτ_ij＝0；nc＝nc+1；

S336、若nc<预定的迭代次数，转步骤S332。

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