留言板

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

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

基于变量投影的结构总体最小二乘算法

吕志鹏 隋立芬

吕志鹏, 隋立芬. 基于变量投影的结构总体最小二乘算法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
引用本文: 吕志鹏, 隋立芬. 基于变量投影的结构总体最小二乘算法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
LÜ Zhipeng, SUI Lifen. Structured Total Least Squares Method Based on Variable Projection[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
Citation: LÜ Zhipeng, SUI Lifen. Structured Total Least Squares Method Based on Variable Projection[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115

基于变量投影的结构总体最小二乘算法

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

国家自然科学基金 41674016

国家自然科学基金 41274016

国家自然科学基金 41904039

详细信息
    作者简介:

    吕志鹏,博士,主要从事空间大地测量数据处理方面的研究。lvzhipeng1989@qq.com

  • 中图分类号: P207

Structured Total Least Squares Method Based on Variable Projection

Funds: 

The National Natural Science Foundation of China 41674016

The National Natural Science Foundation of China 41274016

The National Natural Science Foundation of China 41904039

More Information
    Author Bio:

    LÜ Zhipeng, PhD, specializes in spatial geodetic data processing. E-mail: lvzhipeng1989@qq.com

  • 摘要: 测绘领域诸多实际应用中系数矩阵和观测向量具有结构特征,即系数矩阵和观测向量中包含固定量(甚至固定列) 和随机量,并且不同位置的随机量线性相关。针对这个问题,从变量误差(errors‐in‐variables,EIV)函数模型出发,首先,将系数矩阵和观测向量构成的增广矩阵表示为仿射函数形式,并采用变量投影法对函数模型进行重构;然后,利用拉格朗日法推导出了一种结构总体最小二乘(structured total least squares,STLS) 估计算法。算例分析结果表明,该算法与已有能够解决系数矩阵和观测向量存在结构特征的加权或结构总体最小二乘算法估计结果一致,说明了该算法的有效性,同时阐明了该算法与已有相关算法的关系。
  • 表  1  观测点和对应的权

    Table  1.   Observed Points and Corresponding Weights

    点号 坐标 坐标的权
    u v Wu Wv
    1 0.0 5.9 1 000.0 1.0
    2 0.9 5.4 1 000.0 1.8
    3 1.8 4.4 500.0 4.0
    4 2.6 4.6 800.0 8.0
    5 3.3 3.5 200.0 20.0
    6 4.4 3.7 80.0 20.0
    7 5.2 2.8 60.0 70.0
    8 6.1 2.8 20.0 70.0
    9 6.5 2.4 1.8 100.0
    10 7.4 1.5 1.0 500.0
    下载: 导出CSV

    表  2  不同算法所得直线拟合结果

    Table  2.   Line Fitting Results from Different Methods

    参数 文献[17]算法 文献[18]算法 文献[19]算法 文献[15]算法 本文算法
    $ \widehat{a} $ 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0
    $ \widehat{b} $ -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4
    $ \widehat{\sigma }{}_{0}^{2} $ 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3
    $ \widehat{\sigma }{}_{\widehat{a}}^{2} $ 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0
    $ \widehat{\sigma }{}_{\widehat{b}}^{2} $ 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5
    迭代次数 13 8 7 9 7
    下载: 导出CSV

    表  3  高程数据

    Table  3.   Elevation Data

    点号 高程/mm 点号 高程/mm 点号 高程/mm 点号 高程/mm
    1 26.33 10 25.46 19 28.09 28 26.57
    2 26.27 11 26.12 20 26.78 29 28.36
    3 26.43 12 27.28 21 28.66 30 27.94
    4 25.56 13 26.67 22 26.75 31 26.81
    5 26.82 14 27.95 23 27.24 32 28.50
    6 26.56 15 26.74 24 28.02 33 27.68
    7 25.93 16 27.53 25 26.81 34 26.57
    8 26.43 17 25.31 26 28.50 35 28.36
    9 26.52 18 26.90 27 27.68 36 27.94
    下载: 导出CSV

    表  4  线性自回归模型不同算法所得结果

    Table  4.   Results from Different Methods of Linear Auto-Regression Model

    参数 文献[19]算法 文献[10]算法 文献[21]算法 文献[15]算法 本文算法
    $ {\widehat{\xi }}_{1} $ 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2
    $ {\widehat{\xi }}_{2} $ 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4
    $ {\widehat{\xi }}_{3} $ -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2
    $ \widehat{\sigma }{}_{0}^{2} $ 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1
    $ \widehat{\sigma }{}_{{\widehat{\xi }}_{1}}^{2} $ 0.012 658 067 0 0.012 658 067 0 0.012 658 067 0 0.012 658 067 0
    $ \widehat{\sigma }{}_{{\widehat{\xi }}_{2}}^{2} $ 0.009 481 774 9 0.009 481 774 9 0.009 481 774 9 0.009 481 774 9
    $ \widehat{\sigma }{}_{{\widehat{\xi }}_{3}}^{2} $ 0.009 074 553 9 0.009 074 553 9 0.009 074 553 9 0.009 074 553 9
    迭代次数 43 43 43 60 43
    下载: 导出CSV

    表  5  参数估计表达式比较

    Table  5.   Comparison of Parameter Estimation Expression

    算法 参数估计表达式 结构矩阵提取方法
    文献[10]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \widehat{\mathit{\boldsymbol{B}}}={\widehat{\mathit{\boldsymbol{x}}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}\otimes {\mathit{\boldsymbol{I}}}_{m} $,$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{l}}l}=\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\mathit{\boldsymbol{S}}}^{\mathrm{T}} $;$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵;$ \widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{S}} $为结构矩阵 $ \mathit{\boldsymbol{S}}\cdot \mathit{\boldsymbol{e}}=\left[\begin{array}{cc}{\mathit{\boldsymbol{C}}}_{\mathit{\boldsymbol{A}}}& 0\\ 0& {I}_{m}\end{array}\right]\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $
    文献[21]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,结构矩阵$ {\widehat{\mathit{\boldsymbol{C}}}}_{2}=\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{C}}}_{1} $ $ {\mathit{\boldsymbol{C}}}_{1}\cdot \mathit{\boldsymbol{e}}={\mathit{\boldsymbol{C}}}_{1}\cdot \left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $
    本文算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \mathit{\boldsymbol{\Gamma }}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $,其中,结构矩阵$ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{\psi }} $ $ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{cccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)& \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}\right)\end{array}\right] $通过变量投影法建立$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}} $
    文献[19]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\right) $;$ {\mathit{\boldsymbol{Q}}}_{1} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵
    下载: 导出CSV
  • [1] Fuller W A. Measurement Error Models[M]. New York: Wiley, 1987
    [2] Golub G H, Van Loan C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis, 1980, 17(6): 883-893 doi:  10.1137/0717073
    [3] Van Huffel S, van Dewalle J. The Total Least Squares Problem: Computational Aspects and Analysis[M]. Philadelphia: SIAM, 1991
    [4] 刘经南, 曾文宪, 徐培亮. 整体最小二乘估计的研究进展[J]. 武汉大学学报·信息科学版, 2013, 38(5): 505-512 http://ch.whu.edu.cn/article/id/2641

    Liu Jingnan, Zeng Wenxian, Xu Peiliang. Overview of Total Least Squares Methods[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 505-512 http://ch.whu.edu.cn/article/id/2641
    [5] 王乐洋, 许才军. 总体最小二乘研究进展[J]. 武汉大学学报·信息科学版, 2013, 38(7): 850-856 http://ch.whu.edu.cn/article/id/2703

    Wang Leyang, Xu Caijun. Progress in Total Least Squares[J]. Geomatics and Information Science of Wuhan University, 2013, 38(7): 850-856 http://ch.whu.edu.cn/article/id/2703
    [6] Dunne B E, Williamson G A. QR-Based TLS and Mixed LS-TLS Algorithms with Applications to Adaptive ⅡR Filtering[J]. IEEE Transactions on Signal Processing, 2003, 51(2): 386-394 doi:  10.1109/TSP.2002.806980
    [7] Yan Shijian, Fan Jinyan. The Solution Set of the Mixed LS-TLS Problem[J]. International Journal of Computer Mathematics, 2011, 77(4): 545-561
    [8] 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
    [9] Schaffrin B. A Note on Constrained Total Least-Squares Estimation[J]. Linear Algebra and Its Applications, 2006, 417(1): 245-258 doi:  10.1016/j.laa.2006.03.044
    [10] Fang Xing. A Structured and Constrained Total Least-Squares Solution with Cross-Covariances[J]. Studia Geophysica et Geodaetica, 2014, 58(1), DOI:  10.1007/s11200-012-0671-z
    [11] Xu Peiliang, Liu Jiangnan, Shi Chuang. Total Least Squares Adjustment in Partial Errors-in-Variables Models: Algorithm and Statistical Analysis[J]. Journal of Geodesy, 2012, 86(8): 661-675 doi:  10.1007/s00190-012-0552-9
    [12] Shi Yun, Xu Peiliang, Liu Jingnan, et al. Alternative Formulae for Parameter Estimation in Partial Errors-in-Variables Models[J]. Journal of Geodesy, 2015, 89(1): 13-16 doi:  10.1007/s00190-014-0756-2
    [13] Rosen J B, Park H, Glick J. Total Least Norm Formulation and Solution for Structured Problems[J]. SIAM Journal on Matrix Analysis and Applications, 1996, 17(1): 110-126 doi:  10.1137/S0895479893258802
    [14] Yalamov P Y, Yuan J Y. A Successive Least Squares Method for Structured Total Least Squares[J]. Journal of Computational Mathematics, 2003, 21(4): 463-472
    [15] Markovsky I, Van Huffel S. High-performance Numerical Algorithms and Software for Structured Total Least Squares[J]. Journal of Computational and Applied Mathematics, 2005, 180(2): 311-331 doi:  10.1016/j.cam.2004.11.003
    [16] Neri F, Saitta G, Chiofalo S. An Accurate and Straight Forward Approach to Line Regression Analysis of Error-Affected Experimental Data[J]. Journal of Physics E: Scientific Instruments, 1989, 22(4): 215-217 doi:  10.1088/0022-3735/22/4/002
    [17] Shen Yunzhong, Li Bofeng, Chen Yi. 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
    [18] Snow K. Topics in Total Least-Squares Adjustment within the Errors-in-Variables Model: Singular Cofactor Matrices and Prior Information[D]. Columbus: The Ohio State University, 2012
    [19] Amiri-Simkooei A, Jazaeri S. Weighted Total Least Squares Formulated by Standard Least Squares Theory[J]. Journal of Geodetic Science, 2012, 2(2): 113-124 doi:  10.2478/v10156-011-0036-5
    [20] 王新洲, 陶本藻, 邱卫宁, 等. 高等测量平差[M]. 北京: 测绘出版社, 2006

    Wang Xinzhou, Tao Benzao, Qiu Weining, et al. Higher Surveying Adjustment[M]. Beijing: Surveying and Mapping Press, 2006
    [21] 王乐洋, 许光煜, 温贵森. 一种相关观测的Partial EIV模型求解方法[J]. 测绘学报, 2017, 46(8): 978-987 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201708008.htm

    Wang Leyang, Xu Guangyu, Wen Guisen. A Method for Partial EIV Model with Correlated Observations[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(8): 978-987 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201708008.htm
    [22] Fang Xing. Weighted Total Least Squares Solution for Application in Geodesy[D]. Hanover: Leibniz University Hanover, 2011
    [23] Fang Xing. 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
  • [1] 吕志鹏, 隋立芬.  基于变量投影法的自回归模型方差分量估计 . 武汉大学学报 ● 信息科学版, 2020, 45(2): 205-212. doi: 10.13203/j.whugis20180352
    [2] 吕志鹏, 隋立芬.  基于非线性高斯-赫尔默特模型的结构总体最小二乘法 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1808-1815. doi: 10.13203/j.whugis20180104
    [3] 李思达, 柳林涛, 刘志平, 艾青松.  多变量稳健总体最小二乘平差方法 . 武汉大学学报 ● 信息科学版, 2019, 44(8): 1241-1248. doi: 10.13203/j.whugis20170232
    [4] 吴炜, 沈占锋, 王显威, 吴田军, 王卫红.  一种仿射不变量支持的准稠密控制点匹配算法 . 武汉大学学报 ● 信息科学版, 2018, 43(6): 930-936. doi: 10.13203/j.whugis20160146
    [5] 李霖, 曾祥倚, 段新桥, 朱海红, 张一鸣.  一种基于体绘制的二元变量体可视化方法 . 武汉大学学报 ● 信息科学版, 2017, 42(8): 1082-1089. doi: 10.13203/j.whugis20150288
    [6] 赵俊, 归庆明, 郭飞宵.  基于改进目标函数的partial EIV模型WTLS估计的新算法 . 武汉大学学报 ● 信息科学版, 2017, 42(8): 1179-1184. doi: 10.13203/j.whugis20150180
    [7] 王乐洋, 余航.  总体最小二乘联合平差 . 武汉大学学报 ● 信息科学版, 2016, 41(12): 1683-1689. doi: 10.13203/j.whugis20140670
    [8] 马友青, 刘少创, 魏士俨, 李明磊.  加权总体最小二乘的地面解析摄影测量算法 . 武汉大学学报 ● 信息科学版, 2015, 40(5): 594-598. doi: 10.13203/j.whugis20130387
    [9] 姚宜斌, 黄书华, 孔建, 何军泉.  空间直线拟合的整体最小二乘算法 . 武汉大学学报 ● 信息科学版, 2014, 39(5): 571-574. doi: 10.13203/j.whugis20120104
    [10] 姚宜斌, 孔建.  顾及设计矩阵随机误差的最小二乘组合新解法 . 武汉大学学报 ● 信息科学版, 2014, 39(9): 1028-1032. doi: 10.13203/j.whugis20130030
    [11] 王乐洋, 许才军.  总体最小二乘研究进展 . 武汉大学学报 ● 信息科学版, 2013, 38(7): 850-856.
    [12] 罗志才, 林旭, 周波阳.  自协方差最小二乘噪声估计的改进算法 . 武汉大学学报 ● 信息科学版, 2012, 37(10): 1164-1167.
    [13] 邱卫宁, 齐公玉, 田丰瑞.  整体最小二乘求解线性模型的改进算法 . 武汉大学学报 ● 信息科学版, 2010, 35(6): 708-710.
    [14] 鲁铁定, 周世健.  总体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1351-1354.
    [15] 刘国林, 郝华东, 陶秋香.  卡尔曼滤波相位解缠及其与其他方法的对比分析 . 武汉大学学报 ● 信息科学版, 2010, 35(10): 1174-1178.
    [16] 薛树强, 党亚民, 陈武.  最小二乘估值均方差计算的矩阵体积法 . 武汉大学学报 ● 信息科学版, 2009, 34(9): 1106-1109.
    [17] 邓兴升, 花向红.  动态最小二乘支持向量机学习算法 . 武汉大学学报 ● 信息科学版, 2008, 33(11): 1122-1125.
    [18] 徐天河, 杨元喜.  均方误差意义下正则化解优于最小二乘解的条件 . 武汉大学学报 ● 信息科学版, 2004, 29(3): 223-226.
    [19] 於宗俦.  最小二乘解的一种直接算法 . 武汉大学学报 ● 信息科学版, 1999, 24(1): 68-70.
    [20] 邓德祥, 吴章华, 胡志雄.  一种用于影像镶嵌的快速最小二乘序贯算法 . 武汉大学学报 ● 信息科学版, 1997, 22(1): 29-31.
  • 加载中
计量
  • 文章访问数:  35
  • HTML全文浏览量:  2
  • PDF下载量:  16
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-14
  • 刊出日期:  2021-03-05

基于变量投影的结构总体最小二乘算法

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

    国家自然科学基金 41674016

    国家自然科学基金 41274016

    国家自然科学基金 41904039

    作者简介:

    吕志鹏,博士,主要从事空间大地测量数据处理方面的研究。lvzhipeng1989@qq.com

  • 中图分类号: P207

摘要: 测绘领域诸多实际应用中系数矩阵和观测向量具有结构特征,即系数矩阵和观测向量中包含固定量(甚至固定列) 和随机量,并且不同位置的随机量线性相关。针对这个问题,从变量误差(errors‐in‐variables,EIV)函数模型出发,首先,将系数矩阵和观测向量构成的增广矩阵表示为仿射函数形式,并采用变量投影法对函数模型进行重构;然后,利用拉格朗日法推导出了一种结构总体最小二乘(structured total least squares,STLS) 估计算法。算例分析结果表明,该算法与已有能够解决系数矩阵和观测向量存在结构特征的加权或结构总体最小二乘算法估计结果一致,说明了该算法的有效性,同时阐明了该算法与已有相关算法的关系。

English Abstract

吕志鹏, 隋立芬. 基于变量投影的结构总体最小二乘算法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
引用本文: 吕志鹏, 隋立芬. 基于变量投影的结构总体最小二乘算法[J]. 武汉大学学报 ● 信息科学版, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
LÜ Zhipeng, SUI Lifen. Structured Total Least Squares Method Based on Variable Projection[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
Citation: LÜ Zhipeng, SUI Lifen. Structured Total Least Squares Method Based on Variable Projection[J]. Geomatics and Information Science of Wuhan University, 2021, 46(3): 388-394. doi: 10.13203/j.whugis20190115
  • 变量误差(errors‐in‐variables,EIV)模型[1]是一种同时顾及系数矩阵和观测向量误差的回归模型。对此模型可以采用总体最小二乘(total least squares,TLS)估计进行求解,TLS估计本质上是对最小二乘(least squares,LS)估计的扩展,是LS估计的去正则化处理[2]。在系数矩阵存在误差的情况下,TLS估计优于LS估计,并且具有弱一致性和渐进无偏性[2-3]

    目前,测绘领域对TLS及其各种扩展算法的应用集中于直线拟合、基准变换、自回归分析等问题[4-5],在这些应用中,EIV模型存在结构特征,即系数矩阵和观测向量存在固定量和随机量及不同位置的随机量之间存在线性相关性。可以采用基于广义奇异值分解或者矩阵QR分解的混合最小二乘-总体最小二乘(least squares-total least squares,LS-TLS)算法[3, 6-7]进行求解,这些算法只针对系数矩阵存在固定列的情况;也可以采用观测权阵生成法则[8],将随机量纳入平差的随机模型,通过构造奇异的观测权阵并利用加权总体最小二乘(weighted total least squares,WTLS)估计进行求解,这种算法增加了待估量的个数;甚至可以采用各种约束总体最小二乘(constrained total least squares,CTLS)算法进行求解[9-10],不过这将增加求解的复杂程度。

    EIV模型存在结构特征本质上是结构总体最小二乘(structured total least squares,STLS)问题。目前,已存在诸多算法可以解决这一问题,如文献[11-12]引入伪观测值将随机量纳入平差的函数模型;文献[13]将STLS问题转化为非线性非约束优化问题,然后采用牛顿法进行求解;文献[14]将参数和随机量视为两组待估量,进而将STLS问题转换为两个LS问题;文献[10]分析了线性及二次约束下的STLS问题。上述各种算法仅考虑了系数矩阵的结构特征。

    文献[15]利用变量投影法生成结构矩阵,将STLS问题转化为非线性双重最小化问题,再采用牛顿法进行求解,但是这种算法推导过程相对复杂,迭代效率较低,并且不能给出参数的近似精度估计。对于这一问题,首先,本文顾及系数矩阵和观测向量中随机量的观测权阵,针对系数矩阵和观测向量存在结构特征的情况,利用变量投影法[15]将其展开成仿射函数形式;然后,将STLS问题转化为非线性约束优化问题,并通过拉格朗日法推导出了一种STLS估计算法,进一步对其特征进行了分析;最后,给出了单位权方差和参数协因数阵的一阶近似估计公式。通过两个实例,验证了本文算法的有效性,并分析了本文算法与已有处理STLS问题的相关算法之间的关系。

    • EIV模型如下:

      $$ \mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}=(\mathit{\boldsymbol{A}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}})\mathit{\boldsymbol{x}}, m>n=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{\boldsymbol{A}}\right) $$ (1)

      式中,$ \mathit{\boldsymbol{y}}\in {\mathit{\boldsymbol{R}}}^{m\times 1} $为观测向量;ey为其误差向量;$ \mathit{\boldsymbol{A}}\in {\mathit{\boldsymbol{R}}}^{m\times n} $为系数矩阵;EA为其误差矩阵;$ \mathit{\boldsymbol{x}}\in {\mathit{\boldsymbol{R}}}^{m\times 1} $为参数向量。

      STLS问题的函数模型如式(1)所示,该模型要求系数矩阵的增广矩阵$ \mathit{\boldsymbol{C}}=\left[\mathit{\boldsymbol{A}}\;\;\mathit{\boldsymbol{y}}\right] $与其误差矩阵$ \mathrm{\Delta }\mathit{\boldsymbol{C}}=\left[{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\right] $具有相同的结构。在此假设$ \mathit{\boldsymbol{C}} $仅存在某种仿射结构,由$ \mathit{\boldsymbol{C}} $中独立随机量组成的向量为$ \mathit{\boldsymbol{p}}\in {\mathit{\boldsymbol{R}}}^{{n}_{\mathit{p}}} $,$ {n}_{\mathit{p}} $为独立随机量个数。故根据变量投影法可以将$ \mathit{\boldsymbol{C}} $表示为如下仿射函数[15]

      $$ \mathit{\boldsymbol{C}}=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right)={\mathit{\boldsymbol{S}}}_{0}+\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{S}}}_{i}{\mathit{\boldsymbol{p}}}_{i}, {\mathit{\boldsymbol{S}}}_{0}, {\mathit{\boldsymbol{S}}}_{i}\in {\mathit{\boldsymbol{R}}}^{m\times (n+1)} $$ (2)

      式中,$ {\mathit{\boldsymbol{S}}}_{0} $为固定量仿射结构矩阵,$ {\mathit{\boldsymbol{S}}}_{0} $中非零元素为$ \mathit{\boldsymbol{C}} $中固定量;$ {\mathit{\boldsymbol{S}}}_{i} $为随机量$ {\mathit{\boldsymbol{p}}}_{i} $对应的仿射结构矩阵,一般为稀疏矩阵。STLS问题的数学模型为:

      $$ \left\{\begin{array}{l}\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{p}}-\mathrm{\Delta }\mathit{\boldsymbol{p}}){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=0\\ \mathrm{\Delta }\mathit{\boldsymbol{p}}~N(0, {\sigma }_{0}^{2}\mathit{\boldsymbol{Q}})\end{array}\right. $$ (3)

      式中,$ {\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=\left[\begin{array}{c}\mathit{\boldsymbol{x}}\\ -1\end{array}\right] $为$ \mathit{\boldsymbol{x}} $的增广向量;$ \mathit{\boldsymbol{Q}} $为$ \mathit{\boldsymbol{p}} $的正定协因数阵;$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $为$ \mathit{\boldsymbol{p}} $的观测误差。STLS问题的目标函数为:

      $$ \mathrm{m}\mathrm{i}\mathrm{n}{‖\mathrm{\Delta }\mathit{\boldsymbol{p}}‖}_{{\mathit{\boldsymbol{Q}}}^{-1}}^{2} $$ (4)

      对STLS问题的函数模型,进一步变形为[15]

      $$ \begin{array}{l}\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{p}}-\mathrm{\Delta }\mathit{\boldsymbol{p}}){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=0\iff \\ \;\;\;\;\;\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{S}}}_{i}\mathrm{\Delta }{\mathit{\boldsymbol{p}}}_{i}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\iff \\ \;\;\;\;\;\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{S}}}_{i}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\mathrm{\Delta }{\mathit{\boldsymbol{p}}}_{i}=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\iff \\ \;\;\;\;\;\;\;\;\mathit{\boldsymbol{G}}\left(\mathit{\boldsymbol{x}}\right)\mathrm{\Delta }\mathit{\boldsymbol{p}}=\mathit{\boldsymbol{r}}\left(\mathit{\boldsymbol{x}}\right)\end{array} $$ (5)

      式中,$ \mathit{\boldsymbol{r}}\left(\mathit{\boldsymbol{x}}\right)=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=-(\mathit{\boldsymbol{y}}-\mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}}) $;结构矩阵$ \mathit{\boldsymbol{G}}\left(\mathit{\boldsymbol{x}}\right)=\left[\begin{array}{cccc}{\mathit{\boldsymbol{S}}}_{1}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}& {\mathit{\boldsymbol{S}}}_{2}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}& \cdots & {\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\end{array}\right]=\mathit{\boldsymbol{B}}\mathit{\boldsymbol{\psi }} $,其中,$ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{cccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)& \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}\right)\end{array}\right] $,$ \mathit{\boldsymbol{B}}={\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}\otimes {\mathit{\boldsymbol{I}}}_{m} $,$ \mathrm{v}\mathrm{e}\mathrm{c}(\cdot ) $为矩阵向量化算子,$ \otimes $为克罗内克乘积算子。

      此时,STLS问题可以表示为:$ \underset{\mathit{\boldsymbol{x}}, \mathrm{\Delta }\mathit{\boldsymbol{p}}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathrm{\Delta }\mathit{\boldsymbol{p}}‖}_{{\mathit{\boldsymbol{Q}}}^{-1}}^{2} $。其中,

      $$ \mathit{\boldsymbol{G}}\left(\mathit{\boldsymbol{x}}\right)\mathrm{\Delta }\mathit{\boldsymbol{p}}=\mathit{\boldsymbol{r}}\left(\mathit{\boldsymbol{x}}\right), \mathrm{ }\mathrm{ }\mathrm{\Delta }\mathit{\boldsymbol{p}}~N(0, {\sigma }_{0}^{2}\mathit{\boldsymbol{Q}}) $$ (6)

      式(6)本质上是一个非线性约束优化问题,可以从目标函数出发采用拉格朗日乘数法进行求解。

    • 对于STLS问题,由拉格朗日乘数法构造其目标函数为:

      $$ \mathit{\Phi }(\mathit{\boldsymbol{x}}, \mathrm{\Delta }\mathit{\boldsymbol{p}}, \mathit{\boldsymbol{\lambda }})=\mathrm{\Delta }{\mathit{\boldsymbol{p}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}\mathrm{\Delta }\mathit{\boldsymbol{p}}-2{\mathit{\boldsymbol{\lambda }}}^{\mathrm{T}}\left[\mathit{\boldsymbol{G}}\right(\mathit{\boldsymbol{x}})\mathrm{\Delta }\mathit{\boldsymbol{p}}-\mathit{\boldsymbol{r}}(\mathit{\boldsymbol{x}}\left)\right] $$ (7)

      将式(7)分别对$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $、$ \mathit{\boldsymbol{\lambda }} $、$ \mathit{\boldsymbol{x}} $求偏导并令其等于零得:

      $$ \frac{\partial \mathit{\Phi } }{\partial \mathrm{\Delta }\mathit{\boldsymbol{p}}}=2{\mathit{\boldsymbol{Q}}}^{-1}\mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}-2\mathit{\boldsymbol{G}}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)}^{\mathrm{T}}\widehat{\mathit{\boldsymbol{\lambda }}}=0 $$ (8)
      $$ \frac{\partial \mathit{\Phi } }{\partial \mathit{\boldsymbol{\lambda }}}=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}-\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=0 $$ (9)
      $$ \frac{\partial \mathit{\Phi } }{\partial \mathit{\boldsymbol{x}}}=(\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}{)}^{\mathrm{T}}\widehat{\mathit{\boldsymbol{\lambda }}}=0 $$ (10)

      式中,$ \left(\widehat{\cdot }\right) $表示相应参数的估值,并根据式(2)重构系数矩阵误差矩阵的估值为:

      $$ {\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}=\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}^{\mathit{\boldsymbol{A}}}}{\mathit{\boldsymbol{S}}}_{i}^{\mathit{\boldsymbol{A}}}\left(\mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}_{i}^{\mathit{\boldsymbol{A}}}\right) $$ (11)

      式中,$ \mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}_{i}^{\mathit{\boldsymbol{A}}} $和$ {\mathit{\boldsymbol{S}}}_{i}^{\mathit{\boldsymbol{A}}} $分别为系数矩阵中独立随机量估值和其对应的仿射结构矩阵;$ {n}_{\mathit{\boldsymbol{p}}}^{\mathit{\boldsymbol{A}}} $为系数矩阵中独立随机量的个数。

      假设$ \widehat{\mathit{\boldsymbol{x}}} $已知,式(9)可视为经典平差中的条件平差模型,由条件平差法可得:

      $$ \mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}=\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{\Gamma }}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)}^{-1}\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $$ (12)

      其中,$ \mathit{\boldsymbol{\Gamma }}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $。

      将式(12)代入式(8),并将方程两侧同乘$ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}} $得:

      $$ 2\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)-2\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{\lambda }}}=0 $$ (13)

      $ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $具有凯利逆,故可得:

      $$ \widehat{\mathit{\boldsymbol{\lambda }}}={\mathit{\pmb{\Gamma }}}^{-1}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $$ (14)

      将式(14)代入式(10),变形,得:

      $$ (\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}{)}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-\mathit{\boldsymbol{A}}\widehat{\mathit{\boldsymbol{x}}}+{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}})=0 $$ (15)

      进一步可得:

      $$ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $$ (16)

      式中,$ \widehat{\mathit{\boldsymbol{A}}}=\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}} $。式(16)具有经典平差中参数平差模型最小二乘解的形式。

      算法的迭代过程如下:

      1)初始值:假设$ \mathit{\boldsymbol{A}} $无误差,通过LS估计求得参数的初值$ {\widehat{\mathit{\boldsymbol{x}}}}^{0} $。

      2)由$ {\widehat{\mathit{\boldsymbol{x}}}}^{t} $通过式(12)可得$ \mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}^{t+1} $,进而通过式(11)重构$ {\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1} $,即$ {\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1}=\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}^{\mathit{\boldsymbol{A}}}}{\mathit{\boldsymbol{S}}}_{i}^{\mathit{\boldsymbol{A}}}(\mathrm{\Delta }{\widehat{p}}_{i}^{\mathit{\boldsymbol{A}}}{)}^{t+1} $,进一步得$ {\widehat{\mathit{\boldsymbol{A}}}}^{t+1}=\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1} $,然后得到$ {\mathit{\boldsymbol{\Gamma }}}^{-1}\left({\widehat{\mathit{\boldsymbol{x}}}}^{t}\right) $,最终计算出新的参数估值为:

      $$ {\widehat{\mathit{\boldsymbol{x}}}}^{t+1}=\left[\right({\widehat{\mathit{\boldsymbol{A}}}}^{t+1}{)}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\left({\widehat{\mathit{\boldsymbol{x}}}}^{t}\right){\widehat{\mathit{\boldsymbol{A}}}}^{t+1}{]}^{-1}·\\ \left({\widehat{\mathit{\boldsymbol{A}}}}^{t+1}{)}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\right({\widehat{\mathit{\boldsymbol{x}}}}^{t}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1}{\widehat{\mathit{\boldsymbol{x}}}}^{t}) $$ (17)

      3)重复步骤2),直到满足收敛条件$ ‖{\widehat{\mathit{\boldsymbol{x}}}}^{t+1}-{\widehat{\mathit{\boldsymbol{x}}}}^{t}‖<\epsilon $,其中,$ \epsilon $为收敛阈值。

      对于STLS问题,假设只有$ \mathit{\boldsymbol{y}} $存在误差,可将EIV模型转化为高斯-马尔可夫(Gauss-Markov,GM)模型,即:

      $$ \mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}}=\mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}} $$ (18)

      当$ \mathit{\boldsymbol{A}} $存在误差$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $,将其代入式(18)有:

      $$ (\mathit{\boldsymbol{A}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}})\mathit{\boldsymbol{x}}=(\mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{x}})-({\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{x}})={\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}-{\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}} $$ (19)

      顾及GM模型的假设,式(19)左侧的系数矩阵$ \mathit{\boldsymbol{A}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $无误差。同时,式(16)中ey受到$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $的增益,故$ {\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}} $的协因数阵为$ {\mathit{\boldsymbol{Q}}}_{{\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}}}=\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}=\mathit{\boldsymbol{\Gamma }}\left(\mathit{\boldsymbol{x}}\right) $。进而由参数平差法得:

      $$ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}){\widehat{\mathit{\boldsymbol{y}}}}_{\mathit{\boldsymbol{c}}} $$ (20)

      将式(12)和式(20)相结合用于求解STLS问题。由上述分析可得,TLS问题的求解可以转化为两个LS问题的求解。这一思想在文献[14]中得到了体现,该文献将STLS问题中的待估量分为$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $和$ \mathit{\boldsymbol{x}} $两类,利用式(12)求解$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $,通过直接改正系数矩阵误差求解$ \mathit{\boldsymbol{x}} $,即:

      $$ \widehat{\mathit{\boldsymbol{x}}}=({\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}{\widehat{\mathit{\boldsymbol{A}}})}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}\mathit{\boldsymbol{y}} $$ (21)

      当$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\to 0 $,式(20)趋近于式(21),因而,式(21)可以视为式(20)的简化形式。在建立求解$ \mathit{\boldsymbol{x}} $的GM模型时,式(21)并没有顾及$ \mathit{\boldsymbol{A}} $的误差$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $改正后对GM模型的影响。这一影响不断被$ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}} $吸收,影响迭代效率和迭代稳定性,同时,也无法直接获得TLS估计下$ \mathit{\boldsymbol{y}} $的改正数。

    • 根据TLS原理,上述STLS估计算法的单位权方差可以表示为:

      $$ {\sigma }^{2}=\frac{\mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}\mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}}{d} $$ (22)

      式中,$ d=m-n $为自由度,其数值与观测方程个数和参数个数有关。本文得到的参数$ \mathit{\boldsymbol{x}} $的估计表达式具有LS解的形式,因而,可以采用协因数传播律求解参数的协因数阵。由§1.2的分析可知,$ {\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}} $对应的误差向量$ {\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}} $的协因数阵为$ \mathit{\boldsymbol{\Gamma }}\left(\mathit{\boldsymbol{x}}\right) $,所以根据协因数传播律得:

      $$ {\mathit{\boldsymbol{Q}}}_{\widehat{\mathit{\boldsymbol{x}}}\widehat{\mathit{\boldsymbol{x}}}}=[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1} $$ (23)

      式(23)即为STLS问题参数协因数阵的一阶近似估计公式。

    • 选取的直线拟合算例[16]给出了10组坐标$ ({u}_{i}, {v}_{i}) $用于拟合直线的截距$ a $和斜率$ b $,其中,i=1,2…10。其对应的EIV模型为:

      $$ \left[\begin{array}{c}{v}_{1}\\ {v}_{2}\\ \vdots\\ {v}_{10}\end{array}\right]-\left[\begin{array}{c}{e}_{{v}_{1}}\\ {e}_{{v}_{2}}\\ \vdots \\ {e}_{{v}_{10}}\end{array}\right]=\left(\left[\begin{array}{cc}\begin{array}{c}1\\ 1\\ \vdots\\ 1\end{array}& \begin{array}{c}{u}_{1}\\ {u}_{2}\\ \vdots\\ {u}_{10}\end{array}\end{array}\right]-\left[\begin{array}{cc}\begin{array}{c}0\\ 0\\ \vdots\\ 0\end{array}& \begin{array}{c}{e}_{{u}_{1}}\\ {e}_{{u}_{2}}\\ \vdots\\ {e}_{{u}_{10}}\end{array}\end{array}\right]\right)\cdot \left[\begin{array}{c}a\\ b\end{array}\right] $$ (24)

      式中,$ ({e}_{{u}_{i}}, {e}_{{v}_{i}}) $为10组坐标对应的误差,其中,i=1,2…10。由式(24)可知,截距$ a $对应的系数矩阵元素的误差理论上为零,因而本质上直线拟合问题是STLS问题。10组坐标及其对应的权如表 1所示[16]

      表 1  观测点和对应的权

      Table 1.  Observed Points and Corresponding Weights

      点号 坐标 坐标的权
      u v Wu Wv
      1 0.0 5.9 1 000.0 1.0
      2 0.9 5.4 1 000.0 1.8
      3 1.8 4.4 500.0 4.0
      4 2.6 4.6 800.0 8.0
      5 3.3 3.5 200.0 20.0
      6 4.4 3.7 80.0 20.0
      7 5.2 2.8 60.0 70.0
      8 6.1 2.8 20.0 70.0
      9 6.5 2.4 1.8 100.0
      10 7.4 1.5 1.0 500.0

      表 1中的数据分别采用3种WTLS估计算法(文献[17]算法、文献[18]算法和文献[19]算法)和2种STLS估计算法(文献[15]算法和本文算法)进行求解。各种迭代算法收敛阈值设为$ 1\times {10}^{-10} $,相应的参数估值$ \widehat{a} $和$ \widehat{b} $、验后单位权方差估值$ {\widehat{\sigma }}_{0}^{2} $、参数估计精度$ {\widehat{\sigma }}_{\widehat{a}}^{2} $和$ {\widehat{\sigma }}_{\widehat{b}}^{2} $以及迭代次数如表 2所示。

      表 2  不同算法所得直线拟合结果

      Table 2.  Line Fitting Results from Different Methods

      参数 文献[17]算法 文献[18]算法 文献[19]算法 文献[15]算法 本文算法
      $ \widehat{a} $ 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0
      $ \widehat{b} $ -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4
      $ \widehat{\sigma }{}_{0}^{2} $ 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3
      $ \widehat{\sigma }{}_{\widehat{a}}^{2} $ 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0
      $ \widehat{\sigma }{}_{\widehat{b}}^{2} $ 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5
      迭代次数 13 8 7 9 7

      表 2可知,各种算法均能得到一致的估计结果,即本文算法与现有的WTLS估计算法结果一致,说明本文算法是有效的。它们的收敛次数有差别:文献[17]算法需要13次迭代收敛;文献[15]算法优于文献[17]算法,需要9次迭代收敛;文献[18]算法的迭代收敛减少到8次;文献[19]算法和本文算法迭代收敛次数最少,均为7次;本文算法相对文献[15]算法迭代效率有提高。

      通过构造秩亏的观测权阵[8],各种WTLS估计算法[17-19]可以解决STLS问题,但是WTLS估计需要求解系数矩阵的增广矩阵中的所有元素,因此,其待估量个数和观测权阵维数大于STLS估计。参数估计表达式的差异是决定迭代效率的重要因素,上述各种WTLS估计算法中,文献[19]算法迭代效率最高,并且参数估计表达式具有LS解的形式,本文给出的参数估计表达式结构与之相同,这也是本文算法在上述几种加权或者结构总体最小估计算法中迭代效率最高的原因。

    • 自回归分析是统计上处理时间序列的一种常用方法。文献[20]给出了线性自回归模型的一个测量实例。该实例为对某一建筑物进行了36期沉降观测(如表 3所示[20]),并且认为不同期观测独立等精度。

      表 3  高程数据

      Table 3.  Elevation Data

      点号 高程/mm 点号 高程/mm 点号 高程/mm 点号 高程/mm
      1 26.33 10 25.46 19 28.09 28 26.57
      2 26.27 11 26.12 20 26.78 29 28.36
      3 26.43 12 27.28 21 28.66 30 27.94
      4 25.56 13 26.67 22 26.75 31 26.81
      5 26.82 14 27.95 23 27.24 32 28.50
      6 26.56 15 26.74 24 28.02 33 27.68
      7 25.93 16 27.53 25 26.81 34 26.57
      8 26.43 17 25.31 26 28.50 35 28.36
      9 26.52 18 26.90 27 27.68 36 27.94

      为便于比较,参照文献[20],采用三阶线性自回归模型进行建模。对应的EIV模型如下:

      $$ \left\{\begin{array}{l}{H}_{4}^{\mathrm{*}}={H}_{1}^{\mathrm{*}}{\xi }_{1}+{H}_{2}^{\mathrm{*}}{\xi }_{2}+{H}_{3}^{\mathrm{*}}{\xi }_{3}\\ {H}_{5}^{\mathrm{*}}={H}_{2}^{\mathrm{*}}{\xi }_{1}+{H}_{3}^{\mathrm{*}}{\xi }_{2}+{H}_{4}^{\mathrm{*}}{\xi }_{3}\\ ⋮\\ {H}_{36}^{\mathrm{*}}={H}_{33}^{\mathrm{*}}{\xi }_{1}+{H}_{34}^{\mathrm{*}}{\xi }_{2}+{H}_{35}^{\mathrm{*}}{\xi }_{3}\end{array}\right. $$ (25)

      式中,$ {\xi }_{1} $、$ {\xi }_{2} $、$ {\xi }_{3} $表示自回归模型的参数;$ {H}_{i}^{\mathrm{*}} $表示第i期沉降观测值,其中,i=1,2…36;“*”表示随机量。

      利用本文算法进行解算时,还需要给出$ \mathit{\boldsymbol{\psi }} $,系数矩阵的增广矩阵C可表示为:

      $$ \mathit{\boldsymbol{C}}=\left[\begin{array}{cccc}{H}_{1}^{\mathrm{*}}& {H}_{2}^{\mathrm{*}}& {H}_{3}^{\mathrm{*}}& {H}_{4}^{\mathrm{*}}\\ {H}_{2}^{\mathrm{*}}& {H}_{3}^{\mathrm{*}}& {H}_{4}^{\mathrm{*}}& {H}_{5}^{\mathrm{*}}\\ ⋮& ⋮& ⋮& ⋮\\ {H}_{33}^{\mathrm{*}}& {H}_{34}^{\mathrm{*}}& {H}_{35}^{\mathrm{*}}& {H}_{36}^{\mathrm{*}}\end{array}\right] $$ (26)

      根据式(2)构造$ \mathit{\boldsymbol{C}} $的仿射结构矩阵$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{36} $,$ \mathit{\boldsymbol{C}} $具有Hankel矩阵结构。不难发现,依次保留$ \mathit{\boldsymbol{C}} $各主对角线元素并置为1,即为$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{36} $,进而可得$ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{ccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{36}\right)\end{array}\right] $。

      分别采用1种WTLS估计算法(文献[19]算法)和4种STLS估计算法(文献[10]算法、文献[21]算法、文献[15]算法、本文算法)对表 3中数据进行求解,设置收敛阈值为$ 1\times {10}^{-10} $,相应的参数估值($ {\widehat{\xi }}_{1} $、$ {\widehat{\xi }}_{2} $和$ {\widehat{\xi }}_{3} $)、验后单位权方差估值$ {\widehat{\sigma }}_{0}^{2} $、参数估计精度($ {\widehat{\sigma }}_{{\widehat{\xi }}_{1}}^{2} $、$ {\widehat{\sigma }}_{{\widehat{\xi }}_{2}}^{2} $和$ {\widehat{\sigma }}_{{\widehat{\xi }}_{3}}^{2} $)及迭代次数如表 4所示。

      表 4  线性自回归模型不同算法所得结果

      Table 4.  Results from Different Methods of Linear Auto-Regression Model

      参数 文献[19]算法 文献[10]算法 文献[21]算法 文献[15]算法 本文算法
      $ {\widehat{\xi }}_{1} $ 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2
      $ {\widehat{\xi }}_{2} $ 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4
      $ {\widehat{\xi }}_{3} $ -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2
      $ \widehat{\sigma }{}_{0}^{2} $ 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1
      $ \widehat{\sigma }{}_{{\widehat{\xi }}_{1}}^{2} $ 0.012 658 067 0 0.012 658 067 0 0.012 658 067 0 0.012 658 067 0
      $ \widehat{\sigma }{}_{{\widehat{\xi }}_{2}}^{2} $ 0.009 481 774 9 0.009 481 774 9 0.009 481 774 9 0.009 481 774 9
      $ \widehat{\sigma }{}_{{\widehat{\xi }}_{3}}^{2} $ 0.009 074 553 9 0.009 074 553 9 0.009 074 553 9 0.009 074 553 9
      迭代次数 43 43 43 60 43

      算例1中的直线拟合问题仅是系数矩阵存在固定列,为进一步说明本文算法的有效性,选取了线性自回归问题,该问题中系数矩阵的增广矩阵为Hankel矩阵。由表 4可知,各种算法估计结果具有一致性,验证了本文算法的有效性。在迭代效率方面,文献[15]算法需要60次迭代收敛;其余各种算法优于文献[15]算法,均需要43次迭代收敛,可见,本文给出的STLS估计算法迭代效率同样优于文献[15]算法。

      通过提取独立随机量建立结构矩阵,STLS问题可以利用各种STLS估计算法[10, 15, 21]进行求解。文献[10]算法、文献[21]算法以及本文算法参数估计表达式和结构矩阵提取方式如表 5所示。

      表 5  参数估计表达式比较

      Table 5.  Comparison of Parameter Estimation Expression

      算法 参数估计表达式 结构矩阵提取方法
      文献[10]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \widehat{\mathit{\boldsymbol{B}}}={\widehat{\mathit{\boldsymbol{x}}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}\otimes {\mathit{\boldsymbol{I}}}_{m} $,$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{l}}l}=\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\mathit{\boldsymbol{S}}}^{\mathrm{T}} $;$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵;$ \widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{S}} $为结构矩阵 $ \mathit{\boldsymbol{S}}\cdot \mathit{\boldsymbol{e}}=\left[\begin{array}{cc}{\mathit{\boldsymbol{C}}}_{\mathit{\boldsymbol{A}}}& 0\\ 0& {I}_{m}\end{array}\right]\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $
      文献[21]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,结构矩阵$ {\widehat{\mathit{\boldsymbol{C}}}}_{2}=\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{C}}}_{1} $ $ {\mathit{\boldsymbol{C}}}_{1}\cdot \mathit{\boldsymbol{e}}={\mathit{\boldsymbol{C}}}_{1}\cdot \left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $
      本文算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \mathit{\boldsymbol{\Gamma }}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $,其中,结构矩阵$ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{\psi }} $ $ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{cccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)& \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}\right)\end{array}\right] $通过变量投影法建立$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}} $
      文献[19]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\right) $;$ {\mathit{\boldsymbol{Q}}}_{1} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵

      表 5可知,文献[10]算法、文献[21]算法和本文算法参数估计表达式的结构与文献[19]算法均是相同的,所以它们的迭代效率相同,这种参数估计表达式更早见于文献[22-23]。但是3种方法的结构矩阵提取方法存在差别,其中,文献[10]算法和文献[21]算法结构矩阵提取时,并不考虑系数矩阵和观测向量包含相同随机量的情况,因而,在处理自回归模型时需要利用文献[8]中给出的方法生成秩亏的观测权阵。理论上,文献[21]算法给出的矩阵$ {\mathit{\boldsymbol{C}}}_{1} $形式并不唯一,其中一种形式即为文献[10]算法中的矩阵$ \mathit{\boldsymbol{S}} $。当然各种算法公式推导过程等方面也存在差异,详情参考文献[10, 21]。

    • 文献[15]引入变量投影法构造结构矩阵,这种方法原理简单,规律性强,独立于求解问题。但是文献[15]的算法推导过程相对复杂、计算效率较低,并且不能给出参数的近似精度估计。为此,本文给出了一种将系数矩阵的增广矩阵中随机量纳入随机模型的新思路,通过变量投影法对函数模型进行重构,将系数矩阵的增广矩阵表示为独立随机量的仿射函数形式,采用非线性拉格朗日乘数法得到了一种STLS估计算法,并给出了单位权方差和参数的一阶近似协因数阵估计公式。

      1)两个算例表明,本文算法与已有能够解决系数矩阵和观测向量存在结构特征的加权或者结构总体最小二乘算法估计结果一致,说明了本文算法的有效性,并且相比文献[15]的算法提高了迭代效率。

      2)与通过构造秩亏观测权阵的各种WTLS估计算法相比,本文算法仅需要建立系数矩阵的增广矩阵中独立随机量的正定观测权阵,降低了待估量个数。与其他STLS估计算法相比,本文引入了变量投影法这一结构矩阵构造方法,顾及了EIV模型中系数矩阵和观测向量的整体结构,即考虑了系数矩阵和观测向量包含相同随机量的情况,因而,在处理诸如自回归模型等问题时,无需再建立秩亏的观测权阵。

      3)通过对各种算法的比较可知,STLS估计和WTLS估计分别从从函数模型和随机模型的角度确保TLS估计的统计最优性,WTLS估计和STLS估计具有等价性。

      4)本文给出的STLS算法具有一般性,适用于系数矩阵和观测向量包含随机量、固定量以及固定列等情况,当系数矩阵全为固定值,本文算法转变成加权LS估计算法,当系数矩阵和观测向量仅包含随机量,本文算法等价于WTLS估计算法。

参考文献 (23)

目录

    /

    返回文章
    返回