李杨， 钟宝荣， 许晓宏*， 梁梓君

(1.长江大学地球科学学院, 武汉 430002； 2.长江大学计算机科学学院, 荆州 434100; 3. 中国石油新疆油田公司, 克拉玛依 834000)

在油气勘探开发过程中，三维地质建模能够帮助地质人员了解地质发育特征，辅助油气开采决策[1-2]。空间插值是三维地质建模的关键部分，提高插值的精度可以节省油气勘探开发成本。

目前应用于地质建模的插值方法分为确定性方法和地统计方法。确定性方法通过数学函数，根据相似程度或平滑程度进行插值。地统计方法利用测量点的统计属性进行插值。

在确定性插值方面，赵仕威等[3]利用反距离加权和多重二次曲面函数进行三维地层可视化，在稀疏样本的条件下取得了良好的效果，但是并没有考虑地质发育的各向异性特征。Hillier 等[4]以径向基(radial basis function，RBF)为基础，利用梯度约束的数据驱动特征分析各向异性，能够应用于分散的多元地质结构。庞庆刚等[5]利用紧支撑径向基(compactly supported radial basis functions,CSRBFs)进行三维地质建模，通过样本的法向量来考虑各向异性，实现了在非均匀数据中快速插值。然而基于RBF的地质建模方法普遍过于平滑，在固体矿产方面有优势，但是在油气存储区域不连续且封闭的地质状况中适用性欠佳。另外，RBF中不同的基函数插值效果差异过大。黄牧等[6]利用Voronoi图实现了快速自动化三维地质建模，提高了建模效率，然而Voronoi图是以三角网为基础，与常用的规则六面体相比，必然伴随着数据结构复杂，兼容性差，转换流程耗时等问题。

在地统计插值方面，插值方式主要以克里金(Kriging)为主。张小艳等[7]利用变异过程中修正变异方向的自适应差分进化算法来优化Kriging，进行煤层的三维可视化地质建模。杨阳等[8]通过利用半变异函数理论对样本进行分析，采用普通克里金进行海洋地层厚度的地质建模。蔡文军等[9]利用Kriging方法在含油的灰岩储层中进行三维孔隙压力建模。Li等[10]利用深度神经网络模拟半变异函数进行Kriging插值，减少了理论Kriging方差。然而，Kriging方法具有一定的局限性：首先，Kriging方法应用的前提是样本必须符合特定的随机场(如普通克里金必须要正态分布)；其次，Kriging方法在插值过程需要经过繁琐的人为分析流程，如半变异函数分析，块金、基台、变程等参数的确定。

地理学定律[11-13]认为地理空间分布具有空间相关性和异质性。从空间相关性可以得出地理分布与距离相关，从空间异质性可以得出地理分布具有各向异性特征。在地质统计学中通常利用倾斜角和方位角来计算各向异性特征。由此得出，影响插值的因素包括样本分布的距离、倾斜角、方位角、样本点属性值。相对于其他降维方法，利用多维尺度变换[14-15](multidimemsional scaling，MDS)能够最大限度地减少数据的损失的优势，对影响插值的因素进行分析并降维，提取影响插值的特征，其结果作为训练模型的特征集。相对于人工智能领域中其他的回归算法，利用自适应增强回归[16](adaptive boosting regressor，ABR)具有高精度、强学习器、防止过拟合等优势，对已知样本点的特征进行训练，探索地理分布的空间相关性，建立插值模型。接着，把未知点与样本点的空间相关性信息带入到ABR模型中求得已知点的值。插值模型建立了多个ABR模型，每个ABR模型相对独立，所以最后利用三维均值滤波(mean filtering，MF)进行平滑处理，以此来使结果符合地质属性分布的连续性特征。

在地理学定律和地质统计学的前提下，针对上述文献中的问题，以常见的规则六面体为基础数据结构，提出利用人工智能算法(ABR算法)进行三维地质属性体插值，提高了插值精度，实现自动化插值。本文的研究对于推进地质建模的智能化发展具有一定的意义。

1 技术路线

利用sklearn算法库中的MDS和AdaBoostRegressor(ABR)进行插值。插值方法跟大多数插值方法类似，采用距离未知点最近的n个已知样本点进行未知点的估算。三维空间插值方法概述流程如图1所示，其中P0为未知点P1,P2,…,Pn为距离P0最近的n个已知样本点，其顺序以由近到远排序。

图1 利用MDS、ABR和MF插值流程概述Fig.1 Summary of MDS,ABR and MF interpolation process

2 MDS降维

空间插值中未知点的属性值受其最近邻的若干个已知样本点的距离、已知点的属性值、拓扑关系影响，然而，一些变量对插值结果影响很小。为了提出影响插值的主要因素，减少次要因素的影响，提高插值结果的稳定性，将这些因素作为MDS的降维对象，得到一维数据。

采用未知点(P0)最近邻的n个样本点(P1,P2,P3,…,Pn)来预测未知点的属性值。假设每个已知样本点为未知点，然后根据这些已知样本点建立ABR模型。例如，假设P1为未知点，P2,P3,…,Pn为已知点(共n-1个)。把这些n-1个已知点根据到P1的距离由近到远重新进行排序，排序后的点为Psort1,Psort2,…,Psort(n-1)。对n个已知样本点(P1,P2,P3,…,Pn)依次假设每个样本点为未知点。每次假设就能得到一个训练样本。共产生n个训练样本(每个训练样本的特征有n-1个)，把这些训练样本带入到ABR模型中。当假设P1为未知点时，利用MDS降维的描述如图2所示。

(1)Z表示属性值。Zsort1表示Psort1的属性值，Psort1为距离P1最近的点；Zsort2表示Psort2的属性值，Psort2为距离P1第2近的点；以此类推，Zsort(n-1)表示Psort(n-1)的属性值，Psort(n-1)为距离P1第n-1近的点(距离P1最远的点)。

(2)h表示两个点之间的距离。h1(sort1)表示P1～Psort1的距离；h1(sort2)表示P1～Psort2的距离；以此类推，h1[sort(n-1)]表示P1～Psort(n-1)的距离。

(3) 由于地理现象中的属性特征在不同的方向上呈现出不同的变化速度，所以本文用方位角ɑ表示各向异性特征。α1(sort1)表示P1与Psort1的方位角；α1(sort2)表示P1与Psort2的方位角；以此类推，α1[sort(n-1)]表示P1与Psort(n-1)的方位角。

(4)X1(P1Psort1)表示P1与Psort1的h、Z、α进行降维后的数据。

(5) ABR的训练集中包含n个训练样本，每个训练样本包含n-1个特征。特征根据已知样本点与未知样本点的距离远近获取。特征1表示距离假设的未知点P1最近的样本点Psort1对假设未知点P1的影响，用X1(P1Psort1)表示。特征2表示距离假设的未知点P1第二近的样本点Psort2对假设未知点P1的影响，用X1(P1Psort2)表示。以此类推，共产生n-1个特征。

在图2中需要注意的是，在对所有样本点假设为未知点时，每次假设，在计算特征1时都会产生一个Psort1，这些Psort1并不是代表同一个样本点，而是对假设未知点后，其余的已知点到未知点的距离排序由近到远而产生的。例如，假设P1为未知点时，Psort1表示其余已知点中距离P1最近的点。假设P2为未知点时，Psort1表示其余已知点中距离P2最近的点。

图2 n个样本点MDS降维处理Fig.2 MDS dimension reduction for n sample points

3 ABR空间插值

AdaBoost(adaptive boosting)[16]算法是在多个弱学习模型的基础上建立起来的，把样本带入弱学习模型中进行训练及验证，根据误差对样本权重进行更新，然后带入到新的弱学习模型中训练，最后对这些弱学习模型的结果进行加权求和(或者投票)的方式组合，求得模型的结果。AdaBoost在每次训练弱学习模型时，都对样本进行了验证以及误差计算，所以样本进行了多次验证及误差分析。常用的分类回归算法一般都需要用到训练集及验证集，以此来更新模型权重获取最优模型。AdaBoost算法在计算过程中已经把训练样本当作验证集多次进行验证来获取最优模型，所以额外的验证集并不是必须的。

AdaBoost可以运用于分类(adaptive boosting classifier)和回归(adaptive boosting regressor)，本文中主要运用其回归特性，具体推导过程及应用见文献[16]。

在利用ABR建立模型时，由于样本量较少，很容易产生较好的训练模型，选用的弱学习模型为CART决策树[17]，损失模型为线性模型，最大迭代次数为50，学习率为0.1。

3.1 ABR模型训练

采用ABR方法建立插值模型，一个未知点建立一个ABR模型。若采用n邻域样本点进行插值则训练ABR模型所需的数据集个数为n，训练集的输入中每个样本的特征为n-1个。ABR模型训练时需要确定训练集的输入和输出(图3)。

(1) ABR训练集输入：ABR训练集的输入为经过MDS降维到一维的数据集(图2中降维后的特征，也就是图3中的输入)。

(2) ABR训练集输出：ABR训练集的输出为图2中假设的未知点的真实值(图3中的输出)。

图3 训练ABR插值模型Fig.3 Training ABR interpolation model

所有的训练集的特征之间可能有一定的关联，因为它们都是从P1,P2,…,Pn获取的。然而每个训练样本所表达的意义以及获取的计算方式都是独立的(图2中的特征)。因此，尽管它们具有一定的关联，但是仍然可以用于建立ABR模型。

3.2 未知点预测

通过上述训练得到ABR模型，然后获取未知点P0与其最近的n-1个点之间的相对距离h，方位角α，样本点的属性值Z。P1,P2,…,Pn-1是距离P0最近的n-1个样本点，并且这些样本点是根据到P0的距离，由近到远排序的。如图4所示。

图4 未知点值的预测Fig.4 Prediction of unknown point

(1)Z1表示已知样本点P1的真实值。

(2)h01为未知点P0与已知样本点P1的距离。

(3)α01为未知点P0与已知样本点P1的方位角,用来表示各向异性特征。

(4)X0(P0P1)为未知点P0与已知样本点P1的Z1、h01、α01进行MDS降维后的结果，同理得到X0(P0P2),…,X0(P0Pn-1)。

(5) 把X0(P0P1),X0(P0P2),…,X0(P0Pn-1)代入ABR模型中预测出待估计点P0的属性值Z0。

(6) 特征1, 特征2,…, 特征n-1表示已知样本点距未知点P0由近到远的影响特征。

4 三维MF平滑处理

均值滤波(MF)是一种邻域运算，邻域运算是通过对某一像素的邻域进行加权运算求得新的值代替该像素，加权运算所用到的矩阵叫作局部算子。本文的三维MF采用3×3×3的局部算子，具体计算流程如下。

(1) 设定局部算子为3×3×3，同时在X、Y、Z方向分别扩展一行(列)，新增行的像素值赋值为已知的最近像素值。

(2) 以待处理像素为中心，选取3×3×3的矩阵。

(3) 计算M矩阵中3×3×3个矩阵的均值，并替换掉待处理像素。

(4) 依次遍历其他像素，重复步骤(2)和步骤(3)。

图5中是ABR进行插值的结果(MF处理前)与MF处理后效果对比，可以看出经过MF处理后可以减少多个ABR模型计算的插值结果产生的不连续性，使插值结果符合地质属性分布具有连续性的特征。

图5 均值滤波前后对比Fig.5 Comparison result before and after MF

5 插值结果分析

孔隙度是指岩石孔隙体积与岩石表面体积之比，是评价油气分布的重要指标。以辽河油田某区块孔隙度数据为例，进行了三维空间插值。由于设备或其他原因，某些孔隙度数据可能无效。图6显示了11口井的分布，其中验证样品(3口井，共76个样品)用红色矩形标记。其他9口井用于插值(总共203个样本)。由于钻井位置与含油量、孔隙度、渗透率、饱和度等因素有关，无法保证井的均匀分布。在图6中，左侧的孔隙度较大，这有利于储油、高含油量和更多的井。其他地区孔隙度小，含油量低，井数较少。插值结果的步长在X方向为5 m，Y方向为5 m，Z方向为1 m(每个最小单位网格的插值结果为5×5×1=25 m3)，总共有27 000个网格。工区面积为150×150 m2，地层厚度为30 m(海拔为-1 640～-1 670 m)。

图6 测试样本分布Fig.6 Test sample points distribution

由于反距离加权插值没有考虑地质发育的各向异性特征；Voronoi图插值结果数据结构复杂，可重用性较弱；RBF插值结果过于平滑，只适用于固态矿产的地质体建模，而地质统计学中的克里金方法在含油气层的地质体建模中最为常用。所以将本文方法与地统计学插值方法中的普通克里金(ordinary Kriging，OK)的插值结果进行比较，分析插值效果。OK使用Python中的pykrige包。

当存在多个样本点时，并非所有的样本点都将用于插值。通常选取未知点附近的n个已知样本点建立插值模型，用于计算未知点的值。在本文中，搜索距离未知孔隙度点(图1～图4中的P0)最近的n个已知孔隙度样本点(图1～图4中的P1,P2,…,Pn)。图7显示了MDS-ABR-MF(由MDS、ABR和MF处理的插值结果)和OK插值结果之间的比较，其中n=5、n=10和n=15的插值效果。

图7 不同方法插值效果对比Fig.7 Comparison of interpolation results by different methods

用3口测试井中的76个样本点进行误差分析。通过均方根误差(root mean squared error, RMSE)和平均绝对误差(mean absolute error, MAE)来测试本文方法的插值精度。由表1可以看出MDS-ABR-MF运用本文中的测试数据进行插值比OK插值结果在精度有所提升。

表1 MDS-ABR-MF和OK插值结果与3口测试井之间的RMSE和MAETable 1 The RMSE and MAE between the MDS-ABR-MF and OK interpolation results and 3 test wells

OK与MDS-ABR-MF插值结果的误差产生的原因可能如下。

(1) OK方法假设在标准正态分布的情况下进行插值。然而，无论原始数据如何处理，这些数据都只是近似正态分布，很难达到标准正态分布。

(2) 在各向异性分析和半变异函数分析中，OK插值结果受参数选择的影响。例如，孔隙度范围为0～1，通过测试OK方法，不合理的人为调整可能会产生负值，因为OK插值中计算的权重可能是负值。然而，在ABR插值中，由于弱学习模型为决策树，算法的基本原理决定了未知点的值是已知点的加权计算，所以当ABR算法中决策树为弱学习模型时，决定了权重不能为负，所以最终结果不可能为负值。

(3)图7中，OK插值中的不连续性是因为选取的n个已知样本点的孔隙度值可能会相同。此时，在OK算法中，计算的协方差矩阵无法获得逆矩阵，因此本文取平均孔隙度值，这也是OK算法的一个弊端。

6 结论

利用MDS、ABR和MF进行探索三维空间样本点之间的相关性来实现三维地质属性体插值。摆脱了反距离加权方法中不考虑各向异性特征的弊端。通过对比传统的普通克里金(OK)插值，得出以下结论。

(1) OK插值需要特定的随机场(正态分布)，在非正态分布区域执行OK插值时，插值结果精度可能不高。MDS-ABR-MF方法不受特定随机场的约束。

(2) OK插值需要人工半变异函数分析、各向异性分析、块金、基台和变程等繁琐流程分析，不适当的参数调整也可能产生异常。MDS-ABR-MF插值可以缩短繁琐的参数调整过程，在保证精度的同时实现自动插值。