本文内容来源于《测绘学报》2024年第11期（审图号GS京(2024)2421号）

非线性高斯-赫尔默特模型（nonlinear Gauss-Helmert model, NGHM）是高斯-赫尔默特模型（Gauss-Helmert model, GHM）在非线性平差领域的拓展，可概括表达高斯-马尔可夫模型（Gauss-Markoff model, GMM）[1]和误差变量模型（errors-in-variables, EIV)[2]，已在三维坐标变换[3]、GNSS定位[4]和点云曲面拟合[5-6]等方面得到广泛应用。除显式函数外，NGHM模型函数还可能是隐式函数，如圆曲线拟合模型[7]。显式和隐式非线性函数平差模型的统一表达称为通用NGHM模型。

在测绘领域，已有大量研究人员围绕NGHM模型作了拓展及应用研究[8]。文献[9—10]分别通过将系数矩阵随机变量纳入随机模型的方式和变量投影方法将结构EIV转化为NGHM，实现了模型的概括统一。文献[11]基于NGHM建立了一种适用于大旋转角的三维坐标转换参数求解模型。该模型无须假设条件，可以获得任意旋转角的坐标转换参数。在粗差处理方面，文献[12]用敏感分析法改进了NGHM算法，提升了其粗差探测性能。考虑到方差分量估计的精度会影响抗差估计效果，文献[5]在NGHM模型下研究了融合最优不变二次无偏估计与期望最大化（expectation maximization, EM）算法的粗差探测法。文献[7]针对误差含学生分布的情况，讨论了采用约束最大似然估计和EM算法估计NGHM模型参数的方法。文献[13]讨论了基于EM算法的NGHM抗差估计方法。文献[14]提出了基于等价权原理的NGHM抗差估计算法。文献[6,15]研究了NGHM抗差估计在点云曲面拟合中的应用。文献[16—17]针对海下三维坐标估计问题讨论了NGHM抗差估计法。学者们还将其向序贯平差方面进行了扩展。文献[18]讨论了等式约束下的序贯NGHM模型。文献[19]针对多历元观测数据，提出一种基于广义整体最小二乘原理的NGHM序贯平差法，并以水位流量曲线拟合为例验证了该方法在计算效率方面的优势。此外，文献[20]还研究了病态NGHM的参数估计问题。

NGHM模型的解算通常采用泰勒级数将模型展开，省略掉二次及以上项，用最小二乘推导平差计算公式，根据需要选择迭代计算方法，如高斯-牛顿迭代。在进行迭代计算时，需要提供与真值相差较小的初始值，若相差较大，算法就可能出现不收敛状况。已有研究表明，与传统高斯-牛顿迭代法相比，同伦法收敛的初值范围更大，是大范围收敛方法[21]。

测绘领域的研究人员已成功将同伦方法引入非线性最小二乘平差。文献[22]提出了非线性最小二乘同伦方法，给出了GPS同伦非线性模型和算法，在初值精度较差的情况下，得到了比线性最小二乘求解精度更高的结果。文献[23]提出了一种同时适用于满秩和秩亏网非线性最小二乘平差的统一模型，在初始值精度较低情况下，同伦方法仍能高精度收敛到原方程的估值。文献[24]讨论了粗差影响下和等式约束下的同伦非线性最小二乘平差方法。文献[25]将同伦方法引入坐标转换，研究了基于同伦法的非线性坐标转换模型，扩大了近似值的选取范围。文献[26]将同伦方法用于求解任意旋转角度下三维坐标转换参数。文献[27—28]针对传统同伦法收敛速度慢的问题，讨论了改进切向量及步长的同伦方法，得到了收敛更快、更稳定的算法，并应用于测边网平差。文献[29]针对传统同伦方法局部收敛问题，提出了融合同伦函数与填充函数的非线性最小二乘平差法。文献[30]讨论了同伦方法在相位解缠中的应用，扩大了非线性最小二乘相位解缠的收敛范围。文献[31—32]针对不适定非线性最小二乘问题，讨论了正则同伦方法及算法。

上述研究表明，同伦方法已在以显函数为对象的数据处理中得到广泛应用。本文通过将同伦方法引入通用NGHM模型参数估计中，把研究对象从显式函数扩展至隐式函数，建立起适用于多种情形的统一算法，拓展同伦方法在非线性平差中的应用。

1　通用NGHM的高斯-牛顿解

通用NGHM的误差方程可描述为

(1)

式中，f=[f1f2…fk]T是k×1的向量函数；l是n×1的观测值向量，其真值为；v是n×1观测误差向量；x是m×1的参数向量；m≤k≤n。式（1）包含了显函数和隐函数情况，以及线性和非线性情况，具有较强的通用性。其对应的随机模型为

(2)

式中，E（v）表示观测误差的期望；Var（v）表示观测误差的方差矩阵表示单位权方差；Qll（n×n）表示观测向量的协因数矩阵。

尽管式（1）可以描述线性平差情况，但考虑到仅非线性情况才需要线性化，因此本文主要针对非线性函数下的高斯-赫尔默特模型进行讨论。该模型的解算，通常先用泰勒级数将其线性化，然后采用最小二乘法估计模型参数。

记是v的最佳近似值是x的最佳估计值。l0和分别是和x的近似值，Δl和为对应的改正量，则有，。将式（1）在l0和附件展开，则有线性化后的误差方程

(3)

式中，是k×n的矩阵；是k×m的矩阵；是k×1的向量；和分别表示f对和x的偏导数。

在最小二乘准则vTPllv=min的条件下，由经典平差结论可得估计参数和残差向量的公式

(4)

(5)

式中，NBB=BQllBT。验后单位权方差的估计式为

(6)

式中，权矩阵。估计参数的方差-协方差矩阵计算式为

(7)

当l0和与真值非常接近时，可直接采用附有参数的条件平差估计参数。但当它们与真值相差较大时，需要通过迭代计算来获得更准确的估计值。由式（4）可以看出，其迭代格式与高斯-牛顿迭代法类似，故称其为高斯-牛顿型解法，简称高斯-牛顿法。遗憾的是，当初值与真值相差太远时，该方法存在不收敛的情况。

2　通用NGHM同伦平差解

2.1　通用NGHM同伦解算方程

同伦方法是求解非线性方程组的有效方法，其优点是大范围收敛，对初值精度的限制相对较弱。虽然同伦方法具有更大的收敛范围，但并不代表同伦法一定具有全局收敛性，只有当函数模型具有凸性条件时，才具有全局收敛特性。本文方法在NGHM平差中引入同伦法，以便在高斯-牛顿法不收敛时，仍能获得相对可靠的估计值。

类似显函数中同伦最小二乘平差准则的描述[29]，将通用NGHM同伦平差准则定义为

(8)

式中，a是Rk上的一个已知固定向量，通常取，即参数的初值；0≤t≤1。与显函数中的准则不同，上述准则中v关于参数x的表达式，不一定有直接的显性表达式v=f（x），而可能是隐函数式f（v,x）=0。

将式（8）对参数x求偏导数，可得通用NGHM同伦平差函数

(9)

式中，是n×m的矩阵。由式（9）可见，当t=1时，，此时参数存在唯一解；当t=0时，，此时就是的解。表示关于的函数。将式（9）对参数t求偏导数有[23]

(10)

(11)

(12)

式中，Z是n×m的矩阵，其结构为

(13)

表示关于的二阶偏导数，是一个n×m×m的立体矩阵。考虑到二维矩阵比立体矩阵更符合传统记录方式，也更方便计算，将克罗内克积（Kronecker）和矩阵拉直运算引入上述求导过程，从而避免计算复杂的立体矩阵。

由于是m×1的向量，因此根据矩阵拉直运算公式有如下等式成立

(14)

式中，⊗表示克罗内克积；z=vec(ZT），是mn×1的向量，其结构为

(15)

记，并将z对求偏导数，则有

(16)

显然，G是一个mn×m的矩阵。根据式（16）的结论，式（12）可改写为

(17)

令，，由式（10）可得微分方程

(18)

通过上述转化，将对式（9）的解算转化为在初值下对微分方程式（18）的解算[23]，即

(19)

2.2　固定步长的预测校正公式

通常以式（19）为基础，采用预估-校正法追踪同伦曲线，以获得参数的最佳估计值。预估的方法有很多种，如Runge-Kutta方法和固定步长法等。虽然前者计算的预测点值更准确，但计算量大，而后者的计算相对简单，故本文采用后者来计算预测点值。

记为同伦曲线上一点，i为预测点序号，δt为参数t的衰减步长，为参数的增量。由式（18）可得[33]

(20)

故第i+1号点的预测值为

(21)

由于预测点的坐标值未必一定落在同伦曲线上，因此需要将预测点校正到同伦曲线上。牛顿迭代法是一种有效的校正方法，故本文采用该方法校正预测点。

记，j为第i+1预测点上的第j次校正，则牛顿校正公式可描述为

(22)

牛顿校正是一个迭代计算过程，需要设置一个终止条件，本文采用的校正终止条件为[34]

(23)

式中，εc是给定的阈值；、、表示用第i预测点对应的t值和第j次校正对应的计算的H、J和K。因此，从初始点出发，应用式（21）和式（22），追踪同伦曲线到t=0，就可得到参数的估计值。

2.3　3个关键量的计算公式

由式（11）和式（17）可知，要完成同伦曲线跟踪，需要知道向量，矩阵Z和G的计算公式。在显函数中，可以直接表达为关于估计参数的显性表达式，因此可以直接计算。但在隐函数中，没有关于的显性表达式，无法直接完成计算。这3个量的计算公式推导如下。

若已知估计值，则有成立。将其在处展开，有

(24)

式中，，。显然该表达式同条件平差的函数模型一样。根据条件平差的结论，可得

(25)

虽然C与前面B的结构是相同的，但计算时采用的值略有不同。C的计算采用原始观测值l，而B的计算采用估计值。同样地，L的计算也用原始观测值l，而不用估计值。

记，则可将表示为，它由一系列方程组构成，即

(26)

为方便推导，用f代替fk，并将其对第i个误差vi求偏导，则有

(27)

记，，，则由式（27）可以推得

(28)

将式（28）代入式（13），可得用函数模型偏导数表示的Z，即

(29)

由式（16）可知，G是关于参数的二阶偏导数。将式（27）对第r个参数求偏导数，则有

(30)

记，，。

则式（30）可简化为

(31)

由此可以计算出G中的每个元素。为节约篇幅，这里不再列出其具体结构。在式（27)—式（31）中，i=1，2，…，n；j=1，2，…，m；r=1，2，…，m。

将通过通用同伦平差算法计算的和代入式（6）和式（7），可计算出通用同伦平差的验后单位权方差和参数估计值的方差-协方差矩阵。

在上述推导过程中，有如下几点需要说明。

（1）在推导式（24）时，仅考虑了的一阶项，删除了二阶及以上的高阶项。根据文献[35]可知，当初始参数的精度足够准确时，仅考虑一阶项的方法在通过多次迭代后可以得到与顾及二阶项相当的计算结果。这种处理不会影响算法的收敛性，但是会影响同伦曲线追踪的路径，从而对计算结果产生细微影响。为消除这种影响，在同伦曲线追踪结束后，将所得参数作为高斯-牛顿迭代法的初值进行迭代计算，修正得到更为准确的参数估计值。

（2）在推导式（27）的过程中，隐含有一个观测量只出现在一个误差方程中的假设。但实际中，可能存在一个观测量出现在多个误差方程中的情况，此时式（27）是一个方程组，需要解方程组才能得到式（28)。

（3）在构建G矩阵时，因为是一个关于的显性表达式，用其来推导关于的二阶偏导数，比用f来推导该二阶偏导数更简便直接。

3　非线性高斯-赫尔默特通用同伦平差算法

通用高斯-赫尔默特模型的同伦平差算法（简记通用同伦平差法）的计算步骤描述如下（图1）：①计算衰减后的参数t（i+1）=t（i）-δt，如果t≤0，则终止迭代；②根据给定的参数，采用式（25）计算残差向量；③根据观测量l和参数，采用式（29）和式（31）分别计算矩阵Z和G；④根据前一步计算的Z和G，以及衰减步长δt，采用式（20）和式（21）计算预测值；⑤根据参数和t，采用式（17）和式（9）计算矩阵J和；⑥根据J和，采用式（22）计算校正后的参数值；⑦重复第④、⑤步，直至满足式（23）时停止迭代；⑧令重复第①至⑥步；⑨在上述迭代停止以后，将获得的和代入式（4）和式（5），根据设置的修正次数对参数和进行修正，并用式（6）和式（7）计算和。

图1

图1   通用同伦平差法的计算流程

Fig.1   The flowchart of general homotopy adjustment algorithm

由于该算法可看成是一个循环内嵌一个循环，再外接一个循环的结构，因此其计算的时间复杂度比仅有一个循环的高斯-牛顿迭代法要高很多。

4　算例分析

4.1　算例1　距离定位

本算例的试验数据来自文献[36]，其值见表1。在本算例中，假设坐标（X,Y,Z）不含任何随机误差，仅观测距离d含有随机误差。因此，距离定位的NGHM误差方程为

(32)

式中，l=[d1d2…dk]T；；v=[v1v2…vk]T。在本算例中，设置了两组初值，即和，分别是初值1和初值2。采用前述的高斯-牛顿法和通用同伦平差法对表1的数据进行平差计算。图2—图4描述了算例1中高斯-牛顿法的计算结果，表2描述了通用同伦平差结果。在本节讨论的前3个算例中，高斯-牛顿法的收敛条件设置为；通用同伦平差的固定步长设置为δt=0.000 1；校正终止条件中的εc设置为εc=10-6；第⑨步的修正次数设置为k=3。

表1   算例1的试验数据

Tab.1  Experiment data in example 1

编号	X	Y	Z	d
1	10	20	6	22.180 0
2	8	30	0	30.958 0
3	10	-20	6	23.174 9
4	8	-30	0	30.880 0
5	-10	-20	-6	23.070 0
6	-8	-30	0	30.440 0
7	-10	20	-6	22.700 0
8	-8	30	0	30.720 0

新窗口打开| 下载CSV

图2

图2   高斯-牛顿法估计的参数值（算例1）

Fig.2   Parameter values estimated from Gauss-Newton method (example 1)

图3

图3   高斯-牛顿法估计的验后单位权方差（算例1）

Fig.3   Posteriori unit weight variance from Gauss-Newton method (example 1)

图4

图4   高斯-牛顿法计算的估计参数的方差（算例1）

Fig.4   Variance of the estimated parameters from Gauss-Newton method (example 1)

表2   平差结果（算例1）

Tab.2  The results of adjustment (example 1)

参数	采用初值1		采用初值2
	高斯-牛顿法	通用同伦平差法	高斯-牛顿法	通用同伦平差法
Xu/m	不收敛	-0.214	不收敛	-0.214
Yu/m		0.124		0.124
Zu/m		0.549		0.549
		0.302		0.302
		1.078		1.078
		0.045		0.045
		4.013		4.013

新窗口打开| 下载CSV

由图2—图4可以看出，不论采用初值1还是初值2，高斯-牛顿法估计的参数、验后单位权方差和估计参数的方差都始终在一定范围内摆动，不能收敛于一个固定值。说明该方法在本算例中不收敛。

由表2可以看出，在两种初值条件下，通用同伦平差法的结果都是一致的。这说明通用同伦平差法在两种初值中都能够稳定收敛。

4.2　算例2　GPS伪距定位

本算例的试验数据取自文献[37]。伪距单点定位方程的NGHM误差方程为

(33)

式中，；。与算例1不同，在本算例中，GPS卫星坐标和距离dk都含有随机误差。本例同样设置了两个初值，即初值1：，，，，它们与文献[37]的设置相同；初值2：，，，。

图5列出了高斯-牛顿迭代前100次计算的二范数值。由图5可以看出，该值一直在震荡前进，且始终大于10-10。故当设置的收敛条件小于等于10-10时，算法不会收敛。

图5

图5   迭代计算中的值

Fig.5    in each loop

表3列出了通用同伦平差法结果。由表3可以看出，在不同初值条件下，通用同伦平差法的结果完全一致。值得注意的是，在试验过程中发现，当高斯-牛顿迭代设置的收敛条件大于10-9时，其结果与通用同伦平差的结果完全一致。

表3   通用同伦平差结果（算例2）

Tab.3  The results of general homotopy adjustment (example 2)

参数项	采用初值1		采用初值2
	高斯-牛顿法	通用同伦平差法	高斯-牛顿法	通用同伦平差法
Xu/m		-2 157 555.665		-2 157 555.665
Yu/m		4 380 372.532		4 380 372.532
Zu/m		4 081 041.770		4 081 041.770
		175 322.286		175 322.286
	不收敛	137 879 264.216	不收敛	137 879 264.215
		293 314 174.167		293 314 174.167
		336 342 488.338		336 342 488.338
		266 625 946.277		266 625 946.277
		240 957 927.393		240 957 927.393

新窗口打开| 下载CSV

4.3　算例3　平面圆曲线拟合

与前面的显函数模型不同，本算例采用隐函数模型。平面圆曲线是一个非常典型的隐函数拟合模型，已有大量的文献对其参数估计问题进行了研究[38-45]。圆曲线的标准方程为

(34)

对应的NGHM误差方程为

(35)

式中，l=[X1X2…XkY1Y2…Yk]T；；。假设在某圆曲线上测量了10个坐标点，它们的值见表4。本算例初值1为，单位均为m；初值2为，单位均为m。

表4   圆曲线上的坐标观测值

Tab.4  Coordinate observations on circular curves

点号	X	Y
1	1 651.034 5	2 546.311 1
2	1 147.604 0	2 837.157 6
3	575.052 4	2 736.103 3
4	201.365 7	2 290.801 4
5	201.407 8	1 709.210 3
6	575.008 2	1 263.745 4
7	1 147.459 4	1 162.778 5
8	1 651.067 4	1 453.676 8
9	1 849.899 3	1 999.983 2
10	1 651.360 8	2 546.350 9

新窗口打开| 下载CSV

表5列出了高斯-牛顿法和通用同伦平差法的估计值。由表5可以看出，在初值1条件下，两种方法都能收敛，通用同伦平差法的结果与高斯-牛顿法的结果一致。在初值2条件下，高斯-牛顿法发散，不再收敛，而通用同伦平差法算法仍然收敛，且与初值1条件下高斯-牛顿法的结果一致。

表5   估计结果（算例3）

Tab.5  Estimated results (example 3)

参数项	高斯-牛顿法		通用同伦平差法
	采用初值1	采用初值2	采用初值1	采用初值2
	1 000.009		1 000.009	1 000.009
	1 999.985		1 999.985	1 999.985
	849.993		849.993	849.993
	0.015	发散	0.015	0.015
	0.003		0.003	0.003
	0.003		0.003	0.003
	0.002		0.002	0.002

新窗口打开| 下载CSV

4.4　算例4　测边网平差

本算例的数据来自文献[29]中的测边网数据。设置的两组初值见表6。两种方法估计的坐标值见表7。在本算例中，通用同伦平差法的固定步长设置为δt=0.001；校正终止条件εc设置为εc=10-3；第⑨步的修正次数设置为k=20。高斯-牛顿法的收敛条件同前。

表6   参数初值（算例4）

Tab.6  Initial parameters (example 4)

坐标	采用初值1	采用初值2
XP1	8 990.000	1 000.000
YP1	890.000	2 000.000
XP2	8 500.000	3 000.000
YP2	900.000	4 000.000
XP3	9 000.000	5 000.000
YP3	800.000	6 000.000
XP4	8 700.000	7 000.000
YP4	1 200.000	8 000.000

新窗口打开| 下载CSV

表7   坐标估计值（算例4）

Tab.7  Adjustment results (example 4)

坐标	采用初值1		采用初值2
	高斯-牛顿法	通用同伦平差法	高斯-牛顿法	通用同伦平差法
	9 034.167 03	9 034.167 03		8 551.905 79
	907.528 22	907.528 22		534.879 53
	8 762.945 47	8 762.945 47		8 410.402 64
	1 124.473 48	1 124.473 48	发散	852.059 65
	9 221.056 88	9 221.056 88		8 407.059 58
	1 008.491 04	1 008.491 04		379.506 16
	9 031.113 45	9 031.113 45		8 129.022 37
	1 345.343 85	1 345.343 85		648.287 90
	0.001 87	0.001 87		0.001 87

新窗口打开| 下载CSV

由表7可以看出，在不同的初值条件下，两种方法得到的结果不同，高斯-牛顿法在初值2条件下发散不收敛，而通用同伦平差法是收敛的。尽管通用同伦平差法得到的两组估计值有些差异，但是通过计算两点间的距离可以发现，两种估计值计算的距离几乎相同，这说明在不同初值条件下，控制网整体发生变化，但相对关系没有改变。本算例也说明了通用同伦平差法不一定具有全局收敛性，存在局部收敛的可能。

由算例3和算例4可知，在初值1条件下，两种方法都收敛；在初值2条件下，高斯-牛顿法发散，而通用同伦平差法仍收敛。将初值1条件下同伦平差法在每个步长计数（第几个步长）下的迭代次数与初值2条件下的迭代次数作差，可以描绘出如图6和图7所示的迭代次数差与步长计数的关系图。

图6

图6   两种初值下通用同伦平差法的迭代次数差（算例3）

Fig.6   Difference in the number of iterations of the generalized homotopy adjustment method for initial 1 and 2 (example 3)

图7

图7   两种初值下通用同伦平差法的迭代次数差（算例4）

Fig.7   Difference in the number of iterations of the generalized homotopy adjustment method for initial 1 and 2 (example 4)

由图6、图7可以看出，在两种初值条件下，不同的函数模型、步长和收敛条件，通用同伦平差法的迭代次数差并不完全不同；相同的函数模型、步长和收敛条件，初值不同时，其迭代次数差并不完全等于零，故迭代次数也并不完全相等。

5　结论

将同伦法引入非线性高斯-赫尔墨特模型的解算中，推导了具有通用性的同伦平差解算公式，建立了非线性高斯-赫尔墨特通用同伦平差算法。由于本文讨论的通用同伦平差法是从t=1开始，按照固定步长向前预测一个值，然后采用牛顿迭代将该预测值校正到同伦曲线上，得到一个误差相对小的估计值。在该计算过程中，需要在每一个固定步长下反复迭代计算，存在计算耗时的问题。为提高计算速度，本文设计的通用同伦平差算法中引入了高斯-牛顿法。首先，将步长设置得更长，校正终止条件设置得更大，采用通用同伦平差法估计得到一组参数值（该值的精度相对较低）；然后，将该参数值作为初始值代入高斯-牛顿迭代法中，按照设定的修正次数进行迭代计算，从而得到一个更准确的估计值。从算例来看，这种处理方法具有可行性。通过这样的处理，减少了步长数量，减少了校正迭代的次数，在保证精度的前提下，提高了计算的效率，但是仍然存在计算效率不如高斯-牛顿法的不足之处。通过理论推导和试验分析，得到如下结论。

（1）在高斯-牛顿法收敛时的初值条件下，通用同伦平差法同样收敛且结果完全一致，证明其具有可行性。但通用同伦平差法的计算相对更耗时。

（2）在高斯-牛顿法不收敛时的初值条件下，通用同伦平差法仍可能收敛，并且收敛的初值范围更大。

（3）通用同伦平差法适用于隐函数和显函数，部分变量或全体变量含误差的情况，具有较强的通用性。

初审：侯   琳

复审：宋启凡

终审：金   君