1 引 言

整体最小二乘法(total least squares，TLS)的起源可以追溯到19世纪70年代^[1]，这种方法最初在统计领域开始流行，称为误差变量模型(error in variable，EIV)。 文献^[2]在20世纪80年代初提出“total least squares”这一术语^[2]，实质是从数字分析角度对线性EIV模型重新进行定义。至此以后，TLS法迎来了快速发展期，广泛地应用于系统识别领域^[3]、信号处理领域^[4]、天文学领域^[5]等。TLS平差法是最小二乘(least squares,LS)平差方法的一种自然扩展，在LS基础上考虑系数矩阵受误差干扰的情况，使平差理论更严密，因而越来越受到测绘领域的重视。但是，TLS法引入大地测量数据处理领域的时间相对较晚，1999 年文献^[6]应用其进行曲线拟合。随后，文献^[7]将TLS法引入GPS数据处理中以提高GPS的定位精度。2005年，文献^[8]应用TLS法进行地理统计分析，指出该方法优于加权LS法。文献^[9] 将TLS法应用于变形分析中，但其实质是用TLS法进行坐标变换。目前，TLS法(包含其扩展模型)的应用主要集中在坐标变换方面^{[10, 11, 12, 13, 14, 15, 16, 17, 18]}。TLS平差法还频繁地出现在直线拟合^[19]和线性回归分析^[20]参数估计中。

根据测绘学科自身的特点，TLS平差模型本身也在不断扩展。考虑到估计参数之间满足某种函数关系，文献^{[21, 22, 23]}提出约束TLS平差法;考虑到观测向量为多维情况，文献^{[11, 12]}提出多维TLS平差法；考虑到观测向量和系数矩阵间异方差情况，文献^[20]提出了加权TLS平差方法；针对文献^[20]的解算方法相对复杂的情况，文献^[24]提出了一种简化的加权TLS平差迭代算法。在坐标变换中系数矩阵可能含有重复变量，即结构TLS问题，除采用结构TLS法外也可以通过加权的方式处理^[16]。文献^[25]对含有随机参数的情况进行了详细讨论，并给出加权TLS解存在的必要条件。国内对TLS平差方法的研究也比较多，主要集中在三维坐标变换参数估计^{[26, 27, 28, 29, 30, 31]}、三维激光扫描数据处理^{[28, 32]}、应变参数反演^[33]、自回归参数估计^[34]，对TLS平差算法也进行了一些研究^{[35, 36, 37]}。在模型方面，文献^[38]提出了稳健TLS平差法；文献^[39]研究了广义正则化的TLS平差问题。

然而，上述研究都是在线性函数模型条件下进行的讨论，即系数矩阵的每个元素都是自变量(包括常数)。但是在测量平差中经常出现非线性的情况，例如非线性回归。此时，系数矩阵中的各元素不再是自变量，而是关于自变量的函数，即各元素为函数表达式。因此，有必要对非线性TLS问题进行讨论。尽管非线性EIV模型已经不是一个新问题，但是从测量平差和数据分析角度对非线性整体最小二乘问题(nonlinear total least squares，NTLS)进行讨论并不多见。本文采用拉格朗日极值条件式推导基于牛顿型解法的NTLS平差迭代算法，最后通过非线性回归和曲线拟合算例验证该方法在等精度和非等精度情况下的可行性和有效性。

经典整体最小二乘法的观测方程可以描述为

点击浏览公式

式中，b∈R^m×1是观测向量;e_b∈R^m×1是与观测向量相对应的误差向量；A∈R^m×n是含有误差的系数矩阵；E_A∈R^m×n是系数矩阵相对应的误差矩阵；ξ∈R^n×1是参数向量。尽管在解算TLS问题时其可认为是一非线性问题，但是从系数矩阵元素的结构和类型角度分析其应该是一线性模型，因为系数矩阵中的各元素都是自变量。经典TLS法的平差准则是

点击浏览公式

式中，P_b∈R^m×m和P_A∈R^mn×mn分别是观测向量和系数矩阵元素组成的列向量的权矩阵，在经典TLS平差中它们都是单位矩阵；e_A∈R^mn×1=vec(E_A)，vec(·)表示矩阵的拉直计算，它是将矩阵按列重新排列成一新的列向量。对于经典TLS问题的解算方法，最常用的是SVD分解法，但是当系数矩阵元素之间存在相关性时，SVD分解就不能获得最或然估计结果^[40]。如果上述权矩阵不是单位矩阵，最为常用的解算方法便是基于牛顿型的迭代算法^[20]。值得注意的是，对于按行独立的加权TLS法与附加参数的条件平差法获得的估计结果具有等价性^[41]。
如果系数矩阵元素不再是独立的自变量，而是关于某一自变向量a的函数，此时系数矩阵不再是线性关系，例如

点击浏览公式

此时，系数矩阵A的元素不再是自变量x本身而是关于变量的函数，即系数矩阵的各元素可以由自变矢量a计算得到。将这种TLS平差问题称为非线性整体平差(NTLS)。

根据经典TLS法的描述，NTLS平差的观测方程可以表达为

点击浏览公式

式中，a∈R^q×1是系数矩阵A的输入向量，该向量元素之间可以是相互独立的自变量也可以是相关量；e_a∈R^q×1是输入向量a对应的误差向量。非线性TLS极小条件可以描述为

点击浏览公式

式中，P_a∈R^q×q是输入向量的权矩阵。

NTLS问题可以看成是一个双非线性问题。为能够解算NTLS问题，将式(4)表达成观测误差e_b关于a的函数，即

点击浏览公式

将上式在a⁰处展开并整理，得

点击浏览公式

此时

点击浏览公式

是一个大小为m×q，关于ξ的非线性矩阵，需进一步进行线性化，将式(7)在ξ⁰处展开

点击浏览公式

略去二次极小项乘积整理后，得

点击浏览公式

在上述的公式推导过程中，隐含有如下的等式条件

点击浏览公式

以及。

现在，NTLS就变成在极值条件(5)和方程(10)以及方程(11)的共同条件下求取参数的最佳估值以及观测向量和输入向量的最佳近似值的问题。即

点击浏览公式

式中，θ=[Δξ Δa]^T；O和I分别是(m×n)的零矩阵和(q×q)的单位矩阵。

利用拉格朗日极值条件对方程组(12)构建如下的目标函数

点击浏览公式

对式(13)中各元素分别求导整理后，得

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

式中，“～”表示最佳近似值；“ ^ ”表示最佳估计值。将式(14)和式(15)分别代入式(16)和式(17)中整理后，有

点击浏览公式

点击浏览公式

将求得的拉格朗日参数和代入式(18)中可以求得

点击浏览公式

式中

对式(21)的左边求逆后，可以得到

点击浏览公式

式中，N=Y^TPY。

获得最佳参数估计以后，联合式(19)和式(20)以及式(14)和式(15)可以解得误差量的最佳近似值计算公式

点击浏览公式

点击浏览公式

因此，可以得到方差分量的估计公式

点击浏览公式

式中

点击浏览公式

根据上一节的推导，将NTLS平差的算法流程设计如下。

输入数据：观测向量b，系数矩阵输入向量a、权矩阵P_b和P_a以及迭代终止条件ε(本文各次计算均设置其等于1e-8)。

第1步，令a⁽⁰⁾=a，仅考虑观测向量b的误差，采用LS平差法求得参数ξ的初始估计值ξ⁽⁰⁾，并计算

第2步，采用前一步获得的初始参数值和已知值分别计算

然后代入公式

中计算参数θ的估计值，并将该估计值分解成Δξ⁽ⁱ⁺¹⁾和Δa⁽ⁱ⁺¹⁾两部分，并用公式ξ⁽ⁱ⁺¹⁾=ξ⁽ⁱ⁾+Δξ⁽ⁱ⁺¹⁾和a⁽ⁱ⁺¹⁾=a⁽ⁱ⁾+Δa⁽ⁱ⁺¹⁾计算新的最佳估计量。

第3步，重复第2步直到Δξ⁽ⁱ⁺¹⁾和Δa⁽ⁱ⁺¹⁾的2范数都小于指定的收敛条件ε，并根据式(25)计算验后方差估计值。

输出数据：、a以及验后方差估计值₀²。

本文所讨论的NTLS问题经常出现在非线性回归分析中。多种模型可以通过取对数的方式线性化，然后得到各种新的观测方程。其中，二次曲线回归模型是最为常见的一种非线性回归分析模型。同时，在非线性回归分析当中，各观测量通常都是在相同条件下的观测值，即等权。因此，首先考虑二次曲线回归模型，用其来验证本文给出算法的可行性。

二次曲线回归模型通常可以表达为

点击浏览公式

式中，x和y都是含有误差的观测值；c₁、c₂、c₃是回归系数。

假设在某处通过测量获得如表 1所示的20组观测数据，x和y坐标是等精度观测，中误差为0.05 m。

根据前面的讨论知道，在采用本文给出的NTLS平差迭代算法的过程中，在获得初始参数估计值以后，需要计算Ξ矩阵。在本次的回归系数估计中，该矩阵结构为

点击浏览公式

分别采用LS平差法和NTLS平差法对表 1中的观测数据进行回归系数估计，将估计得到的参数结果列于表 2中。

需要注意的是，在采用LS平差法的过程中，因为认为系数矩阵元素不含有任何误差，故直接采用x的值和元素表达式计算得到的结果作为系数据矩阵元素值。

从表 2获得的结果可以发现，NTLS法获得的估计值与真值之间的差值比LS平差法获得的估计值与真值之间的差值小很多，特别是截断参数的估计。这说明，本文给出的算法能够获得比LS平差法更好的估计结果。

上面的试验是在单独一次估计中获得的结果，为进一步说明NTLS法比LS法优越，将上面的试验进行改进，设计如下：

回归系数的真值仍然采用表 2中描述的值，对真值坐标附加上从0.001 m开始，步长为0.003 m，到0.1 m结束的中误差，在每个中误差情况下模拟500次。不同中误差情况下获得的参数估计值的平均值与真值之差同先验中误差值的关系如图 1所示。

图 1 参数估计值与真实值的差值 Fig. 1 Difference of the parameters between the estimated values and the true ones

图选项

从图中可以看出，当误差比较小时，LS平差法和NTLS平差法获得的参数估计结果非常接近。在本例中，当中误差小于0.013 m时，它们的估计结果几乎完全相同。同时发现高次项对应的系数，不管在什么误差情况下两种方法都能够获得几乎与真值完成相同的结果。随着误差的增加以及幂的降低，系数的估计结果与真值的差距越来越大，但是NTLS方法获得的结果总比LS法获得的结果更接近于真实值。

将模拟500次后获得的验后中误差平均值与先验的中误差值的差值随误差增长的变化关系描述在图 2中。

图 2 估计中误差与先验值的差值 Fig. 2 Difference of the MSE between the prior value and the estimated ones

图选项

从图 2的描述来看，LS平差法很难获得正确的验后中误差估计，相反，NTLS平差法能获得与先验值几乎一致的估计结果。因为从图上看，其估计结果与先验值的差一直都在零附近波动。这些都说明，在等权情况下，NTLS平差法处理系数矩阵非线性情况下的平差问题时比LS平差法更好。

前面对等权情况下的NTLS平差法进行了试验讨论，然而在测量数据处理中，非等精度观测的情况经常出现，有必要对非等权情况下的NTLS平差问题进行讨论。假设有如下的对数曲线模型

点击浏览公式

式中，a和b是待估计系数；x和y都是观测量。假设它们的观测精度分别为0.08 m和0.02 m。同样分别采用NTLS平差法和LS平差法进行估计。在NTLS平差的过程中，Ξ矩阵的结构应该为

点击浏览公式

两种方法都模拟500次，将获得的参数平均值以及平均值与设计真值的差值列于表 3。

表 3 两种平差方法结果比较 Tab. 3 Comparison of the estimated parameters between the classical LS and the NTLS adjustment

参数	LS平差法			NTLS平差法		真值
	估计值	绝对差		估计值	绝对差
â	2.341 91	0.001 91		2.340 98	0.000 98	2.34
	3.099 38	-0.000 62		3.099 68	-0.000 32	3.10

表选项

从表 3不难发现，NTLS平差法获得的参数估计值比LS平差法获得的估计值更接近于真值，这说明NTLS平差法比LS平差法在参数估计部分更有效。尽管在此取得平均值，但在试验过程中发现它同等权情况时一样，每个点的估计值都比LS平差法更接近真实值，因此，避免重复，此处就没有再将其列出。

在定权的过程中，先验单位权方差取为1 m²。LS平差法采用

点击浏览公式

计算验后单位权方差，而TNLS平差法采用式(25)计算。将模拟500次获得的单位权方差估计值作为y轴，模拟次数作为x轴绘制成图 3。

图 3 两种方法的估计单位权方差 Fig. 3 Estimated variance components of the classical LS and the NTLS algorithm

图选项

图 3中的横线表示两种方法获得的单位权方差估计平均值。此图说明，NTLS平差获得单位权方差估计值与输入的先验单位权方差值几乎完全一致；而LS平差获得的单位权方差估计值出现较大偏差，这主要是因为其没有考虑系数矩阵中的误差。这说明本文给出的算法是有效并合理的。

本文对NTLS问题进行了研究，并给出了一种迭代计算方法。可以得出以下几点结论：

(1) NTLS方法保证了系数矩阵中不同位置的非线性元素的同一自变量获得相同的误差改正，常数元素项不获得任何改正值，其平差理论相对严密。

(2) 本文给出的算法，对于初始值的要求相对较低，不需要非线性LS法来获取初始值，直接采用线性最小二乘结果即可。

(3) 不管是等权或是非等权情况，NTLS法获得参数估计结果都更接近真实值，方差估计结果也更接近先验值。

文章仅就一种算法进行了讨论，而对于NTLS问题的其他方面没有涉及，例如模型在测绘科学技术中更广泛的应用等，这是笔者需要进一步努力研究的方向。