留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

二维直角建筑物规则化的加权总体最小二乘平差方法

李林林 周拥军 周瑜

李林林, 周拥军, 周瑜. 二维直角建筑物规则化的加权总体最小二乘平差方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
引用本文: 李林林, 周拥军, 周瑜. 二维直角建筑物规则化的加权总体最小二乘平差方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
LI Linlin, ZHOU Yongjun, ZHOU Yu. A Weighted Total Least Squares Adjustment Method for Building Regularization[J]. Geomatics and Information Science of Wuhan University, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
Citation: LI Linlin, ZHOU Yongjun, ZHOU Yu. A Weighted Total Least Squares Adjustment Method for Building Regularization[J]. Geomatics and Information Science of Wuhan University, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510

二维直角建筑物规则化的加权总体最小二乘平差方法

doi: 10.13203/j.whugis20160510
基金项目: 

国家自然科学基金 41274012

详细信息
    作者简介:

    李林林, 博士生, 主要从事加权总体最小二乘方法的研究。lilinlin@sjtu.edu.cn

    通讯作者: 周拥军, 博士, 讲师。yjzhou@sjtu.edu.cn
  • 中图分类号: P207;P237

A Weighted Total Least Squares Adjustment Method for Building Regularization

Funds: 

The National Natural Science Foundation of China 41274012

More Information
    Author Bio:

    LI Linlin, PhD candidate, specializes in weighted total least squares. E-mail: lilinlin@sjtu.edu.cn

    Corresponding author: ZHOU Yongjun, PhD, lecturer. E-mail:yjzhou@sjtu.edu.cn
  • 摘要: 针对基于遥感数据的二维建筑物的直角化问题,以建筑物边界点的坐标为观测值,以顾及边界正交限制条件的直线斜率和截距为参数,建立附有限制条件的变量误差(errors-in-variables,EIV)模型。考虑观测向量和设计矩阵相关的情况,给出了增广设计矩阵的协方差阵的计算方法,推导了附限制条件的通用加权总体最小二乘(weighted total least squares,WTLS)平差算法,以及近似精度评定算法和仅含二次型限制条件的WTLS平差方法。理论和算例分析表明,在建筑物重建问题中,附有限制条件的EIV模型比经典附有限制条件的Gauss-Helmert模型易于构建,所提的WTLS算法快速收敛速度快,对拓展WTLS平差方法的应用具有理论与实践意义。
  • 图  1  有限制条件和无限制条件的WTLS平差方法的矩形房屋重建效果

    Figure  1.  Performance of Constrained and Unconstrained WTLS Algorithms for Building Reconstruction

    表  1  坐标观测值及其误差

    Table  1.   Measured Coordinates and TheirCovariance Matrices

    线段 编号 x/m y/m σx/m σy/m ρ
    AB 1 11.829 2 11.085 1 0.2 0.2 0.158 6
    2 13.257 0 11.898 5 0.2 0.2 0.363 6
    3 14.89 61 12.669 0 0.2 0.2 0.067 6
    4 16.643 2 13.889 6 0.2 0.2 0.480 5
    5 18.133 2 14.782 9 0.2 0.2 0.291 8
    6 19.569 3 15.874 3 0.2 0.2 -0.347 4
    7 20.889 9 16.326 2 0.2 0.2 0.333 0
    8 22.476 9 17.493 4 0.2 0.2 -0.308 1
    9 24.000 3 17.961 9 0.2 0.2 0.139 0
    10 25.810 7 18.986 1 0.2 0.2 0.169 0
    BC 1 26.873 0 21.542 0 0.4 0.4 0.272 1
    2 25.429 1 23.148 4 0.4 0.4 -0.120 2
    3 25.855 6 24.298 1 0.4 0.4 -0.058 4
    4 23.367 4 25.526 9 0.4 0.4 -0.016 9
    5 23.210 0 27.384 0 0.4 0.4 0.108 1
    CD 1 21.082 6 27.340 3 0.6 0.6 -0.324 0
    2 19.306 1 27.328 6 0.6 0.6 -0.498 0
    3 18.642 1 26.241 6 0.6 0.6 0.290 2
    4 16.225 0 25.370 7 0.6 0.6 0.013 6
    5 14.494 1 25.455 8 0.6 0.6 -0.286 8
    6 12.768 9 23.845 0 0.6 0.6 -0.396 6
    7 11.334 2 22.192 1 0.6 0.6 -0.342 7
    8 10.277 9 21.232 1 0.6 0.6 -0.092 5
    9 8.124 3 21.55 10 0.6 0.6 -0.092 2
    10 6.785 1 19.495 8 0.6 0.6 -0.447 3
    DA 1 5.636 7 16.556 6 0.8 0.8 0.441 8
    2 6.677 8 16.086 5 0.8 0.8 -0.350 0
    3 6.860 9 14.313 9 0.8 0.8 -0.115 6
    4 7.999 9 14.213 5 0.8 0.8 -0.188 9
    5 10.213 6 9.386 0 0.8 0.8 -0.331 5
    下载: 导出CSV

    表  2  附限制条件的WTLS方法与边界单独平差结果的比较

    Table  2.   Comparison of Constrained WTLS Method and Adjustment Individually

    线段 附二次型限制条件的WTLS法 无限制条件的WTLS法 真值
    $\mathit{\hat a}\ $ $\mathit{\hat b}\ $ Δa Δb $\mathit{\hat a}\ $ $\mathit{\hat b}\ $ Δa Δb a b
    AB 0.575 6 4.288 4 -0.001 8 0.061 9 0.579 6 4.245 3 0.002 2 0.018 8 0.577 4 4.226 5
    BC -1.737 4 67.705 1 -0.005 3 0.384 6 -1.431 4 60.118 4 0.300 6 -7.202 1 -1.732 1 67.320 5
    CD 0.575 6 15.976 9 -0.001 8 0.203 4 0.549 5 16.342 5 -0.027 9 0.569 0 0.577 4 15.773 5
    DA -1.737 4 27.201 0 -0.005 3 -0.119 6 -1.632 1 26.325 5 0.100 0 -0.995 0 -1.732 1 27.320 5
    下载: 导出CSV
  • [1] Markovsky I, van Huffel S. Overview of Total Least Squares Methods[J]. Signal Process, 2007, 87(10):2283-2302 doi:  10.1016/j.sigpro.2007.04.004
    [2] Markovsky I, Rastello M, Premoli P, et al. The Element-Wise Weighted Total Least Squares Problem[J]. Computational Statistics & Data Analysis, 2006, 50(1):181-209 http://dl.acm.org/citation.cfm?id=1648009
    [3] Schaffrin B, Wieser A. On Weighted Total Least-Squares Adjustment for Linear Regression[J]. Journal of Geodesy, 2008, 82(7):415-421 doi:  10.1007/s00190-007-0190-9
    [4] Shen Y Z, Li B F, Chen Y. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85(4):229-238 doi:  10.1007/s00190-010-0431-1
    [5] Mahboub V. On Weighted Total Least-Squares for Geodetic Transformations[J]. Journal of Geodesy, 2012, 86(5):359-367 doi:  10.1007/s00190-011-0524-5
    [6] Snow K. Topics in Total Least-Squares Adjustment Within the Errors-in-Variables Model: Singular Cofactor Matrices and Prior Information[D]. Ohio: Ohio State University, 2012
    [7] Fang X. Weighted Total Least Squares:Necessary and Sufficient Conditions, Fixed and Random Parameters[J]. Journal of Geodesy, 2013, 87(8):733-749 doi:  10.1007/s00190-013-0643-2
    [8] Zhou Yongjun, Kou Xinjian, Zhu Jianjun. A Newton Algorithm for Weighted Total Least-Squares Solution to a Specific Errors-in-Variables Model with Correlated Measurements[J]. Stud Geophys Geod, 2014, 58(3):349-375 doi:  10.1007/s11200-013-0254-7
    [9] van Huffel S, Lemmerling P. Total Least Squares and Errors-in-Variables Modeling:Analysis, Algorithms and Applications[M].Dordrecht:Kluwer Academic Publishers, 2002
    [10] Schaffrin B.A Note on Constrained Total Least-Squares Estimation[J].Linear Algebra Appl, 2006, 417(1):245-258 http://www.sciencedirect.com/science/article/pii/S0024379506002023
    [11] Golub G H, Hansen P C, O'Leary D P. Tikhonov Regularization and Total Least Squares[J].SIAM J Matrix Anal Appl, 1999, 21(1):185-194 doi:  10.1137/S0895479897326432
    [12] Beck A, Ben-Tal A.On the Solution of the Tikhonov Regularization of Total Least-Squares[J].SIAM J Optim, 2006, 17(1):98-118 doi:  10.1137/050624418
    [13] Schaffrin B, Felus A Y. An Algorithmic Approach to the Total Least-Squares Problem with Linear and Quadratic Constraints[J]. Stud Geophys Geod, 2009, 53(1):1-16 doi:  10.1007/s11200-009-0001-2
    [14] Mahboub V, Sharifi M A. On Weighted Total Least-Squares with Linear and Quadratic Constraints[J]. J Geod, 2013, 87(3):279-286 doi:  10.1007/s00190-012-0598-8
    [15] Zhang S L, Tong X H, Zhang K. A Solution to EIV Model with Inequality Constraints and Its Geodetic Applications[J]. J Geod, 2013, 87(1):23-28 doi:  10.1007/s00190-012-0575-2
    [16] Zeng W, Liu J, Yao Y. On Partial Errors-in-Variables Models with Inequality Constraints of Parameters and Variables[J]. J Geod, 2015, 89(2):111-119 doi:  10.1007/s00190-014-0775-z
    [17] Fang X. Weighted Total Least-Squares with Constraints:A Universal Formula for Geodetic Symmetrical Transformations[J]. J Geod, 2015, 89(2):1-11 doi:  10.1007/s00190-015-0790-8
    [18] 童小华, 邓愫愫, 史文中.数字地图合并的平差原理与方法[J].武汉大学学报·信息科学版, 2007, 32(7):621-625 http://ch.whu.edu.cn/CN/abstract/abstract1932.shtml

    Tong Xiaohua, Deng Susu, Shi Wenzhong.A New Least Squares Adjustment Method for Map Conflation[J]. Geomatics and Information Science of Wuhan University, 2007, 32(7):621-625 http://ch.whu.edu.cn/CN/abstract/abstract1932.shtml
    [19] Sampath A, Shan J. Building Boundary Tracing and Regularization from Airborne LiDAR Point Clouds[J]. Photogrammetric Engineering & Remote Sensing, 2007, 73(7):805-812 http://d.old.wanfangdata.com.cn/NSTLQK/NSTL_QKJJ022178597/
    [20] Xu P L, Liu J N, Shi C. Total Least Squares Adjustment in Partial Errors-in-Variables Models:Algorithm and Statistical Analysis[J]. J Geod, 2012, 86(8):661-675 doi:  10.1007/s00190-012-0552-9
    [21] Nocedal J, Wright S. Numerical Optimization[M]. Berlin:Springer, 2006
  • [1] 王乐洋, 孙坚强.  总体最小二乘回归预测模型的方差分量估计 . 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
    [2] 吕志鹏, 隋立芬.  基于非线性高斯-赫尔默特模型的结构总体最小二乘法 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1808-1815. doi: 10.13203/j.whugis20180104
    [3] 王乐洋, 余航.  附有相对权比的加权总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2019, 44(8): 1233-1240. doi: 10.13203/j.whugis20170265
    [4] 王乐洋, 余航.  火山Mogi模型反演的总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
    [5] 王乐洋, 李海燕, 陈晓勇.  拟牛顿修正法解算不等式约束加权总体最小二乘问题 . 武汉大学学报 ● 信息科学版, 2018, 43(1): 127-132. doi: 10.13203/j.whugis20150333
    [6] 王乐洋, 许光煜, 陈晓勇.  附有相对权比的PEIV模型总体最小二乘平差 . 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
    [7] 王乐洋, 许光煜.  系数矩阵误差对地壳应变参数反演的影响 . 武汉大学学报 ● 信息科学版, 2017, 42(10): 1453-1460. doi: 10.13203/j.whugis20160010
    [8] 王乐洋, 吴飞, 吴良才.  GPS高程转换的总体最小二乘拟合推估模型 . 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
    [9] 汪奇生, 杨德宏, 杨腾飞.  EIV模型参数估计的新方法 . 武汉大学学报 ● 信息科学版, 2016, 41(3): 356-360. doi: 10.13203/j.whugis20140182
    [10] 王乐洋, 许光煜.  加权总体最小二乘平差随机模型的验后估计 . 武汉大学学报 ● 信息科学版, 2016, 41(2): 255-261. doi: 10.13203/j.whugis20140275
    [11] 马友青, 刘少创, 魏士俨, 李明磊.  加权总体最小二乘的地面解析摄影测量算法 . 武汉大学学报 ● 信息科学版, 2015, 40(5): 594-598. doi: 10.13203/j.whugis20130387
    [12] 姚宜斌, 黄书华, 张良, 胡羽丰, 李国平.  求解三维坐标转换参数的整体最小二乘新方法 ! . 武汉大学学报 ● 信息科学版, 2015, 40(7): 853-857. doi: 10.13203/j.whugis20130302
    [13] 王乐洋, 许才军.  总体最小二乘研究进展 . 武汉大学学报 ● 信息科学版, 2013, 38(7): 850-856.
    [14] 刘经南, 曾文宪, 徐培亮.  整体最小二乘估计的研究进展 . 武汉大学学报 ● 信息科学版, 2013, 38(5): 505-512.
    [15] 周拥军, 邓才华.  加权和不加权TLS方法及其在不等精度坐标变换中的应用 . 武汉大学学报 ● 信息科学版, 2012, 37(8): 976-979.
    [16] 鲁铁定, 周世健.  总体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1351-1354.
    [17] 王乐洋, 许才军, 鲁铁定.  病态加权总体最小二乘平差的岭估计解法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1346-1350.
    [18] 徐天河, 杨元喜.  均方误差意义下正则化解优于最小二乘解的条件 . 武汉大学学报 ● 信息科学版, 2004, 29(3): 223-226.
    [19] 张祖勋, 张剑清, 胡翔云.  基于物方空间几何约束最小二乘匹配的建筑物半自动提取方法 . 武汉大学学报 ● 信息科学版, 2001, 26(4): 290-295.
    [20] 高士纯.  附有限制条件的间接分组平差模型与公式 . 武汉大学学报 ● 信息科学版, 1996, 21(1): 36-40.
  • 加载中
图(1) / 表(2)
计量
  • 文章访问数:  703
  • HTML全文浏览量:  61
  • PDF下载量:  252
  • 被引次数: 0
出版历程
  • 收稿日期:  2017-04-21
  • 刊出日期:  2019-03-05

二维直角建筑物规则化的加权总体最小二乘平差方法

doi: 10.13203/j.whugis20160510
    基金项目:

    国家自然科学基金 41274012

    作者简介:

    李林林, 博士生, 主要从事加权总体最小二乘方法的研究。lilinlin@sjtu.edu.cn

    通讯作者: 周拥军, 博士, 讲师。yjzhou@sjtu.edu.cn
  • 中图分类号: P207;P237

摘要: 针对基于遥感数据的二维建筑物的直角化问题,以建筑物边界点的坐标为观测值,以顾及边界正交限制条件的直线斜率和截距为参数,建立附有限制条件的变量误差(errors-in-variables,EIV)模型。考虑观测向量和设计矩阵相关的情况,给出了增广设计矩阵的协方差阵的计算方法,推导了附限制条件的通用加权总体最小二乘(weighted total least squares,WTLS)平差算法,以及近似精度评定算法和仅含二次型限制条件的WTLS平差方法。理论和算例分析表明,在建筑物重建问题中,附有限制条件的EIV模型比经典附有限制条件的Gauss-Helmert模型易于构建,所提的WTLS算法快速收敛速度快,对拓展WTLS平差方法的应用具有理论与实践意义。

English Abstract

李林林, 周拥军, 周瑜. 二维直角建筑物规则化的加权总体最小二乘平差方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
引用本文: 李林林, 周拥军, 周瑜. 二维直角建筑物规则化的加权总体最小二乘平差方法[J]. 武汉大学学报 ● 信息科学版, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
LI Linlin, ZHOU Yongjun, ZHOU Yu. A Weighted Total Least Squares Adjustment Method for Building Regularization[J]. Geomatics and Information Science of Wuhan University, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
Citation: LI Linlin, ZHOU Yongjun, ZHOU Yu. A Weighted Total Least Squares Adjustment Method for Building Regularization[J]. Geomatics and Information Science of Wuhan University, 2019, 44(3): 422-428. doi: 10.13203/j.whugis20160510
  • 近年来,总体最小二乘方法(total least squares, TLS),特别是加权总体最小二乘方法(weighted TLS, WTLS)在测绘领域受到了广泛的关注,已成为测量平差领域的热点研究问题之一。现有的WTLS方法以形如y形式的变量误差模型(errors-in-variables, EIV)为基础,以设计矩阵和观测向量所有元素的加权最小二乘为估计准则,是一种有别于经典平差的方法。

    针对无限制条件的EIV模型的TLS平差方法,文献[1-2]根据系数矩阵的权结构将TLS问题分为简单TLS、广义TLS、按元素独立的WTLS、附有限制条件的WTLS问题等。测绘领域研究最多的是WTLS问题,基于无限制条件的EIV模型,文献[3]针对设计矩阵的协因数阵QA可分解两个非奇异矩阵的Kronecker积的情形,提出了一种高效的WTLS方法;文献[4]针对WTLS问题提出了一种非线性迭代解法;文献[5]对一般形式的协因数阵QA提出了改进的WTLS方法,但未考虑Ay的相关性;文献[6-8]研究了Ay相关情况的通用WTLS方法。

    在实际问题中,观测系统中往往包含一些必要的限制条件如等式约束或不等式约束等。这些约束或附加必要的几何条件以避免系统只有零解的情况,或者减弱偏差以保证估值的稳定性。对于线性EIV模型,文献[9-10]研究了附有线性约束的简单TLS问题;文献[11-12]研究了附有二次型约束的简单TLS问题;文献[13]研究了附二次型限制条件的EIV模型的拉格朗日乘数迭代法。但以上方法均不能处理设计矩阵和观测向量相关的情况。文献[14]提出了一种同时附有线性约束和二次约束的WTLS方法,该方法未考虑观测值和设计矩阵相关的情况,且在某些情况下无法收敛;文献[15-16]研究了附不等式约束条件的EIV问题;文献[17]研究了附有限制条件的WTLS方法,并将其应用于三维坐标变换,但对参数初值的获取和精度评价不够完善。

    建筑物的直角化是GIS数据分析、三维激光扫描数据处理、三维重建中的基本问题,常用的方法是选择房屋角点坐标为未知参数,以直角约束条件为条件方程按经典平差方法求解[18],由于需要将坐标观测数据转换为角度条件并线性化,导致建模和解算较复杂,在观测数据和参数较多的情况下,其弱点尤为明显。文献[19]选择边界的直线方程参数作为未知数,将建筑物边界分别拟合,然后考虑边界的正交条件沿用传统的平差方法,该方法虽然较传统方法简洁,但未考虑观测值不等精度的问题,理论上也不严密。本文针对二维建筑物直角化问题,建立统一的附限制条件的EIV模型,并考虑观测值和设计矩阵相关的情况,推导出附限制条件的WTLS平差方法以及相应的近似精度评定方法,并通过一个矩形房屋重建算例进行验证。本文旨在为以激光探测与测量(light detection and ranging,LiDAR)或高分遥感为主要数据来源的二维建筑物高精度重建问题提供EIV建模和WTLS平差方法,拓展TLS平差方法在摄影测量与遥感领域的应用。

    • 设有如下的线性系统:

      $$ \mathit{\boldsymbol{\tilde y}} = A\left( {\tilde x} \right)\xi $$ (1)

      式中,$\mathit{\boldsymbol{\tilde x}}\mathit{\boldsymbol{、\tilde y}} $分别表示q×1、m×1阶的实数输入和输出向量的真值;$\mathit{\boldsymbol{A}}\left( {\mathit{\boldsymbol{\tilde x}}} \right) $表示m×p阶映射矩阵,即其元素为常数或向量${\mathit{\boldsymbol{\tilde x}}} $的函数;ξp×1阶的参数向量。在实际问题中,往往只能得到含误差的输入输出变量,设观测值满足以下关系:

      $$ \mathit{\boldsymbol{\tilde y}} = \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}},\mathit{\boldsymbol{\tilde x}} = \mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}} $$ (2)

      式中,xy分别表示输入和输出向量的观测值;exey分别表示相应的误差向量,并满足如下的随机模型:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{l}}}: \sim N\left[ {{\bf{0}},\sigma _0^2{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{l}}}} \right]}\\ {\mathit{\boldsymbol{l}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{x}}\\ \mathit{\boldsymbol{y}} \end{array}} \right],{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{l}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{x}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \end{array}} \right],{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{l}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{x}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yx}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{xy}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}} \end{array}} \right]} \end{array} $$ (3)

      式中,N[μ, D]表示数学期望为μ、方差为D的正态分布。很显然,A(x)的元素可以是观测值的线性或非线性形式,本文仅考虑A(x)为观测值向量的线性组合的情况。目前TLS方法应用中的直线拟合、仿射变换、Bursa变换与自回归模型参数估计等A矩阵中含误差的元素均为观测值的线性函数。在测量平差问题中,由于xy向量均含有误差的模型属于EIV模型,宜采用WTLS方法求解,现有的WTLS平差方法可表示为求解以下EIV问题:

      $$ \left. \begin{array}{l} \min :\frac{1}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}^ + \mathit{\boldsymbol{e}}\\ {\rm{s}}.\;{\rm{t}}.:\left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{\xi }} = \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} \end{array} \right\} $$ (4)

      式中,A表示含有随机误差EAm×p阶设计矩阵;QB+表示协因数矩阵QB的广义逆矩阵(由于B矩阵常包含常数或重复元素,因此QB多为奇异矩阵)。A中含有常数的EIV模型也被称为部分EIV模型[20]。在大地测量中,含有误差的输入输出向量可等价为观测值,其先验误差可以近似或精确给定,根据误差传播原理,当设计矩阵中的随机元素是观测向量的线性组合时,可得到如下的等价随机模型:

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{B = }}\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{A}}&\mathit{\boldsymbol{y}} \end{array}} \right];{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{B}}} \sim N\left[ {{\bf{0}},\sigma _0^2{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}} \right]}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{B}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{A}}}}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\rm{vec}}\left( {{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right)}\\ {{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \end{array}} \right],{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{A}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Ay}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yA}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}} \end{array}} \right]} \end{array} $$ (5)

      QB表示增广设计矩阵的协因数阵,可以通过以下方法计算[8]

      $$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{A}}}}&{{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{Ay}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yA}}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{J}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{J}}^{\rm{T}}}}&{\mathit{\boldsymbol{J}}{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{xy}}}}}\\ {{\mathit{\boldsymbol{Q}}_{\mathit{\boldsymbol{yx}}}}{\mathit{\boldsymbol{J}}^{\rm{T}}}}&{{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}}} \end{array}} \right] $$ (6)

      式中,$\mathit{\boldsymbol{J = }}{\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{J}}_{11}^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{J}}_{m1}^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{J}}_{1p}^{\rm{T}}}& \cdots &{\mathit{\boldsymbol{J}}_{mp}^{\rm{T}}} \end{array}} \right]^{\rm{T}}} \in {{\bf{R}}^{mp \times q}} $,q表示观测值向量的维数;Jij表示矩阵元素Aij对应于向量x的Jacobi矩阵,即:

      $$ {\mathit{\boldsymbol{J}}_{ij}} = \frac{{\partial {\mathit{\boldsymbol{A}}_{ij}}\left( \mathit{\boldsymbol{x}} \right)}}{{\partial \mathit{\boldsymbol{x}}}}\left( {1 \le i \le m,1 \le j \le p} \right) $$ (7)

      显然,当A的随机元素可表示为x的线性形式时,J为常数矩阵。若向量xy相关,即Qxy0Ay也相关。

      在测量平差问题中,很多情况下,参数向量ξ需满足一定的限制条件,限制条件用来提供观测系统的先验误差信息,或表示参数间的物理或几何约束,或用来减弱观测系统的病态,通常可以表示为等式条件或不等式约束形式[13]。本文仅考虑等式约束,设约束条件的一般形式可表示为:

      $$ \mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{\xi }} \right) = {\bf{0}} $$ (8)

      式中, c(ξ)表示s×1的函数向量,表示约束条件由s个独立函数组成。得到附限制条件的WTLS平差的通用模型:

      $$ \left. \begin{array}{l} \min :\frac{1}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}^ + \mathit{\boldsymbol{e}}\\ {\rm{s}}.\;{\rm{t}}.\;:\;\left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{\xi }} = \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}\\ \mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{\xi }} \right) = {\bf{0}} \end{array} \right\} $$ (9)
    • 根据拉格朗日乘数法原理,构造拉格朗日函数:

      $$ \begin{array}{*{20}{c}} {\min :l\left( {\mathit{\boldsymbol{\xi }},\mathit{\boldsymbol{\lambda }},\mathit{\boldsymbol{\mu }}} \right) = }\\ {\frac{1}{2}{\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}^ + \mathit{\boldsymbol{e}} + {\mathit{\boldsymbol{\lambda }}^{\rm{T}}}\left[ {\left( {\mathit{\boldsymbol{A}} - {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right)\mathit{\boldsymbol{\xi }} - \mathit{\boldsymbol{y}} + {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \right] + {\mathit{\boldsymbol{\mu }}^{\rm{T}}}\mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{\xi }} \right)} \end{array} $$ (10)

      式中,λμ分别表示m×1、s×1阶拉格朗日乘数向量。为方便表达,令$ \mathit{\boldsymbol{\eta = }}{\left( {{\bf{ \pmb{\mathit{ ξ}} }} - {\bf{1}}} \right)^{\rm{T}}} \in {{\bf{R}}^{\mathit{p + }{\rm{1}}}}$,不考虑限制条件的WTLS问题等价于以下无约束条件的极值问题[8]

      $$ \min :f\left( \mathit{\boldsymbol{\xi }} \right) = \frac{1}{2}{\mathit{\boldsymbol{\eta }}^{\rm{T}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{\left[ {\left( {{\mathit{\boldsymbol{\eta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_m}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}\left( {\mathit{\boldsymbol{\eta }} \otimes {\mathit{\boldsymbol{I}}_m}} \right)} \right]^{ - 1}}\mathit{\boldsymbol{B\eta }} $$ (11)

      其中,⊗表示Kronecker积,得到考虑限制条件的目标函数为:

      $$ \min :\mathit{\boldsymbol{l}}\left( {\mathit{\boldsymbol{\xi }},\mathit{\boldsymbol{\mu }}} \right) = f\left( \mathit{\boldsymbol{\xi }} \right) + {\mathit{\boldsymbol{\mu }}^{\rm{T}}}\left( {\mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{\xi }} \right)} \right) $$ (12)

      式(12)是一个无约束条件的非线性函数的极值问题,可以用常用的数值解算方法,如梯度法、拟牛顿法、牛顿法、Levenberg-Marquardt法求解[21]。由于牛顿法具有二阶收敛速度,本文采用牛顿迭代法,为此需要对式(12)求二阶导数。为方便表达,分别对f(ξ)和c(ξ)求导数,得到:

      $$ {\mathit{\boldsymbol{J}}_f} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial f}}{{\partial {\xi _1}}}}\\ \vdots \\ {\frac{{\partial f}}{{\partial {\xi _p}}}} \end{array}} \right],{\mathit{\boldsymbol{H}}_f} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}f}}{{\partial {\xi _1}\partial {\xi _1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {\xi _1}\partial {\xi _p}}}}\\ \vdots &{}& \vdots \\ {\frac{{{\partial ^2}f}}{{\partial {\xi _p}\partial {\xi _1}}}}& \cdots &{\frac{{{\partial ^2}f}}{{\partial {\xi _p}\partial {\xi _p}}}} \end{array}} \right] $$ (13)

      式中,JfRp×1HfRp×p分别表示函数f(ξ)的Jacobi向量和Hessian矩阵,其解析表达式在文献[8]中已给出。限制条件是一个由若干函数组成的向量,利用矩阵求导公式:

      $$ \frac{{\partial \left( {\mathit{\boldsymbol{RS}}} \right)}}{{\partial \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}}} = \frac{{\partial \mathit{\boldsymbol{R}}}}{{\partial \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}}}\left( {{\mathit{\boldsymbol{I}}_t} \otimes \mathit{\boldsymbol{S}}} \right) + \left( {{\mathit{\boldsymbol{I}}_s} \otimes \mathit{\boldsymbol{S}}} \right)\frac{{\partial \mathit{\boldsymbol{S}}}}{{\partial \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}}} $$ (14)

      式中,RS分别表示含有未知参数Γ的矩阵;Γ表示s×t的参数矩阵。由此得到限制条件的一阶和二阶导数分别为:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}} = \frac{{\partial \mathit{\boldsymbol{c}}}}{{\partial {\mathit{\boldsymbol{\xi }}^{\rm{T}}}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_{{c_1}}}}& \cdots &{{\mathit{\boldsymbol{J}}_{{c_s}}}} \end{array}} \right],{\mathit{\boldsymbol{H}}_\mathit{\boldsymbol{c}}} = \frac{{{\partial ^2}\mathit{\boldsymbol{c}}}}{{\partial \mathit{\boldsymbol{\xi }}\partial {\mathit{\boldsymbol{\xi }}^{\rm{T}}}}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_{{c_1}}}}& \cdots &{{\mathit{\boldsymbol{H}}_{{c_s}}}} \end{array}} \right]}\\ {{\mathit{\boldsymbol{J}}_{{c_i}}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial {c_i}}}{{\partial {\xi _1}}}}\\ \vdots \\ {\frac{{\partial {c_i}}}{{\partial {\xi _p}}}} \end{array}} \right],{\mathit{\boldsymbol{H}}_{{c_i}}} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}{c_i}}}{{\partial {\xi _1}\partial {\xi _1}}}}& \cdots &{\frac{{{\partial ^2}{c_i}}}{{\partial {\xi _p}\partial {\xi _1}}}}\\ \vdots &{}& \vdots \\ {\frac{{{\partial ^2}{c_i}}}{{\partial {\xi _1}\partial {\xi _p}}}}& \cdots &{\frac{{{\partial ^2}{c_i}}}{{\partial {\xi _p}\partial {\xi _p}}}} \end{array}} \right]\left( {i = 1 \cdots s} \right)} \end{array} $$ (15)

      式中,JcRp×sHcRp×ps分别表示函数c(ξ)的一阶、二阶导数。若用JH分别表示式(12)中的函数l(ξ, μ)相对于参数[ξT    μT]T向量的梯度向量和Hessian矩阵,则:

      $$ \mathit{\boldsymbol{\bar J}} = \left[ {\begin{array}{*{20}{c}} {\frac{{\partial \mathit{\boldsymbol{l}}\left( \mathit{\boldsymbol{\xi }} \right)}}{{\partial \mathit{\boldsymbol{\xi }}}}}\\ {\frac{{\partial \mathit{\boldsymbol{l}}\left( \mathit{\boldsymbol{\xi }} \right)}}{{\partial \mathit{\boldsymbol{\mu }}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_f} + {\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}}{\mathit{\boldsymbol{\mu }}^{\rm{T}}}}\\ {\mathit{\boldsymbol{c}}\left( \mathit{\boldsymbol{\xi }} \right)} \end{array}} \right] $$ (16)
      $$ \mathit{\boldsymbol{\bar H}} = \left[ {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}\mathit{\boldsymbol{l}}\left( \mathit{\boldsymbol{\xi }} \right)}}{{\partial {\mathit{\boldsymbol{\xi }}^2}}}}&{\frac{{{\partial ^2}\mathit{\boldsymbol{l}}\left( \mathit{\boldsymbol{\xi }} \right)}}{{\partial \mathit{\boldsymbol{\xi }}\partial \mathit{\boldsymbol{\mu }}}}}\\ {\frac{{{\partial ^2}\mathit{\boldsymbol{l}}\left( \mathit{\boldsymbol{\xi }} \right)}}{{\partial \mathit{\boldsymbol{\xi }}\partial \mathit{\boldsymbol{\mu }}}}}&{\frac{{{\partial ^2}\mathit{\boldsymbol{l}}\left( \mathit{\boldsymbol{\xi }} \right)}}{{\partial {\mathit{\boldsymbol{\mu }}^2}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_f} + {\mathit{\boldsymbol{H}}_\mathit{\boldsymbol{c}}}\left[ {{\mathit{\boldsymbol{I}}_p} \otimes \mathit{\boldsymbol{\mu }}} \right]}&{{\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}}}\\ {\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}^{\rm{T}}}&{\bf{0}} \end{array}} \right] $$ (17)

      根据牛顿迭代法的基本原理,设广义参数向量为$\mathit{\boldsymbol{\overline \xi }} = {\left[ {{\mathit{\boldsymbol{\xi }}^{\rm{T}}}{\mathit{\boldsymbol{\mu }}^{\rm{T}}}} \right]^{\rm{T}}} $,在已知参数当前值的情况下,得到参数的改正数为:

      $$ \delta \mathit{\boldsymbol{\bar \xi }} = - {{\mathit{\boldsymbol{\bar H}}}^{ - 1}}\mathit{\boldsymbol{\bar J}} $$ (18)

      以上给出了算法的解析表达式,需要经数值迭代实现。附有限制条件的EIV模型的WTLS平差方法归纳如下:

      1) 计算参数的初值${\mathit{\boldsymbol{\hat \xi }}^{\left( 0 \right)}}、{\mathit{\boldsymbol{\hat \mu }}^{\left( 0 \right)}} $,由式(6)计算增广设计矩阵的协方差阵QB, 并设初始化迭代次数i=0。

      2) 按以下顺序计算:$\mathit{\boldsymbol{\eta = }}{\left[ {{\mathit{\boldsymbol{\xi }}^{\left( i \right)}} - {\bf{1}}} \right]^{\rm{T}}} $,$\mathit{\boldsymbol{R = }}\left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}}{\mathit{\boldsymbol{I}}_m}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}{\left( {\mathit{\boldsymbol{\hat \eta }} \otimes {\mathit{\boldsymbol{I}}_\mathit{m}}} \right)^{ - 1}} $,$\mathit{\boldsymbol{u = RB\hat \eta }} $,${\mathit{\boldsymbol{N}}_\mathit{i}} = \left( {\mathit{\boldsymbol{e}}_i^{\rm{T}}{\mathit{\boldsymbol{I}}_\mathit{m}}} \right)\mathit{\boldsymbol{QB}}\left( {\mathit{\boldsymbol{\hat \eta }} \otimes {\mathit{\boldsymbol{I}}_\mathit{m}}} \right) + \left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}}{\mathit{\boldsymbol{I}}_\mathit{m}}} \right)\mathit{\boldsymbol{QB}}\left( {{\mathit{\boldsymbol{e}}_i} \otimes {\mathit{\boldsymbol{I}}_\mathit{m}}} \right)$(eip+1的单位向量,除第i个元素为1外,其他的元素均为零)。

      $$ \left. \begin{array}{l} {\left[ {{\mathit{\boldsymbol{J}}_f}} \right]_i} = \mathit{\boldsymbol{B}}_i^{\rm{T}}\mathit{\boldsymbol{u}} - \frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{u}}\\ {\left[ {{\mathit{\boldsymbol{H}}_f}} \right]_{ij}} = \left\{ \begin{array}{l} \mathit{\boldsymbol{B}}_j^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{B}}_i} - \mathit{\boldsymbol{B}}_j^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{B}}_i^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_j}\mathit{\boldsymbol{u}} + \frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_j}\mathit{\boldsymbol{u}} + \frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{N}}_j}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{u}},i \ne j\\ \mathit{\boldsymbol{B}}_j^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{B}}_i} - \mathit{\boldsymbol{B}}_j^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{B}}_i^{\rm{T}}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_j}\mathit{\boldsymbol{u}} + \frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_j}\mathit{\boldsymbol{u}} + \frac{1}{2}{\mathit{\boldsymbol{u}}^{\rm{T}}}{\mathit{\boldsymbol{N}}_j}\mathit{\boldsymbol{R}}{\mathit{\boldsymbol{N}}_i}\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^{\rm{T}}}\left( {\mathit{\boldsymbol{e}}_i^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}\left( {\mathit{\boldsymbol{e}}_i^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{u}},i = j \end{array} \right. \end{array} \right\} $$ (19)

      式中,i=1…mj=1…p;[Jf]i表示向量Jf的第i行元素;[Hf]ij表示矩阵Hf的第i行、第j列元素;由式(15)计算限制条件的JcHc,并进一步计算参数的改正数:

      $$ {\rm{ \mathsf{ δ} }}\hat {\bar \xi} = \left[ {\begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat \xi }}}\\ {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat \mu }}} \end{array}} \right] = - {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_f} + {\mathit{\boldsymbol{H}}_\mathit{\boldsymbol{c}}}\left[ {{\mathit{\boldsymbol{I}}_p} \otimes \mathit{\boldsymbol{\mu }}} \right]}&{{\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}}}\\ {\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}^{\rm{T}}}&{\bf{0}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_f} + {\mathit{\boldsymbol{\mu }}^{\rm{T}}}\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}^{\rm{T}}}\\ {\mathit{\boldsymbol{c}}\left( \xi \right)} \end{array}} \right] $$ (20)

      3) 判定$\frac{{\left\| {\delta \mathit{\boldsymbol{\hat {\overline \xi} }} } \right\|}}{{\left\| {\mathit{\boldsymbol{\hat {\overline \xi} }} } \right\|}} \le \varepsilon $(ε为给定的阈值)是否成立,若成立则结束计算;否则更新参数${\mathit{\boldsymbol{\hat \xi }}^{\left( {i + 1} \right)}} = {\mathit{\boldsymbol{\hat \xi }}^{\left( i \right)}} + \mathit{\delta }\mathit{\boldsymbol{\hat \xi }} $,${\mathit{\boldsymbol{\hat \mu }}^{\left( {i + 1} \right)}} = {\mathit{\boldsymbol{\hat \mu }}^{\left( i \right)}} + \mathit{\delta }\mathit{\boldsymbol{\hat \mu }} $,并令迭代次数i:=i+1,继续执行步骤2)。

      上文给出了附有限制条件的EIV模型的WTLS平差方法,由于需要求目标函数的二阶导数,导致算法的解析表达式较为复杂,但算法清晰地给出了一阶、二阶偏导的解析表达式,易于编程实现,而且牛顿算法具有较高的收敛速度。本文算法可解算与设计矩阵和观测向量相关的问题,适用于任何形式的限制条件或无限制条件的EIV模型的平差问题,是一种通用的WTLS平差方法。

    • 相关研究已证明总体最小二乘目标函数是一个非凸函数,初值的选取将影响收敛速度,甚至决定能否准确收敛。由于测量平差问题属于超定问题而且误差通常较小,因此可以利用部分观测值确定参数$\mathit{\boldsymbol{\hat \xi }} $初值,然后通过将附有限制条件的EIV模型线性化为附有限制条件的条件问题确定联系数$\mathit{\boldsymbol{\hat \mu }} $的初值,具体的方法在有关测量平差的教材中已有详细推导,本文不再赘述。

    • 精度评定是测量平差的必要步骤,诸多文献已证明WTLS平差方法是有偏的[4, 7, 20]。在小误差情况下,参数的估值偏差极小可以忽略。文献[8]将EIV模型线性化为附有参数的条件平差法问题来进行近似的精度评定。由于本文采用的牛顿迭代算法在计算过程中已经得到了目标函数的一阶、二阶导数,因此可以在精度评定中加以利用。本文提出一种新的近似精度评定方法,适用于小误差情况。首先计算单位权方差的估值:

      $$ \begin{array}{*{20}{c}} {\hat \sigma _0^2 = \frac{{\tilde e_\mathit{\boldsymbol{B}}^{\rm{T}}{\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{B}}}\tilde e_\mathit{\boldsymbol{B}}^{\rm{T}}}}{{m - p + s}} = }\\ {\frac{{{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{{\left[ {\left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_m}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}\left( {\mathit{\boldsymbol{\hat \eta }} \otimes {\mathit{\boldsymbol{I}}_m}} \right)} \right]}^{ - 1}}\mathit{\boldsymbol{B\hat \eta }}}}{{m - p + s}}} \end{array} $$ (21)

      式中,m-p+s表示多余观测数,也称为平差系统的自由度。在EIV模型中,参数的误差是由[Ay]的误差引起的,若不考虑参数的误差,式(20)中只有Jf受观测值误差的影响:

      $$ \Delta {\mathit{\boldsymbol{J}}_f} = \left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}}{\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{m}}}} \right)\Delta \mathit{\boldsymbol{B}} $$ (22)

      式中,ΔJf表示ΔBJf产生的影响量。将式(22)代入式(20)中,并运用误差传播定律得到:

      $$ {\mathit{\boldsymbol{Q}}_{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat {\bar \xi }}}}} = {{\mathit{\boldsymbol{\bar H}}}^{ - 1}}\mathit{\boldsymbol{\bar U}}{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{B}}}{{\mathit{\boldsymbol{\bar U}}}^{\rm{T}}}{{\mathit{\boldsymbol{\bar H}}}^{ - 1}} $$ (23)

      其中,

      $$ \mathit{\boldsymbol{\bar H}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_f} + {\mathit{\boldsymbol{H}}_\mathit{\boldsymbol{c}}}\left[ {{\mathit{\boldsymbol{I}}_p} \otimes \mathit{\boldsymbol{\mu }}} \right]}&{{\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}}}\\ {\mathit{\boldsymbol{J}}_\mathit{\boldsymbol{c}}^{\rm{T}}}&{\bf{0}} \end{array}} \right],\mathit{\boldsymbol{\bar U}} = \left[ {\begin{array}{*{20}{c}} {\left( {{{\mathit{\boldsymbol{\hat \eta }}}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_\mathit{\boldsymbol{m}}}} \right)}\\ {\bf{0}} \end{array}} \right] $$ (24)
    • 本文提出的附有限制条件的平差方法适用于设计矩阵A是观测值的线性组合、Ay相关的EIV模型,而且限制条件可以是任意非线性函数的情况。函数模型、随机假设、限制条件均可以拓展或简化以满足特定应用。如在房屋重建问题中,限制条件通常只包含一个二次型限制条件,文献[14]中附有线性限制条件和二次型条件也可以简化为二次型限制条件,设二次型限制条件的通式为:

      $$ {\mathit{\boldsymbol{\xi }}^{\rm{T}}}\mathit{\boldsymbol{M\xi }} = 1 $$ (25)

      式中,M是一个p×p的方阵,按照本文方法可以得到以下的迭代公式:

      $$ \left[ {\begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat \xi }}}\\ {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{\hat \mu }}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{H}}_f} + \mathit{\boldsymbol{\mu M}}}&{2\mathit{\boldsymbol{M\xi }}}\\ {2{\mathit{\boldsymbol{\xi }}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}}&{\bf{0}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_f} + 2\mathit{\boldsymbol{\mu M\xi }}}\\ {{\mathit{\boldsymbol{\xi }}^{\rm{T}}}\mathit{\boldsymbol{M\xi }} - 1} \end{array}} \right] $$ (26)

      类似地可以得到相应的数值计算方法,通过模型和算法的简化,可以拓展EIV模型和WTLS方法的应用。

    • 本文算例模拟某矩形房屋的重建问题,设房屋长10 m、宽5 m, AB边与Y轴正向的夹角为π/6,A点的坐标为(10, 10),沿长度方向测量10个点,沿宽度方向测量5个点,其坐标的理论值分别沿长度和宽度方向等距离选取。为模拟观测值的不等精度问题,ABBCCDDA边的xy坐标的误差分别为±0.2、±0.4、±0.6、±0.8 m,同一点的xy坐标的相关系数在[-0.5,0.5]之间随机产生,用于模拟设计矩阵和观测向量相关的情况。将模拟得到的测量误差和坐标理论值相加得到观测值。设直线ABBC的斜率分别为kukv,4条边的截距分别为b1b2b3b4。为保证房屋的两条邻边互相垂直,斜率需满足kukv=-1的限制条件。由此得到附有二次型限制条件的EIV模型,其中各矩阵可表示为:

      $$ \mathit{\boldsymbol{y}} = \left( {\begin{array}{*{20}{c}} {y_{AB}^1}\\ \vdots \\ {y_{AB}^{10}}\\ \vdots \\ {y_{DA}^1}\\ \vdots \\ {y_{DA}^5} \end{array}} \right),\mathit{\boldsymbol{A}} = \left( {\begin{array}{*{20}{c}} 1&0&0&0&{x_{AB}^1}&0\\ \vdots &{}&{}&{}&{}& \vdots \\ 1&0&0&0&{x_{AB}^{10}}&0\\ \vdots &{}&{}&{}&{}& \vdots \\ 0&0&0&1&0&{x_{DA}^1}\\ \vdots &{}&{}&{}&{}& \vdots \\ 0&0&0&1&0&{x_{DA}^5} \end{array}} \right),\mathit{\boldsymbol{\xi }} = \left( {\begin{array}{*{20}{c}} {{b_1}}\\ {{b_2}}\\ {{b_3}}\\ {{b_4}}\\ {{k_u}}\\ {{k_v}} \end{array}} \right),\mathit{\boldsymbol{M}} = \left( {\begin{array}{*{20}{c}} 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&0\\ 0&0&0&0&0&{ - 1/2}\\ 0&0&0&0&{ - 1/2}&0 \end{array}} \right) $$ (27)

      式中,(xABi, yABi)表示线段AB上第i个点的坐标值,其他边上的测量点的表达方法类似。至此,二维建筑物重建问题已经被表示为附有限制条件的EIV模型,和文献[18]的传统平差方法相比,本文方法构建简单,计算时无需对非线性的条件方程线性化。表 1给出了一组模拟数据及其误差,将其作为观测值组成附有限制条件的EIV模型。很显然,上述问题属于附有线性和二次型限制条件的WTLS方法的特殊形式,将其用文献[14]方法求解,不能准确收敛。由此可见,文献[14]方法在某些情况下不能适用,而用本文方法在经过5次迭代后能准确收敛,从而说明本文方法是一种有效且通用的方法。

      表 1  坐标观测值及其误差

      Table 1.  Measured Coordinates and TheirCovariance Matrices

      线段 编号 x/m y/m σx/m σy/m ρ
      AB 1 11.829 2 11.085 1 0.2 0.2 0.158 6
      2 13.257 0 11.898 5 0.2 0.2 0.363 6
      3 14.89 61 12.669 0 0.2 0.2 0.067 6
      4 16.643 2 13.889 6 0.2 0.2 0.480 5
      5 18.133 2 14.782 9 0.2 0.2 0.291 8
      6 19.569 3 15.874 3 0.2 0.2 -0.347 4
      7 20.889 9 16.326 2 0.2 0.2 0.333 0
      8 22.476 9 17.493 4 0.2 0.2 -0.308 1
      9 24.000 3 17.961 9 0.2 0.2 0.139 0
      10 25.810 7 18.986 1 0.2 0.2 0.169 0
      BC 1 26.873 0 21.542 0 0.4 0.4 0.272 1
      2 25.429 1 23.148 4 0.4 0.4 -0.120 2
      3 25.855 6 24.298 1 0.4 0.4 -0.058 4
      4 23.367 4 25.526 9 0.4 0.4 -0.016 9
      5 23.210 0 27.384 0 0.4 0.4 0.108 1
      CD 1 21.082 6 27.340 3 0.6 0.6 -0.324 0
      2 19.306 1 27.328 6 0.6 0.6 -0.498 0
      3 18.642 1 26.241 6 0.6 0.6 0.290 2
      4 16.225 0 25.370 7 0.6 0.6 0.013 6
      5 14.494 1 25.455 8 0.6 0.6 -0.286 8
      6 12.768 9 23.845 0 0.6 0.6 -0.396 6
      7 11.334 2 22.192 1 0.6 0.6 -0.342 7
      8 10.277 9 21.232 1 0.6 0.6 -0.092 5
      9 8.124 3 21.55 10 0.6 0.6 -0.092 2
      10 6.785 1 19.495 8 0.6 0.6 -0.447 3
      DA 1 5.636 7 16.556 6 0.8 0.8 0.441 8
      2 6.677 8 16.086 5 0.8 0.8 -0.350 0
      3 6.860 9 14.313 9 0.8 0.8 -0.115 6
      4 7.999 9 14.213 5 0.8 0.8 -0.188 9
      5 10.213 6 9.386 0 0.8 0.8 -0.331 5

      考虑二次型限制条件的WTLS平差方法与不考虑限制条件,采用WTLS方法单独拟合4条边界的结果见表 2,其中$\mathit{\hat a}、\mathit{\hat b} $分别表示各边界的斜率和截距,Δa、Δb表示估计参数与相应的真值间的差值。然后通过求出直线交点得到房屋角点坐标,上述两种方法的重建效果见图 1。从表 2的结果可以看出,附二次型限制条件的估计结果与真值的最大差值在0.4 m之内,而不考虑正交限制条件的单独拟合方法与真值的最大差值达7.2 m,说明前者的重建精度明显高于后者。增加模拟次数,比较参数的均方根误差也可以得出相同的结论,其主要原因是:附限制条件的WTLS方法顾及了矩形参数间的几何约束条件,并将所有的点整体参与平差,从而增加了整个平差系统的多余观测数,因而提高了重建精度。从图 1中可以看出,附限制条件的重建方法能保持矩形的形状,而无限制条件的WTLS方法无法保持,说明了附有限制条件的WTLS方法的应用价值。

      表 2  附限制条件的WTLS方法与边界单独平差结果的比较

      Table 2.  Comparison of Constrained WTLS Method and Adjustment Individually

      线段 附二次型限制条件的WTLS法 无限制条件的WTLS法 真值
      $\mathit{\hat a}\ $ $\mathit{\hat b}\ $ Δa Δb $\mathit{\hat a}\ $ $\mathit{\hat b}\ $ Δa Δb a b
      AB 0.575 6 4.288 4 -0.001 8 0.061 9 0.579 6 4.245 3 0.002 2 0.018 8 0.577 4 4.226 5
      BC -1.737 4 67.705 1 -0.005 3 0.384 6 -1.431 4 60.118 4 0.300 6 -7.202 1 -1.732 1 67.320 5
      CD 0.575 6 15.976 9 -0.001 8 0.203 4 0.549 5 16.342 5 -0.027 9 0.569 0 0.577 4 15.773 5
      DA -1.737 4 27.201 0 -0.005 3 -0.119 6 -1.632 1 26.325 5 0.100 0 -0.995 0 -1.732 1 27.320 5
    • 理论研究方面,本文针对基于遥感数据的二维建筑物重建问题,考虑观测向量和设计矩阵相关的情况,推导了附限制条件的EIV模型的WTLS平差方法,并将二次型限制条件作为限制条件的特例给出了相应的算法和近似精度评定方法。本文采用的牛顿迭代算法具有二阶收敛速度,研究成果对于完善WTLS平差方法、拓展总体最小二乘方法的应用有一定的理论价值。本文仅考虑设计矩阵是观测值的线性组合的情况,而且采用的牛顿算法需要计算目标函数的二阶导数,从而增加了解析表达式和相应算法的复杂程度,所采用的精度评定方法是近似方法,只适用于小误差的情况,以上问题尚需进一步研究和完善。

      在应用方面,本文以基于LiDAR或高分遥感数据的二维直角建筑物高精度重建为应用背景,研究了EIV模型构建及相应的WTLS平差算法。和现有平差方法相比,本文方法具有建模简洁、解算精度高的优势,可进一步拓展至复杂建筑的三维重建问题。本文研究成果在测绘技术从以确定点的坐标和高程为主的单点测量模式逐步转换为以测量二维或三维几何模型为主的信息测绘模式、遥感技术大规模应用的技术背景下,具有重要的应用价值。

      图  1  有限制条件和无限制条件的WTLS平差方法的矩形房屋重建效果

      Figure 1.  Performance of Constrained and Unconstrained WTLS Algorithms for Building Reconstruction

参考文献 (21)

目录

    /

    返回文章
    返回