随着GPS技术在测绘领域的发展和完善，采用部分GPS水准重合点实现高程拟合越来越受到测绘工作者的青睐^[1]。获取高精度的高程异常是高程拟合的核心内容^[2]。对于地势起伏平缓的研究区域，常规高程拟合方法能够满足实际拟合精度的要求。但当地势起伏较大、地形结构较复杂时，拟合结果与实际高程异常往往存在较大的偏差。多面函数拟合法能以任意精度逼近复杂曲面，但中心节点的选取较为繁琐。本文提出将蚁群算法引入多面函数建模，并考虑粗差的影响，对中心节点的数量和点位分布进行合理的选取，采用稳健估计剔除粗差对拟合模型的影响，保证所选取数据可靠的情况下提高拟合模型精度。

1 多面函数 1.1 多面函数基本原理

多面函数是一种将若干简单曲面叠加形成复杂连续曲面的方法^[3]。任何一个圆滑曲面都可以通过多面函数法以任意精度逼近，该方法非常适合于在地势起伏较大且地形条件复杂区域的GPS高程拟合。该方法在每个已知的数据点建立曲面，在方向上将各旋转曲面按照一定比例叠加成连续曲面，使其严格过各数据点，从而达到逼近圆滑曲面的效果^[4]。

设已知曲面函数为f(x, y)，若函数φ(x, y)满足方程：

$ \sum {{{\left( {f\left( {{x_i},{y_i}} \right) - \varphi \left( {{x_i},{y_i}} \right)} \right)}^2}} = \min $

(1)

则称φ(x, y)为曲面函数f(x, y)的逼近函数。将逼近函数φ(x, y)表示为：

$ \varphi \left( {x,y} \right) = \sum\limits_{j = 1}^u {{\beta _j}{Q_j}\left( {\left( {x,y} \right),\left( {{x_j},{y_j}} \right)} \right)} $

(2)

式中，β为待定系数，Q为核函数，u为核函数的数量，(x_j, y_j)为选取的中心节点坐标。

核函数Q可选择任意的简单函数，为了计算方便通常选取为对称型或距离型核函数。有3种常用核函数：

倒双曲函数

$ Q = {\left( {{{\left( {x - {x_j}} \right)}^2} + {{\left( {y - {y_j}} \right)}^2} + \delta } \right)^{ - 1/2}} $

(3)

正双曲函数

$ Q = {\left( {{{\left( {x - {x_j}} \right)}^2} + {{\left( {y - {y_j}} \right)}^2} + \delta } \right)^{1/2}} $

(4)

三次曲面型函数

$ Q = {\left( {{{\left( {x - {x_j}} \right)}^2} + {{\left( {y - {y_j}} \right)}^2}} \right)^{3/2}} + \delta $

(5)

式中，δ为核函数平滑因子，可调节核函数的形状，是影响函数模型精度的因素之一。

设有n组观测值(x_i, y_i, ζ_i)，i=1，2，…，n，建立误差方程：

$ \begin{array}{l} {v_1} = \sum\limits_{j = 1}^u {{{\hat \beta }_j}{Q_j}\left( {\left( {{x_1},{y_1}} \right),\left( {{x_{0j}},{y_{0j}}} \right)} \right)} - {\zeta _1}\\ {v_2} = \sum\limits_{j = 1}^u {{{\hat \beta }_j}{Q_j}\left( {\left( {{x_2},{y_2}} \right),\left( {{x_{0j}},{y_{0j}}} \right)} \right)} - {\zeta _2}\\ \;\;\;\;\;\; \cdots \\ {v_n} = \sum\limits_{j = 1}^u {{{\hat \beta }_j}{Q_j}\left( {\left( {{x_n},{y_n}} \right),\left( {{x_{0j}},{y_{0j}}} \right)} \right)} - {\zeta _n} \end{array} $

(6)

式中，(x_0j, y_0j)为中心节点的坐标，ζ为坐标对应的观测值。用矩阵表示为：

$ \mathit{\boldsymbol{V}} = \mathit{\boldsymbol{A\hat \beta }} - \mathit{\boldsymbol{\zeta }} $

(7)

式中，

$ \begin{array}{*{20}{c}} {\mathop {\mathit{\boldsymbol{A}}}\limits_{n \times u} = }\\ {\left[ \begin{array}{l} {Q_1}\left( {\left( {{x_1},{y_1}} \right),\left( {{x_{01}},{y_{01}}} \right)} \right), \cdots ,{Q_u}\left( {\left( {{x_1},{y_1}} \right),\left( {{x_{0u}},{y_{0u}}} \right)} \right)\\ {Q_1}\left( {\left( {{x_2},{y_2}} \right),\left( {{x_{01}},{y_{01}}} \right)} \right), \cdots ,{Q_u}\left( {\left( {{x_2},{y_2}} \right),\left( {{x_{0u}},{y_{0u}}} \right)} \right)\\ \;\;\;\;\;\; \cdots \\ {Q_1}\left( {\left( {{x_n},{y_n}} \right),\left( {{x_{01}},{y_{01}}} \right)} \right), \cdots ,{Q_u}\left( {\left( {{x_n},{y_n}} \right),\left( {{x_{0u}},{y_{0u}}} \right)} \right) \end{array} \right]} \end{array} $

(8)

利用最小二乘原理求得系数项：

$ \mathit{\boldsymbol{\hat \beta }} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{\zeta }} $

(9)

将$\mathit{\boldsymbol{\hat \beta }}$代入式(2)求得多面函数拟合模型，再根据拟合模型代入已知点位数据，即可完成高程异常拟合。

1.2 稳健权改进的多面函数

稳健估计是在数据存在粗差时，选用某种适合的估计方法使参数的估计值尽可能少受粗差的影响，得出最佳估计值^[5]。本文将稳健估计中的选权迭代法与多面函数结合，用于数据中含有粗差时，求解多面函数拟合参数的最佳估值，从而有效提高拟合模型的精度。

选权迭代法的关键是选择ρ函数或者φ函数，然后利用w_i=φ_i/v_i和P_i=p_iw_i关系式构造权因子w和等价权P，最后根据等价权进行迭代求解^[6]。构造稳健估计的权函数，ρ函数表示为：

$ \rho \left( v \right) = \left\{ \begin{array}{l} {v^2}/2,\left| v \right| < 1.5\sigma \\ \left| v \right|,1.5\sigma < \left| v \right| < 2.5\sigma \\ d,\left| v \right| > 2.5\sigma \end{array} \right. $

(10)

式中，v为残差值，d为任意常数，σ为单位权中误差，权函数ρ对应的权因子可表示为：

$ w\left( v \right) = \left\{ \begin{array}{l} 1,\left| v \right| < 1.5\sigma \\ 1/\left( {\left| v \right| + {k_0}} \right),1.5\sigma < \left| v \right| < 2.5\sigma \\ 0,\left| v \right| > 2.5\sigma \end{array} \right. $

(11)

式中，w(v)为取任意v时的权值大小，k₀通常取较小值。

在计算过程中，初始权矩阵设为单位权，再根据多面函数模型的误差方程式(7)，按照最小二乘法求解参数估值及其残差表达式：

$ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( 1 \right)}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{PA}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P\zeta }} $

(12)

$ {\mathit{\boldsymbol{V}}^{\left( 1 \right)}} = \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \beta }}}^{\left( 1 \right)}} - \mathit{\boldsymbol{\zeta }} $

(13)

根据式(13)并结合式(10)和式(11)，可求得等价权矩阵P。利用式(12)进行迭代计算，选定某一微小量ε，当第k次与第k-1次迭代所得的参数估值之差的绝对值小于给定的阈值ε时，结束迭代，即

$ \left| {{{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}} - {{\mathit{\boldsymbol{\hat \beta }}}^{\left( {k - 1} \right)}}} \right| \le \varepsilon $

(14)

满足条件后停止迭代，得出最后的估值结果：

$ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}{{\mathit{\boldsymbol{\bar P}}}^{\left( {k - 1} \right)}}\mathit{\boldsymbol{A}}} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{{\mathit{\boldsymbol{\bar P}}}^{\left( {k - 1} \right)}}\mathit{\boldsymbol{\zeta }} $

(15)

$ {\mathit{\boldsymbol{V}}^{\left( k \right)}} = \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat \beta }}}^{\left( k \right)}} - \mathit{\boldsymbol{\zeta }} $

(16)

通过式(15)的估值结果可知，迭代后的参数估值剔除了粗差对模型的影响。多面函数法的模型求解过程中，核函数的选择直接关系到拟合的效果，通常根据观测数据自身较强的相关性来选择核函数^[7]。平滑因子取不同的值会影响核函数的细节变化。徐南等^[8]将核函数加权系数和光滑因子共同参与遗传算法搜索，对光滑因子进行合理选取。中心节点位置及数量的选择是多面函数中的关键影响因素，点位选取的质量决定模型的精度。针对多面函数的特征点位选取比较繁琐和人工选取存在的缺陷，本文提出基于蚁群算法寻找多面函数中心节点的方法。

2 基于蚁群算法的多面函数

蚁群算法是通过模拟蚂蚁在寻找食物源过程中的正反馈原理寻找最短路径的仿生算法^[9]。蚂蚁寻找食物时，在经过的路径上释放易挥发的物质信息素(pheromone)，蚂蚁之间通过信息素进行相互通信^[10]。蚁群算法寻优过程分为3步：找出初始路径，构成完整规则；对构造的规则进行剪枝，去掉与实际条件不相关的节点；更新信息素浓度，为后续蚂蚁选择路径提供依据。蚁群算法寻优原理是蚂蚁在寻找完整路径时通过重复选择属性节点来实现的，属性节点的选取可以用随机方式进行。具体算法见文献[9-10]。

蚁群多面函数法是将蚁群算法用于多面函数建模时中心节点的选取，使所选择的特征点能够充分反映实际地形地貌特征。该方法的基本思想是：在已知数据覆盖的研究区随机分配一定数量的蚂蚁，给定相同的初始信息素浓度，设置信息素浓度挥发系数；然后让每只蚂蚁根据数据情况分别寻找地形复杂的山脊线和山谷线附近区域，通过构造初始规则寻找最优路径向山脊线或山谷线处靠拢；再根据信息素挥发后的浓度更新路径，确定地形复杂的特征点位置；最后汇集特征点。为避免局部区域的过度拟合，通过设置局部范围筛选密集特征点，选定地形复杂区域内均匀分布的中心节点，选取中心节点后构建拟合模型，得出最优估值。

3 算例分析

采用文献[2]中的算例数据构建模型，该地区观测数据覆盖面积约100 km²，研究区域地势起伏明显，海拔最高3 602 m、最低2 781 m。本文从收集的GPS水准重合数据中选择400个用于构建拟合模型，采用蚁群算法从建模数据中搜索最优中心节点，剩余部分数据用于检核。

3.1 蚁群多面函数法拟合GPS高程

对于地形复杂、起伏比较明显的地区，利用多面函数法求解拟合模型前需要找出研究区内最佳特征点。根据已知拟合数据数量，随机分配200只蚂蚁，设置迭代次数为80次，分别寻找地势起伏较大的特征点。通过蚁群寻优路径原理，蚂蚁快速寻找具有代表性的山脊线或山谷线上的点，通过迭代循环找出适用于构建拟合模型的最佳特征点。蚂蚁寻找特征点的前后位置分布见图 1。

由图 1看出，首先随机分配蚂蚁的初始位置，通过算法迭代及正反馈原理寻找特征点位后，蚂蚁向山脊线或山谷线靠拢，进一步缩小特征点的分布范围。通过设置一定阈值控制特征点的密集程度，找出山脊线和山谷线附近的分布点。为了使拟合效果更好，需要合理选取少量非特征点，采用均匀格网法补充选取变化平缓区域的中心节点，共同作为中心节点参与建模。

利用蚁群算法结合均匀格网找到研究区内的最优中心节点，用于多面函数构建拟合模型。因数据的坐标单位和高程异常的单位不统一，在参与建模之前需要对数据进行规格化处理，即完成坐标数据的预处理后可提高计算效率，同时也避免系数矩阵数值因过大而超界，降低计算效率。给定初始计算的权阵为单位阵，采用最小二乘原理解出系数矩阵，利用稳健估计的多面函数拟合模型迭代最优参数估值$\mathit{\boldsymbol{\hat \beta }}$后，代入式(2)反求已知坐标点位的高程异常，并与已知高程异常值进行比对。最后根据拟合所得高程异常值与已知拟合点的高程异常对比，分别绘出三维效果(图 2)。

蚁群算法结合稳健估计的多面函数拟合模型在剔除粗差影响情况下求得拟合模型的内符合精度σ₀=±0.468 7 mm，从图 2也可直观地看出拟合模型精度较高。完成模型构建以后，用高精度拟合模型拟合其余80个检测点，并对拟合结果进行精度分析。

3.2 结果分析

采用蚁群算法寻找多面函数拟合模型的中心节点，利用含稳健权的多面函数拟合模型求解系数矩阵$\mathit{\boldsymbol{\hat \beta }}$，评定拟合模型的内符合精度。为证明蚁群算法结合稳健估计的优势，将求得结果与相同条件下分别采用均匀格网法加入稳健估计前后拟合结果进行比对。表 1为采用不同方法分别统计10次的内符合精度并计算其平均值的对比结果。

由表 1可知，采用蚁群算法构建拟合模型所得结果精度高于均匀格网，加入稳健估计后所得拟合精度进一步提高。表明蚁群算法用于含稳健权的多面函数拟合模型有效可行，且模型精度高。

评定模型的内符合精度后，以80个已知点为检核数据，将该模型用于拟合检测点的GPS高程异常值，从检测点的拟合结果中随机抽取10组残差进行对比(表 2)。

由表 2可知，蚁群算法拟合结果波动较小，拟合精度趋于稳定；格网法拟合结果精度不太稳定，在选择多面函数中心节点的同时存在弊端，不能合理选取具有代表性的中心节点，从而很难提供稳定可靠的模型精度。

为合理评估蚁群算法结合稳健估计用于多面函数构建拟合模型的精度，求解80个检测点外符合精度，在大量实验的基础上，随机选取10次拟合计算结果与均匀格网结果进行对比(表 3)。

由表 3可知，蚁群算法结合稳健估计拟合结果的精度为±7.1 mm，用均匀格网加入稳健估计拟合结果精度为±9.6 mm，采用蚁群算法的总体精度比格网法提高26%。对于已知数据中含有粗差的情况，在保证拟合数据可靠的同时，也提高了拟合模型的精度及可靠性。为了更加形象直观地体现蚁群算法的有效性及稳定性，分别绘制采用蚁群算法和格网法拟合检测点的残差对比图(图 3)。

从图 3看出，格网法拟合结果残差波动性较大，稳定性较差；加入蚁群算法后拟合结果趋于稳定，且总体波动区间也有缩小，表明该方法对于精度的提高有效可行，并保证了拟合模型具有较好的稳定性。

4 结语

考虑到多面函数法构建拟合模型参数难以选择以及剔除粗差的影响，本文将蚁群算法与稳健估计同时引入多面函数拟合法。该方法能够对地势起伏较大的区域合理选取特征点，结合格网选取少量非特征点作为中心节点参与多面函数模型构建。分别采用蚁群算法和格网法加入稳健估计前后对比分析，得到不同情况下的拟合模型结果。分析表明，将蚁群算法和稳健估计同时用于多面函数求解模型的参数更加合理，且拟合模型的精度得到很大提高，同时保证了拟合模型的稳定性。