空间插值在测氡探测隐蔽火源中的应用比较

Comparison of Spatial Interpolations Applied to Radon Measurement for Hidden Fire Location

摘要

为了选择适用于局部观测尺度下测氡探测隐蔽火源的最优空间插值方法， 本文以北岩煤矿火区样本数据为研究对象， 基于Arc GIS地统计分析模块， 通过交叉检验定量筛选出4种插值方法下误差最小的插值模型， 定性比较了不同模型的空间分布特征， 确定了最优插值方法并进行了钻孔验证。结果表明， 4种方法中误差最小的插值模型分别为反距离权重法选择模型幂指数p=1， 径向基函数法为薄板样条函数模型， 泛克里金法与普通克里金法均为高斯函数模型； 综合空间分布图内插、 外推表面分布特征及钻孔验证结果， 插值方法排序为普通克里金法（高斯函数）＞泛克里金法（高斯函数）＞径向基函数法（薄板样条函数）＞反距离权重法（p=1）， 优先选择普通克里金法（高斯函数）进行插值来获取煤自燃火区地表氡分布差异信息。

Abstract

In order to select the optimal spatial interpolation applied to radon measurement for hidden fire location at the local scale， it took the sample data from the fire area of Beiyan Coal Mine as the study object based on the Arc GIS geostatistic wizard， then quantitatively screened out the most accurate interpolation model among four interpolations by cross-checking， qualitatively discerned the spatial distribution characteristics of different models， and ascertained the optimal interpolation by drilling verification. The results show that the most accurate model among four interpolations： the Inverse Distance Weight selects the Power p=1 model， the Radial Basis Function uses the Thin Plate Spline model， and both the Universal Kriging and the Ordinary Kriging are a Gaussian model； with the comprehensive deliberation on the distribution characteristics both in interpolating and extrapolating surface， as well as drilling verification， the interpolations are ordered as Ordinary Kriging（Gaussian）>Universal Kriging（Gaussian）>Radial Basis Function （Thin Plate Spline）>Inverse Distance Weighting（p=1）， the Ordinary Kriging（Gaussian） is preferred to obtaining surface radon distribution difference in coal spontaneous combustion fire areas.

Graphical abstract

关键词

隐蔽火源探测 / 空间插值 / 氡分布 / 普通克里金

Key words

hidden fire detection / spatial interpolation / radon distribution / ordinary Kriging

0 引　言

非接触隐蔽火源探测技术是煤矿安全领域应重点突破的关键技术之一， 大量灭火工程实践表明测氡法有望为隐蔽火源精准探测提供可行的解决途径。该方法通过获取煤自燃火区地表氡分布差异信息， 从而圈定隐蔽火源位置、 范围及发展趋势^［1］。实际探测过程中， 氡分布差异信息的获取主要分为两个阶段： 1） 采用α杯法获得等间距的氡计数散点值； 2） 将散点值转换为连续的空间分布图。但基于有限数量的氡计数， 无法生成高精度的空间分布图， 需要使用空间插值方法以提高散点数据密度^［2‐4］。

目前， 插值方法可分为确定性插值法和地统计插值法， 前者常用的有反距离权重法（Inverse Distance Weight， IDW）、 径向基函数法（Radial Basis Function， RBF）， 后者常用的有普通克里金法（Ordinary Kriging， OK）、 泛克里金法（Universal Kriging， UK）。针对不同研究对象， 所适用的最优插值方法有所不同^［5-8］。李华等^［9］对比了OK和IDW在核辐射剂量场分布的插值效果， 认为多数情况下OK优于IDW， 但在数据极少的区域OK的误差较大。解恒燕等^［10］认为IDW在获得美国Upper Sangamon流域月降水量数据时插值精度要优于OK与RBF。王长虹等^［11］对比了OK与UK在粉质黏土层黏聚力和内摩擦角等参数的预测精度， 认为UK优于OK。姜晓剑等^［12］认为RBF在对559个气象站点的逐日基本气象要素插值时精度要高于OK与IDW。

此外， 即使研究对象相同， 不同观测尺度下所适用的最优插值方法也有所不同。何 等^［13］对市区日平均降水量插值时， 认为反距离权重法插值精度高于普通克里金法， 但在自治区级尺度下普通克里金法更为理想。Chen等^［14］将观测尺度扩大至全国， 发现普通克里金法插值效果最优。

不合理的插值方法会降低空间分布图的精度， 进而影响对火源位置的判断^［15-16］。测氡法以局部观测尺度下的地表氡计数为研究对象， 目前尚缺乏对其所适配最优插值方法的应用比较。本文以北岩煤矿火区样本数据为研究对象， 基于Arc GIS地统计分析模块进行IDW、 RBF、 OK及UK插值， 通过交叉检验定量筛选4种方法下误差最小的插值模型， 定性比较不同模型的空间分布特征， 确定最优插值方法并对其分布图所呈现的火源位置、 外推表面特征进行钻孔验证， 以期为精准定位隐蔽火源提供依据。

1 研究区域与研究方法

1.1 研究区域概况与数据采集

研究区域位置及测点分布如图 1 所示， 测场方位为NS， 起点为西翼A17^#点， 向东延伸270 m至R17^#点， 向南最远延伸255 m。

测场西南及东南方向存在基岩裸露区域， 无法布置采样点， 整体呈倒“凸”形网格状分布。测点布置于网格中心， 间距为15 m×15 m， 探测区域面积为52 650 m²。

1.2 空间插值方法

1.2.1 反距离权重法

IDW以待插值点与若干邻近样本之间距离的倒数为权重进行加权平均， 其计算公式为

Z={\displaystyle \sum _{i=\mathrm{1}}^n}\frac{\mathrm{1}/d_i^p}{{\displaystyle \sum _{i=\mathrm{1}}^n}\left(\mathrm{1}/d_i^p\right)}{Z}_i,（1）

式中： Z为某待插点的值； Z_i 为第i（i=1，2，…，n）个已知采样点的值； d_i 为第i个已知采样点与待插点的距离； p为幂指数， 当p一定时， 随d_i 的增加， 采样点所占权重呈指数递减。

1.2.2 径向基函数法

RBF是用不同的基函数来完成逼近实值函数F=F（x）某点x函数值的方法， 其核心是构造一个具有下述形式的逼近函数S（x）， 其计算公式为

S\left(x\right)={\displaystyle \sum _{i=\mathrm{1}}^n}{\alpha }_i\phi \left(‖x-x_i‖\right)，x\in R^d,（2）

式中： \phi \left(t\right)为基函数（t≥0）； ‖x-x_i‖为待估计点x与基函数的逼近节点x_i （i=1，2，…，n）之间的欧氏距离； α_i 为第i个基函数对应的待定系数。

1.2.3 普通克里金法与泛克里金法

OK与UK通过变异函数分析样本数据反映的区域化变量的结构信息， 进而根据待插值点与有限邻域内样本点的空间位置关系、 样本点间的相互位置关系， 对待插值点进行无偏最优估计， 其计算公式为

Z^{*}\left(x\right)={\displaystyle \sum _{i=\mathrm{1}}^n}{\lambda }_{i}Z\left(x_i\right),（3）

式中： {\lambda }_i为权重系数； Z^{*}\left(x\right)为待插值点的估计值； Z\left(x_i\right)为邻近样本点的已知值。

两种方法均以变异函数\gamma \left(h\right)为基础构建插值模型， \gamma \left(h\right)的计算公式为

\gamma \left(h\right)=\frac{\mathrm{1}}{\mathrm{2}N\left(\mathit{h}\right)}{\displaystyle \sum _{i=\mathrm{1}}^{N(\mathit{h})}}{\left[Z\left(x_i\right)-Z\left(x_i+\mathit{h}\right)\right]}^{\mathrm{2}},（4）

式中： h 为两样本点间的向量； Z（x_i ）， Z（x_i +h）为x_i 和x_i+h处的已知测量值（i=1，2，…，N\left(h\right)）； N（ h ）为相距h的样本点对个数。

两种方法的不同之处在于有限邻域内变量是否平稳。OK适用于满足二阶平稳假设的变量， 其数学期望为未知常数m。UK则为非平稳变量， 数学期望为多个单项式的线性组合m\left(x\right)。

1.3 定量检验方法

采用交叉检验定量比较不同方法的插值精度， 其原理为首先移除单个样本点的氡测量值， 然后使用插值方法获得同一样本点的预测值， 依次计算各已知样本点处实际测量值与预测值之间的误差。利用平均误差（Mean Error， ME）、 均方根误差（Root Mean Square Error， RMSE）、 平均值标准化误差（Mean Standardized Error， MSE）、 标准化均方根误差（Root Mean Square Standardized Error， RMSSE）作为交叉检验指标， 计算公式为

ME=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i=\mathrm{1}}^N}\left[Z^{*}\left(x_i\right)-Z\left(x_i\right)\right],（5）

RMSE=\sqrt[]{\frac{\mathrm{1}}{N}\bullet {{\displaystyle \sum _{i=\mathrm{1}}^N}\left[Z^{*}\left(x_i\right)-Z\left(x_i\right)\right]}^{\mathrm{2}}},（6）

MSE=\frac{{\displaystyle \sum _{i=\mathrm{1}}^N}\left\{\left[Z^{*}\left(x_i\right)-Z\left(x_i\right)\right]/{\sigma }^{*}\left(x_i\right)\right\}}{N},（7）

RMSSE=\sqrt[]{\frac{{{\displaystyle \sum _{i=\mathrm{1}}^N}\left\{\left[Z^{\mathrm{*}}\left(x_i\right)-Z\left(x_i\right)\right]/{\sigma }^{\mathrm{*}}\left(x_i\right)\right\}}^{\mathrm{2}}}{N}},（8）

式中： N为已知测点的数量； {\sigma }^{\mathrm{*}}\left(x_i\right)为估计值的标准差。

交叉检验指标的评价标准为： 1） ME的绝对值接近于0； 2） MSE接近于0； 3） RMSE越小越好； 4） RMSSE应接近于1。

2 结果与讨论

2.1 地统计插值法的交叉检验

通过有限的样本数据来拟合实验变异函数曲线{\gamma }^{\mathrm{*}}\left(h\right)， 进而推断出理论变异函数， 并配以相应的理论模型。所以， 最优插值法的实质是最优理论模型的选择过程。

分别利用球面、 指数、 高斯和三角函数等4种常用的理论模型进行OK和UK插值， 其交叉检验结果如图 2（a）、 图 2（b） 所示。

OK中， 高斯函数模型的|ME|以及MSE最近接近于0， RMSE最小， RMSSE最接近于1。因此， 高斯函数模型可作为OK的最优理论模型。

UK中， 高斯函数模型的|ME|以及MSE最近接近于0， RMSE最小， RMSSE较接近1。因此， 选择高斯函数模型作为UK的最优理论模型。

2.2 确定性插值法的交叉检验

IDW中， 幂指数p能够控制权重值的变化速度， 所以p的选择会影响插值精度。通常p＞5时， 只有紧邻待插值点的个别样本点会对插值结果产生影响， 插值失去意义， 所以本文中p的取值介于1~5之间， 交叉检验的结果如图 2（c） 所示。

RMSE随p的增大而减小， 其中最大值RMSE _p=1仅是最小值RMSE _p=5的1.25倍， 说明RMSE相差不大； |ME|随p的增大而增大， 其中最大值|ME| _p=5是最小值|ME| _p=1的3.96倍， 所以p=1为最优权重。

RBF利用最小化表面总曲率的基函数来估值， 基函数的选择决定了插值表面通过所有已知样本点的方式， 进而影响插值精度。Arc GIS地统计分析模块提供了规则样条函数、 张力样条函数、 高次曲面函数、 反高次曲面函数和薄板样条函数等5种常用基函数， RBF交叉检验的结果如图 2（d） 所示。

RMSE中最大的薄板样条函数仅是最小规则样条函数的1.18倍， 说明各基函数的RMSE相差不大， 而薄板样条函数的|ME|明显最接近于0， 所以薄板样条函数为最优基函数。

2.3 定性分析与讨论

通过交叉检验定量筛选出4种插值方法下的最优模型， 其空间分布图如图 3 所示。此外还需要结合以下实际情况对其空间分布图进行定性分析：

1） 氡测值在火区与非火区间变化剧烈， 最优插值方法应能适用于较大幅度的数值波动；

2） 已被划定为疑似火区的范围内存在无法布置采样点的区域（Ⅰ、 Ⅱ表示基岩裸露区）， 最优插值方法应有合理的外推精度。

4种模型所呈现的空间分布均显示了3处火区A、 B、 C， 但不同模型在表面平滑程度、 局部精度以及外推区域分布特点上存在一定差异。

图 3（a） 中， IDW的插值表面较为粗糙， 其中火区C内的“牛眼”影响对高温火源位置的判断。此外， 3处火区均出现相邻氡计数区间边界呈“齿轮”状交叉的情况， 不符合煤自燃温度演化规律。这是由于IDW的计算方式， 无法充分考虑样本点的空间分布。当待插值点附近的样本点数据很大或很小时， 待插值点的结果容易受极值的影响而产生明显的“牛眼”现象。图 3（b） 中， RBF的插值表面数值波动较大， 3处火区的周边及外推表面Ⅱ均存在氡计数为负值的系统误差区域。这是由于作为一种确定性插值方法， 插值表面必须通过每一个已知样点， 同时又要保证表面的总曲率最小。

火区范围内已知样点的氡测值较大， 插值表面通过时向上凸起的程度加大， 总曲率增大； 相邻的非火区内， 已知样点的氡测值较小， 为降低总曲率， 在已知样点邻近区域内形成氡计数为负值的插值点集， 插值表面向下凹陷拟合， 总曲率减小， 形成“火区高估， 非火区低估”的系统误差区域。

图 3（c） 中， UK的插值表面整体较为平滑， 但其外推部分存在异常分布区域。表面Ⅰ内存在呈半椭圆形的无效数据点集（阴影部分）； 表面Ⅱ内除1处无效数据点集外， 整个外推区域的氡计数远小于本底值， 与邻近的内插区域割裂， 显然和实际情况不符。这是由于测场内探测数据集呈现“整体漂移， 局部平稳”的趋势， 即氡计数由非火区向火区A、 B、 C逐渐增大，称为“整体漂移”， 但单个测点所在15 m×15 m的范围内氡计数波动较小， 测量值可代表局部均值， 称为“局部平稳”。这与泛克里金法所适用的“整体平稳， 局部漂移”相反。

图 3（d） 中， OK能够规避上述3种方法存在的问题， 插值表面光滑且能够明显展示高温火源位置， 充分考虑样本点的空间分布， 外推区域估计值分布在合理区间， 所以高斯函数模型的OK是最优插值方法。

2.4 钻孔验证

由地面施工钻孔至煤层位置， 利用钻孔温度测量仪测量孔内最高温度， 进而对定性分析结果进行验证， 结果如表 1 所示。

图 3（a） 中易将火区C内的“牛眼”判定为火源位置， 但C₁、 C₂并不处于“牛眼”区域， 其钻孔验证结果却仍显示明显的火源特征， 说明“牛眼”会干扰火源位置的判断。

钻孔Ⅱ₁位于图3（b）中氡计数为负值的系统误差区域； 钻孔Ⅰ₁、 Ⅱ₁位于图 3（c） 中的无效数据区域， 均与钻孔验证结果明显不符。

图 3（d） 中所呈现的火源位置信息、 空间分布特征， 与各钻孔验证结果基本一致。

3 结　论

1） 交叉检验结果表明4种方法中误差最小的插值模型分别为： 反距离权重法（幂指数p=1）， |ME|为0.210 1， RMSE为128.270 5； 径向基函数法（薄板样条函数）， |ME|为0.119 9， RMSE为151.279 5； 泛克里金法（高斯函数）， |ME|为8.5471， RMSE为165.9394； 普通克里金法（高斯函数）， |ME|为8.2240， RMSE为162.1121。

2） 综合空间分布图内插、 外推表面分布特征及钻孔验证结果， 插值方法排序为： 普通克里金法（高斯函数）＞泛克里金法（高斯函数）＞径向基函数法（薄板样条函数）＞反距离权重法（幂指数p=1）。因此， 优先选择普通克里金法（高斯函数）进行插值来获取煤自燃火区地表氡分布差异信息。

参考文献

[1] 王俊峰.煤地下自燃时覆岩中氡气运移规律及应用研究［D］.太原： 太原理工大学， 2010.

[2] QI X Y， CHEN L Z， ZHANG L B， et al. In situ FTIR study on real-time changes of active groups during lignite reaction under low oxygen concentration conditions［J］. Journal of the Energy Institute， 2019， 92（5）： 1557-1566.

[3] 周斌， 董智宇， 王俊峰， 等.氧化升温过程中煤岩介质体氡析出特性实验研究［J］.煤炭学报， 2020， 45（S2）： 859-866.

[4] ZHOU Bin， DONG Zhiyu， WANG Junfeng， et al. Experimental study on radon exhalation characteristics of coal and rock during oxidation and heating［J］.Journal of China Coal Society， 2020， 45（S2）： 859-866. （in Chinese）

[5] 朱政江， 张雪英.基于BP和RBF神经网络的煤与瓦斯突出模型研究［J］.中北大学学报（自然科学版）， 2014， 35（5）： 571-575.

[6] ZHU Zhengjiang， ZHANG Xueying. Research on coal and gas bursting models based on BP and RBF neural networks［J］.Journal of North University of China （Natural Science Edition）， 2014， 35（5）： 571-575.（in Chinese）

[7] MARIANO C， MÓNICA B. A random forest-based algorithm for data-intensive spatial interpolation in crop yield mapping［J］. Computers and Electronics in Agriculture， 2021， 184： 106094.

[8] YAO X L， FU B J， LÜ Y H， et al.Comparison of four spatial interpolation methods for estimating soil moisture in a complex terrain catchment［J］.PLoS One， 2013， 8（1）： e54660.

[9] CHOI H， KIM H， YEOM S， et al. An indoor environmental quality distribution map based on spatial interpolation methods［J］. Building and Environment， 2022， 213： 108880.

[10] 智协飞， 吴柏莹， 罗忠红， 等.华东地区地面和高空风场的多模式集成精细化预报研究［J］. 大气科学学报， 2023， 46（6）： 917-927.

[11] ZHI Xiefei， WU Baiying， LUO Zhonghong， et al. Research on multi-model integrated refined forecasting of surface and upper-altitude wind fields in east China［J］. Transactions of Atmospheric Sciences， 2023： 46（6）： 917-927. （in Chinese）

[12] 李华， 赵原， 刘立业， 等. 基于克里金理论的辐射剂量场插值重构与可视化研究［J］. 辐射防护， 2019， 39（6）： 475-482.

[13] LI Hua， ZHAO Yuan， LIU Liye， et al. Research on Interpolation reconstruction and visualization of radiation does field based on kriging theory［J］. Radiation Protection， 2019， 39（6）： 475-482.（in Chinese）

[14] 解恒燕， 张深远， 侯善策， 等. 降水量空间插值方法在小样本区域的比较研究［J］.水土保持研究， 2018， 25（3）： 117-121.

[15] XIE Hengyan， ZHANG Shenyuan， HOU Shance， et al. Comparison research on rainfall interpolation methods for small sample areas［J］.Research of Soil and Water Conservation， 2018， 25（3）： 117-121.（in Chinese）

[16] 王长虹， 朱合华， 钱七虎.克里金算法与多重分形理论在岩土参数随机场分析中的应用［J］.岩土力学， 2014， 35（S2）： 386-392.

[17] WANG Changhong， ZHU Hehua， QIAN Qihu. Application of Kriging methods and multi-fractal theory to estimate of geotechnical parameters spatial distribution［J］.Rock and Soil Mechanics， 2014， 35（S2）： 386-392. （in Chinese）

[18] 姜晓剑， 刘小军， 黄芬， 等. 逐日气象要素空间插值方法的比较［J］. 应用生态学报， 2010， 21（3）： 624- 630.

[19] JIANG Xiaojian， LIU Xiaojun， HUANG Fen， et al. Comparison of spatial interpolation methods for daily meteorological element［J］. Chinese Journal of Applied Ecology， 2010， 21（3）： 624-630. （in Chinese）

[20] 何䶮， 傅德平， 赵志敏， 等.基于GIS的新疆降水空间插值方法分析［J］. 水土保持研究， 2008， 15（6）： 35-37.

[21] HE Yan， FU Deping， ZHAO Zhimin， et al. Analysis of spatial interpolation methods to precipitation based on GIS in Xinjiang［J］.Research of Soil and Water Conservation， 2008， 15（6）： 35-37.（in Chinese）

[22] CHEN D， OU T， GONG L， et al. Spatial interpolation of daily precipitation in China： 1951-2005［J］. Advances in Atmospheric Sciences， 2010， 27（6）： 1221-1232.

[23] 高歌， 龚乐冰， 赵珊珊， 等. 日降水量空间插值方法研究［J］. 应用气象学报， 2007， 18（5）： 732-736.

[24] GAO Ge， GONG Lebing， ZHAO Shanshan， et al. Spatial interpolation methods of daily precipitation［J］. Journal of Applied Meteorological Science， 2007， 18（5）： 732-736.（in Chinese）

[25] ZHANG G， SU X， AYANTOBO O O， et al. Spatial interpolation of daily precipitation based on modified ADW method for gauge-scarce mountainous regions： A case study in the Shiyang River Basin［J］. Atmospheric Research， 2021， 247： 105167.

基金资助

国家自然科学基金资助项目(5207040968)

山西省青年科技研究基金(202203021212311)