本文以某工业园区项目的GPS水准测量数据，利用地质统计学中的Ordinary Kriging方法，进行了GPS高程异常拟合的不同方案试验研究，结果发现普通克立格方法合理拟合方案下的估计精度可满足测量规范要求，推荐应用于GPS高程异常拟合工程。
GPS高程转换实质上主要是通过求取地面点的高程异常值，然后将外业GPS测量的大地高减去内业获取的高程异常值即可得到实用的正常高。而高程异常是地球重力场的重要参数之一，从理论上讲，实现GPS大地高向正常高转换的最好方法是综合利用GPS测量数据、重力测量数据、地球重力场模型和其他数学方法，而对于一般单位而言，在无法获取必要的重力资料的情况下，数学拟合方法仍然是目前进行GPS高程转换的首选方案。
关于GPS高程的转换问题，测绘界的许多专家学者提出了多种解决办法。这些方法归纳起来主要有：重力法、GPS水准法(数学模型拟合法)、数学模型抗差估计法、数学模型优化方法、GPS三角高程法、平差转换法、整体平差法、神经网络法等。如常用的移动加权平均法、曲线拟合法、曲面拟合法均属于数学模型拟合法。以上诸多方法在应用过程中均取得了一定的成效，但各自也都存在一些缺点。主要由于受地球区域密度分布异常、地表地形的复杂程度等因素的影响，有的方法理论精度很高，但与实际情况的吻合程度不高;而有的方法计算过程比较复杂，且实际效果也不够理想。因此，如何来进行高精度的GPS高程拟合与转换仍然是测量工程实践中需要不断研究和探讨的问题。
本文在介绍地质统计学基本理论及普通克立格Ordinary Kriging方法的原理及应用基础上，以某测区的实测GPS高程异常数据为例，进行了基于普通克立格方法的GPS高程拟合的相关试验分析与研究，按照不同已知高程异常点数据个数设计了四种方案，并对四种方案的拟合结果进行对比分析，证实了普通克立格方法用于GPS高程异常拟合的可行性，并探讨了普通克立格方法高程异常拟
合的精度，以及已知高程异常点数量、分布与高程拟合精度的关系。
1 基本原理
地质统计学是在经典统计方法的基础上，充分考虑地质变量的空间变化特征-相关性和随机性，并以反映地质现象区域化的随机函数-变差函数作为基本工具，来研究地质和采矿等工作中的各种问题。在结构分析的基础上采用各种Kriging法，来估计实际问题的。克立格法(Kriging)是一种推求最优、线性、无偏估计量的方法。即Kriging法是在考虑信息样品的形状、大小及其与待估块段间的空间分布位置等几何特征以及品位的空间结构之后，为了达到线性、无偏和最小估计方差的估计，而对每一样品值分别赋予一定的权系数，最后进行加权平均来估计块段品位的方法。
1.1变差函数
变差函数既能描述区域化变量的结构性变化，又能捕述其随机性变化，且它的计算是许多其他地质统计学计算的基础。
变差函数的定义如下：变差函数为区域化变量Z(x)和Z(x+h)的增量平方的数学期望，即区域化变量增量的方差。当区域化变量满足二阶平稳假设、本征假设或(准)二阶平稳假设、(准)本征假设时，变差函数只依赖于步长h，与x取值无关，就可以得到一维实验变差函数的计算公式：
式(1)即为Matheron推荐的传统计算公式，式中：N(h)是x轴上相隔h的点的对数，Z(xi)和Z(xi+h)是观测值Z(x)和Z(x+h)的N(h)对现实。如果Z(x)是定义在二维、三维空间的区域化变量，则x是二维、三维空间中的点，h是二维、三维空间中的向量。
变差函数不仅是Kriging估值、条件模拟等计算的基础，更重要的是它反映和刻画了区域化变量的许多性质，是分析区域化变量空间变异性的重要工具。一般情况下，代表高程异常的变差函数事先未知，可选用实用变差函数模型进行描述，然后根据实测资料进行拟合。设随机场是二阶平稳且为统计各向同性的，则常用的理论变差函数与交叉协方差函数理论球状模型如下：
式(2)中：C0为块金常数(Nugget Effect)，块金常数反映了变差函数在空间上的突变性，产生这种现象的原因可能是由于空间不相关的观测误差以及小于取样间距尺度上的空间变异引起的;C0+C为基台值(Sill)，C为拱高。变差函数基台值的大小可反映变量在该方向上变化幅度或总的变异程度的大小。块金常数大小可反映区域化变量随机性大小。块金值与基台值之比C0/(C0+C)反映块金方差占总空间异质性变异的大小。
1.2普通克立格法
将所研究的区域化变量，如高程异常记为U(x)，且U(x)满足二阶平稳假设或本征假设，E[U(x)]=m为未知常数，对区域内任一点x0处的变量U(x0)进行最优、无偏、线性估计的方法称为普通克立格法，且x0处的估计量可表示为：
对高程异常的变异函数和结构分析结果表明高程异常存在空间相关性，可视为区域化变量，则可利用普通克立格方法进行内插或外推。该方法的实质是根据高程异常这一区域化变量的原始数据和变异函数的结构特点，利用已知样点对未知样点的高程异常进行线性无偏、最优估计。估计时考虑了已知样点与未知高程异常样点的相互空间位置关系，以及变异函数提供的结构信息。
2 案例分析
从某测区的静态GPS控制网获得100个GPS高程异常数据点，分别选用20、30、40、50个样本点作为已知数据集，其余为未知点，用于外符合精度的检验。本文借助Arcgis9.0软件成熟的地质统计学模块，计算并拟合变差函数，并进行高程异常的空间性相关分析;最后利用普通克立格方法拟合高程异常，并进行精度分析与方案对比。探讨普通克立格方法在GPS高程异常拟合中的问题。
2.1Arcgis9.0软件的地质统计学模块简介
以20个样本点作为已知数据集的计算为例，首先启动Arcgis9.0软件的地质统计学模块，点击Geostatistical Analyst，选择Geostatistical Wizard(地统计分析向导)弹出如下对话框：
Method栏中选择Kriging方法，Dataset1中Inputdata选择事先做好的已知高程异常Excel数据表，Attribute选择“高程异常”，点击Next，
Geostatistical Method栏中选择Ordinary Kriging方法下的Prediction Map，其他默认后点击Next，弹出
选择变异函数模型为球状模型，不考虑带状各向异性，调节Number of Lags和Lag Size，获得左上角的绿色变异函数拟合曲线，用于克立格估计。20个已知点数据变异函数块金值0，基台值0.13476，变程0.24143。获得理想的变异函数模型后，点击Next，进行交叉验证，
20个高程异常点的预测值与已知值的分布曲线接近于绿色理论曲线，可由预测值与已知值的差值计算分析估计误差，除有一点高程异常拟合误差接近0.40m外，其余21个点拟合误差均小于0.10m，与普通水准测量精度相当。点击Finish后，完成克立格估计，生成测区高程异常差值。
2.2四种计算方案对比
四种方案变差函数图计较，
变差函数的计算采用球状模型拟合，四种方案的变差函数均具有较强的空间的结构性。变差函数块金值C0均为零，表明高程异常测量误差可以忽略不计，高程异常这一区域化变量的随机性较弱。C0+C为基台值，变差函数基台值的大小可反映变量在该方向上变化幅度或总的变异程度的大小，表1反映出高程异常具有很强的空间相关性，其空间相关半径(变程值)约在0.277附近。
显示四种方案的高程异常拟合外符合精度0.017m-0.047m，当已知点数量由20%增加至50%时，外符合精度逐渐提高，均为厘米水平，达到了四等水准测量的精度。而内符合精度0.104m-0.044m，当已知点数量由20%增加至50%时，内符合精度逐渐提高，除20%已知点方案外，其余均为厘米水平;值得注意的是当已知点为40%时，再增加已知点数量对精度的提高并无实际意义。
3 结论
本文建立了基于Kriging方法的高程异常变差函数模型，利用普通克立格公式进行插值计算，最后利用实测数据进行了交叉验证与精度分析。借助ArcGIS地质统计学模块中克立格插值工具，采用普通克立格(Ordinary Kriging)方法对本例数据进行高程异常拟合的试验计算，得到以下结论：
(1)高程异常数据的统计分布接近正态分布，高程异常的变差函数具有空间较强的空间结构性，可将其视为区域化变量，利用最佳拟合的变差函数，对未知高程异常点进行普通克立格插值。
(2)从四种方案的高程异常拟合图以及交叉误差表，可直观地分析拟合成果以及误差分布情况，进一步统计分析后，发现普通克立格方法具有较高的内、外符合精度，其外符合精度均达到了厘米级水平，与四等水准测量精度相当。
(3)随着已知数据的增多，高程异常估计精度逐渐提高，但当已知点为总采样点数的40%时，再增加已知点数量对高程异常拟合精度的提高已无实际意义。建议做高程拟合的已知点数量控制在总点数的40%以内即可。
作者：贾建峰 来源：西部资源
