地质雷达数据中克里金插值采样数据选择算法

王振武， 卜异亚， 马键强

(中国矿业大学(北京) 机电与信息工程学院，北京 100083)



地质雷达数据中克里金插值采样数据选择算法

王振武， 卜异亚， 马键强

(中国矿业大学(北京) 机电与信息工程学院，北京 100083)

采样数据选择的合理性是克里金算法插值结果是否准确的前提。针对如何从地质雷达数据中选择恰当的克里金插值采样数据的问题，本文提出了两种克里金插值采样数据选择算法，即动态三维球体覆盖(D3DGC)选点算法和双次反距离选点(DRD)算法。两种算法在基于粒子群优化的普通克里金算法下验证了采样数据选择的合理性，实验结果分析表明，D3DGC算法在最佳选点个数下误差率比传统选点方法降低了6.7%，在其他选点个数下误差率(<10%)也明显低于传统选点方法，DRD算法在参数e和c为4以及最佳选点个数的条件下误差率也达到了6.2%，在其他情况下误差率也低于传统选点方法。

克里金算法；插值；地质雷达；采样数据选择算法；动态三维球体覆盖选点算法；双次反距离选点算法；粒子群优化；误差率

地质雷达(ground penetrating radar，GPR)技术是一种基于电磁波反射原理，用于浅层地质构造探测和工程治理检测的地球物理方法[1]。为了深入勘查或检测探测区域的地质构造或地下结构，需要对测线数据做插值处理来获取探测区域的整体数据。在常见的空间插值方法中以克里金(Kriging)插值方法为代表的地质统计方法提供了揭示空间变量结构相关性与随机性的理论，是地质领域应用最广泛的插值方法之一[2]。在对克里金算法的研究中，有些学者根据不同的应用情况，将不同类型的克里金算法组合使用。牛文杰等将同位置协同克里金(Collocated CoKriging)和贝叶斯克里金方法相结合，推导出一种新的贝叶斯同位置协同克里金估值方法，综合了上述两种方法的优点[3]。李晓军等针对不连续地层，提出了一种采用指示克里金和普通克里金相结合的方法，即先用指示克里金法估计地层分布范围，再用普通克里金方法估计分布范围内的地层厚度，更多的学者将克里金方法和其他方法结合使用[4]。王辉赞等针对变异函数难以刻画实际数据分布的不足，采用最小二乘支持向量机从实际资料场拟合重构变异函数，表现出较好的针对性和客观性[5]。贾雨等将粒子群算法和克里金算法相结合，在粒子群优化过程中，通过高斯变异、样本点权重系数设定、搜索范围约束等方式提高了插值精度[6]。颜华等针对三维温度场重建问题，首先采用最小二乘法获得一个由少量像素描述的温度场，然后运用克里金内插和外推运算获得整个被测温度场的细致描述，这种方法可实现复杂的三维温度场的高精度重建[7]。徐驰等将高光谱反演与协同克里金方法相结合，将土壤表层含水率反演结果作为协同克里金插值的协变量来辅助进行协同克里金插值，提高了预测效率[8]。ZHANG等先采用克里金算法进行大型球磨机的功能适配，然后利用遗传算法优化了齿轮传送的可靠性问题[9]。FANG等将序列克里金算法和人工蜂群算法结合来优化点焊接头以改善其疲劳寿命，该方法用克里金模型来近似表达点焊接头设计模型的目标函数，然后采用人工蜂群算法进行迭代优化[10]。

在上述的研究中，对克里金插值方法的研究多集中在对不同类型克里金算法的组合[3-4]以及与其他算法理论的结合[5-11]中，鲜有学者研究克里金插值方法中采样数据的选择问题。一般情况下，按照空间相关性原理，地质数据应该选取距离插值点较近的数据点作为采样数据。文献[12]针对某储层数据，研究了空间相关分析的主要影响因素与克里金估计结构的关系，指出影响插值结构的因素和参数包括采样点的位置分布特征、变差函数模型的块金效应、变差函数模型的变程、变差函数模型的类型、空间结构的各向异性以及空间搜索范围等。文献[11]进一步指出，距离估计点较近的采样数据对估计结构的影响较大，且大多数情况下应该选用球状模型或指数模型。但文献[11]也只是给出了一般性的选点原则，对具体数据而言，应该如何选择采样数据并没有进行深入分析。文献[2]针对传统的克里金方法采用欧式距离测度描述空间结构相关性的缺点，基于最小方差准则，提出了综合曼氏、欧式和切氏三种距离测度的多测度加权克里金方法，并通过实验表明该方法提高了插值精度和可靠性，但是该方法也仅仅是对距离测度公式做了扩展，没有深入讨论采样数据如何选择的问题，而选择合适的采样数据是提高克里金插值准确性的前提。论文针对这一问题提出了动态三维球体覆盖(dynamic three-dimensional globe cover，D3DGC)选点算法和双次反距离(double reverse distance，DRD)选点算法，并基于粒子群优化的普通克里金算法验证了两种算法的有效性。

1 克里金插值采样数据选择算法描述

1.1 基本概念

本文在已获取地质雷达测线数据的基础上研究克里金插值的采样数据选择问题。为了对所提出的算法进行描述，首先给出相关概念。

定义1 地质数据点δ=〈ζ，γ〉：δ是一个二元组，其中ζ=〈x，y，z〉表示三维坐标值，γ表示地质点的属性值，如物质对电磁波的反射值等。

(1)

式中：δij表示测线数据上的第i行j列的地质数据点，且1≤i≤m，1≤j≤n。为表述方便，设一共p条测线数据，第l条测线表示为l，将测线上第h行的地质数据点集合定义为δh，h∈[1，m]，l∈[1，p]。

定义3 测线数据集Γ：Γ表示对同一探测区域的同一次探测活动中获取的测线数据总和，即

Γ={1，2，…，lp}

(2)

定义4 测线数据凸包Δ：Δ表示当前测线数据集所描述的最大三维探测区域。

Δ={δ|δ.ζ.x∈[Xmin，Xmax]，δ.ζ.y∈[Ymin，Ymax]，
δ.ζ.z∈[Zmin，Zmax]}

Xmin=min{k.δij.ζ.x}，Xmax=max{k.δij.ζ.x}
Ymin=min{k.δij.ζ.y}，Ymax=max{k.δij.ζ.y}
Zmin=min{k.δij.ζ.z}，Zmax=max{k.δij.ζ.z}

(3)

式中：k∈[1，p]，i∈[1，mk]，j∈[1，nk]，Xmin和Xmax分别表示x坐标的最大和最小值，min{ }和max{ }表示取集合的最小值和最大值操作。

定义5 待插值数据点集ψ：ψ=Δ-Γ，克里金插值的目的就是从Γ中选择合适的采样点对ψ进行插值。

定义6 凸包水平分割层η：η表示Δ按照测线平行方向的层次划分，即

(4)

Ψ={φ1，φ2，…,φu}

(5)

式中：φr表示第r层的待插值数据点的集合， 1≤r≤u且存在关系测线l中的δi满足δi⊆l.δi，l∈[1，p]，i∈[1，m]。

定义7 待插值点的采样数据点集φr[ε]→π[μ]：假设选择μ个采样点进行克里金插值，则φr中第ε个插值点的μ个采样点的集合表示为φr[ε]→π[μ]。

定义8 距离计算公式τ：论文采用欧几里得距离作为距离测度公式。设ψ中的某待插值点为θ(a，b，c)，测线上某点l中δ的ζ值为(x，y，z)，则距离τ定义为

(6)

1.2 采样数据选择算法原理

地质数据具有空间相关性和随机性的特点，针对不同数据的特点，本文提出了D3DGC选点算法和DRD选点算法。

1.2.1 D3DGC 算法

根据地质数据的相关性原理，一般而言选取距离插值点较近的数据点作为采样点，针对地质雷达数据特点，论文提出了D3DGC选点算法，算法原理如图1所示。

图1 D3DGC算法原理图Fig.1 Theory graph of D3DGC

定义9 球与直线交点方程：设待插值点为θ(a，b，c)， 球的半径为R，l上的δh与球的交点φ(x，y，z)计算公式为

(7)

设测线垂方向采样点步长为d，每个待插值点选取μ个采样点进行插值。基于如上定义，D3DGC算法的流程如图2所示。

图2 D3DGC算法流程图Fig.2 The flow graph of D3DGC

对于每一个待插值点θ(a，b，c)，球的半径以d为单位递增，计算l的δh与球体交截后包含在球体内的地质数据点个数kl，这里直到条件满足时为止。通过逐层计算ηi上的φi，直到所有待插值点的采样点选择完毕为止。

如果θ(a，b，c)的插值采样点只从ηi中的δi选择，其中δi属于l，1≤≤p则D3DGC算法就简化为二维模式，如图3所示，称之为动态圆覆盖(dynamic circle cover，DCC)选点算法。

准时制（JIT)，可概括为“任何时间、任何地点、任何事情”，都可以及时解决和处理，以达到零时间、零距离的高效作业模式。审批平台可以引入此理论在审批业务流程的设计和模型的架构当中，以实现审批的高效、便捷。具体JIT模式的优势见图1。

1.2.2 DRD算法

图3 DCC算法原理图Fig.3 Theory graph of DCC algorithm

定义10 点到直线的距离ϑ。设φi中当前的待插值点为θ(a，b，c)，l上的δh可由式(7)的方程2描述，若将方程2简记为则θ(a，b，c)到l上的δh的距离可定义为

(8)

设待插值点选取μ个样本点进行插值，.δh的点数为n(为插值准确性和效率考虑，通常μ≪n)，按照式(8)计算距离后，分配到第l条测线上的采样点数由式(9)计算：

(9)

(10)

图4 DRD算法原理图Fig.4 Theory graph of DRD

图5 DRD算法流程图Fig.5 The flow graph of DRD

2 算法实验分析

由于主观设置的变差函数模型参数对克里金算法的插值结果影响较大，为了客观准确地验证两种采样点选择算法，论文采用了粒子群优化(particle swarm optimization，PSO)算法对普通克里金算法的变差函数模型参数进行优化。

PSO算法将优化问题的候选解假设为一个无体积和质量的粒子，在d维空间中粒子通过自身和其他粒子的经验来调整自己的位置。设有n个粒子组成的粒子群，其中第i个粒子在d维空间中的速度为vi=[vi1vi2…vid]T，位置为xi=[xi1xi2…xid]T，粒子i的当前个体最优位置为pid，整个粒子群的当前全局最优位置为gd，粒子i通过式(11)的迭代来更新自己的速度和位置。

(11)

定义11 适应度函数F:本文给定的PSO算法的适应度函数表示为

(12)

式中：F(δo) 为第o个粒子(地质数据点)的适应度函数值，V*(δo)表示粒子的克里金估计量，V(δo)表示例子的实际值。论文采用普通Kriging算法进行实验结果的验证，变差函数采用球状模型，拱高f、块金值g、变程s，参数采用PSO算法进行优化，粒子的速度和位置用向量〈f,g,s〉表示，粒子群的最大迭代次数α=1 000，适应度函数阈值β=10-6。PSO优化的克里金(PSO-Kriging)算法的流程如图6所示。

图6 PSO-Kriging算法流程图Fig.6 The flow graph of PSO-Kriging

本文基于地质雷达探测数据来验证所提出的两种采样点选择算法的有效性。地质雷达仪为中国矿业大学(北京)自主研发的MTGR-4F车载地质雷达，雷达参数为：探测深度5m，时间窗口100ns，采样点数512，天线主频200MHz。实验数据探测路段为北京市朝阳区东三环中路，探测目的有3个：地底的分层结构、有无地下管线以及有无空洞。实验数据取值范围为[-23 679，15 421]；测线条数为3，采样道数分别为18 632、18 088和18 838，实验所用数据的采集工作照如图7所示。

图7 实验数据的获取Fig.7 The acquirement of experiment data

实验硬件和软件环境信息如下：CPU为Pentium(R)Dual-CoreT4300 2.10GHz，内存为RAM2GB，操作系统为Windows7 32位，实验平台为VS.NET2012，开发语言为C++。

传统的采样点选择方式是对同一ηi，1≤i≤u上所有l的δi数据按照式(6)进行排序，选择距离最近的某些地质数据点作为克里金插值的采样点。本文将所提出的两种算法和这种传统的选点方法进行了比较，为准确起见，误差率统计的为N0=2 000个点的平均误差率，误差率ρ由式(13)计算：

(13)

实验探测数据的最优采样数据点数与采样道数、采样点数等具体的地质雷达仪参数有关。考虑到插值算法的运算效率和插值结果精度等因素，本实验选取的数据点个数μ∈[160，300]，依次按10个点的数量递增。如图8所示，在不同的选点个数下，D3DGC算法的误差率明显低于传统选点算法。当选点个数为210时，两种算法的误差率都为最小，D3DGC算法的误差率仅为1.3%，而传统选点算法的误差率却为8%，而且当选点个数在180～270时，D3DGC算法的误差率都小于10%，而传统选点算法只有在200～210误差才小于10%，且都大于8%，在其他范围内的最大误差达到了31.1%。与传统选点算法相比，D3DGC算法效果好的原因在于，该算法在三维空间中考虑相关性，按照球状半径递增的原理逐步寻找采样数据点，比传统选点算法更为准确合理。

图8 D3DGC算法实验结果Fig.8 The experiment results of D3DGC

表1 DRD算法实验结果

3 结论

1)D3DGC算法是在三维空间中按照球状半径递增原理选取的采样点，因此其克里金插值误差率明显地比传统选点方法好；2)针对数据量大且相邻数据值较为接近的地质雷达数据，DRD算法综合考虑了地质数据的空间相关性和随机性，在参数e和c取值为4的条件下也取得了不错的效果。

D3DGC和DRD算法均已在道路塌陷的检测问题中得到了应用，使得克里金插值的精度得到了提高。对于DRD算法，针对不同的地质雷达数据，如何实现参数e和c的自适应取值是将来需要进一步研究的问题。

[1]李晋平，邵丕彦，谷牧.地质雷达在铁路隧道工程质量检测中的应用[J].中国铁道科学， 2006， 27(2)： 56-59.

LI Jinping， SHAO Piyan， GU Mu. Application and analysis of GPR in railway tunnel engineering quality inspection[J]. China railway science， 2006，27(2)： 56-59.

[2]刘志平，何秀凤，张淑辉.多测度加权克里金法在高边坡变形稳定性分析中的应用[J].水利学报， 2009， 40(6)： 709-715.

LIU Zhiping， HE Xiufeng， ZHANG Shuhui. Multi-distance measures weighted Kriging method for deformation stability analysis of steep slopes[J]. Shuili xuebao， 2009，40(6)： 709-715.

[3]牛文杰，朱大培，陈其明.贝叶斯同位置协同克里金(Bayesian Collocated CoKriging)估值法的研究[J]计算机辅助设计与图形学学报， 2002， 14(4)： 343-347.

NIU Wenjie， ZHU Dapei， CHEN Qiming. Bayesian collocated CoKriging[J]. Journal of computer-aided design & computer graphics， 2002， 14(4)： 343-347.

[4]李晓军，张振远.基于指示和普通克里金的不连续地层厚度估计方法[J].岩土力学， 2014， 35(10)： 2881-2887.

LI Xiaojun， ZHANG Zhenyuan. A combined indicator-ordinary Kriging method for estimating thicknesses of discontinuous geological strata[J]. Rock and soil mechanics， 2014， 35(10)： 2881-2887.

[5]王辉赞，张韧，刘巍，等.支持向量机优化的克里金插值算法及其海洋资料对比试验[J].大气科学学报， 2011， 34(5)： 567-573.

WANG Huizan， ZHANG Ren， LIU Wei， et al. Kriging interpolation method optimized by support vector machine and its application in oceanic data[J]. Transactions of atmospheric sciences， 2011， 34(5)： 567-573.

[6]贾雨，邓世武，姚兴苗，等.基于约束粒子群优化的克里金插值算法[J].成都理工大学学报：自然科学版， 2015， 42(1)： 104-109.

JIA Yu， DENG Shiwu， YAO Xingmiao， et al. Kriging interpolation algorithm based on constraint particle swarm optimization[J]. Journal of Chengdu University of Technology：science & technology edition， 2015， 42(1)： 104-109.

[7]颜华，李欣，王善辉.基于最小二乘法和克里金插值的三维温度场重建[J].沈阳工业大学学报， 2014， 36(3)： 303-307.

YAN Hua， LI Xin， WANG Shanhui. Reconstruction of three-dimensional temperature field based on least-square method and kriging interpolation[J]. Journal of Shenyang University of Technology， 2014， 36(3)： 303-307.

[8]徐驰，曾文治，黄介生，等.基于高光谱与系统克里金的土壤耕作层含水率反演[J].农业工程学报， 2014， 30(13)： 94-103.

XU Chi， ZENG Wenzhi， HUANG Jiesheng， et al. Modelling soil water content in plow layer using co-kriging interpolation based on hyperspectral data [J]. Transactions of the Chinese society of agricultural engineering， 2014，30(13)： 94-103.

[9]ZHANG Guanyu， WANG Guoqiang， LI Xuefei， et al. Global optimization of reliability design for large ball mill gear transmission based on the Kriging model and genetic algorithm[J]. Mechanism and machine theory， 2013， 69:321-336.

[10]FANG Jianguang， GAO Yunkai， SUN Guangyong， et al. Optimization of spot-welded joints combined artificial bee Colony algorithm with sequential sriging optimization[J]. Advances in mechanical engineering， 2014: 1-10.

[11]DENG Mingguo， LI Wenchang， LI Bo， et al. Application of log kriging on estimated reserves of the 10-9 ore body of lutangba in the gejiu tin deposits[J]. Journal of China university of mining & technology， 2007， 17(2)： 286-289.

[12]刘永社，印兴耀，贺维胜.空间相关分析因素对储层建模中克里金估计结果的影响[J].石油大学学报：自然科学版， 2004， 28(2)： 24-29.

LIU Yongshe， YIN Xingyao， HE Weisheng. Effects of the factors for spatial correlation analysis on Kriging estimation in reservoir modeling[J]. Journal of the University of Petroleum： science & technology edition， 2004， 28(2)： 24-29.

本文引用格式：

王振武， 卜异亚， 马键强. 地质雷达数据中克里金插值采样数据选择算法[J]. 哈尔滨工程大学学报， 2017， 38(5)： 784-790.

WANG Zhenwu，BU Yiya，MA Jianqiang. Sampling data selection algorithm for Kriging interpolation based on ground penetrating radar data[J]. Journal of Harbin Engineering University， 2017， 38(5)： 784-790.

Sampling data selection algorithm for Kriging interpolation based on ground penetrating radar data

WANG Zhenwu，BU Yiya，MA Jianqiang

(School of Mechanical Electronic & Information Engineering， China University of Mining & Technology (Beijing)， Beijing 100083， China)

The rationality of the data sampling selection method is a prerequisite for judging whether the Kriging interpolation results are accurate. Considering the problem of how to properly choose sampling data from the probing data of ground penetrating radar， in this paper， we propose two sampling selection algorithms for Kriging interpolation: the dynamic three-dimensional globe cover (D3DGC) algorithm and the double reverse distance (DRD) algorithm. We verified the rationality of these two methods with respect to the ordinary Kriging algorithm based on particle swarm optimization. The experimental results show the D3DGC algorithm to reduce the error rate to 6.7% compared with the traditional method for an optimal number of sampling points. With respect to other sampling numbers， D3DGC also demonstrated obvious advantages (error rate<10%). The DRD algorithm error rate was 6.2% for the optimal number of sampling points when the parameter e and c values were 4， respectively. For other sampling numbers， the DRD algorithm also obtained better results than the traditional method.

Kriging algorithm; interpolation; ground penetrating radar; sampling selection algorithm; dynamic three-dimensional globe cover; double reverse distance; particle swarm optimization(PSO); error

2016-01-07.

日期：2017-04-26.

国家高技术研究发展计划重大专项(2012AA12A308)；核设施退役及放射性废物治理科研项目(FZ1402-08)；国家自然科学基金项目(61302157)；北京市高等学校青年英才计划项目(YETP0939).

王振武(1978-)， 男， 副教授，博士.

王振武，E-mail: wangzhenwu@126.com.

10.11990/jheu.201601022

TP39

A

1006-7043(2017)05-0784-07

网络出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20170426.1152.046.html