利用熵函数模型解算GPS高程异常

邓捷利,陈天伟,杨茜茜,郑旭东

(桂林理工大学 a.广西空间信息与测绘重点实验室;b.测绘地理信息学院,广西 桂林 541006)

GPS定位技术被广泛应用于测量工作之中,GPS测量获取的是大地高程,但工程中通用的高程系统为正常高,可采用水准测量方法获得。大地高与正常高之间存在差值,称为高程异常,即H=h+ξ, 其中，H表示大地高,h为正常高,ξ为高程异常。若能获取某点的高程异常就能将该点的GPS高程应用于工程之中,所以确定高程异常在实际工程应用中尤为重要。

获取高程异常有多种方法,其中一种是在进行GPS测量的同时,再对部分点进行水准测量,就可得知这部分点的高程异常,通过这些已知高程异常的点,就能插值或拟合出测区内其他点的高程异常[1]。不同拟合方法适用条件不同,控制点呈直线排列的情形时(如公路、 铁路等),使用曲线拟合方法;当控制点呈面状分布且高程异常变化不明显才适合使用平面拟合;地形复杂则使用曲面拟合；还有一些常用拟合模型,加权平均值法(多用于山区地形)、曲面样条法、神经网络法等[2],本文主要讨论通过加权均值法给出权系数初值再使用最大熵法插值。

待定点所使用的已知点是离散的、随机的,所以插值得到高程异常变量也有着同样的性质。为使插值模型使用的数据充分离散,让模型更稳健,使用最大熵法。依据最大熵原理[3],当系统熵值最大,模型误差将会最小。

1 加权均值法原理

加权均值法是先对已知值进行定权，然后通过加权平均来求得期望值：

(1)

该插值方法的核函数是具有地理意义的,表达了点间距对高程异常插值的影响。核函数具有地理意义的还有克里金法、杨赤中法、分形插值等,这些插值方法的权系数也适合作为最大熵法的权系数初值。

2 最大熵原理

最大熵理论的数学模型为

(2)

其中:S[p]为系统的熵值;pi为该方法求得的权系数，pi≥0;gj表示各阶统计矩函数； E(gj)表示统计矩的期望;i=1,2,…n;j=1,2,…,m。

对于待定点而言,联测点(已知高程异常)的分布是离散的、随机的,为了能反映出空间属性,则该最大熵模型的约束条件应该包括待定点高程异常点的数学期望值、方差。

模型系统的熵值:

S[p]=-(p1lnp1+p2lnp2+…+pnlnpn)/ln 2;

(3)

约束条件为

(4)

以上约束条件从上到下分别为0、1、2阶矩约束。E1=E(ξ)是高程异常的期望值,因为随机变量的原点矩等于相应的样本矩[5],可由加权均值法算得。E2可以通过高程异常样本值和E1按计算方差估值的公式算得,此式即为二阶矩约束条件。本文的目的是在约束条件下使得系统的熵值S[p]趋向于最大,使得模型最大程度地依赖样本,则模型更稳健。

待定点的高程异常值ξ的估计量为

(5)

3 解算方法

最大熵模型是个附有限制条件的极值问题, 由于该函数是非凸函数, 难以使用拉格朗日乘数法求解, 多用数值分析进行求解。 因为传统的优化算法存在一些缺点, 对于某些目标函数不易求得最优解。

遗传算法是一种新兴算法,该算法解算快并且可以求出全局最优解。经过多年的研究,该方法因其多种优点得到了广泛的认同。遗传算法(GA)借鉴遗传学中自然选择进化过程来求解答案,该方法大致可以分为4个步骤:编码、初始群体的生成、适应度价值评估检查、操作算子[6]。下文具体论述初始群体生成中的初值确定和适应度价值评估检查中如何处理约束的问题。

3.1 初值的确定

遗传算法的寻优性能会受到初值设置的影响,初值设置十分重要。本文熵函数待求值pi是样本点(已知点)的权系数,而插值的方法多种多样,权系数的确定方法也复杂多样,基于上文加权均值法的优点,采用加权均值法所得到的结果作为一阶约束的样本矩，以反距离定权作为初值。

3.2 约束处理

遗传算法通常解决无约束问题,而的最大熵模型则包括了多个限制条件,因为没有限制条件的最大熵是没有意义的[7]。工程中常使用罚函数法来处理约束问题。该方法是通过一个因子把目标函数与约束条件链接成新的无约束的目标函数,构造方法一般有两种:一种是采用加法形式val(x)=f(x)+c(x); 另一种是乘法形式val(x)=f(x)·c(x), 其中c(x)是惩罚项。当解在非可行域,就会进行惩戒,目标函数就会增大；当解可行,惩罚值就低,目标函数就变小,通过迭代使得目标函数最小,从而求出最优解。

4 实例分析

一般来说, 地形平缓区域的建模精度要优于地形复杂区域, 故实例采用地形复杂的数据。 本文实验数据来源于黎剑[8]对区域GPS高程异常拟合及建模方法的研究, 对象为云南某地区, 该地区为山区地形,共27个GPS水准联测点(大地高和正常高已知), 剔除3个点, 将16个点作为样本点, 8个点作为待定点。 点分布如图1所示。

大地高和正常高相减可得到8个待定的高程异常, 记录在表1的第2列。 参照文献[8]中使用二次曲面拟合以及多面函数拟合得到的值记录在表1的3、 4列。 通过第1节所描述的方法即可得到某个待定点周围已知点的权系数,通过周围已知点权系数以及高程异常即可获取该待定点高程异常的插值。在搜索方式采用邻域搜索圆,搜索点数为4个的条件下使用加权均值法,得到的插值记录在表1第5列。

图1 24个点点位分布图Fig.1 Distribution view of the 24 points

表1 高程异常以及4种方法的插值结果

Table 1 Elevation anomalies and interpolation results of four methods m

点号真值二次曲面多面函数加权均值最大熵法431.97231.88631.79831.94131.927732.31232.28632.29232.22732.2691032.64932.43232.51932.55332.5331232.08832.20432.15132.12332.1051632.27832.2732.29832.22832.2571832.02932.10932.03632.08532.082132.47232.43732.49432.45232.4582332.17232.20232.20432.16632.153

在Matlab平台利用遗传算法函数进行编程,对最大熵模型进行解算,对该模型的解算本质就是附有限制条件的极值问题,因为函数本身性质,只能通过数值分析的方法进行求解。解算大致分为3个步骤:首先,将式(3)作为目标函数,将该函数编制函数文件sangobj.m;然后,将式(4)的约束条件编为罚函数文件sangcstr.m,其中边界约束pi>0在主程序中列出;最后,设置初值,如本文第3节描述,使用加权均值法得到的权系数作为种群初值。由以上步骤解算得各pi值为权系数。通过以上最大熵法获得的权系数和对应已知点的高程异常则可得到待定点的插值,结果见表1第6列。通过真值和拟合值之差可以得到残差,使用残差进行精度分析，结果见表2。

表2 各模型精度分析

Table 2 Precision analysis of each model cm

二次曲面多面函数加权均值法最大熵法最大残差21.717.49.611.6最小残差0.820.6-1.7中误差9.88.25.34.9

加权均值法与前两种方法对比得出:由于加权均值法中反距离加权法的核函数具有地理意义,更能反映地形变化对高程异常的影响,并使用搜索半径进行分段拟合,所以加权均值法在地形较为复杂的地区比起前两种方法拟合精度要更高一些。加权均值法在7号点拟合表现差于前两种方法,是因为7号点位于图幅边缘,该点并未被已知点包围,已知点分布不均匀,所以导致插值精度降低。16号点因为该点附近的样本点高程异常都比较相近,可推测出16号附近的高程异常变化不大,所以各种拟合方法差别不大。21号点各拟合方法精度也比较接近,与16号点相似,该点周围高程异常差变化不明显,且该点周围的4个已知点接近对称分布,对该点的拟合如同对相近的数值求均值,所以各种拟合方法差别很小。

对比加权均值法和最大熵法得出:使用最大熵模型结合遗传算法得出的插值数据精度大多高于加权均值法,虽然在最大残差、最小残差上最大熵法并未表现出优越性,但最大熵法中误差较低,比起加权均值法,其整体更接近真值,只有在少部分点表现出较差的结果,如23号点改进后比改进前要差,因为该点附近的已知点分布不均匀,多位于该点上方,该点的搜索范围下方缺少已知点信息，由于对样本依赖程度较大,插值结果偏向下方已知点的信息，精度反而变低。

5 结束语

当已知点分布较为均匀时,与核函数没有地理意义的二次曲面拟合法以及多面函数法等相比,采用加权平均值法,即使山区地形大地水准面变化较大,也能较好地解决高程异常拟合问题。

最大熵法最大程度地排除了主观因素的影响,让模型最大程度地依赖样本数据,使得模型更加稳健。与加权平均值法相比，最大熵法拟合出的数据整体上更为靠近真值,证明在地形复杂的区域使用最大熵法进行拟合有着良好的表现。