地质雷达(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〉表示三维坐标值，γ表示地质点的属性值，如物质对电磁波的反射值等。

定义2 测线数据ℓ:ℓ是地质雷达仪沿测线方向移动所探测到的数据的集合。为讨论方便，本文假设同一探测区域的测线具有相同的采样点数(m)和采样道数(n)，即

$\ell \text{=}\left[ \begin{matrix} {{\delta }_{11}} & \cdots & {{\delta }_{1n}} \\ \vdots & \vdots & \vdots \\ {{\delta }_{m1}} & \cdots & {{\delta }_{mn}} \\ \end{matrix} \right]\text{=}\left[ \begin{matrix} {{\delta }_{1}} \\ \vdots \\ {{\delta }_{m}} \\ \end{matrix} \right]$

(1)

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

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

$\mathit{\Gamma }=\left\{ {{\ell }^{1}},{{\ell }^{2}},\cdots ,{{\mathit{l}}^{p}} \right\}$

(2)

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

$\begin{align} & \mathit{\Delta }=\left\{ \delta |\delta .\zeta .x\in \left[ {{X}_{\rm{min}}},{{X}_{\rm{max}}} \right],\delta .\zeta .y\in \left[ {{Y}_{\rm{min}}},{{Y}_{\rm{max}}} \right], \right. \\ & \quad \quad \quad \quad \left. \delta .\zeta .\mathit{z}\in \left[ {{Z}_{\rm{min}}},{{\mathit{Z}}_{\rm{max}}} \right] \right\} \\ & {{X}_{\rm{min}}}=\rm{min}\left\{ {{\ell }^{\mathit{k}}}.{{\delta }_{\mathit{ij}}}.\zeta .\mathit{x} \right\},{{X}_{\rm{max}}}=\max \left\{ {{\ell }^{\mathit{k}}}.{{\delta }_{\mathit{ij}}}.\zeta .\mathit{x} \right\} \\ & {{Y}_{\rm{min}}}=\rm{min}\left\{ {{\ell }^{\mathit{k}}}.{{\delta }_{\mathit{ij}}}.\zeta .\mathit{y} \right\},{{Y}_{\rm{max}}}=\max \left\{ {{\ell }^{\mathit{k}}}.{{\delta }_{\mathit{ij}}}.\zeta .\mathit{y} \right\} \\ & {{Z}_{\rm{min}}}=\rm{min}\left\{ {{\ell }^{\mathit{k}}}.{{\delta }_{\mathit{ij}}}.\zeta .\mathit{z} \right\},{{Z}_{\rm{max}}}=\max \left\{ {{\ell }^{\mathit{k}}}.{{\delta }_{\mathit{ij}}}.\zeta .\mathit{z} \right\} \\ \end{align}$

(3)

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

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

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

${{\mathit{\boldsymbol{ }}\!\!\eta\!\!\rm{ }}_{t}}=\left[ \begin{matrix} {{\delta }_{11}} & \cdots & {{\delta }_{1v}} \\ \vdots & \vdots & \vdots \\ {{\delta }_{u1}} & \cdots & {{\delta }_{uv}} \\ \end{matrix} \right]$

(4)

式中：显然u≥m，v≥n，且Δ=$\sum\limits_{t=1}^{u}{{{\eta }_{t}}}$。ψ按测线平行方向分割可以定义为

$\mathit{\Psi }=\{{{\varphi }_{1}},{{\varphi }_{2}},\cdots ,\rm{ }{{\varphi }_{u}}\}$

(5)

式中：φ_r表示第r层的待插值数据点的集合，1≤r≤u且存在关系测线ℓ^l中的δ_i满足δ_i⊆η_i，η_i=φ_i+${{\sum\limits_{i=1}^{p}{\ell }}^{l}}$.δ_i，l∈[1，p]，i∈[1，m]。

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

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

$\tau ={{\left( x-a \right)}^{2}}+{{\left( y-b \right)}^{2}}+{{\left( z-c \right)}^{2}}$

(6)

1.2 采样数据选择算法原理

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

1.2.1 D3DGC算法

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

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

$\left\{ \begin{array}{l} {\left( {\mathit{x} - \mathit{a}} \right)^2} + {\left( {\mathit{y} - \mathit{b}} \right)^2} + {\left( {\mathit{z} - \mathit{c}} \right)^2} = {R^2}\\ \left\{ \begin{array}{l} \frac{{\mathit{x} - {\ell ^l}.{\delta _{i1}}.\zeta .\mathit{x}}}{{{\ell ^l}.{\delta _{i1}}.\zeta .x - {\ell ^l}.{\delta _{in}}.\zeta .\mathit{x}}} = \\ \frac{{y - {\ell ^l}.{\delta _{i1}}.\zeta .y}}{{{\ell ^l}.{\delta _{i1}}.\zeta .y - {\ell ^l}.{\delta _{in}}.\zeta .\mathit{x}}} = \\ \frac{{z - {\ell ^l}.{\delta _{i1}}.\zeta .z}}{{{\ell ^l}.{\delta _{i1}}.\zeta .z - {\ell ^l}.{\delta _{in}}.\zeta .\mathit{x}}} \end{array} \right. \end{array} \right.$

(7)

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

对于每一个待插值点θ(a, b, c)，球的半径以d为单位递增，计算ℓ^l的δ_h与球体交截后包含在球体内的地质数据点个数k_l，这里h=h±v，v∈[1，R/d]，直到条件$\sum\limits_{l = 1}^p {{k_l}}$≥μ满足时为止。通过逐层计算η_i上的φ_i，直到所有待插值点的采样点选择完毕为止。

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

1.2.2 DRD算法

当ℓ中δ的γ较为相似且测线数据量较大时，只按照相关性原理进行采样点选择，插值效果可能不理想，这时需要将随机性考虑进去，基于此论文提出了DRD选点算法。

定义10 点到直线的距离ϑ。设φ_i中当前的待插值点为θ(a, b, c)，ℓ^l上的δ_h可由式(7) 的方程2描述，若将方程2简记为$\frac{{x - {x_0}}}{X} = \frac{{y - {y_0}}}{Y} = \frac{{z - {z_0}}}{Z}$则θ(a, b, c)到ℓ^l上的δ_h的距离可定义为

${\vartheta _{lh}} = \frac{{{{\sqrt {{{\left| \begin{array}{l} b - {y_0}\quad c - {z_0}\\ \quad Y\quad \quad \;Z \end{array} \right|}^2} + {{\left| \begin{array}{l} c - {z_0}\quad a - x\\ \quad Z\quad \quad \,X \end{array} \right|}^2} + \left| \begin{array}{l} a - {x_0}\quad b - {y_0}\\ \quad X\quad \quad \,X \end{array} \right|} }^2}}}{{\sqrt {{X^2} + {Y^2} + {Z^2}} }}$

(8)

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

${\pi _l} = \left\lfloor {\mu \times \frac{{\frac{1}{{{{({\vartheta _{lh}})}^e}}}}}{{\sum\limits_{k = 1}^p {\frac{1}{{{{({\vartheta _{lh}})}^e}}}} }}} \right\rfloor $

(9)

且满足$\sum\limits_{l = 1}^p {{\pi _l}} $=μ，参数e是常数，表示反距离比例调整因子，此值越大，则距离待插值点越近的测线分配的采样点越多。

接下来的问题是如何确定分配到第l条测线上的采样点数π_l的分布。如图 4(a)所示，设第l条测线到待插值点p距离最近的点为A_l，把测线l等间距划分为π_l个区间，区间i内的点数为n/π_l」，i∈[1，π_l]。设A_l落于第j(j≤π_l)个区间内，则落于区间i内的采样点数按照式(10) 计算。

${\sigma ^l}_i = \left\lfloor {\frac{{{{\left( {\frac{1}{{\left| {i - j} \right| + 1}}} \right)}^e}}}{{{{\sum\limits_{i = 1}^{{\pi _l}} {\left( {\frac{1}{{\left| {i - j} \right| + 1}}} \right)} }^e}}} \times {\pi _l}} \right\rfloor $

(10)

对于分配于区间i内的σ_i^l个采样点，按等间距原则从区间i中的n/π_l」个地质数据点中选择采样点。这里c是常数，表示反距离比例调整因子，此值越大，则距离A_l越近的点越容易被选择为采样点，如图 4(b)所示。DRD算法的流程如图 5所示。

2 算法实验分析

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

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

$\left\{ \begin{array}{l} {\mathit{\boldsymbol{v}}_i}^{(t + 1)} = \mathit{\boldsymbol{\omega }}{\mathit{\boldsymbol{v}}^t}_i + {c_1}{r_1}\left( {{p^t}_{id} - {x^t}_i} \right) + {c_2}{r_2}\left( {{g^t}_d - {x^t}_i} \right)\\ {\mathit{\boldsymbol{x}}^{(t + 1)}}_i = {\mathit{\boldsymbol{x}}^t}_i + {\mathit{\boldsymbol{v}}_i}^{(t + 1)} \end{array} \right.$

(11)

式中：v_i^(t+1)和v_i^t分别为粒子i下一时刻和当前时刻的速度，x_i^(t+1)和x_i^t分别为粒子下一时刻和当前时刻的位置，ω为惯性权重，ω较大时算法具有较强的全局搜索能力，ω较小则算法倾向于局部搜索。c₁和c₂是学习因子，为非负常数，r₁和r₂是两个[0, 1]的随机数。

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

$F({{\delta }_{o}})=\left| {{V}^{}}^{*}({{\delta }_{o}})-V({{\delta }_{o}}) \right|$

(12)

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

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

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

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

$\rho =\frac{100}{{{N}_{0}}}\sum\limits_{i=1}^{{{N}_{0}}}{\left| \frac{{{v}_{i}}-{{v}^{}}{{_{i}}^{*}}}{{{v}_{i}}} \right|}\times 100%$

(13)

式中：v_i和v_i^*分别表示第i个插值点的实际值和计算值。考虑地质数据的空间相关性，可以运用D3DGC算法进行采样点的计算。一般来说，应该选择距离插值点较近的点作为采样数据点，但采样数据点选择越多，它们距离插值点就越远，最终的插值效果未必最好，因此如何选择采样点个数对于克里金插值结果至关重要。

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

当ℓ数据量大且相邻数据值较为接近时，在考虑地质数据的空间相关性的同时也可以把随机性考虑进来，DRD算法是兼顾相关性和随机性的一种算法，算法得到的实验结果如表 1所示。如表 1所示，实验结果分别给出了e和c的值为3、4、5的情况，当e=c=4时，DRD算法的误差率低于传统选点算法，当选点个数在210附近时误差率最小(6.2%)。所以当整体数据相似性过高而需要加权离散选点时，同时考虑空间相关性和随机性是一种可行的选择。参数e和c的选择与ℓ数据量大小以及相邻数据的相似程度有关，一般来说，按照经验e和c选取[3, 5]之间的值。DRD算法在ℓ^l，l∈[1，p]以及ℓ^l的不同区间上数据点的分配随距离呈几何关系衰减，对于实验所用数据，当e和c的值为3和5时，实验结果与其值为4时的结果相差不大。

传统选点算法的时间复杂度与ℓ^l，l∈[1，p]的采样道数n、选点个数μ以及ℓ的条数p有关，即O(n×p×μ)，由于p≤μ≤n，所以O(n×p×μ)≈O(n)。而D3DGC算法由于在三维空间中考虑相关性，因此时间复杂度与采样点数m、采样道数n、ℓ的条数p以及选点个数μ都相关，即O(m×n×p×μ)，同样地，因为p≤μ≤n，所以O(m×n×p×μ)≈O(m×n)。对于DRD算法而言，它是三种算法时间复杂度最小的，只与选点个数μ和ℓ的条数p相关，即O(p×μ)。

3 结论

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

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