属性建模中基于变差函数理论的数据分析方法

2020-12-14王大成刘义坤付青春曲国辉白军辉李文龙

成都理工大学学报（自然科学版） 2020年6期

王大成, 刘义坤, 付青春, 曲国辉, 白军辉, 陈奋, 李文龙

(1.东北石油大学提高油气采收率教育部重点实验室，黑龙江大庆 163000；2.大庆油田有限责任公司第三采油厂，黑龙江大庆 163113; 3.中国石油大庆油田有限责任公司,黑龙江大庆 163113)

在地质建模的数据分析过程中，变差函数是最基本与最重要的模拟工具，它用于描述数据值的空间相关性，其结果直接影响模型中各种属性的计算和分布，是属性建模的计算依据与核心[1-2]。近年来关于变差函数的理论描述和应用研究较多，大部分都集中在概念和理论方法的描述、操作流程设计等方面，还有一部分开展了变差函数在沉积微相自动识别中的应用研究。这些成果的取得，使地质建模人员在变差函数的基本理解、认识和应用方面有了较大地提高[3-7]。但是综合国内外研究成果发现，由于其理论应用难度大、相关参数较多且相互关联，实际应用中调整过程复杂，相关研究依然缺少结合实际模拟过程数据分析[8-10]；尤其是数据分析过程“假拟合”及“无效搜素”现象尚无专家学者发现、解释及评价。由于数据分析结果直接影响属性预测实际计算精准程度[11]，笔者在变差函数理论基础上，结合区域地质情况，通过变程对属性计算影响实验，对变差函数的参数设置界限和范围进行研究，给出“假拟合”和“无效搜索”现象合理解释，设计既符合变差函数原理又结合实际认识的技术流程，明确参数设置方法及范围，实现属性分布精准模拟。

1 变差函数原理

变差函数是地质统计学基本工具，克里金插值和随机模拟中基于象元算法大部分都需要用到变差函数，它能够表征区域化变量的空间变化特征，特别是透过随机性反映区域化变量的空间结构性[12]。

1.1 变差函数定义与邻域参数

设Z(x)为满足二阶平稳假设的区域化随机变量，h为2个样本点空间分隔距离，Z(xi)和Z(xi+h)分别是Z(x)在空间位置xi和xi+h处的实测值，那么把在某一α方向上相距h的2个区域化变量增量的方差之半称为变差函数[13]，计算公式为

(1)

其中N(h)是间距为h的所有样本点对数量；i=1,2,…,N(h)。在本征假设条件下，变差函数仅依赖于分隔距离和方向，而与位置无关，因此它是一个关于矢量距离的连续函数。

变差函数在实际求取过程中，通常样本数据有限且多数不规则分布，如果按照精确条件不可能获得足够的点对数进行变差函数计算[14]。为了获得可靠的变差函数，在数据分析过程中，引入经典的截断楔形来定义样品邻域，方向和步长分别给出容差范围。同时由于角度容差和步长容差的存在，使搜索域偏离主方向的范围越来越大，为此使用偏离主方向线的一个固定带宽来限制[15]。只要点对在上述条件范围内，就认为该点对可以参与到变差函数值计算中来。各参数含义见图1。

1.2 函数计算与曲线拟合

已知一个研究区内有N个数据点，为了拟合某一方向变差函数曲线，需要计算该方向上若干不同距离变差函数值，其中选取的基本距离h称为基本滞后距[16]。首先从某点开始根据给定方向、基本滞后距等条件筛选点对。由于函数计算精度与数据对个数成正比关系[17]，因此为了保证有足够多的点对参与计算，将所有满足邻域参数设定范围的数据点对进行筛选，再进行变差函数值计算(图2-A)。以此类推，分别依次计算2h、3h、…、mh倍滞后距变差函数值γ(h)，将其标注在坐标轴上，就得到变差函数散点图，用曲线将γ(h)数据点连接就得到变差函数趋势图(图2-B)。其中sill是基台值，表示变量在空间的总变异程度；C0是块金值，表示在很小间距时的空间变异程度；拱高Ct表示取得有效数据尺度上可观测到的变异程度；a是变程，反映区域变量相关性范围，变程范围内数据具有相关性。

在实际地质建模过程中，通常都是利用主流大型软件的数据分析模块，通过手动拟合的方式，充分考虑地质因素的基础上，以sill、C0、Ct数值为参考依据进行不断调整拟合，求得主方向、变程值等参数进行后续克里金插值和属性计算。

2 “假拟合”与“无效搜索”现象

按照数据分析过程和序贯高斯参数引入条件，通常在确定主方向情况下，通过调整各类楔形参数范围，求取的变程值是唯一的确定的参数。但是在数据分析过程中，数据、模型、计算和分析方法均相同情况下，通过大量的变差函数曲线拟合实践发现了“假拟合”和“无效搜索”现象。

“假拟合”现象是指在同一沉积单元、同一种属性、同一方向，只是楔形参数不同的情况下，出现多套实际数据与理论模型拟合，从而计算得出多个不同变程值(图3)。如果按照传统方法及曲线拟合标准，多个变程值均符合要求；但是由于序贯高斯模拟只能引入一套变程参数，多变程结果误导数据分析计算结果，这种求解的不确定性直接影响属性预测准确性。

“无效搜索”现象是指曲线拟合过程中其他参数不变的情况下，当搜索半径达到一定数值后，增加滞后距数目和搜索半径，拟合曲线均不发生明显变化(图4)，极大地消耗工时却没有效果。

上述两种现象是在大量数据分析精细研究过程及变差函数曲线拟合实践中发现的。由于变程是属性计算核心，因此笔者通过变程对属性计算影响对比实验进行研究。

3 变程对属性计算影响对比实验

3.1 实验原理

序贯高斯模拟是应用最广泛的连续性变量随机模拟方法，该方法是被认为模拟连续变量的首选方法，其对变差函数的计算结果具有明显的依赖性和敏感性[18]。同理在设定基本条件的情况下，序贯高斯的属性模拟结果也可直接反映变差函数计算效果。因此，为了对比变程对属性计算结果的影响，结合序贯高斯模拟算法要求，对研究区按照模型、网格系统、算法等同一性，因素改变、评价指标对比唯一性，对比数据以计算结果为准客观性3个原则，分别计算不同变程模型渗透率总值、均值，通过数据对比及规律分析，明确变程对属性计算的影响程度。

3.2 实验方法

首先优选高密度井网研究区的高质量地质模型作为实验模型，以满足变差函数参数求取以及序贯高斯模拟过程点对数据数量和质量的首要需求，保证基础资料同一性；第二步对测井解释渗透率数据网格离散化，通过数据分布的前后端截断滤掉奇异值，利用对数变换减少渗透率属性的离散性，同时按照累计概率分布曲线完成正态变换，保证了数据输入端完全符合变差函数计算和序贯高斯模拟精准计算的要求；第三步是参数设定方面。搜索角度为南北向-5°，随机种子数设定为5 771，保证计算过程随机路径统一。同时为满足点对采样需求、数据多样性和实验结果可靠性，笔者结合研究区井距实际加密实验频次：平面主变程按照50 m间距等差递增；次变程充分考虑各向异性影响，按照主次变程10∶9的比例递增；垂向变程根据储层平均厚度设定为5 m，其他参数按照软件默认设定。这样既满足随机模拟过程参数获取规则，又严格保证变程值成为实验唯一改变因素，模型参数见表1。最后为了保障具有足够数据的统计学规律，开展了67组不同变程的属性计算实验，采用序贯高斯球状模型算法计算得出不同变程条件下属性体总值，再通过计算结果多方面对比研究变程对属性计算结果的影响。

3.3 实验结果

按照上述实验原理与设计方法，将不同变程值与模型渗透率体计算结果制成图5(为表示清晰，横坐标用主变程代表主次变程，下同)。从整体看渗透率体总值较高，但是未出现各类数据间无规律突变，说明实验在同一和唯一方面设计具有合理性。具体来看属性体在50～150 m变程阶段，渗透率体总值呈上升趋势，150 m阶段达到最大值,为1.01×105μm2。从200 m变程开始总值迅速降低后逐渐平缓下降，在 3 350 m变程时渗透率体总值达到最低,为0.77×105μm2。

3.4 实验分析

3.4.1 整体规律分析

由于研究区平均井距为250 m，根据变差函数原理认为当变程小于数据点对距离时，将会忽视重要的结构信息，从而使地质统计学基于空间相关性进行估计的优势不复存在[19]。为了进一步证实该理论，笔者以研究区某储层为例建立50～150 m变程条件下计算模型(图6-A、B、C)，通过与图6-D的沉积微相解剖成果对比可以看出属性分布非常离散，几乎体现不出储层发育的连续性和趋势性，说明变差函数和序贯高斯模拟失去本质意义，由此判断变程从50～150 m计算结果的短暂突增阶段不具有参考性。因此根据图5实验结果发现规律：从200 m变程值开始，随着变程值的不断增大，属性体计算总值呈现先陡后平的一致性递减趋势，从某一变程值开始，随着变程增加计算结果趋于稳定，变程数据不再具有明显计算意义，笔者将这个值称为“变程阈值”。

表1 地质模型及属性计算基本参数Table 1 Basic parameters of geological model and attribute calculation

实践表明，变差函数曲线拟合十分繁琐，如果能够确定变程阈值，将会很大程度避免“无效搜索”现象，减少数据分析工作量，提高工作效率。更为重要的是由于阈值对前期数据分析过程变程值求取的范围限定具有重要指导作用，使变程值有了较为明确的量化范围，可有效避免盲目试算和经验判定带来的不确定性，从而进一步提高属性分布预测的精度。

3.4.2 确定变程阈值

研究表明搜索半径通常小于最远两个数据点的1/2，根据变差函数理论，变程值一定小于搜索半径，因此可以得出变程值小于研究区最远数据点1/2，即小于 1 600 m。为了进一步明确范围，根据前述实验计算结果，假设每个不同变程值对应模型属性体值为γ(h)，将间隔d距离的值相减即为值差γ(h)d。以值差后的变程值为横坐标，以属性体差值为纵坐标就可以得到每相差d距离条件下影响趋势图(图7)。为了让值差影响程度表现得更直观和明显，将每个网格分别赋渗透率值(10、20、30)×10-3μm2计算其总值γ(h)*，用γ(h)d/γ(h)*的倍数关系对比相差d距离变程对属性计算结果的影响程度。按照上述方法，根据实验数据绘制不同变程变化50 m时属性体差值对比图(图7)。

图7曲线整体趋势表明原始变程值越小，相差50 m情况下对其属性计算结果影响越大，变程值越大影响程度越小。实验数据显示，影响程度最大的是变程值200 m时，相当于出现18个均值渗透率10×10-3μm2、9个均值渗透率20×10-3μm2、6个均值渗透率30×10-3μm2储层误差，最小相差倍数趋近于0。这是由于变程阈值后属性体总值趋于稳定差值，接近于0，导致影响程度低；但是这不说明阈值后的数据误差更小，反而进一步证明了前面论述的阈值后属性体总值数据没有意义以及确定阈值对属性影响精准评估和计算的重要性。进一步统计发现变程在1 100 m以后属性体差值在0.02×105μm2数值上连续保持稳定的6个数据点后逐渐趋于0值，次坐标轴的属性相对影响曲线显示，在1 100 m以后20×10-3μm2、30×10-3μm2相对影响曲线接近重合且相对影响储层个数小于1，说明在该变程值后属性值变化趋于稳定，对属性计算结果影响很小。再根据实验初期设置的容差参数±50%计算，实验结果与变差函数理论计算最大变程值范围吻合，说明1 100 m阈值分析结果是可靠的，因此综合分析认为1 100 m为该模型的变程阈值。

3.4.3 变程值对属性计算影响

笔者按照变程值从200 m到1 100 m阈值间的影响储层相对个数平均值对影响程度评价发现，变程值变化50 m时对储层属性差值影响相当于均值30×10-3μm2储层相差2.5个，说明变程值的变化对属性计算影响较大。为了完善这项研究成果，进一步证明50 m变程差的相关研究不是个例，对研究区分别进行变程值相差50 m、100 m、150 m、…、900 m的18组属性体影响情况数据统计(图8)。

从原始变程值与渗透率体值差关系曲线发现，所有变程值差均呈现随着原始变程值增大，渗透率体值差变小的一致性趋势，说明变程值变化对较小变程属性计算影响较大，对较大变程属性计算影响较小。从实际统计数据看，以均值30×10-3μm2储层为对比标准，变化50 m是渗透率体差值最小约2.5个储层；随着变程值变化幅度的增加，影响储层相对数量不断增加，当变化幅度达到900 m时，相对影响储层数量达到28个，进一步说明变程值求取大小和精准程度对属性计算结果影响较大。此外，根据统计数据进行了趋势线关系公式的拟合，精准度达到 0.998 5，因此在实际属性建模过程中，结合误差允许范围，可以用关系公式进行变差值对属性影响评估。同时根据实际误差允许范围需要实现变程差值的反算，这对进一步明确变程范围具有重要意义。

上述一系列实验、数据统计及分析结果表明变程值的变化对属性计算影响具有较强趋势性和规律性，但是由于对计算结果影响偏差较大，因此属性建模过程要基于变差函数理论及序贯高斯模型，充分结合储层发育实际及各参数调整范围进行数据分析，减少变程求解过程“假拟合”和“无效搜索”等不确定因素，从而提高属性预测精度。

4 “假拟合”与“无效搜索”的解释

两种现象都给数据分析过程造成了很大的困扰，一方面是无法保障精准计算；另一方面是在基于变差函数理论的数据分析过程中，由于工作繁琐、难度较大且不容易理解，初学者往往在拟合过程耗费大量时间，一直无法求解变程，经验者往往对数据分析过程不够重视或忽略，盲目设定导致计算结果偏差较大。笔者尝试从点对搜索域对变差函数影响的角度进行解释，为高效精准的数据分析提供理论支撑。

“假拟合”现象是在其他参数确定条件下通过调整搜索半径和滞后距数目过程中发现的，滞后距基本保持不变。根据数据点对搜索策略，搜索方向和滞后距保持不变说明搜索数据点对条件一致，假设原始滞后距数目为m个，拓展搜索半径滞后距数目增多到(m+n)个，那么拓展后增加的点对数据全部来自于搜索半径变大产生n倍滞后距内点对。由于变差函数是依赖点对数据进行计算的矢量值，因此两者前期m个滞后距的拟合效果必然相同。根据上述原理，如果拓展范围较小，拓展半径对变差函数计算及模型拟合趋势并不会有明显影响；但是如果出现较大幅度拓展(如图3的倍数级增加)，点对数目大幅度增长就可能出现不同的拟合趋势，从而形成“假拟合”现象。需要注意的是，增加搜索半径毕竟增加了点对数目，增加了数据计算结果更多的可能性，因此“假拟合”现象的发现表明仅仅按照传统方法和标准完成函数曲线拟合，得到的计算结果可能不是最终确定值。根据前述实验数据和变差函数原理可知，若数据分析实践中出现“假拟合”现象，要充分考察实际区域地质、参数设置和调整范围的合理性，同等条件按照“比较选用较大变程值”原则有效避免“假拟合”影响，减少属性计算不确定性。

“无效搜索”产生的关键就是对搜索半径和滞后距数目同时较小范围调整，这个过程的实质是对于计算核心的基本滞后距并未明显改变，在变差函数拟合过程中只是发挥扩大搜索半径和增加点对数目，这一点与前面多重拟合现象类似；另外一个原因就是搜索半径等参数设置已经覆盖整个研究区域，所有点对数据已经全部参与计算，那么此时再增加搜索半径和滞后距数目与最大覆盖条件下所形成的邻域范围基本一致，因此必然产生“无效搜索”现象。

综合上述的理论解释和实践研究表明：变差函数值是滞后距的函数，它是以基本滞后距为基础，通过楔形邻域控制，获取点对再进行数值计算的过程。因此滞后距是变差函数计算的核心，滞后距调整优先级要高于搜索半径。同时，数据分析过程中参数设置和范围调整至关重要，其调整与设定具有限定范围，合理设定及量化求解是属性精准计算的重要基础。

5 技术流程设计及参数设置

为了能够通过参数设置获取更精确的建模效果，提高工作效率，笔者基于变差函数的基本理论，在前人研究成果基础上设计了新的数据分析技术流程，综合本文实验分析和理论解释成果进一步明确参数调整范围。

5.1 技术流程设计及参数范围

整个数据分析技术流程共分为8个步骤，每个步骤明确参数范围设定(图9)。

①准备及处理数据：通过截断变换滤掉不合理奇异值，使数据成正态分布，再根据建模算法需要对数据进行变换，如序贯高斯模拟算法对渗透率建模时，就需要对数据做对数和标准正态分布变换。

②选择拟合模型：常用模型主要有指数、球状、高斯3种模型，各类模型的适用条件[20]如表2。

表2 变差函数拟合模型类别Table 2 Categories of variogram fitting models

③给定搜索方向：主方向是变差函数拟合和序贯高斯模拟引入方向，是第一优先设置参数，需结合储层精细描述成果，通常与沉积物源方向一致。

④设置基本滞后距：根据前面发现的两个现象以及所带来的启示，技术流程优先级上，滞后距设置要高于搜索半径，这是与其他很多传统流程不一致的地方。根据实验分析结果，为保证数据点对都能参与到计算中来，基本滞后距至少应该大于平均井距2/3，通常设定为沿某一方向平均井距。

⑤明确搜索半径：搜索半径是滞后距个数与基本滞后距乘积，由于上一步已经确定基本滞后距，因此可通过调整滞后距数目确定搜索半径。结合区域实际大小，根据数据点对实际需求及搜索半径的邻域法则，基本滞后距≤搜索半径≤最远数据点距离1/2。

⑥微调滞后容差：由于滞后容差和角度容差等参数在变差函数计算中具有意义，就是通过增加邻域提高点对数目，从而提高与模型曲线拟合精度，因此应是设定主参数之后进行滞后容差和角度容差调整。主流软件已经给出滞后容差最小值50%，但是最大限定值并没有明确。根据楔形邻域法则，滞后容差是为了满足点对数目条件而设定的搜索条件域之一，每一个范围内点对数量是由其搜索域决定，如果相邻搜索域重叠，其域内必然出现点对的重复搜索和计算，这也是“假拟合”现象出现的原因之一。这就要求滞后容差既要保证点对数量设定邻域，还要减少搜索域重叠带来的影响。由于滞后容差是按照参数正负比确定范围，因此根据点对搜索策略，在两个相邻滞后距搜索域内，只要前一滞后距最大值小于等于相邻滞后距的最小值，即可满足上述条件要求。其原理和范围限定可表述如下：假设滞后容差增加百分比为x，基本滞后距为h，任意滞后距为基本滞后距m倍，表达为mh，相邻滞后距为(m+1)h，根据前述理论分析，则应有

(mh+0.5mh)(1+x)≤[(m+1)h+

0.5(m+1)h](1-x)

(2)

对公式(2)整理得公式(3)

x≤1/(2m+1)

(3)

其中公式(2)和公式(3)中m是≥1的正整数。

从公式(3)可以看出，滞后容差参数增加百分比与滞后距倍数有关，随着滞后距倍数增加，滞后容差参数增加百分比逐渐变小。由于搜索半径包括多个倍数滞后距，上述公式求得的x是不同倍数滞后距下的动态数值，因此若想保障搜索域全部不重叠，m值应是上述第⑤步确定的最大滞后距个数，将其代入公式(3)即可得到增加容差最小百分比，即接近于0。由于m的最小值为1，根据公式(3)计算x最大值为1/3，因此滞后容差增加调整幅度是0%～33.3%。

⑦调整角度容差：根据行业认可带宽设定方法，将带宽设置为2倍滞后距；但是当容差角小于某一角度时，带宽设定将失去明显意义。因此，根据带宽参数设定意义，假设容差角为α，搜索半径为nh，带宽为2h，根据三角函数关系有tanα=2h/nh=2/n，α=arctan(2/n)，因此角度容差调整范围在α=[arctan(2/n)]°～90°之间。

⑧曲线拟合求取变程：通过参数调整，按照块金值尽可能小(接近于0)、基台值接近1、变程值尽可能大的原则拟合模型曲线，保证计算函数曲线与模型曲线基本重合来求取变程值。这里需要说明的是，变程值尽可能大是为了有效避免“假拟合”现象的误导而提出，并非完全无限制条件取极大值。因此，实际建模中要结合“变程阈值”范围约束同时避免“无效搜索”现象的存在。

综上所述，各参数间虽然有相关性，但是参数设置却有先后顺序及范围限定，参数调整可以按照一定步骤、一定参考值范围进行。这样不仅符合变差函数的基本理论，实现了变程值全量化参数求取，保证属性预测的精准程度，更是便于设计技术流程，极大提高工作效率。

5.2 建模效果

研究区储层整体划分为3种类型，分别选取Ⅰ、Ⅱ、Ⅲ类储层的典型层，按照新设计流程和参数设置方法进行数据分析，获得3个典型层的主次变程值分别是818/706 m、664.5/615.1 m、408/370 m，再利用序贯高斯算法进行属性分布计算，图10-A、B、C是3种类型典型层的渗透率预测结果。

为验证新方法属性模拟的精准程度，开展井网抽稀验证，结果表明模型精度达到91.4%，较该区域原始模型精度提高4%以上。此外，由于沉积微相精细描述成果能够直接有效地反映属性分布特征，因此为了进一步验证和体现建模整体效果，笔者再用沉积微相精细描述成果(图10-D、E、F)对渗透率的属性分布模型进行比对验证。前人研究已经表明研究区整体属于大型河道或相对稳定的三角洲沉积环境，河道砂体在沉积过程中对储层发育具有良好的稳定性和连续性。从图10-A、B、C属性预测结果看，渗透率分布整体连续性较好地反映了上述储层发育机理，代表渗透率预测值较高的红色区域的外部形态、方向和条带趋势，均很好地反映了沉积微相精细描述成果中的主河道整体走向、弯曲形态以及分流河道的窄小特征，这在3种类型储层中均有所体现。同理，对于表外和尖灭发育较多的薄差部位，在模型中也通过渗透率的极低值有所反映，这在Ⅱ类和Ⅲ类储层的效果对比中体现较为明显(图10中B、E和C、F对比)。因此综合Ⅰ、Ⅱ、Ⅲ类储层模拟对比发现，该方法不仅能够对好的储层分布模拟，还能对发育不好、连续性较差的储层进行较好模拟。