2. 江西省数字国土重点实验室，江西 南昌 330013
2. Jiangxi Province Key Laboratory for Digital Land,Nanchang 330013,China
1 引 言
大地测量数据处理领域中的病态问题是广泛存在的,例如大地测量反演[1 ，2]、控制网平差[3]、GPS快速定位[4]、航空重力向下延拓[5]等方面。当处理病态问题同时需要顾及系数矩阵误差时，即病态总体最小二乘问题的解算，是目前测量数据处理领域研究的热点问题之一[6, 7, 8, 9, 10, 11, 12, 13, 14, 15]；处理病态总体最小二乘问题的常用方法主要有截断奇异值解法[6, 15]、岭估计法[12, 13, 15]、Tikhonov正则化方法[7, 8, 9, 10, 11]、广义正则化解法等[14]。在解决病态总体最小二乘问题的这些方法中，主要是基于数学理论上的，虽然有些文献是从测量平差角度出发，但是在大地测量实际问题中的物理意义却不是那么的清晰，即存在以下的不足[16]：①偏重数学上算法的引用，对算法在大地测量实际问题中的物理意义的合理性未加重视；②过于强调如何缩小病态观测方程的系数阵(或法矩阵)的条件数，忽略了病态原因的分析和从机制上采取有效措施。虚拟观测法在最小二乘病态问题及其他准则带参数的模型中得到了广泛研究和应用[16, 17, 18, 19, 20, 21, 22]；准则参数表示的是实际观测方差与虚拟观测方差的比值，具有重要的物理意义；虚拟观测法在数据处理中的应用都可以得到与传统方法等价或较优的计算结果[22]。基于虚拟观测法的思想，本文提出了一种基于虚拟观测的总体最小二乘病态问题的岭估计法，该方法的公式推导中明确了岭参数的物理意义，通过数值算例表明虚拟观测法在处理病态总体最小二乘问题时是非常有效的。
2 病态总体最小二乘问题的虚拟观测法 2.1 病态总体最小二乘问题的数学模型设有平差线性模型
式中，L1为观测值向量;e1为观测值的噪声;A1是维数为n×m的系数矩阵;EA1为系数矩阵的噪声;X为m个未知参数，且
,
为EA1按列拉直得到的列向量。将式(1)改写为
求解上述方程的总体最小二乘方法可以表示为约束优化问题
约束条件L1+e1∈Range(A1+EA1) 式中，
是矩阵的Fronenius范数。
当系数矩阵的法方程N = A1TA1的条件数很大，即它的行列式接近零时，矩阵出现病态。病态性影响了模型估计结果的好坏，使得解的情况很不稳定。
2.2 病态总体最小二乘问题的虚拟观测解法病态的原因是设计矩阵存在复共线性关系,即观测值之间相关而产生秩亏或奇异现象。若参数之间是相互独立的，用观测方程表示如下，作为第2类观测
式中，L2=[0 0 … 0]T为虚拟观测；A2为单位阵；X为参数向量。e2 ～N(0,σ22Im)为虚拟观测的噪声。
将式(1)改写为
上式可以表示为
令
因此，式(1)可以表示为
根据式(8)，由协因数传播定律容易得到e的协因数阵为
令e的平差值为V1，e2的平差值为V2，X的平差值为
，Qe的平差值为
，将实际观测式(8)与虚拟观测式(4)联立得
在总体最小二乘准则下
由于e～N(0,σ12(1+XTX)In)，e2～N(0,σ22Im)，
，根据文献[22]知“准则参数表示的是实际观测方差与虚拟观测方差的比值”，所以
令
则有
由式(11)根据求函数自由极值的方法得
将式(4)、式(8)、式(9)代入式(14)得
即
令
则有 式(18)整理得
式(19)通过多次迭代运算可以得到较准确的参数解，该式与总体最小二乘岭估计式是一致的。
单位权方差估值可以按照下式计算
由式(11)、式(13)和式(17)得
式中，λL称为准则参数或岭参数;λ称为准则子参数。
通过上述岭估计虚拟观测法的推导可以看出，准则子参数λ是实际观测方差与虚拟观测方差的比值，所以通过虚拟观测法给参数λ赋予了恰当的物理意义。虽然，本文中应用的方法是将先验信息作为一项约束条件，与附有先验约束的间接平差理论上是等价的，但是，在附有约束的间接平差中，是直接运用拉格朗日求极值方法构造目标函数，而本文中是将约束条件与观测方程联立在原有估计准则的基础上得出新的准则，即式(21)。该准则式(21)与根据正则化理论处理病态总体最小二乘问题的估计准则是一致的。
病态总体最小二乘问题的虚拟观测法迭代计算步骤为:
(1) 采用最小二乘法求得待定参数的初始值
(2) 
(3)
(4) 重复第(2)、(3)步，当
时，迭代终止；ε可选一较小正值作为阈值。
在上述迭代算法中，准则子参数λ是实际观测方差与虚拟观测方差的比值，若两类观测值方差的先验信息已知，则可以根据式(13)求得λ，即验前单位权方差法[23]。当两类观测值方差的先验信息未知时，可以采用下面的岭迹法来确定该参数。
基本思想为：取0到一个相对大的数作为参数λ的取值区间，如本文算例1中的区间[0, 20]，选定一个尽量小的步长，如本文中的步长为0.001，让参数λ以该步长取遍区间内的所有值，将每个λ值代入迭代公式中进行迭代运算，迭代终止时每个λ值都对应一组参数解
，将的各个分量的岭迹画在同一张图上，即横坐标为λ值，纵坐标为的各个分量值，选取使各分量的岭迹大都稳定的那个点所对应的λ值作为准则子参数，此时的即为最优解。
采用文献[15]中病态设计矩阵为
未知参数有5个，X=[x1,x2,x3,x4,x5]T，真值为
。加入随机噪声，对于观测值L，其观测噪声e～N(0,σ02I)，σ0=0.1,I为单位矩阵。设计矩阵A中的元素与观测值之间的元素相互独立，且其误差eA=vec(EA)～
。随机误差由Matlab随机数发生器产生。加入随机误差后的设计矩阵和观测向量变为[15]
矩阵N=ATA的条件数为2.0837×104，病态性严重。为了比较不同情况下解的差异，分别用总体最小二乘法、TLS岭估计的L曲线法、最小二乘法、最小二乘岭估计的L曲线法以及本文提出的虚拟观测法进行解算，并比较计算所得的估值以及估计值与真值
差值的范数
的大小，计算结果如表 1所示；虚拟观测法的岭迹见图 1，图 2绘制了相应准则子参数λ所对应的差值范。
| TLS法 | L曲线法 (TLS)[15] | LS法 | L曲线法 (LS) | 虚拟 观测法 | 真值 | |
|
3.3041 | 1.1899 | 1.3944 | 1.2157 | 1.1567 | 1 |
| -2.8033 | 0.3903 | 0.1223 | 0.3728 | 0.4028 | 1 | |
| 0.0600 | 0.8122 | 0.7791 | 0.8280 | 0.7771 | 1 | |
| 3.5874 | 0.6109 | 0.2628 | 0.5983 | 0.6028 | 1 | |
| 2.9027 | 1.3060 | 1.4414 | 1.3157 | 1.2980 | 1 | |
| 参数λ | — | 0.1500 | — | 0.2991 | 9.7210 | — |
 | 6.7322 | 0.8290 | 1.3088 | 0.8547 | 0.8231 | — |
|
| 图 1 虚拟观测法的岭迹图 Fig. 1 The ridge mark of virtual observation method |
|
| 图 2 准则子参数对应的差值范数图 Fig. 2 The norm of difference with sub-criterion-parameter |
上述算例中仅已知实际观测的先验信息，而未知虚拟观测的先验信息值。从式(4)中可以知道，虚拟观测向量的方差等于参数的方差，为此，利用虚拟观测法得出的结果代入式(20)中得到参数单位权方差的估值，亦即虚拟观测值的单位权方差。根据式(13)可以得到参数λ的估值，再利用迭代公式求解参数，结果如表 2所示。
| 方法 | 参数X的估值 | 参数λ |
|
| 先验值未知的虚拟观测法 | [1.1567 0.4028 0.7771 0.6028 1.2980] | 9.7210 | 0.8231 |
| 先验值已知的虚拟观测法 | [1.1640 0.4004 0.7854 0.6005 1.2997] | 8.4936 | 0.8234 |
采用病态测边网的算例[3, 15, 24](图 3)。模拟一个空间测边网。P1、P2、…、P9为9个已知点，其坐标与其到两个未知点P10、P11的观测距离如表 3所示，并假设未知点真值分别为(0，0，0)和(7，10，-5)，两个未知点之间的观测距离为d10,11=13.1078m，各距离为等精度观测，中误差为0.001m。要求根据19个观测距离确定两个未知点的坐标[24]。
| m | ||||||
| 点 号 | 坐标 | 观测距离 | ||||
| X | Y | Z | di,10 | di,11 | ||
| P1 | 23.000 | 10.000 | 0.010 | 25.07869 | 16.76517 | |
| P2 | -10.000 | 9.990 | 0.000 | 14.13451 | 17.71965 | |
| P3 | 35.000 | 10.010 | -0.010 | 36.41588 | 28.44294 | |
| P4 | 100.000 | 19.990 | 0.005 | 101.47943 | 93.16839 | |
| P5 | -36.000 | 10.005 | 0.000 | 37.36422 | 43.29905 | |
| P6 | 0.000 | 10.010 | -0.005 | 10.01004 | 8.60060 | |
| P7 | 56.000 | 9.995 | 0.010 | 56.99606 | 49.25618 | |
| P8 | -15.000 | 10.015 | -0.010 | 18.03590 | 22.55966 | |
| P9 | -1.7000 | 10.008 | 0.015 | 10.15063 | 10.04382 | |
算例中N=ATA的条件数为4.5851×103，严重病态。在计算过程中，未知点的坐标取近似值为[3, 15, 24](0.01，-0.01，0.02)和(7.01，9.99，-5.01)。同样，为了比较不同情况下解的差异，分别用总体最小二乘法、TLS岭估计的L曲线法、最小二乘法、最小二乘岭估计的L曲线法以及本文提出的虚拟观测法进行解算，并比较计算所得的估计值以及估计值与真值差值范数的大小，计算结果如表 4所示；虚拟观测法的岭迹见图 4，图 5绘制了相应准则子参数λ所对应的差值范数。
|
| 图 4 虚拟观测法的岭迹图 Fig. 4 The ridge mark of virtual observation method |
|
| 图 5 准则子参数对应的差值范数图 Fig. 5 The norm of difference with sub-criterion-parameter |
(1) 从表 1和表 4的解算结果可以看出，虚拟观测法优于L曲线(TLS)法，算例1中差值范数Δ分别为0.8231和0.8290；算例2中差值范数分别为0.0280和0.9170。
| TLS法 | L曲线法 (TLS)[15] | LS法 | L曲线法 (LS) | 虚拟 观测法 | 真值 | |
| -0.0560 | -0.0408 | -0.0481 | -0.0487 | 0.0051 | 0 |
| 0.0033 | 0.0418 | -0.0379 | -0.0317 | -0.0109 | 0 | |
| 19.8651 | 0.8534 | 8.9686 | -0.0199 | 0.0193 | 0 | |
| 7.1006 | 6.9837 | 7.0186 | 6.9549 | 7.0015 | 7 | |
| -0.0701 | 9.7598 | 5.2929 | 9.8471 | 9.9874 | 10 | |
| -4.8405 | -4.7736 | -5.0026 | -5.0722 | -5.0103 | -5 | |
| 参数λ | — | 0.2000 | — | 0.4353 | 46.9580 | — |
| 22.2726 | 0.9170 | 10.1290 | 0.1855 | 0.0280 | — |
(2) 虚拟观测法在获取较优结果的同时，具有更明确的物理含义，即参数λ是实际观测方差与虚拟观测方差的比值；参数λL较好地均衡了准则式(21)中实际观测和虚拟观测两部分的拟合残差；由式(19)可以看出，虚拟观测法通过包含参数λ的准则参数(岭参数)λL，较好地改善了法方程系数矩阵的病态性。
(3) 从上面两个算例的结果来看，L曲线(TLS)与虚拟观测法得出的参数估计的结果相差较小，算例1更为相近，但得到的参数λ的估值却相差很大。文献[22]指出：关于准则参数的确定目前已经开展了大量研究，存在一个共同的不足就是准则参数只是表示估计准则中观测相关部分及非观测部分之间平衡的一个模糊量，在客观上没有真实的值(真值)和物理(几何)含义，各种方法所确定的准则参数的数值表示的并不是同一个量的估计值，相互之间没有可比性[22]。
(4) 综合两个算例来看，虚拟观测法在求解病态总体最小二乘问题中不但能得到比较理想的结果，而且它的计算过程相比较而言更加简单快捷。
(5) 处理病态问题的总体最小二乘岭估计法和最小二乘岭估计法都比相应的普通总体最小二乘法和普通最小二乘法更好地降低了病态性对解的影响。
4 结 论本文提出了一种解决总体最小二乘病态性问题的方法——虚拟观测法。该方法是利用平差中参数之间的相互独立性作为先验约束条件，并利用虚拟观测方程形式表示出来当做第2类观测量与实际观测模型进行联立求解。在公式的推导中得到了准则参数的物理意义，并且通过算例表明了利用虚拟观测法求解病态总体最小二乘问题可以得到与传统方法等价(算例1)或较优的计算结果(算例2)，显示了虚拟观测法是一种十分有效的方法。参数λ的选取是虚拟观测法的关键，如何更加有效合理地选取该参数还需要作进一步研究。
| [1] | WANG Leyang, XU Caijun, WANG Jianjun. Research on Equality Constraint Inversion with Ill-posed Constraint Matrix[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5): 397-401. (王乐洋, 许才军, 汪建军. 附有病态约束矩阵的等式约束反演问题研究[J]. 测绘学报, 2009, 38(5): 397-401.) |
| [2] | WANG Leyang, XU Caijun. Ridge Estimation Method in Inversion Problem with Ill-posed Constraint[J]. Geomatics and Information Science of Wuhan University, 2011, 36(5): 612-616. (王乐洋, 许才军. 附有病态约束的反演问题的岭估计法[J]. 武汉大学学报:信息科学版, 2011, 36(5): 612-616.) |
| [3] | GUO Jianfeng. The Diagnosis and Process of Ill-conditioned Adjustment System. [D]. Zhengzhou: Information Engineering University, 2002: 45-47. (郭建锋. 测量平差系统病态性的诊断与处理[D]. 郑州:信息工程大学, 2002: 45-47.) |
| [4] | WANG Zhenjie, OU Jikun, LIU Lintao. Investigation on Solutions of Ill-conditioned Problems in Rapid Positioning Using Single Frequency GPS Receivers[J]. Acta Geodaetica et Cartographica Sinica, 2005, 34(3): 197-201. (王振杰, 欧吉坤, 柳林涛. 单频GPS快速定位中病态问题的解法研究[J]. 测绘学报, 2005, 34(3): 197-201. ) |
| [5] | JIANG Tao, LI Jiancheng, WANG Zhengtao, et al. Solution of Ill-posed Problem in Downward Continuation of Airborne Gravity[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(6): 684-689. (蒋涛, 李建成, 王正涛, 等. 航空重力向下延拓病态问题的求解[J]. 测绘学报, 2011, 40(6): 684-689.) |
| [6] | FIERRO R D, GOLUB G H, HANSEN P C, et al. Regularization by Truncated Total Least Squares[J]. SIAM Journal on Scientific and Statistical Computing, 1997, 18(1): 1223-1241. |
| [7] | GOLUB G H, HANSEN P C, O’LEARY D P. Tikhonov Regularization and Total Least Squares[J]. SIAM Journal on Matrix Analysis and Applications, 1999, 21(1): 185-194. |
| [8] | SIMA D M, HUFFEL S V, GOLUB G H. Regularized Total Least Squares Based on Quadratic Eigenvalue Problem Solvers[J]. BIT Numerical Mathematics, 2004, 44: 793-812. |
| [9] | GUO H, RENAUT R. A Regularized Total Least Squares Algorithm[C]//Proceedings of Total Least Squares and Errors-in-Variables Modeling. Amsterdam: Kluwer Academic Publishers, 2002: 57-66. |
| [10] | BECK A, BEN T A. On the Solution of the Tikhonov Regularization of Regularized Total Least Squares Problem[J]. SIAM Journal on Optimization, 2006, 17(1): 98-118. |
| [11] | YUAN Zhenchao, SHEN Yunzhong, ZHOU Zebo. Regularized Total Least-squares Solution to Ill-posed Error-in-variable Model[J]. Journal of Geodesy and Geodynamics, 2009, 29(2): 131-134.(袁振超, 沈云中, 周泽波. 病态总体最小二乘模型的正则化算法[J]. 大地测量与地球动力学, 2009, 29(2): 131-134.) |
| [12] | WANG Leyang, XU Caijun, LU Tieding. Ridge Estimation Method in Ill-posed Weighted Total Least Squares Adjustment[J]. Geomatics and Information Science of Wuhan University, 2010, 35(11): 1346-1350. (王乐洋, 许才军, 鲁铁定. 病态加权总体最小二乘平差的岭估计解法[J]. 武汉大学学报:信息科学版, 2010, 35(11): 1346-1350.) |
| [13] | WANG Leyang. Research on Theory and Application of Total Least Squares in Geodetic Inversion[D]. Wuhan: Wuhan University, 2011(王乐洋. 基于总体最小二乘的大地测量反演理论及应用研究[D]. 武汉: 武汉大学, 2011.) |
| [14] | GE Xuming, WU Jicang. Generalized Regularization to Ill-posed Total Least Squares Problem[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 372-377. (葛旭明, 伍吉仓. 病态总体最小二乘问题的广义正则化[J]. 测绘学报, 2012, 41(3): 372-377.) |
| [15] | LU Tieding, NING Jinsheng. Total Least Squares Adjustment Theory and Its Applications[J]. Beijing: China Science and Technology Press, 2011: 89-102. (鲁铁定, 宁津生. 总体最小二乘平差理论及其应用[M]. 北京:中国科学技术出版社, 2011: 89-102.) |
| [16] | FENG Guangcai, DAI Wujiao, ZHU Jianjun, et al. A Virtual Observational Approach to Ill-conditioned Problem[J]. Science of Surveying and Mapping, 2007, 32(2): 38-39. (冯光财, 戴吾蛟, 朱建军, 等. 基于虚拟观测的病态问题解法[J]. 测绘科学, 2007, 32(2): 38-39.) |
| [17] | OUYANG Wensen, ZHU Jianjun. Expanding of Classical Surveying Adjustment Model[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(1):12-15. (欧阳文森, 朱建军. 经典平差模型的扩展[J]. 测绘学报, 2009, 38(1):12-15.) |
| [18] | ZHU Jianjun, FENG Guangcai, DAI Wujiao. A Quasi Observation Approach for Semi-parameter Regression[J]. Geotechnical Investigation and Surveying, 2006, (9): 54-57. (朱建军, 冯光财, 戴吾蛟. 半参数模型解算的一种虚拟观测法[J]. 工程勘察, 2006, (9): 55-57.) |
| [19] | ZHAO Haitao, GUO Guangli, CHA Jianfeng. The Virtual Observation Algorithm with Restriction Condition Parameter Adjustment[J]. Science of Surveying and Mapping, 2008,33 (2):33-34. (赵海涛,郭广礼,查剑锋. 附有限制条件间接平差的虚拟观测算法[J]. 测绘科学,2008,33 (2):33-34.) |
| [20] | LE Kejun, QIU Bin, ZHU Jianjun. Semi-parametric Regression Model Constrained with Priori Information and Its Virtual Observation Solution[J]. Journal of Geodesy and Geodynamics, 2008, 28(6): 107-111. (乐科军,邱斌,朱建军.具有先验信息约束的半参数回归模型及其虚拟观测值解法[J].大地测量与动力学, 2008, 28(6): 107-111.) |
| [21] | DAI Wujiao, FENG Guangcai, ZHU Jianjun. A Method for Selecting Ridge Parameter Based on Helmert Variance Components Estimation[J]. Journal of Geodesy and Geodynamics, 2006, 26(4): 30-33. (戴吾蛟, 冯光财, 朱建军. 一种基于Helmert方差分量估计的岭参数确定方法[J]. 大地测量与动力学, 2006, 26(4): 30-33.) |
| [22] | ZHU Jianjun, TIAN Yumiao, TAO Xiaojing. United Expression and Solution of Adjustment Criteria with Parameters[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(1): 8-13. (朱建军, 田玉淼, 陶肖静. 带准则参数的平差准则及其统一与解算[J]. 测绘学报, 2012, 41(1): 8-13.) |
| [23] | WANG Leyang, XU Caijun. Total Least-squares Adjustment with Weighting Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2011, 36(8): 887-890. (王乐洋,许才军. 附有相对权比的总体最小二乘平差[J].武汉大学学报:信息科学版, 2011, 36(8): 887-890.) |
| [24] | GUI Qingming, ZHANG Jianjun, GUO Jianfeng. Shrunken Type-Robust Estimation[J]. Acta Geodaetica et Cartographica Sinica, 2000, 29(3): 224-228. (归庆明, 张建军, 郭建锋. 压缩型抗差估计[J]. 测绘学报, 2000, 29(3): 224-228.) |