传统点质量模型中，虚拟点通常等间隔地分布在某一球面上(单层模型)或多个球面上(多层模型)^[1-8]。传统方法对数据格网化要求较高，当输入值离散分布时无法充分顾及输入值的空间分布特性，影响解算精度^[9]，而且当输入数据之间距离过近时，会使得设计矩阵中系数相似程度极高，影响方程稳定求解。针对该问题，Lehmann提出了一种自由点质量模型^[10]，该方法将扰动质点的位置、埋深以及大小看作待求量，该算法本质上是一种非线性模型，理论上其解是最优的。但是此种方法求解困难，而且最后的解受初始值影响较大^[11-12]。吴怿昊等^[13]在格网化的基础上将虚拟质点的埋深和分辨率看作变量，通过“标靶”法确定模型的最优结构，但质点模型格网化带来的系统性问题仍未解决。

针对以上问题，本文建立了一种简单半约束点质量模型，该模型可以利用少量离散数据进行建模。此外利用迭代过程中的残差特性，建立了固定加权点半约束点质量模型和自适应半约束点质量模型。并利用模拟数据对以上3种模型在局部(似)大地水准面拟合领域中的适应性进行了验证。

1 传统点质量模型

假设扰动质点M_j(j=1…m)均匀分布在半径为R_j的Bjerhanmmar球面上^[2], 各变量之间关系如图 1所示。

图 1中，P点为待求点，P₀, P₁…P_n为已知重力场元点，O代表地球球心，ρ为地球半径，l_pj为质点到计算点距离。

待求点P的扰动位T_P为：

$ T_{P}=G \sum\limits_{j=1}^{m} \frac{M_{j}}{l_{p j}} $

(1)

高程异常ζ可由布隆斯公式求得：

$ \zeta \approx \frac{T_{P}}{\gamma} $

(2)

式中，γ为P点平均正常重力。联立式(1)和式(2)可得高程异常(为表示方便，将公式中的约等于写为等于)：

$ \zeta=G \sum\limits_{j=1}^{m} \frac{M_{j}}{l_{p j} \gamma} $

(3)

2 简单半约束点质量模型

半约束点质量模型是将扰动质点固定于实测GPS/水准点的下方，即固定扰动质点的平面位置，其埋深作为待求量^[11]。模型如图 2所示。

图 2中, G_i和G_j为地表两个实测GPS/水准点，设定G_i下方有一个扰动质点，D_ij为两点之间的平面距离，R_im为G_i点下方的第m个质点M_im的埋深(m=1, 2, 3…M_iN)，M_iN为G_i下的扰动质点的总数。

将式(3)改化为：

$ \zeta_{j}=G \sum\limits_{i=1}^{G_{N}} \sum\limits_{m=1}^{M_{i N}} \frac{M_{i m}}{\sqrt{R_{i m}^{2}+D_{i j}^{2} \gamma}} \quad\left(j=1, 2, 3 \cdots G_{N}\right) $

(4)

式中，G_N为拟合点个数；ζ_j为所有质点对G_j点影响(j=1, 2, 3…G_N)。

由式(4)可得，当质点位于计算点下方时影响最大，利用这一特性展开迭代。同时在迭代过程中采用了Coredell解法，即认为相邻最近点之间的关系最为可信^[14]。

质点M_im对G_j的影响ζ′_j为：

$ \zeta_{j}^{\prime}=\left\{\begin{array}{l}{G \frac{M_{i m}}{\sqrt{R_{im}^{2}+D_{i j}^{2} \gamma}}, j \neq i} \\ {G \frac{M_{i m}}{\sqrt{R_{i m}^{2} \gamma}}, j=i}\end{array}\right. $

(5)

根据式(5)给出如下迭代步骤：

1) 优先对最大值点进行拟合：遍历观测数据，寻找绝对值最大的点(如图 2中G_i，记为遍历点)，设定其正下方某一深度处有一扰动质点M_im；

2) 确定距离该点距离最近的点(如图 2中G_j)。根据最近点置信度最大原则，假设G_i和G_j点的高程异常是由M_im引起的，根据式(6)可得参数k:

$ k=\frac{\zeta_{j}^{\prime}}{\zeta_{i}^{\prime}}=\frac{R_{i m}}{\sqrt{R_{i m}^{2}+D_{i j}^{2}}}(0<k<1) $

(6)

3) 求解R_im：

$ R_{i m}=\frac{D_{i j} k}{\sqrt{1-k^{2}}} $

(7)

4) 根据式(5)求定M_im：

$ M_{i m}=\xi_{i} \sqrt{R_{im}^{2}+D_{i j}^{2}} \gamma / G $

(8)

5) 根据式(5)求定M_im对每一个实测点高程异常的影响；

6) 将实测值与步骤5)中的结果作差作为新值，设定限差进行迭代。

3 加权半约束点质量模型

本实验将置信点扩展到最近的Q个，利用反距离加权法^[15]对Q个置信点赋权。在此过程中，加权点个数Q的选择尤为重要，Q过小会导致拟合信息不足影响拟合效果，反之则可能会造成过度拟合，引起“龙格现象”^[16-17]。

3.1 固定加权点个数

利用“标靶”思想，步进Q并分析实验结果，进而确定区域最优的加权点。

1) 同简单半约束点质量模型中第1)步；

2) 确定距离该点距离最近的Q个点，求取k₁, k₂…k_Q，根据每一点距质点的距离设定权值p₁, p₂…p_Q；

3) 根据式(6)求解D₁, D₂…D_Q，根据式(7)确定Q个质点埋深；

4) 根据式(8)分别求得G_i正下方的Q个质点大小；

5) 根据式(5)并结合步骤2)中确定的p₁, p₂…p_Q，求定M_im对每一个已知点的影响，利用差值作为步骤1)中的已知值循环上述步骤进行迭代, 直到满足限差要求，求得加权点个数为Q时的拟合均方差δN_Q；

6) 设定加权点个数为Q′=Q+1，得到δN_Q′，比较δN_Q与δN_Q′的大小，选取拟合均方差最小的值作为最优加权点个数。

3.2 自适应确定加权点个数

在迭代过程中遍历点外部环境不同，将每一步迭代中的加权点个数都设定为定值并不是最优的方法^[18-19]。随着迭代的深入，高程异常残差会逐渐减小并趋于平滑，这时需要的加权点个数相对于迭代初期要少，此时固定加权点个数会造成计算冗余。

将§3.1每一步迭代过程中Q设定为变量。拟合点距离遍历点距离越远，该点权越小，与该拟合点相关的变量k、埋深R_im、高程异常ζ以及质量M_im绝对值越小。所以对每一个遍历点而言，在最优加权点个数之前，随着Q的增大，拟合互差呈现收敛的趋势。在§3.1中步骤2)用如下步骤替换(示意图见图 3)：

1) 将Q_i设为最小值2，并根据§3.1中计算步骤求出Q_i个质点的计算值ΔN_Qi；

2) 将Q_i+1=Q_i+1，并根据§3.1求出Q_i+1个质点的计算值ΔN_Qi+1；

3) 将ΔN_Qi与ΔN_Qi+1作差，得到δN；

4) 判断δN的标准差与限差m的大小，当标准差小于m时，可认定增加加权点不能提高模型精度，所以该遍历点的最优加权点个数为Q_i；反之则Q_i=Q_i+1继续计算，直到确定符合限差的Q；

5) 按照§3.1中步骤2)~5)过程和限差进行迭代。

4 数值实验

实验采用NGA公布的2 160阶EGM2008模型进行模拟实验(下同)^[20]。选取两块实验区进行实验(如图 3(a)和3(b)所示)，区域一高程落差120 m，平均高程133 m。模拟该区域GPS/水准点3 600个，169个作为建模数据(图 3(a)中圆点)，3 431个点作为模型检核数据。区域二地形变化剧烈，高程落差1 120 m，平均高程973 m，数据生成方法同上。

作为对比实验，传统点质量模型利用2 160/36阶EGM2008组合场模拟残差重力异常进行建模，模型采用经典模型参数^[7]。第1层模型为A模型，第1、2层组合为B模型，第1~3层组合为C模型，第1~4层组合为D模型。

4.1 简单半约束法效果分析

A~D 4种模型和简单半约束点质量(简称E)模型区域一计算结果见表 1。

表 1中，A~D模型存在一定的系统性偏差，这主要由传统点质量模型的截断误差造成。本文§2中提出的简单半约束点质量的计算差值均值较小，因为模型摒弃了传统质点层的概念，迭代过程中生成不同埋深质点，理论上将拟合因子拓展为无限个，利用其线性组合实现重力场元的全频域逼近。

对A~D模型进行系统差处理，检核差值如图 4所示，将部分检核点B~E模型的拟合差值扩大显示(图 4右上部)。表 1、图 4中，传统分层组合在一定的分辨率之内(5′)，拟合均方差随着质点层数的增多而变小，D模型相比于C模型精度提升十分有限，因为拓展到2 160阶的EGM2008对应的模型空间分辨率为5′，更密的点本质上是对其格网数据的平滑，故精度影响较小，而且D模型使得模型计算量扩大了近一倍，故此采用C模型进行对比实验。

区域一内相比于模型C, 本文提出的简单半约束点质量模型的拟合均方差减小了1.1 cm。表 1中，本文方法所需基础数据和生成质点数仅为模型C的3.53%和2.38%，这使得建模工作量减少了数十倍，计算效率得到数量级的提升。

利用相同方法对区域二的数据进行处理，得到结果见表 2。

相比于区域一，区域二重力场信息变化更为复杂，A~D模型拟合结果中存在最小3.34 cm的系统差。为方便观察计算结果聚散程度，去系统差后拟合差值如图 5所示。

分析表 2信息，区域二中半约束点质量法生成的质点数达到870个，约为区域一质点数的8倍。因为区域二内重力场变化更为剧烈，需要更多的质点来表征其特性。

E模型的均方差较C模型的均方差大1.4 cm, 说明在C模型消除系统差之后精度优于本模型。分析发现C模型中格网数据保证了检核点高程异常的精度，使得每个检核点周边都有数个已知点对其进行精度控制。E(红色)模型大部分检核点的拟合差值较蓝色更加集中，E模型中一部分检核点周围没有拟合点进行精度保障，造成图 5中存在部分红色长“毛刺”，误差可达20 cm，使得结果震荡剧烈。

重新优化区域二中拟合点分布，使其按照均匀格网排列(如图 6所示)。得到实验结果见表 3、图 7(模型记为F)。

从图 7和表 3可以看出，F模型均方差相比于E、C分别减小了约1.6 cm、0.3 cm，这说明基础数据格网化后建模效果变好，验证了图 5、表 3中关于拟合精度不佳原因的分析。比较图 7中的E、F模型，F模型中的“毛刺”点大幅减小，拟合精度较为平稳，这说明图 6中数据结构对检核点的控制能力比图 3(b)强，采用格网化输入数据建模精度优于离散化数据建模精度。而且，相比于E，F模型生成的质点数缩减近一倍，说明建模数据的合理分布还能达到简化模型的效果(以下实验区域二采用图 6数据分布形式)。

4.2 加权半约束法效果分析 4.2.1 固定加权点个数

根据§4.1中的思路，将两个区域的初始加权点个数各设为2，依次增多加权点进行实验，得到不同加权点下拟合效果和生成的质点总数，如表 4所示。

由表 4可知，随着加权点的增多，生成的质点数逐渐增多，计算效率受到一定的损失。区域一、二中均方差最小的加权点分别为4和3个，差值均方差相比于简单半约束点质量方法分别减小了30.87%和47.31%，这说明加权半约束点质量模型是可行的，且拟合精度有了一定的提升。此方法中加权处理在收敛中起到了平滑作用，减少了单个质点在逼近过程中可能出现的震荡现象。

4.2.2 自适应加权点个数

按照§3.2步骤，将每一步迭代过程中的加权点个数设定为变量，两区域内计算结果如图 8、图 9和表 5所示，其中固定加权点模型采用两个区域内最优加权点个数进行实验。

表 5中自适应加权法在区域一内精度(均方差)与表 5中最优加权点模型下相当，总质点数减少一倍有余；区域二中自适应加权法精度获得了近3 mm的改善，质点数仅为原来一半。这是由于随着迭代的深入，拟合点高程异常整体趋于平滑，每一步迭代中所需加权点个数逐渐变少，此时仍固定个数势必会导致过度拟合和计算冗余，自适应方法将每一步迭代过程中的加权点设为自适应变量，结合实验结果说明自适应方法是合理的、可行的，而且可以简化模型。

5 结语

本文对传统的点质量模型进行改化，提出简单、固定加权点和自适应半约束点质量模型，并采用不同条件的模拟数据进行实验，得到如下结论：

1) 本文提出的简单半约束点质量模型可以用于局部(似)大地水准面拟合中，相比于传统的点质量方法具有无系统性偏差、建模数据少、计算效率高等优点；

2) 固定加权点半约束点质量模型在精度方面优于简单半约束点质量模型方法，但计算速度略有损失；

3) 自适应加权半约束点质量模型在保证精度的前提下改善了固定加权点半约束点质量模型的计算效率；

4) 本文提出的3种半约束点质量模型利用迭代的方法代替了传统模型中的矩阵求逆过程，避免了设计求逆过程中的不适定问题，对重力场其他领域研究有一定的参考意义。