针对地质云钻孔数据的空间插值方法选择*

张 灏，王 娇,2，郑新奇

(1. 中国地质大学(北京)信息工程学院，北京 100083；2. 中国科学院地理科学与资源研究所 资源与环境信息系统国家重点实验室，北京 100101)

9月20日习近平主席在郑州座谈会中首次将黄河流域生态保护和高质量发展上升为重大国家战略。但现今黄河流域中很多省份的矿区因采矿活动造成采空塌陷、地下水疏干、地质地貌景观破坏等问题，严重危害矿区人民正常生产生活。为了正确模拟和反映矿区地下岩层的分布形态，指导矿区可持续性发展[1]，陈健峰等学者基于地质云钻孔数据利用Civil 3d软件进行建模分析[2]。然而由离散地质云钻孔数据所创建的曲面很难真实表达岩层分布状况[3]，需要确定精确的空间插值方法。

鉴于此本研究针对黄河流域内某矿区地质云钻孔数据，结合GIS空间分析地统计学方法和相关软件，分析出数据特征，利用克里金插值无偏、最优的特性[4]选取多模型参数进行插值分析，得到最优插值方式，以辅助相关研究分析和决策。

1 研究区概况与数据获取

1.1 研究区概况

选取矿区位于宁夏省吴忠市，地处宁夏回族自治区中部，地貌形态为山地、丘陵、黄河冲积平原和库区等，黄河干流流经吴忠69 km，山区海拔1 300～1 900 m。境内矿产资源丰富，主要有石油、煤炭、矿石、天然气等30多种矿产资源，随着大规模的矿井开采，矿区污染状况日益严重，地质环境保护工作成为当地所面临的重要问题。

1.2 数据获取

选取矿井样本点数据共80个，来源于全国重要地质钻孔数据库服务平台[5]提供的实地钻孔报表。研究区域为石记场矿区，每一份钻孔报表中含有约10种岩层的数据，选取其中一层石膏岩数据作为分析对象。用每个矿井的高程减去对应石膏岩的换层深度，得到这些岩层的高程，整理成excel表格数据。ArcMap10.2提供了添加x，y数据工具，可以将excel数据转换为属性表。全国重要地质钻孔数据库服务平台记录数据经纬度采用WGS1984坐标系，使用ArcMap10.2定义投影工具给数据定义相应的坐标系[6]，得到钻孔数据空间分布如图1所示。

图1 矿区采样点分布

2 矿区钻孔数据空间特征分析

2.1 地质统计学基础

克里金法(Kriging)是地统计学的主要内容之一，根据待插值点与临近实测高程点的空间位置, 对待插值点的高程值进行线性无偏最优估计, 通过生成一个关于高程的克里格插值图来表达研究区域的原始地形。总的公式是:

式中,Z(x0)表示未知样点的值，Z(xi)表示未知样点周围的已知样本点的值，n为已知样本点的个数;λi为第i个样本点的权重。它的确定是通过半方差图分析获取的, 根据统计学上无偏和最优的要求, 利用拉格朗日极小化原理, 可推导权重值和半方差之间的公式[7]。

2.2 空间数据分析

在获得样本点数据后, 首先需要对数据进行分析, 检验数据的分布, 分析数据的分布趋势等。ArcMap10.2地统计向导提供了直方图法、QQ图法、趋势分析和半变异/协方差函数云等工具来检查分析数据[8]。

2.2.1 矿井岩层数据分布检验

ArcMap10.2 地统计学向导提供了两种验证方法：直方图法和正态Q-Q图法。通过地统计分析工具对生成80个样本点进行可视化分析，比较Log变化前后正态分布情况。发现数据在进行Log变换后更符合正态分布，从直方图和Q-Q图结果中可以看到样本点基本沿正态直线分布。

图2 矿区地质云数据直方图分布图

图3 矿区地质云数据Q-Q图分布图

2.2.2 趋势效应

运用趋势分布工具对数据进行分析如图4所示，结果显示数据在X方向呈二阶分布，在Y方向上呈轻微的一阶分布，总体分布趋势需通过试验验证。

图4 矿区地质云数据趋势分布图

2.2.3 空间自相关

运用半变异函数协方差云图对数据进行探索性分析，从半变异函数结果看数据集中在左下方，符合区域化变量距离越近，相似性越大的特点，表明数据可用于克里金插值。

图5 半变异函数云图

3 空间插值

3.1 数据插值结果验证方法

由于训练样本数据不足，采用交叉验证方式评价[9]，参数选取平均误差(ME)、均方根误差(RMSE)、平均标准误差(ASE)和标准化均方根误差(RMSSE)评估插值结果。判断插值精度可按以下标准综合进行：平均误差和标准化平均误差的绝对值最接近于0；均方根误差越小越好；平均标准误差与均方根误差最接近[6]。

3.2 岩层高程分布趋势剔除

为比较数据与一阶和二阶趋势面符合情况，采取普通克里金、简单克里金和泛克里金三种插值方法，其余变量参数统一设置为步长12，高斯变异函数模型。得到各差值误差结果如下，从结果看一阶预测误差值较小。表明数据区域总体上呈现一阶分布，因此之后的分析模型都以一阶趋势面为基础。

表1 趋势面阶数误差结果

3.3 克里金不同插值方法插值比较

3.3.1 变异函数参数选择

不同半方差函数模型和步长会对插值结果产生影响[6]，为了获取最优的模型参数，本文进行多模型多搜索半径的插值比较。在对数据进行Log变换，一阶趋势面剔除和不考虑各向异性基础上，选取应用广泛的高斯、球面和指数模型[8]；设置初始搜索半径步长值为6，终止步长值为18，间隔为3，插值结果如图6所示。

图6 不同模型参数空间插值结果

3.3.2 普通克里金法插值结果分析

普通克里金法要求区域化变量满足二阶平稳或本征假设，至少满足准二阶平稳或准本征的假设。此条件下，在估计邻域内要求期望值为常数。插值结果如表格所示，从结果中分析，指数函数模型相对误差较大，选取高斯球面模型的平均误差ME和标准化平均误差MSE的绝对值最接近于0,球面函数模型的平均标准误差ASE与均方根误差RMSSE最接近。相比之下高斯球面模型各误差值都较低，故选取它作为最优参数。

表2 普通克里金误差统计

为确定最优步长数，单独对高斯模型不同步长数结果进行展开。要求RMSSE值越小越好，依曲线图所示，选取高斯函数模型步长数为9时RMSSE值最低，故确定步长数9为最优参数。

3.3.3 泛克里金法插值结果分析

泛克里金法，就是在漂移的形式和非平稳随机函数的协方差函数或变异函已知的条件下，一种考虑有漂移的无偏线性估计量的地统计学方法[10]，本

图7 RMSSE误差曲线1

区域数据满足正态分布，符合泛克里金插值条件。通过泛克里金法进行插值得到误差统计结果如表3所示。

表3 泛克里金插值误差统计

上述通过逐一比较误差进行精度评价的方式复杂且含有不确定性因素，因此定义3个评价指数s，t，h，计算公式为：

(1)

(2)

(3)

其中,n为不同步长数的个数，本文统一为6，RME为平均误差，RRMSE为均方根误差，RASE为平均标准误差，s，t，h三个值越小代表插值精度越高。根据公式(1)(2)(3)计算新的评价指数如表4所示。

表4 泛克里金误差评价指数

从误差结果看高斯函数各误差值最小，故选取高斯函数为最优模型，对其RMSSE误差值随步长数增加的情况进行可视化如图，从结果来看选取步长数为9时为最优。

图8 RMSSE误差曲线2

3.3.4 简单克里金法插值分析

简单克里金法可以使用半变异函数或协方差和变换，并且允许测量误差[10]。使用该方法进行插值并计算s，t，h评价指数的值如下表5所示。

表5 简单克里金误差评价指数

从结果看球面函数各项指标值较低，误差最小，故选取球面函数为最优模型。RMSSE基于步长变化的误差曲线如图，选取步长数为6时误差相对最小。

图9 RMSSE误差曲线3

3.3.5 三种方法插值结果比较

从插值结果来看，选取高斯函数模型，步长数为9时三种克里金插值方法交叉验证误差都较小。故以此为统一模型参数比较不同插值的误差大小如表6所示。

表6 高斯模型步长数9误差统计

从表中看选择普通克里金方法插值效果最好，RMSE、MSE和ME误差值都为最低。故对于本文研究区域选取普通克里金方法，参数选取高斯球面函数，步长数9为最优。

3.4 三维可视化展示

在ArcMap10.2中导出插值栅格地图，利用ArcScene软件高度拉伸工具[6]对插值结果进行三维可视化结果如图。观察区域的分布特征，整个矿区石膏岩层高程从西向东降低，在矿区中部岩层高程发生剧烈的变化，西北和东南地区高程较高。

图8 普通克里金插值三维可视化

4 结 论

本研究采用克吕格插值多模型插值方法对80个地质云钻孔数据进行空间插值，从趋势面分析和预测误差大小来看，克里格插值方法的一阶趋势预测效果要好于二阶趋，表明该区域数据总体上更接近一阶分布趋势。从多模型参数交叉验证结果来看，采用普通克里格方法插值，高斯函数，步长9进行空间插值精度最高，能最准确地反映该矿区石膏岩地层的空间分布形态。不足之处是该区域地质云钻孔数据数量有限，如果增加样本点个数进行进一步误差检验得到更精确的误差分析结果。