1 引 言

TLS算法最初由文献^[1]运用于数值计算领域，文献^[2]将该算法推广，并进一步将TLS理论及算法运用于科学与工程领域。另一方面，在处理线性模型Ax≈b时可能面临设计矩阵存在病态的情况，这将导致因观测量的微小波动而造成模型参数求解的巨大波动。在LS算法基础上，文献^[3]提出的正则化理论以及文献^[4]提出的截断奇异值(TSVD)的办法通常被用来解决线性估计下LS的病态问题。根据文献^[1]及文献^[5]推导得到的TLS算法可知，TLS的求解过程为一个降正则化的过程，因此在TLS解算中，设计矩阵的病态不仅将导致求解参数的不稳定，并且随着病态性的加强，其求解的不稳定性也将逐步增大。文献^[6]利用广义奇异值分解(GSVD)的办法得到了TLS方法下的截断奇异值法(T-TLS)。基于Tikhonov正则化理论，文献^[7]推导了总体最小二乘正则化方法(R-TLS)。随后文献^{[8, 9, 10]}也相继推导了多种形式的R-TLS方法。在Tikhonov正则化方法中，求解参数不仅应满足残差最小的约束目标函数，还应满足解的二次范数不等式条件，从而避免求解参数因观测值的微小变动而造成的巨大波动。在该方法的约束条件中需人为引入一个正常数，进而使得解的二次范数满足不等式约束，但是正常数的引入并没有绝对的标准，也没有一个客观的参考值，因此其值的引入存在较强的主观性。文献^[11]首先考虑了利用误差限制的条件处理TLS的病态问题，并利用牛顿迭代法推导和证明了该算法的可靠性。国内已有不少学者将TLS方法引入到测量数据的处理中，文献^[12]利用加权总体最小二乘算法(WTLS)对平面靶标及球靶标进行拟合；文献^[13]利用TLS方法建立了更加可靠的坐标转换算法。但对病态总体最小二乘的分析及其在测量数据处理中的应用研究还不多，相关文献也较少，其解算方法也主要是基于文献^[1]推导的R-TLS方法及文献^[14]提出的L曲线法。文献^[15]推导了病态总体最小二乘模型的正则化迭代算法；文献^[16]基于平差模型推导了加权病态总体最小二乘的岭估计算法；文献^[17]推导了一种基于L曲线的广义病态总体最小二乘正则化算法。

基于Tikhonov正则化理论及文献^[8]和文献^[11]等的分析结果，本文探讨了利用误差限制下的TLS正则化方法，并结合文献^[18]提出的最小-最大(max-min,MM)残差法，推导出一种基于误差限的TLS正则化算法(max-min residual regularization TLS,MMR-RTLS)，以下简称R-RTLS。

2 总体最小二乘解及正则化算法 2.1 总体最小二乘方法

本文所采用的线性估计基础模型为

式中，A为设计矩阵且列满秩；b为观测值向量；x为待求参数。总体最小二乘求解约束条件可表达为 式中，A^#、b^#分别为设计矩阵及观测值向量真值；‖‖_F表示F范数。根据约束条件(2)，基础模型的TLS解为^{[1, 2, 5]} 式中，σ_n+1为增广矩阵[A b]的第n+1个奇异值；I为单位阵。如果A为病态矩阵，则上述问题称之为病态总体最小二乘问题。由式(3)可以看出，TLS的求解过程为一个降正则化过程。所以，对于病态总体最小二乘问题而言，其解的不稳定程度将超过对应的最小二乘解。

2.2 Tikhonov正则化

Tikhonov正则化理论已广泛运用于处理病态问题。在LS的基础上，Tikhonov正则化方法需满足式(4)的约束^[3]

式中，δ为给定正常数;L通常为基于一阶或二阶偏导数的约束矩阵。在式(4)中不等式约束条件有效的情况下，通过建立拉格朗日多项式可以得到Tikhonov正则化解x_λ应满足的更一般的最小约束表达 式中，α为正则化参数。因此Tikhonov正则化解x_α为 在式(6)的求解中，正则化参数α可由文献^[14]提出的L曲线法确定。这里式(6)的求解是基于式(5)的约束条件，而式(5)是在约束条件式(4)中不等式约束有效且正常数δ足够小的情况下才成立。由此可以看到，在Tikhonov正则化方法中正常数δ起到了重要的控制作用。

2.3 总体最小二乘的正则化

文献^[7]等利用Tikhonov正则化理论推导并证明了总体最小二乘的正则化方法(R-TLS)。基于Tikhonov正则化方法的约束条件，建立R-TLS的约束条件式(7)

通过文献^[7]的推导可知，R-TLS下的解x_RTLS应满足式(8)

式中

式中，μ为拉格朗日因子。

文献^[9]在Golub推导基础上，运用特征值的办法推导并证明了R-TLS下的解x_RTLS也同时满足式(9)

式中，λ_L、λ_I、δ的定义同式(8)。

文献^[10]证明了R-TLS下类似于式(5)的约束式

2.4 误差限的正则化方法

在式(4)与式(7)中，δ均可视为不等式约束条件的上限，因此在处理R-TLS问题时正常数δ的给定十分重要。但在实际问题处理中，δ通常是由经验假设给出，而真实的上限是不知道的，因此在实际求解中经常会遇到因δ给定不适所造成的求解过度平滑或未平滑的情况。但另一方面，由于观测量总是满足一定的精度，观测值误差导致的设计矩阵和观测向量的误差上限η和η_b是可以根据观测值精度预先确定的，基于式(2)可将误差限制表示为

将R-TLS的目标函数与不等式约束条件对换，在新的目标函数下将式(7)转换为 从式(12)可以看到，在该约束条件下去除了式(7)中所含有的正常数δ，取而代之的是2个约束条件不等式的约束上限，它们可由观测值精度先验确定，从而克服正则化方法中正常数不易正确给出的缺点。

3 误差约束 3.1 MM约束条件

根据式(12)中约束条件的表达可知，在观测精度一定时设计矩阵与观测值向量的误差范围可表达为

在实际的测量中，设计矩阵和观测值向量的取值可表达为图 1所示圆周范围内的某一个值，即圆周表示理论真值，半径为客观的观测精度误差最大值。用数学表达式表达为a_ij^#=a_ij+η_ij,b_i^#=b_i+η_bi，其中a_ij^#为理论设计矩阵的第i行j列的元素，a_ij为实际设计矩阵的第i行j列的元素，而η_ij为对应设计矩阵元素的观测误差，b_i^#、b_i、η_bi的含义同设计矩阵元素。本文只考虑等精度观测情况，不等精度观测通过加权等处理可转换为等精度问题。

对于基础方程式(1)的每一个解，都可以得到一系列的残差

在所有的解中筛选出一个解，该解使得式(14)的最大值(max)达到最小(min)，即

由式(15)可得 则式(15)的约束条件可转换为：

3.2 MMR-RTLS方法

在2.4节所讨论的基于测量误差的R-TLS方法中，将解的范数定为目标函数，利用观测值误差范围作为约束条件，这样可以有效地避免因正常数δ的不适所造成的影响；另一方面，在3.1节的分析中，基于观测值误差范围，利用MM约束条件的方法可以有效地控制求解值的稳定性。

同样在基于观测值误差范围的前提下，建立如下约束

式(18)可以近似转换为 建立目标函数式 显然式(20)为一凹函数，对x求导并且令其等于零，可以得到与式(8)一致的法方程，此时

注意上述求导过程中，η_b对x求导等于0，所以解中不再含η_b，但η_b决定了式(19)中的目标函数极值。设计迭代算法如下：

(1) 利用LS方法得到x_LS作为初始值。

(2) 根据观测量的精度及设计矩阵形式确定设计矩阵误差上限η。在实际运用中，由于观测受到外界不确定因素的影响，观测量的精度往往达不到仪器标称精度，所以在确定误差上限时可以适当考虑在标称精度的基础上略有加大。

(3) 利用初始值及误差上限设定值，根据式(21)分别求解λ_L和λ_I的初始值。

(4) 根据式(8)，求解正则化值。

(5) 判断‖Δx‖₂＜ε，若不成立则将得到值作为初始值回到第3步，直至收敛(Δx表示每次迭代中对x的改正数，ε表示一个小正数，比如取10^－6)。

4 算 例

使用文献^[19]“Regularization Tools 2007”中的病态试验数据来测试本文所探讨的病态总体最小二乘的正则化方法。在算例中同时对设计矩阵及观测值向量加入不同程度的随机扰动误差用于模拟观测误差，并分别使用上述几种方法求解未知参数。

算例1：采用文献^[19]“Regularization Tools 2007”病态模拟数据“phillips”,设计矩阵条件数的数量级为10⁵。设计矩阵的元素方差为0.01，观测值向量方差为0.005。图 2(a)给出Golub算法中给定不同正常数所得到的结果和本文方法得到的结果，图中R值即为上文中的正常数δ。图 2(b)给出了本文方法的结果与L曲线下R-LS及R-TLS的结果。

为对比不同方法所得到的结果相对真值的偏差及波动大小，表 1分别列出了几种方法求解的偏差范数‖x－x^#‖₂与相对偏差范数。表 1中R-RTLS为本文所提出方法结果，R-TLS为Golub方法在给定不同正常数R情况下的结果，R-TLS(L)与R-LS(L)为利用L曲线所得到的结果。

算例2：采用文献^[19]“Regularization Tools 2007”病态模拟数据“show”。设计矩阵条件数的数量级为10¹⁹。设计矩阵的元素方差为0.01，观测值向量方差为0.005。图 3和表 2的注释和表头说明同图 2和表 1。

算例3：采用文献^[19]“Regularization Tools 2007”病态模拟数据“show”。设计矩阵的元素方差为0.1，观测值向量方差为0.01。图 4和表 3中的注释和表头说明同图 2和表 1。

从以上3组算例可以看到，本文所提出的方法可以有效地处理总体最小二乘的病态问题。Golub提出的R-TLS算法，在正常数给定恰当时可以得到较好的解算结果，但当正常数给定超过或不足理论最佳值时，解算结果会出现过度平滑或未平滑状态，并且波动值较大，参见图 2(a)～4(a)。本文提出的算法在无需假设正常数的情况下，通过设定观测值误差上限的办法从而稳定求解，因Golub算法中正常数的真值并不知道，所以对比给定较为恰当的正常数值时，本文方法较之Golub算法的结果精度相当，并略好于该算法。算例1中正常数在0.6处达到可靠值，该解算结果较之文本方法所得结果有0.077的偏差。算例2中正常数恰当时，其结果与本文方法结果有0.788 5的偏差，而算例3中该偏差为0.672 6。在病态性程度加强且误差增大时，本文算法相对于Golub算法其结果提高了近40%和24.4%。其次，在L曲线法较为理想时，其值也与本文提出算法所得结果精度相当。最后，算例2与算例3对比可以发现，在误差值较大的情况下，本文所提出算法的稳定性更强，相比其他方法可以更好地保持原有数据的总体形态及趋势，其稳定性较好。

算例4：采用文献^[20]模拟空间测边网算例。P₁、P₂、…、P₁₀为10个已知点，其坐标具体数据略去。10个已知点到3个未知点P₁₁、P₁₂、P₁₃(假设模拟坐标真值分别为(0,0,0)，(68，-26,9)和(14,41，-11))的距离，以及3个未知点间的距离假定已通过测量得到。设各距离为等精度观测，中误差为±0.01 m。根据33个距离观测值确定3个未知点坐标。计算中3个未知点坐标近似值分别取(0.03 m,-0.025 m,0.01 m),(68.03 m,-25.97 m,8.98 m)和(14.04 m,40.97 m,-11.04 m)。该测边网所建立观测方程的系数阵A为严重病态，法矩阵条件数的数量级为10⁵。表 4列出各种算法结果及误差范数等。

从表 4中可以看到，本文提出的算法在实际边角网平差解算中效果较好。例4中给定的初始值并不理想，LS方法得到的解是完全不可靠的。R-LS(L)与R-TLS(L)方法所得到结果均有一定程度的提高，相对误差系数也达到相对可靠的水平，但运用新方法所得到解算结果明显好于该两组解算结果，其偏差水平达到了观测误差给定水平(厘米级)，相对偏差也有明显的提高。其次，Golub方法的解算结果表明，当正常数给定较为合适时，其解算结果与本文提出方法的解算结果相当，仅存在千分之几的偏差，这也证明了新方法的可靠性。

通过采用不同初始值实验可以发现(限于篇幅，具体数据略去)，在初始值给定较为理想的情况下，LS的偏差明显减小，L曲线下的R-LS及合适R 值下的R-TLS的解算精度也有明显的提高。本文提出的新方法的解算精度同样也有一定程度的提高，并且明显好于其他方法。在正常数给定适当的情况下，Golub方法得到的解算结果与本文新方法的解算结果相当，也仅有万分之几的区别。

5 结 论

本文通过分析Golub所提出的R-TLS方法(下称原方法)，并结合文献^[11]关于误差限分析的结果，提出了一种更为可靠的病态问题正则化求解方法(R-RTLS，下称新方法)。通过本文的公式推导和算例分析，可以得到以下几点结论：

(1) 新方法通过“观测量的精度”这一客观标准给定误差上限值，克服因病态性所造成的总体最小二乘解的不稳定。该方法避免了原有方法中需主观给定正常数的过程，使其解算过程更加合理。

(2) 在给定误差上限时，新方法充分考虑了因外界因素对观测值精度的影响，即在标称精度或设计精度上略有增加，使该方法更具可操作性，同时也避免了原有方法在给定正常数不适当时所造成解的剧烈波动。

(3) 新方法对初值精度的要求不高，其抗差性较原有方法有一定的提高，扩大了该方法的使用范围。