Kriging的两步极小法在GPS高程拟合中的应用*

王东东，郭 戬，赵旭坤，刘 兵，王 静

(1.浙江建设职业技术学院，浙江 杭州 311231; 2.陕西铁路工程职业技术学院，陕西 渭南 714000)

工程实践过程中，利用GPS手段求取正常高，关键要先得到高程异常值。高程异常的获取可以通过地球重力场模型求取，重力场方面的数据由于保密和工程成本的原因，很少使用；通过函数拟合的方法求取高程异常方法简单，可操作性比较强，因此，这种方法有广阔的应用前景。由于随机场往往包含趋势性和随机性，因此采用两步极小法进行处理更加合理[1]。Kriging方法是法国学者G.Materon于1962年在“应用地质统计学论”中首次提出和发展的。它是地质统计学的核心，是对局部区域随机过程的一种预测方法,是建立在变异函数和结构分析的基础之上，具有求线性、无偏和最优的内插估计值，统计特性是方差最小和平滑效应。普通Kriging法是以待测点为中心的邻域内一些已测点数据、已测点之间的距离和变异函数提供的结构信息，得到对每个已测点的权系数，然后进行加权平均计算未测点的估值[2-5]。所谓两步极小法，即首先利用二次多项式处理趋势性部分[6-8]，然后Kriging处理随机性部分，该方法兼顾了趋势和随机两个部分，使构建模型更加符合实际情况。在变异函数的拟合过程中，对拟合点的选取很少考虑拟合点对拟合的贡献大小，采用Cook距离筛选拟合点，找出一些对拟合贡献较大的点，同时剔除一些对拟合产生干扰的点，可以提高变异函数的参数求解精度。

1 二次多项式拟合法

采用二次多项式函数作为曲面拟合的函数模型：

(1)

式中，ζi为高程异常值；a0、a1、a2、a3、a4、a5为待求参数；xi、yi为平面坐标；Δi为误差。利用该函数模型拟合一个曲面，作为大地水准面，对待求点进行插值。该函数模型的矩阵形式为：

(2)

即

ζ=AX+Δ

(3)

误差方程为：

v=ζ-AX

(4)

通过n(n≥6)个水准重合点，利用最小二乘法解算待求参数，带入式(3)，解算出插值点的高程异常。

2 普通Kriging法原理

Kriging是基于变异函数的局部区域的最优无偏估计。假设区域变量满足本征假设或二阶平稳，其属性值为Z(xi)。未测点的估计值是由该区域已测点的观测值加权求和得到的，即为：

(5)

式中，Z*(x0)为未测点的估计值；λi为已测点属性值的权系数。

由无偏性和方差最小可得：

(6)

由二阶平稳和本征假设进一步可得：

(7)

式中，j=1,2,3,…,n;γ(xi-xj)为变异函数；μ为拉格朗日乘子。上式的矩阵形式为：

(8)

式中，γij为已测点xi和xj之间的变异函数；γi0为已测点xi和未测点x0之间的变异函数。由变异函数及式(8)解得λi和μ。变异函数的模型有指数模型、高斯模型和球状模型[4]等，GPS高程拟合一般采用球面模型，其模型为：

(9)

式中，C0为块金值；C为拱高；a为变程。

实际应用中，变异函数可由下式计算：

(10)

式中，N为已测点中距离为hm的点的对数。从模型上可以看出，它是一个距离的函数。

若要求取球面模型的参数，需要利用式(10)计算各个间隔上的变异函数值。计算出所有已测点之间的距离并按间距d进行分组。由于距离比较离散，点对之间的距离若在我们设定的步长范围之内时，就归为一组，此范围内所有点对距离的平均值为本组变异函数的距离h，步长值可设置为d/2。分组模型为：

(11)

式中，maxh和minh为所有已测点之间距离的最大值和最小值；Nh为所有已测点之间距离的分组数。将求得的变异函数的值带入式(9)，利用最小二乘法求解未知参数C0、C和a。由于球面模型是非线性函数，可将其线性化为：

γ(hm)=C0+k1h+k2h3

(12)

3 Cook距离的原理

1977年，Cook距离从参数置信区域提出，经过多年的发展已经应用在各种统计模型中。本文主要应用Cook距离筛选变异函数的拟合点，提高参数求解的精度，进而提高Kriging插值的精度。一般线性模型为：

(13)

(14)

(15)

Cook距离越小，表明剔除i行数据后，其对回归参数的影响比较小，该点适合参与变异函数的拟合；Cook距离越大，表明剔除i行数据后，其对回归参数的影响比较大，则该点不利于变异函数的参数的求解，拟合时应该将其剔除，从而达到提高变异函数参数求解的精度。

4 算例分析

为了验证Cook距离对Kriging法中变异函数的参数求解中的有效性，选取了某地区地形平缓，无粗差、同精度的GPS点和水准重合点17个，其GPS网按照国家B级网施测，正常高按照二等水准测量施测。高程异常通过大地高和正常高容易得到，相关数据如表1所示。

表1 已知点坐标和高程异常[11]/m

本文设置三组实验进行计算，三组实验都采用两步极小法进行处理。先进行二次多项式拟合，对数据的趋势项部分进行处理；然后利用普通Kriging法对随机部分进行处理。实验过程中，采用球面模型作为变异函数。

方案Ⅰ：令点5、26、4、28为检核点，其余13个点作为变异函数的拟合点。

方案Ⅱ：令点5、26、4、28、22、6、9为检核点，其余10个点作为变异函数的拟合点。

方案Ⅲ：从17个点中选取4个Cook距离最大的值作为检核点，如表2所示，分别为26、22、6、9，其余13个点作为变异函数的拟合点。

表 2 已知数据的Cook距离

从表2可以看出，6、9、22、26点的距离数值比较大，对回归参数的影响较大，不利于模型的构建。

通过对三个实验方案进行解算，得到二次多项式法和基于二次多项式的Kriging法的外符合精度，如表3所示。

表3 三个方案的外符合精度/mm

由实验方案和表2可知，在拟合点的选择上，方案Ⅰ是随机的，方案Ⅱ和方案Ⅲ是经过Cook距离选择的。三种方案的解算结果如表3所示，从计算结果比较中可以得出：

(1)方案Ⅲ中采用二次多项式法和基于二次多项式的Kriging法的外符合精度分别为4.2 mm、3.2 mm，均优于方案Ⅰ，表明经过Cook距离筛选拟合点是可行的，可以剔除部分强扰动点，进一步提高模型精度。

(2)方案Ⅰ、方案Ⅱ和方案Ⅲ中，基于二次多项式法的外符合精度分别为6.0 mm、3.1 mm、4.2 mm；基于二次多项式的Kriging法的外符合精度分别为4.9 mm、2.9 mm、3.2 mm。可以得出，基于二次多项式的Kriging法的外符合精度相比二次多项式法的外符合精度有一定的改善，验证了两步极小法能够提高GPS高程拟合的精度。

(3)由实验方案可知，方案Ⅱ采用的拟合点数量要比方案Ⅰ、方案Ⅲ少，方案Ⅱ采用二次多项式法和基于二次多项式的Kriging法的外符合精度均优于方案Ⅰ、方案Ⅲ，这是由于方案Ⅱ采用的拟合点是Cook距离较小的点，提高了模型解算精度。可以得出：在模型构建的过程中，不仅要顾及拟合点的数量，也要考虑拟合点扰动性的强度。

5 结 论

Kriging方法的关键之处是变异函数选择和参数的求解，GPS高程拟合往往采用球面模型，此时只有得到高精度的参数解，才能得到很好的插值精度。变异函数参数求解过程中，由于一些点的强扰动性影响了参数求解的精度。利用Cook距离剔除了一些对求解过程中会产生较大干扰的点，可以提高变异函数参数求解的精度。Cook距离较大时表明该点离趋势面较远，属于离群数据，选择拟合点时要剔除这些离群数据，进一步提高建模的精度。两步极小法可以兼顾随机场的趋势性和随机性，实验也验证了这一理论的正确性。由于本文采用的数据点比较少，下一步将会对更大范围的数据进行相关研究。