留言板

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

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

附有相对权比的PEIV模型总体最小二乘平差

王乐洋 许光煜 陈晓勇

王乐洋, 许光煜, 陈晓勇. 附有相对权比的PEIV模型总体最小二乘平差[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
引用本文: 王乐洋, 许光煜, 陈晓勇. 附有相对权比的PEIV模型总体最小二乘平差[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
WANG Leyang, XU Guangyu, CHEN Xiaoyong. Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
Citation: WANG Leyang, XU Guangyu, CHEN Xiaoyong. Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001

附有相对权比的PEIV模型总体最小二乘平差

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

国家自然科学基金 41664001

国家自然科学基金 41204003

江西省杰出青-人才资助计划项目 20162BCB23050

国家重点研发计划 2016YFB0501405

测绘地理信息公益性行业科研专项 201512026

江西省教育厅科技项目 GJJ150595

流域生态与地理环境监测国家测绘地理信息局重点实验室项目 WE2015005

对地观测技术国家测绘地理信息局重点实验室项目 K201502

东华理工大学博士科研启动金 DHBK201113

详细信息
    作者简介:

    王乐洋, 博士, 副教授, 主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P207

Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor

Funds: 

The National Natural Science Foundation of China 41664001

The National Natural Science Foundation of China 41204003

the Outstanding Youth Talents in Jiangxi Province 20162BCB23050

the National Key Research and Development Program 2016YFB0501405

the National Department Public Benefit Research Foundation (Surveying, Mapping and Geoinformation) 201512026

the Science and Technology Project of the Education Department of Jiangxi Province GJJ150595

the Key Laboratory of Watershed Ecology and Geographical Environment Monitoring WE2015005

the Project of Key Laboratory of Mapping from Space, NASG K201502

the Scientific Research Foundation of ECIT DHBK201113

More Information
    Author Bio:

    WANG Leyang, PhD, associate professor, specializs in geodetic inversion and geodetic data processing. E-mail:wleyang@163.com

  • 摘要: 针对观测向量和系数矩阵权分配不合理、验前随机模型不准确的情况,以部分误差变量(partial errors-in-variables,PEIV)模型为基础,推导了附有相对权比的总体最小二乘平差算法;通过在平差准则中加入相对权比,自适应调整观测向量和系数矩阵随机元素对模型参数估计的贡献,给出了确定相对权比的验前单位权方差法和判别函数最小化迭代算法,该算法普遍适用于一般性的系数矩阵和权矩阵。通过直线拟合和坐标转换模拟算例的比较分析,发现当观测值和系数矩阵的验前单位权方差已知,且较准确时,验前单位权方差法确定相对权比和参数估计的效果较好;而以${{\overline{\mathit{{\mathit{\Phi}}}}}_{1}}\left( \hat{\varepsilon },{{{\hat{\varepsilon }}}_{a}} \right)={{\hat{\varepsilon }}^{\text{T}}}\hat{\varepsilon }+\hat{\varepsilon }_{a}^{\text{T}}{{\hat{\varepsilon }}_{a}} $作为判别函数是判别函数最小化迭代算法中效果最好的。
  • 图  1  判别函数最小化法流程图

    Figure  1.  Flowchart of Minimum Discriminate Function Method

    图  2  直线拟合算例各方案判别函数与相对权比之间的关系图

    Figure  2.  Relationship Between Discriminate Function and Weight Scaling Factor in Straight Line Fitting

    图  3  坐标转换算例各方案判别函数与相对权比之间的关系图

    Figure  3.  Relationship Between Discriminate Function and Weight Scaling Factor in Coordinate Transformation

    表  1  各方案所用方法

    Table  1.   Method of Each Project

    方案 方法
    1 最小二乘方法(即λ=1)
    2 加权总体最小二乘方法(即λ=0.5)
    3 验前单位权方差法λ=σ022/(σ012+σ022),即λ=0.6000
    4 ${{\overline {\mathit{\Phi}} }_1}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) = {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}\mathit{\boldsymbol{\hat \varepsilon }} + \mathit{\boldsymbol{\hat \varepsilon }}_a^{\rm{T}}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}} $为判别函数
    5 $ {{\overline {\mathit{\Phi}} }_2}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) = {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon }} + \mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_a}$为判别函数
    6 ${{\overline {\mathit{\Phi}} }_3}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = \lambda {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon }} + \left( {1-\lambda } \right)\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_a} $为判别函数
    7 $ {{\overline {\mathit{\Phi}} }_4}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = \sum\limits_{i = 1}^n {\left| {{{\mathit{\boldsymbol{\hat \varepsilon }}}_i}} \right|} + \sum\limits_{i = 1}^t {\left| {{{\mathit{\boldsymbol{\hat \varepsilon }}}_{{\mathit{\boldsymbol{a}}_i}}}} \right|} $为判别函数
    下载: 导出CSV

    表  2  模型参数的计算结果

    Table  2.   Results of Model Parameter

    方案 $ {{\hat \beta }_1}$ ${{\hat \beta }_2} $ λ $\left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\| $
    真值 2.000 0 1.500 0 / /
    1 1.864 9 2.597 1 1.000 0 1.784 4
    2 2.006 8 1.523 1 0.500 0 0.980 9
    3 2.016 9 1.465 3 0.600 0 0.965 1
    4 2.018 8 1.455 0 0.620 0 0.963 2
    5 2.006 8 1.523 1 0.500 0 0.980 9
    6 2.043 9 1.323 1 0.990 0 1.003 1
    7 2.043 9 1.323 1 0.990 0 1.003 1
    下载: 导出CSV

    表  3  原始坐标系和目标坐标系中13个点的坐标值

    Table  3.   Simulated Points in the Start System and Target System for the Affine Transformation

    点号 1 2 3 4 5 6 7 8 9 10 11 12 13
    原始坐标
    Xs 1.0 1.0 1.0 1.5 2.0 2.5 3.0 3.0 3.0 2.5 2.0 1.5 1.0
    Ys 1.0 2.0 3.0 3.5 4.0 3.5 3.0 2.0 1.0 0.5 0.0 0.5 1.0
    目标坐标
    Xt 1.1 0.3 -0.5 -0.45 -0.4 0.45 1.3 2.1 2.9 2.85 2.8 1.95 1.1
    Yt 6.3 7.0 7.7 8.35 9.0 8.95 8.9 8.2 7.5 6.85 6.2 6.25 6.3
    下载: 导出CSV

    表  4  模型参数的计算结果

    Table  4.   Results of Model Parameter

    方案 $ {{\hat a}_1}$ $ {{\hat b}_1}$ ${{\hat c}_1} $ $ {{\hat a}_2}$ $ {{\hat b}_2}$ $ {{\hat c}_2}$ λ $\left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\| $
    真值 0.900 0 -0.800 0 1.000 0 0.600 0 0.700 0 5.000 0 / /
    1 0.849 6 -0.776 1 1.046 3 0.553 6 0.697 3 5.093 9 1.000 0 0.341 7
    2 0.894 0 -0.793 6 1.002 1 0.591 8 0.700 5 5.016 1 0.500 0 0.298 3
    3 0.899 3 -0.795 7 0.996 4 0.597 8 0.701 9 5.002 1 0.600 0 0.297 2
    4 0.895 1 -0.794 0 1.000 9 0.593 0 0.700 7 5.013 3 0.520 0 0.297 9
    5 0.894 0 -0.793 6 1.002 1 0.591 8 0.700 5 5.016 1 0.500 0 0.298 3
    6 0.921 2 -0.802 7 0.965 4 0.627 2 0.715 2 4.920 5 0.990 0 0.317 0
    7 0.921 2 -0.802 7 0.965 4 0.627 2 0.715 2 4.920 5 0.990 0 0.317 0
    下载: 导出CSV
  • [1] Adcock R J. Note on the Method of Least Squares[J]. The Analyst, 1877, 4(6):183-184 doi:  10.2307/2635777
    [2] Pearson K. On Line and Planes of Closest Fit to Systems of Points in Space[J]. Philosophical Magazine Series 6, 1901, 2(11):559-572 doi:  10.1080/14786440109462720
    [3] 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
    [4] 王乐洋.总体最小二乘解性质研究[J].大地测量与地球动力学, 2012, 32(5):48-52 http://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201205011.htm

    Wang Leyang. Research on Properties of Total Least Squares Estimation[J]. Journal of Geodesy and Geodynamics, 2012, 32(5):48-52 http://www.cnki.com.cn/Article/CJFDTOTAL-DKXB201205011.htm
    [5] 王乐洋, 许才军.总体最小二乘研究进展[J].武汉大学学报·信息科学版, 2013, 38(7):850-856 http://ch.whu.edu.cn/CN/abstract/abstract2703.shtml

    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/CN/abstract/abstract2703.shtml
    [6] 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
    [7] 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
    [8] 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
    [9] Xu Peiliang, Liu Jingnan, 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
    [10] Amiri-Simkooei A R. Application of Least Squares Variance Component Estimation to Errors-in-Variables Models[J]. Journal of Geodesy, 2013, 87(10-12):935-944 doi:  10.1007/s00190-013-0658-8
    [11] Xu Peiliang, Liu Jingnan. Variance Components in Errors-in-Variables Models:Estimability, Stability and Bias Analysis[J].Journal of Geodesy, 2014, 88(8):719-734 doi:  10.1007/s00190-014-0717-9
    [12] 王乐洋, 许才军.附有相对权比的总体最小二乘平差[J].武汉大学学报·信息科学版, 2011, 36(8):887-890 http://ch.whu.edu.cn/CN/abstract/abstract629.shtml

    Wang Leyang, Xu Caijun. Total Least Squares Adjustment with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2011, 36(8):887-890 http://ch.whu.edu.cn/CN/abstract/abstract629.shtml
    [13] 杨元喜, 景一帆, 曾安敏.自适应参数估计与内外部精度的关系[J].测绘学报, 2014, 43(5):441-445 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201405002.htm

    Yang Yuanxi, Jing Yifan, Zeng Anmin. Adaptive Parameter Estimation and Inner and External Precision[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(5):441-445 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201405002.htm
    [14] 王乐洋. 基于总体最小二乘的大地测量反演理论及应用研究[D]. 武汉: 武汉大学, 2011

    Wang Leyang. Research on Theory and Application of Total Least Squares in Geodetic Inversion[D]. Wuhan:Wuhan University, 2011
    [15] 王乐洋, 许才军, 张朝玉.一种确定联合反演中相对权比的两步法[J].测绘学报, 2012, 41(1):19-24 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201201006.htm

    Wang Leyang, Xu Caijun, Zhang Chaoyu. A Two-step Method to Determine Relative Weight Ratio Factors in Joint Inversion[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(1):19-24 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201201006.htm
    [16] Mahboub V. On Weighted Total Least-squares for Geodetic Transformation[J]. Journal of Geod-esy, 2012, 86(5):359-367 doi:  10.1007/s00190-011-0524-5
  • [1] 王乐洋, 谷旺旺, 赵雄, 许光煜, 高华.  利用偏差改正的方差分量估计方法确定联合反演相对权比 . 武汉大学学报 ● 信息科学版, 2022, 47(4): 508-516. doi: 10.13203/j.whugis20200216
    [2] 曾安敏, 明锋, 吴富梅.  坐标参考框架长期解内约束模型 . 武汉大学学报 ● 信息科学版, 2022, 47(9): 1447-1451. doi: 10.13203/j.whugis20190453
    [3] 王乐洋, 邹传义.  PEIV模型参数估计理论及其应用研究进展 . 武汉大学学报 ● 信息科学版, 2021, 46(9): 1273-1283, 1297. doi: 10.13203/j.whugis20200312
    [4] 王乐洋, 孙坚强.  总体最小二乘回归预测模型的方差分量估计 . 武汉大学学报 ● 信息科学版, 2021, 46(2): 280-288. doi: 10.13203/j.whugis20180450
    [5] 谢建, 龙四春, 周璀.  不等式约束PEIV模型的最优性条件及SQP算法 . 武汉大学学报 ● 信息科学版, 2020, 45(7): 1002-1007. doi: 10.13203/j.whugis20180297
    [6] 赵俊, 郭飞霄, 李琦.  PEIV模型WTLS估计的Fisher-Score算法 . 武汉大学学报 ● 信息科学版, 2019, 44(2): 214-220. doi: 10.13203/j.whugis20170061
    [7] 吕志鹏, 隋立芬.  基于非线性高斯-赫尔默特模型的结构总体最小二乘法 . 武汉大学学报 ● 信息科学版, 2019, 44(12): 1808-1815. doi: 10.13203/j.whugis20180104
    [8] 王乐洋, 余航.  附有相对权比的加权总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2019, 44(8): 1233-1240. doi: 10.13203/j.whugis20170265
    [9] 王乐洋, 李海燕, 陈晓勇.  拟牛顿修正法解算不等式约束加权总体最小二乘问题 . 武汉大学学报 ● 信息科学版, 2018, 43(1): 127-132. doi: 10.13203/j.whugis20150333
    [10] 王乐洋, 余航.  火山Mogi模型反演的总体最小二乘联合平差方法 . 武汉大学学报 ● 信息科学版, 2018, 43(9): 1333-1341. doi: 10.13203/j.whugis20160469
    [11] 王乐洋, 吴飞, 吴良才.  GPS高程转换的总体最小二乘拟合推估模型 . 武汉大学学报 ● 信息科学版, 2016, 41(9): 1259-1264. doi: 10.13203/j.whugis20140421
    [12] 王乐洋, 余航.  总体最小二乘联合平差 . 武汉大学学报 ● 信息科学版, 2016, 41(12): 1683-1689. doi: 10.13203/j.whugis20140670
    [13] 龚循强, 李志林.  一种利用IGGII方案的稳健混合总体最小二乘方法 . 武汉大学学报 ● 信息科学版, 2014, 39(4): 462-466. doi: 10.13203/j.whugis20120015
    [14] 黄令勇, 吕志平, 任雅奇, 陈正生, 王宇谱.  多元总体最小二乘在三维坐标转换中的应用 . 武汉大学学报 ● 信息科学版, 2014, 39(7): 793-798.
    [15] 王乐洋, 许才军.  总体最小二乘研究进展 . 武汉大学学报 ● 信息科学版, 2013, 38(7): 850-856.
    [16] 王乐洋, 许才军.  附有相对权比的总体最小二乘平差 . 武汉大学学报 ● 信息科学版, 2011, 36(8): 887-890.
    [17] 孔建, 姚宜斌, 吴寒.  整体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(6): 711-714.
    [18] 鲁铁定, 周世健.  总体最小二乘的迭代解法 . 武汉大学学报 ● 信息科学版, 2010, 35(11): 1351-1354.
    [19] 吴波, 汪小钦, 张良培.  端元光谱自动提取的总体最小二乘迭代分解 . 武汉大学学报 ● 信息科学版, 2008, 33(5): 457-460.
    [20] 徐天河, 杨元喜.  坐标转换模型尺度参数的假设检验 . 武汉大学学报 ● 信息科学版, 2001, 26(1): 70-74.
  • 加载中
图(3) / 表(4)
计量
  • 文章访问数:  1684
  • HTML全文浏览量:  168
  • PDF下载量:  295
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-02-04
  • 刊出日期:  2017-06-05

附有相对权比的PEIV模型总体最小二乘平差

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

    国家自然科学基金 41664001

    国家自然科学基金 41204003

    江西省杰出青-人才资助计划项目 20162BCB23050

    国家重点研发计划 2016YFB0501405

    测绘地理信息公益性行业科研专项 201512026

    江西省教育厅科技项目 GJJ150595

    流域生态与地理环境监测国家测绘地理信息局重点实验室项目 WE2015005

    对地观测技术国家测绘地理信息局重点实验室项目 K201502

    东华理工大学博士科研启动金 DHBK201113

    作者简介:

    王乐洋, 博士, 副教授, 主要研究方向为大地测量反演及大地测量数据处理。wleyang@163.com

  • 中图分类号: P207

摘要: 针对观测向量和系数矩阵权分配不合理、验前随机模型不准确的情况,以部分误差变量(partial errors-in-variables,PEIV)模型为基础,推导了附有相对权比的总体最小二乘平差算法;通过在平差准则中加入相对权比,自适应调整观测向量和系数矩阵随机元素对模型参数估计的贡献,给出了确定相对权比的验前单位权方差法和判别函数最小化迭代算法,该算法普遍适用于一般性的系数矩阵和权矩阵。通过直线拟合和坐标转换模拟算例的比较分析,发现当观测值和系数矩阵的验前单位权方差已知,且较准确时,验前单位权方差法确定相对权比和参数估计的效果较好;而以${{\overline{\mathit{{\mathit{\Phi}}}}}_{1}}\left( \hat{\varepsilon },{{{\hat{\varepsilon }}}_{a}} \right)={{\hat{\varepsilon }}^{\text{T}}}\hat{\varepsilon }+\hat{\varepsilon }_{a}^{\text{T}}{{\hat{\varepsilon }}_{a}} $作为判别函数是判别函数最小化迭代算法中效果最好的。

English Abstract

王乐洋, 许光煜, 陈晓勇. 附有相对权比的PEIV模型总体最小二乘平差[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
引用本文: 王乐洋, 许光煜, 陈晓勇. 附有相对权比的PEIV模型总体最小二乘平差[J]. 武汉大学学报 ● 信息科学版, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
WANG Leyang, XU Guangyu, CHEN Xiaoyong. Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
Citation: WANG Leyang, XU Guangyu, CHEN Xiaoyong. Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor[J]. Geomatics and Information Science of Wuhan University, 2017, 42(6): 857-863. doi: 10.13203/j.whugis20150001
  • 总体最小二乘平差理论及应用研究是目前测绘数据处理领域的热点研究问题,总体最小二乘准则最早由Adcock[1]提出,Pearson[2]推导得出了总体最小二乘数值解法。Golub和Van Loan[3]提出了著名的SVD算法,并正式命名总体最小二乘估计方法(TLS)。普通总体最小二乘方法能够有效解决误差变量(errors-in-variables,EIV)模型问题,相对经典的Gauss-Markov模型,EIV模型考虑了系数矩阵也含有随机误差的情况[4, 5];针对EIV模型中普遍存在观测值精度不同且相关的情况,Schaffrin等提出了加权总体最小二乘算法(weighted total least squares,WTLS)[6, 7],Fang研究了一般性权矩阵下的WTLS算法[8];Xu等提出了基于部分变量误差(partial errors-in-variables,PEIV)模型的WTLS算法,形成了一个EIV问题的统一求解模式[9]。PEIV模型是EIV模型的一种扩展形式,主要适用于系数矩阵中同时存在常数元素和随机元素的情况。在WTLS算法中,均认为观测向量和系数矩阵的验前随机模型具有相同的单位权方差,但在实际问题中,观测向量和系数矩阵元素往往来自不同的观测和平差,二者的验前单位权方差也是不相同的,此时观测值和系数矩阵的权阵分配并不合理。针对这种情况,文献[10, 11]在加权总体最小二乘法中采用方差分量估计方法,直接求解出观测向量和系数矩阵的单位权方差估值来确定观测向量和系数矩阵的合理权阵;文献[12]提出附有相对权比的总体最小二乘平差方法,通过在平差准则中加入相对权比,自适应调整观测向量和系数矩阵随机元素对模型参数估计的贡献[13];文献[14]提出了基于总体最小二乘的联合反演方法,并用于长白山天池火山区垂直位移和水平位移的联合反演,为了简化计算对两类数据附加相同的权比。

    本文在文献[12]的基础上,对函数模型进行推广,以PEIV模型为基础,提出了附有相对权比的PEIV模型总体最小二乘平差算法,该算法普遍适用于一般性的系数矩阵和权矩阵。同时推导了确定相对权比的验前单位权方差法和判别函数最小化迭代算法;在判别函数最小化迭代算法中,加入多种确定相对权比的判别函数进行比较和分析。直线拟合和坐标转换的模拟算例表明,当观测值和系数矩阵的验前单位权方差已知且较准确时,验前单位权方差法效果较好;判别函数最小化迭代算法中在确定相对权比时选取方案4(${{\overline{\mathit{{\mathit{\Phi}}} }}_{1}}\left( \boldsymbol{\hat{\varepsilon }},{{{\boldsymbol{\hat{\varepsilon }}}}_{a}} \right)={{\boldsymbol{\hat{\varepsilon }}}^{\text{T}}}\boldsymbol{\hat{\varepsilon }}+\boldsymbol{\hat{\varepsilon }}_{a}^{\text{T}}{{\boldsymbol{\hat{\varepsilon }}}_{a}} $)作为判别函数,算法效果显著。

    • 线性平差模型[4-6]

      $$ \mathit{\boldsymbol{y = \bar A\beta }} $$ (1)

      式中,$\boldsymbol{\bar y} $为n×1观测向量y的真值;${\boldsymbol{\bar A}} $为n×m系数矩阵A的真值;βm×1的待估模型参数。

      PEIV函数模型可以表示为[9]:

      $$ \begin{gathered} \left\{ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{y}} = \left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\left( {\mathit{\boldsymbol{h + B\bar a}}} \right) + \mathit{\boldsymbol{\varepsilon }}} \\ {\mathit{\boldsymbol{a = \bar a + }}{\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;} \end{array}} \right. \hfill \\ \;\;\;\;\;\;\mathit{\boldsymbol{h + B\bar a}} = {\rm{vec}}\left( {\mathit{\boldsymbol{\bar A}}} \right) \hfill \\ \end{gathered} $$ (2)

      式中,$\mathit{\boldsymbol{y}} = \mathit{\boldsymbol{\bar y}} + \mathit{\boldsymbol{\varepsilon }} $,$\mathit{\boldsymbol{A = \bar A}} + {\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{A}}} $,εεA分别为观测向量和系数矩阵的随机误差向量和随机误差矩阵;h是系数矩阵$ {\mathit{\boldsymbol{\bar A}}}$中非随机元素组成的nm×1列向量;Bnm×t的固定矩阵,其阶数取决于系数矩阵${\mathit{\boldsymbol{\bar A}}} $中随机元素的数目(此处假设为t);a是系数矩阵A中随机元素组成的t×1列向量,其真值用${\mathit{\boldsymbol{\bar a}}} $表示,εa为列向量a的随机误差向量,vec(·)为矩阵拉直运算,⊗表示Kronecker直积。

      其平差的随机模型为:

      $$ \left\{ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{D}}\left( \mathit{\boldsymbol{\varepsilon }} \right) = \sigma _{01}^2{\mathit{\boldsymbol{Q}}_{\varepsilon \varepsilon }}} \\ {\mathit{\boldsymbol{D}}\left( {{\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}}} \right) = \sigma _{02}^2{\mathit{\boldsymbol{Q}}_{{\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}}{\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}}}}} \end{array}} \right. $$ (3)

      式中,D表示方差-协方差矩阵;σ012为观测向量的单位权方差;Qε ε=P1-1为观测向量的协因数阵;σ022为系数矩阵元素的单位权方差;Qεaεa=P2-1为随机误差向量εa的协因数阵;通常σ012σ022

      假设观测值与系数矩阵元素是相互独立的:

      $$ {\rm{Cov}}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = 0 $$ (4)

      通过观测向量y和系数矩阵元素a估计PEIV模型中$ {\mathit{\boldsymbol{\bar a}}}$和β的附有相对权比的PEIV模型总体最小二乘平差准则为:

      $$ \varphi = \lambda {\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\varepsilon }} + \left( {1-\lambda } \right)\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}} = \min $$ (5)

      式中,λ为相对权比[12],0 < λ < 1。

      将式(2) 代入式(5) 得附有相对权比的PEIV模型总体最小二乘平差准则为:

      $$ \begin{gathered} \varphi = \lambda {\left\{ {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\left( {\mathit{\boldsymbol{h + B\bar a}}} \right)-\mathit{\boldsymbol{y}}} \right\}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\left\{ {\left( {{\mathit{\boldsymbol{\beta }}^{\rm{T}}} \otimes {\mathit{\boldsymbol{I}}_n}} \right)\left( {\mathit{\boldsymbol{h + B\bar a}}} \right)-\mathit{\boldsymbol{y}}} \right\} \hfill \\ + \left( {1-\lambda } \right){\left( {\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{P}}_2}\left( {\mathit{\boldsymbol{\bar a}} - \mathit{\boldsymbol{a}}} \right) = \min \hfill \\ \end{gathered} $$ (6)

      设$ {\mathit{\boldsymbol{\bar a}}}$和β的加权总体最小二乘估值分别为${\mathit{\boldsymbol{\hat{\bar{a}}}}} $和$ {\boldsymbol{\hat \beta }}$,则φ对$ {\mathit{\boldsymbol{\bar a}}}$和β求偏导数并令其等于零得

      $$ \begin{align} &\frac{\partial \varphi }{\partial \boldsymbol{\bar{a}}}=2\lambda {{\boldsymbol{B}}^{\text{T}}}\left( {\boldsymbol{\hat \beta }}\otimes {{\boldsymbol{I}}_{\text{n}}} \right){{\boldsymbol{P}}_{1}}\left\{ \left( {\boldsymbol{\hat \beta }}{{}^{\text{T}}}\otimes {{\boldsymbol{I}}_{n}} \right)\left( \boldsymbol{h}+\boldsymbol{B\hat{\bar{a}}} \right)-\boldsymbol{y} \right\} \\ &\text{+2}\left( \text{1-}\lambda \right){{\boldsymbol{P}}_{\text{2}}}\left( \boldsymbol{\hat{\bar{a}}}\text{-}\boldsymbol{a} \right)\text{=0} \\ \end{align} $$ (7)
      $$ \begin{align} & \frac{\partial \varphi }{\partial \boldsymbol{ \beta }}\rm{=}2\lambda \left\{ \left[\begin{matrix} \mathit{\boldsymbol{h}}_{1}^{\rm{T}} \\ \mathit{\boldsymbol{h}}_{2}^{\rm{T}} \\ \vdots \\ \mathit{\boldsymbol{h}}_{m}^{\rm{T}} \\ \end{matrix} \right]\rm{+}\left[\begin{matrix} {{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{1}^{\rm{T}} \\ {{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{2}^{\rm{T}} \\ \vdots \\ {{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{m}^{\rm{T}} \\ \end{matrix} \right] \right\} \\ & {{\mathit{\boldsymbol{P}}}_{1}}\left\{ \left( {{{\hat{\beta }}}^{\rm{T}}}\otimes {{\mathit{\boldsymbol{I}}}_{n}} \right)\cdot \left( \mathit{\boldsymbol{h+B\hat{\bar{a}}}} \right)-\mathit{\boldsymbol{y}} \right\}=0 \\ \end{align} $$ (8)

      式中,h=[h1T   h2T  …  hmT]Thin×1的列向量;B=[B1T  B2T  …  BmT]TBin×t阶矩阵。

      由式(7) 和式(8) 得:

      $$ \mathit{\boldsymbol{\hat{\bar{a}}}}={{\left\{ \lambda \mathit{\boldsymbol{S}}_\boldsymbol{\beta }^{\rm{T}}{{\mathit{\boldsymbol{P}}}_{1}}{{\mathit{\boldsymbol{S}}}_\boldsymbol{\beta }}+\left( 1-\lambda \right){{\mathit{\boldsymbol{P}}}_{2}} \right\}}^{-1}}\left\{ \left( 1-\lambda \right){{\mathit{\boldsymbol{P}}}_{2}}\mathit{\boldsymbol{a}}-\lambda \mathit{\boldsymbol{S}}_\boldsymbol{\beta }^{\rm{T}}{{\mathit{\boldsymbol{P}}}_{1}}{{\mathit{\boldsymbol{S}}}_{\mathit{\boldsymbol{h}}}}+\lambda \mathit{\boldsymbol{S}}_\boldsymbol{\beta }^{\rm{T}}{{\mathit{\boldsymbol{P}}}_{1}}\mathit{\boldsymbol{y}} \right\} $$ (9)
      $$ \mathit{\boldsymbol{\hat \beta = }}{\left( {{\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{h}}} + {\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{B}}} + {\mathit{\boldsymbol{N}}_{\mathit{\boldsymbol{Bh}}}} + {\mathit{\boldsymbol{N}}_\mathit{\boldsymbol{h}}}_\mathit{\boldsymbol{B}}} \right)^{-1}}\left( {{\mathit{\boldsymbol{u}}_h} + {\mathit{\boldsymbol{u}}_B}} \right) $$ (10)

      式中,Sβ=($\mathit{\boldsymbol{\hat \beta }} $TIn)BSh=($\mathit{\boldsymbol{\hat \beta }} $TIn)hNhNBNBhNhB都是m×m的矩阵,矩阵元素分别为:Nh(i, j)=hiTP1hjNB(i, j)= ${{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{i}^{\rm{T}}{{\mathit{\boldsymbol{P}}}_{1}}{{\mathit{\boldsymbol{B}}}_{j}}\mathit{\boldsymbol{\hat{\bar{a}}}} $,NBh(i, j)= ${{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{i}^{\rm{T}}{{\mathit{\boldsymbol{P}}}_{1}}{{\mathit{\boldsymbol{h}}}_{j}} $,NhB(i, j)=hiTP1Bj ${\mathit{\boldsymbol{\hat{\bar{a}}}}} $(i, j=1, 2, …, m)。uhuB分别为:

      $$ {{\mathit{\boldsymbol{u}}}_{\mathit{\boldsymbol{h}}}}=\left[\begin{matrix} \mathit{\boldsymbol{h}}_{1}^{\rm{T}} \\ \mathit{\boldsymbol{h}}_{1}^{\rm{T}} \\ \vdots \\ \mathit{\boldsymbol{h}}_{m}^{\rm{T}} \\ \end{matrix} \right]{{\mathit{\boldsymbol{P}}}_{1}}\mathit{\boldsymbol{y}}, {{\mathit{\boldsymbol{u}}}_{\mathit{\boldsymbol{B}}}}=\left[\begin{matrix} {{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{\rm{1}}^{\rm{T}} \\ {{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{2}^{\rm{T}} \\ \vdots \\ {{{\mathit{\boldsymbol{\hat{\bar{a}}}}}}^{\rm{T}}}\mathit{\boldsymbol{B}}_{m}^{\rm{T}} \\ \end{matrix} \right]{{\mathit{\boldsymbol{P}}}_{1}}\mathit{\boldsymbol{y}} $$

      单位权方差估值为:

      $$ \hat \sigma _0^2 = \frac{{\lambda {\mathit{\boldsymbol{ \hat \varepsilon }}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon }} + \left( {1-\lambda } \right)\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}}}{{n-m}} $$ (11)

      h=0,B=Inm时,PEIV模型退化为普通EIV模型,其解可以表示为:

      $$ \mathit{\boldsymbol{\hat{\bar{A}}=A+}}{{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{A}}}}\left( {\boldsymbol{\hat \beta }}\otimes {{\mathit{\boldsymbol{I}}}_{n}} \right)\mathit{\boldsymbol{Q}}_\boldsymbol{\varepsilon \varepsilon }^{-1}\left[\mathit{\boldsymbol{y-}}\left( \mathit{\boldsymbol{A}}-{{{\mathit{\boldsymbol{\hat{E}}}}}_{\mathit{\boldsymbol{A}}}} \right){\boldsymbol{\hat \beta }} \right] $$ (12)
      $$ {\boldsymbol{\hat \beta }}={{\left( \mathit{\boldsymbol{\hat{\bar{A}}}}{{\mathit{\boldsymbol{P}}}_{1}}\mathit{\boldsymbol{\hat{\bar{A}}}} \right)}^{-1}}\mathit{\boldsymbol{\hat{\bar{A}}}}{{\mathit{\boldsymbol{P}}}_{1}}y $$ (13)

      式中,QA=BQεaεaBT; $ {{\mathit{\boldsymbol{\hat E}}}_\mathit{\boldsymbol{A}}} = {\rm{ve}}{{\rm{c}}^{-1}}\left( {\mathit{\boldsymbol{h + B}}{{\mathit{\boldsymbol{\hat \varepsilon }}}_a}} \right)$; $ \boldsymbol{\hat{\bar{A}}=}\text{ve}{{\text{c}}^{-1}}\left( \boldsymbol{h+B\hat{\bar{a}}} \right)$, vec-1(·)表示拉直运算的逆运算, “^”表示估值。

      退化情况的条件设定为h=0,B=Inm,即表示此时系数矩阵A中不存在常数元素,且不存在随机元素重复出现的情况。普通EIV模型问题中一般需要求解的未知量包括观测向量yn个改正数,系数矩阵An×m个改正数,m个未知参数β,若在求解过程中构造了拉格朗日极值函数,则还需求解n×1的拉格朗日乘子;相对于Gauss-Markov模型的n+m个未知量,所需求解的未知量过多,也是制约EIV模型求解方法广泛应用于实际问题的一个重要因素,特别是对于海量数据处理问题。

      PEIV模型的构建在很大程度上减少了模型所需求解的未知量个数,从式(9) 和式(12) 的对比中可以看出,PEIV模型在求解过程中只需求解系数矩阵中t个随机元素的改正数;而EIV模型需要考虑系数矩阵中所有元素的改正数,同时在求解过程中需要采取一定方式固定系数矩阵中的常数元素,保证常数元素不发生改正。

      对式(10) 进行简单的变形可以得到:

      $$ \begin{align} & {{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{h}}}}+{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{B}}}}+{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{Bh}}}}+{{\mathit{\boldsymbol{N}}}_{\mathit{\boldsymbol{hB}}}}={{{\mathit{\boldsymbol{\hat{\bar{A}}}}}}^{\rm{T}}}{{\mathit{\boldsymbol{P}}}_{1}}\mathit{\boldsymbol{\hat{\bar{A}}}} \\ & \ \ \ \ \ \ {{\mathit{\boldsymbol{u}}}_{\mathit{\boldsymbol{h}}}}+{{\mathit{\boldsymbol{u}}}_{\mathit{\boldsymbol{B}}}}={{{\mathit{\boldsymbol{\hat{\bar{A}}}}}}^{\rm{T}}}{{\mathit{\boldsymbol{P}}}_{1}}\mathit{\boldsymbol{y}} \\ \end{align} $$

      结合式(13) 可以看出,在求解参数β的过程中,采用PEIV模型和EIV模型并无区别。

    • 将式(5) 改写为:

      $$ \begin{gathered} \varphi = \mathit{\boldsymbol{\varepsilon }}\left( {\lambda \cdot {\mathit{\boldsymbol{P}}_1}} \right)\mathit{\boldsymbol{\varepsilon }} + \mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}\left[{\left( {1-\lambda } \right){\mathit{\boldsymbol{P}}_2}} \right] \cdot {\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}} \hfill \\ = {\mathit{\boldsymbol{\varepsilon }}^{\rm{T}}}{{\mathit{\boldsymbol{\hat P}}}_1}\mathit{\boldsymbol{\varepsilon + \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{\mathit{\boldsymbol{\varepsilon }}_\mathit{\boldsymbol{a}}} \hfill \\ \end{gathered} $$ (14)

      式中,

      $$ {{\mathit{\boldsymbol{\hat P}}}_1} = \lambda \cdot {\mathit{\boldsymbol{P}}_1}, {{\mathit{\boldsymbol{\hat P}}}_2} = \left( {1-\lambda } \right) \cdot {\mathit{\boldsymbol{P}}_2} $$ (15)

      ${{\mathit{\boldsymbol{\hat P}}}_1} $和$ {{\mathit{\boldsymbol{\hat P}}}_2}$分别表示通过相对权比调整后的观测向量和系数矩阵中随机元素的权阵。

      通过简单公式推导可得:

      $$ \lambda = \sigma _{02}^2/(\sigma _{01}^2 + \sigma _{02}^2)\;, 0 < \lambda < 1 $$ (16)

      式(16) 即为文献[12]中验前单位权方差法确定相对权比的公式。

      由上面的推导可以看出,通过附加相对权比λ可以自适应调整观测向量和系数矩阵随机元素对模型参数估计的贡献。从式(16) 可知,当一个平差问题中,确定的相对权比0 < λ < 0.5时,可得到σ022 < σ012,即在给系数矩阵定权时所选取的单位权方差偏小,相对于观测向量中元素所分配到的权值,系数矩阵中元素所分配到的权值过大,从而使得在迭代计算过程中,计算得到的系数矩阵中元素的改正数小于其实际值;当确定的相对权比0.5 < λ 1时,则反之。当相对权比满足式(16) 时,调整后观测向量和系数矩阵随机元素的权阵${{\mathit{\boldsymbol{\hat P}}}_1} $和${{\mathit{\boldsymbol{\hat P}}}_2} $更为合理。

    • 在数据处理过程中,观测向量和系数矩阵的验前单位权方差σ012σ022往往是无法准确已知的,因此很难通过式(16) 直接求得相对权比的准确数值。此时可采用图 1所示的迭代算法求解相对权比以及相应的参数估值。

      图  1  判别函数最小化法流程图

      Figure 1.  Flowchart of Minimum Discriminate Function Method

      在判别函数最小化法中,相对权比λ从靠近零处开始,以固定步长d,遍历整个取值区间0 < λ < 1,每一个确定的λ可以得到相对应的观测向量的改正数$ {{\boldsymbol{\hat \varepsilon }}}$和系数矩阵中随机元素的改正数${{{\boldsymbol{\hat \varepsilon }}}_a}$,以及未知参数估值${\mathit{\boldsymbol{\hat \beta }}} $,通过$ {\mathit{\boldsymbol{\hat \varepsilon }}}$和${{{\mathit{\boldsymbol{\hat \varepsilon }}}_a}} $可以计算相应的判别函数${\overline {\mathit{\Phi}} ^{\left( j \right)}}\left( {{\boldsymbol{\hat \varepsilon }},{{{\boldsymbol{\hat \varepsilon }}}_\boldsymbol{a}}} \right) $的大小,在算法迭代过程中,λ的取值为d、2d、…、1-d,相应的未知参数估值和判别函数为${{\mathit{\boldsymbol{\hat \beta }}}^{\left( 1 \right)}} $、${{\mathit{\boldsymbol{\hat \beta }}}^{\left( 2 \right)}} $、…、$ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( {\frac{1}{d}-1} \right)}}$和$ {{\overline {\mathit{\Phi}} }^{\left( 1 \right)}}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right)$,$ {{\overline {\mathit{\Phi}} }^{\left( 2 \right)}}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right)$, …,${{\overline {\mathit{\Phi}} }^{\left( {\frac{1}{d}-1} \right)}}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) $。本文依据判别函数$ {{\overline {\mathit{\Phi}} }^{\left( j \right)}}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right)$的大小来搜索合适的相对权比λ的取值,把其中$ {{\overline {\mathit{\Phi}} }^{\left( j \right)}}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right)$的最小值对应的参数向量和此时的相对权比作为最终结果。

      ${{\overline {\mathit{\Phi}} }^{\left( j \right)}}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) $是$ {\mathit{\boldsymbol{\hat \varepsilon }}}$和${{{\mathit{\boldsymbol{\hat \varepsilon }}}_a}} $的函数用来确定相对权比,常用$\overline {\mathit{\Phi}} \left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) $函数主要有$ {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}\mathit{\boldsymbol{\hat \varepsilon }} + \mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}, {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon + \hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}、\lambda {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon }} + \left( {1-\lambda } \right)\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}、\sum\limits_{i = 1}^n {\left| {{{\mathit{\boldsymbol{\hat \varepsilon }}}_i}} \right|} + \sum\limits_{i = 1}^t {\left| {{{\mathit{\boldsymbol{\hat \varepsilon }}}_{{\mathit{\boldsymbol{a}}_i}}}} \right|} $等[15]

      (1) 判别函数$\overline {\mathit{\Phi}} _1^{(j)}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) = {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\left( j \right){\rm{T}}}}{{\mathit{\boldsymbol{\hat \varepsilon }}}^{\left( j \right)}} + \mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\left( j \right)T}\mathit{\boldsymbol{\hat \varepsilon }}_a^{\left( j \right)} $以观测数据偏离其估计值的距离平方和的大小为判别标准,等价于标准差这一通用评判标准,更具合理性。

      (2) 判别函数$\overline {\mathit{\Phi}} _2^{(j)}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) = {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\left( j \right){\rm{T}}}}{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat \varepsilon }}}^{\left( j \right)}} + \mathit{\boldsymbol{\hat \varepsilon }}_a^{\left( j \right)T}{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\left( j \right)} $与$\overline {\mathit{\Phi}} _1^{(j)}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) $相比增加了观测向量和系数矩阵中随机元素相对应的初始权阵,由于在第二步参数估计过程中已通过相对权比对观测向量和系数矩阵随机元素的改正数进行了合理分配,因此在判别函数中加入权阵则不太合理。

      (3) 判别函数$\overline {\mathit{\Phi}} _3^{(j)}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = \lambda {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\left( j \right){\rm{T}}}}{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat \varepsilon }}}^{\left( j \right)}} + \left( {1-\lambda } \right)\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\left( j \right){\rm{T}}}{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\left( j \right)} $与自由度的比值为单位权方差,而单位权方差的大小并不能反映平差结果的优劣,因此采用该判别函数并不能有效确定相对权比。

      (4) 判别函数$\overline {\mathit{\Phi}} _4^{(j)}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = \sum\limits_{i = 1}^n {\left| {\mathit{\boldsymbol{\hat \varepsilon }}_i^j} \right|} + \sum\limits_{i = 1}^t {\left| {\mathit{\boldsymbol{\hat \varepsilon }}_{{\mathit{\boldsymbol{a}}_i}}^j} \right|} $为观测数据偏离其估计值距离的绝对值之和,当观测数据中含有粗差时,是一种较为合适的确定相对权比的判别函数,当观测数据中不含粗差时,选取该判别函数效果不明显[15]

    • 模拟一个普通直线拟合问题,有y=β1x+β2。待估的未知参数有2个,分别为直线的斜率β1和截距β2,真值为β1=2,β2=1.5。自变量x从2到10之间每隔0.8模拟一个观测点, 并计算相对应的因变量y的值。在观测值xy的坐标上分别加上均值为0,方差为σx2=σ022·Px-1σy2=σ012·Py-1的随机误差,其中σ012=1,σ022=1.5,相应的权Px=diag(10.0, 3.0, 6.0, 1.0, 1.0, 8.0, 4.0, 3.0, 6.0, 1.0)、Py=diag(1.0, 2.0, 3.0, 1.0, 5.0, 4.0, 2.0, 7.0, 2.0, 10.0),其中diag(p1, p2,…, pn)表示的是对角线元素为p1, p2, …, pn的对角矩阵。做100次随机实验,分别采用表 1中7种方案进行参数求解,并比较计算所得的估计值$ {\mathit{\boldsymbol{\hat \beta }}}$与真值${\mathit{\boldsymbol{\tilde \beta }}} $的差值范数$ \left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\|$的大小,计算结果如表 2所示,表中${\mathit{\boldsymbol{\hat \beta }}} $、λ、$ \left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\|$都取100次结果的平均值。图 2表示方案4、方案5、方案6、方案7的判别函数与相对权比之间的关系图。本文将迭代的阈值δ0设为10-8

      表 1  各方案所用方法

      Table 1.  Method of Each Project

      方案 方法
      1 最小二乘方法(即λ=1)
      2 加权总体最小二乘方法(即λ=0.5)
      3 验前单位权方差法λ=σ022/(σ012+σ022),即λ=0.6000
      4 ${{\overline {\mathit{\Phi}} }_1}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) = {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}\mathit{\boldsymbol{\hat \varepsilon }} + \mathit{\boldsymbol{\hat \varepsilon }}_a^{\rm{T}}{{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}} $为判别函数
      5 $ {{\overline {\mathit{\Phi}} }_2}\left( {\mathit{\boldsymbol{\hat \varepsilon }}, {{\mathit{\boldsymbol{\hat \varepsilon }}}_\mathit{\boldsymbol{a}}}} \right) = {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon }} + \mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_a}$为判别函数
      6 ${{\overline {\mathit{\Phi}} }_3}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = \lambda {{\mathit{\boldsymbol{\hat \varepsilon }}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat \varepsilon }} + \left( {1-\lambda } \right)\mathit{\boldsymbol{\hat \varepsilon }}_\mathit{\boldsymbol{a}}^{\rm{T}}{\mathit{\boldsymbol{P}}_2}{{\mathit{\boldsymbol{\hat \varepsilon }}}_a} $为判别函数
      7 $ {{\overline {\mathit{\Phi}} }_4}\left( {{\bf{ \pmb{\mathsf{\hat ε}} }}, {{{\bf{ \pmb{\mathsf{\hat ε}} }}}_\mathit{\boldsymbol{a}}}} \right) = \sum\limits_{i = 1}^n {\left| {{{\mathit{\boldsymbol{\hat \varepsilon }}}_i}} \right|} + \sum\limits_{i = 1}^t {\left| {{{\mathit{\boldsymbol{\hat \varepsilon }}}_{{\mathit{\boldsymbol{a}}_i}}}} \right|} $为判别函数

      表 2  模型参数的计算结果

      Table 2.  Results of Model Parameter

      方案 $ {{\hat \beta }_1}$ ${{\hat \beta }_2} $ λ $\left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\| $
      真值 2.000 0 1.500 0 / /
      1 1.864 9 2.597 1 1.000 0 1.784 4
      2 2.006 8 1.523 1 0.500 0 0.980 9
      3 2.016 9 1.465 3 0.600 0 0.965 1
      4 2.018 8 1.455 0 0.620 0 0.963 2
      5 2.006 8 1.523 1 0.500 0 0.980 9
      6 2.043 9 1.323 1 0.990 0 1.003 1
      7 2.043 9 1.323 1 0.990 0 1.003 1

      图  2  直线拟合算例各方案判别函数与相对权比之间的关系图

      Figure 2.  Relationship Between Discriminate Function and Weight Scaling Factor in Straight Line Fitting

      假设直线方程为:

      $$ \mathit{\boldsymbol{y = }}{\beta _1}\mathit{\boldsymbol{x + }}{\beta _2} $$ (17)

      直线拟合的PEIV模型为式(2)。

      式中,$\underset{10\times 1}{\mathop{\mathit{\boldsymbol{y}}}}\, =\left[\begin{matrix} {{y}_{1}} \\ {{y}_{2}} \\ \vdots \\ {{y}_{9}} \\ {{y}_{10}} \\ \end{matrix} \right];\underset{10\times 1}{\mathop{\mathit{\boldsymbol{a}}}}\, =\left[\begin{matrix} {{x}_{1}} \\ {{x}_{2}} \\ \vdots \\ {{x}_{9}} \\ {{x}_{10}} \\ \end{matrix} \right] $;固定矩阵$ \mathit{\boldsymbol{B}}=\left[\begin{matrix} 1 \\ 0 \\ \end{matrix} \right]\otimes {{\mathit{\boldsymbol{I}}}_{10}};\mathit{\boldsymbol{h}}=\underbrace{\left[0, 0, \cdots, 0, \right.}_{10}{{\underbrace{\left. 1, 1, \cdots 1 \right]}_{10}}^{\rm{T}}};\underset{10\times 1}{{\boldsymbol{\hat \beta }}}\, =\left[\begin{matrix} {{{\hat{\beta }}}_{1}} \\ {{{\hat{\beta }}}_{2}} \\ \end{matrix} \right];\underset{10\times 1}{\mathop\ \boldsymbol{\varepsilon }}\, =\left[\begin{matrix} {{\varepsilon }_{{{y}_{1}}}} \\ {{\varepsilon }_{{{y}_{1}}}} \\ \vdots \\ {{\varepsilon }_{{{y}_{9}}}} \\ {{\varepsilon }_{{{y}_{10}}}} \\ \end{matrix} \right];\underset{10\times 1}{\mathop{{\boldsymbol{\varepsilon }_{\mathit{\boldsymbol{a}}}}}}\, =\left[\begin{matrix} {{\varepsilon }_{{{x}_{1}}}} \\ {{\varepsilon }_{{{x}_{1}}}} \\ \vdots \\ {{\varepsilon }_{{{x}_{9}}}} \\ {{\varepsilon }_{{{x}_{10}}}} \\ \end{matrix} \right]$。

    • 六参数坐标转换模型可写成[16]

      $$ \begin{gathered} \left[{\begin{array}{*{20}{c}} {{X_t}} \\ {{Y_t}} \end{array}} \right] \approx \left[{\begin{array}{*{20}{c}} {\cos \alpha }&{-\sin \alpha } \\ {\sin \alpha }&{\cos \alpha } \end{array}} \right]\left[{\begin{array}{*{20}{c}} 1&0 \\ \delta &1 \end{array}} \right] \cdot \left[{\begin{array}{*{20}{c}} {{s_X}{X_s}} \\ {{s_Y}{Y_s}} \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {{c_1}} \\ {{c_2}} \end{array}} \right] \hfill \\ = \left[{\begin{array}{*{20}{c}} {{a_1}}&{{b_1}} \\ {{a_2}}&{{b_2}} \end{array}} \right] \cdot \left[{\begin{array}{*{20}{c}} {{X_s}} \\ {{Y_s}} \end{array}} \right] + \left[{\begin{array}{*{20}{c}} {{c_1}} \\ {{c_2}} \end{array}} \right] \hfill \\ \end{gathered} $$

      当控制点个数为k(k≥3) 时,六参数转换模型可以表示为:

      $$ \left[{\begin{array}{*{20}{c}} {{X_{t1}}} \\ {{Y_{t1}}} \\ \cdots \\ {{X_{tk}}} \\ {{Y_{tk}}} \end{array}} \right]{\rm{ = }}\left[{\begin{array}{*{20}{c}} {{X_{s1}}}&{{Y_{s1}}}&1&0&0&0 \\ 0&0&0&{{X_{s1}}}&{{Y_{s1}}}&1 \\ \cdots &{}&{}&{}&{}& \cdots \\ {{X_{sk}}}&{{Y_{sk}}}&1&0&0&0 \\ 0&0&0&{{X_{sk}}}&{{Y_{sk}}}&1 \end{array}} \right] \cdot \left[{\begin{array}{*{20}{c}} {{a_1}} \\ {{b_1}} \\ {{c_1}} \\ {{a_2}} \\ {{b_2}} \\ {{c_2}} \end{array}} \right] $$ (18)

      式中,(Xs, Ys)、(Xt, Yt)分别为原始坐标系和目标坐标系中控制点的坐标值;sXsY为沿X轴和Y轴的缩放比例,α为旋转参数;c1c2为沿X轴和Y轴的平移参数;δ为关于坐标系两个坐标轴旋转的尺度因子。

      在文献[16]算例的基础上,通过Matlab软件进行数值模拟实验。选取13个点的二维坐标,设转换参数a1=0.9,b1=-0.8,c1=1,a2=0.6,b2=0.7,c2=5,按六参数转换模型得到新坐标系下的坐标值(Xt, Yt),实验数据如表 3。在两套坐标上均加上均值为0,方差阵为DSDT的随机误差,其中,

      表 3  原始坐标系和目标坐标系中13个点的坐标值

      Table 3.  Simulated Points in the Start System and Target System for the Affine Transformation

      点号 1 2 3 4 5 6 7 8 9 10 11 12 13
      原始坐标
      Xs 1.0 1.0 1.0 1.5 2.0 2.5 3.0 3.0 3.0 2.5 2.0 1.5 1.0
      Ys 1.0 2.0 3.0 3.5 4.0 3.5 3.0 2.0 1.0 0.5 0.0 0.5 1.0
      目标坐标
      Xt 1.1 0.3 -0.5 -0.45 -0.4 0.45 1.3 2.1 2.9 2.85 2.8 1.95 1.1
      Yt 6.3 7.0 7.7 8.35 9.0 8.95 8.9 8.2 7.5 6.85 6.2 6.25 6.3

      DS=0.015·diag(1, 2, 3, 1, 5, 4, 2, 7, 2, 1, 8, 3, 6),

      DT=0.010·diag(1, 3, 6, 1, 1, 8, 4, 3, 6, 5, 4, 5, 2)

      做100次随机实验,分别采用表 1中所列方案进行参数求解,并比较计算所得的估计值${\mathit{\boldsymbol{\hat \beta }}} $与真值${\mathit{\boldsymbol{\tilde \beta }}} $的差值范数$ \left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\|$的大小(βa1b1c1a2b2c2组成的参数列向量),计算结果如表 4所示,同样表中${\mathit{\boldsymbol{\hat \beta }}} $、λ、$ \left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\|$都取100次结果的平均值。图 3表示方案4~方案7的判别函数与相对权比之间的关系图。

      表 4  模型参数的计算结果

      Table 4.  Results of Model Parameter

      方案 $ {{\hat a}_1}$ $ {{\hat b}_1}$ ${{\hat c}_1} $ $ {{\hat a}_2}$ $ {{\hat b}_2}$ $ {{\hat c}_2}$ λ $\left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\| $
      真值 0.900 0 -0.800 0 1.000 0 0.600 0 0.700 0 5.000 0 / /
      1 0.849 6 -0.776 1 1.046 3 0.553 6 0.697 3 5.093 9 1.000 0 0.341 7
      2 0.894 0 -0.793 6 1.002 1 0.591 8 0.700 5 5.016 1 0.500 0 0.298 3
      3 0.899 3 -0.795 7 0.996 4 0.597 8 0.701 9 5.002 1 0.600 0 0.297 2
      4 0.895 1 -0.794 0 1.000 9 0.593 0 0.700 7 5.013 3 0.520 0 0.297 9
      5 0.894 0 -0.793 6 1.002 1 0.591 8 0.700 5 5.016 1 0.500 0 0.298 3
      6 0.921 2 -0.802 7 0.965 4 0.627 2 0.715 2 4.920 5 0.990 0 0.317 0
      7 0.921 2 -0.802 7 0.965 4 0.627 2 0.715 2 4.920 5 0.990 0 0.317 0

      图  3  坐标转换算例各方案判别函数与相对权比之间的关系图

      Figure 3.  Relationship Between Discriminate Function and Weight Scaling Factor in Coordinate Transformation

    • (1) 从表 2表 4的解算结果可以看出,由于最小二乘方法(方案1) 没有顾及系数矩阵的误差,得到的结果最差;加权总体最小二乘方法(方案2) 考虑了系数矩阵的误差,以及观测值不等精度、相关的情况,因此方案2优于方案1,但在该方案中未将观测向量和系数矩阵权分配不合理、验前随机模型不准确的情况考虑在内,导致平差结果出现偏差,解算结果与方案3、方案4相比较差。

      (2) 直线拟合和坐标转换算例中验前单位权方差法(方案3) 效果都比较好,这与理论相符,因为给定的权比即为实际的权比。然而这种方法不具操作性,因为在数据处理之前,观测向量和系数矩阵的验前单位权方差是很难准确已知的。

      (3) 综合上面两个算例,以参数结果与真值的差值范数$ \left\| {\Delta \mathit{\boldsymbol{\hat \beta }}} \right\|$作为评判标准,判别函数最小化迭代算法中在确定相对权比时选取方案4中的判别函数,算法效果是最好的,且确定的相对权比与方案3中的最为接近,这与§3中对该判别函数进行的分析一致。方案6和方案7的解算结果相同,但从图 2图 3可以看出,两种方案并不等价,该两种方案确定的相对权比均为0.99,此时系数矩阵中随机元素对模型参数估计的贡献基本趋于零,显然不太合理。上面两个算例数据中均不含粗差,因此采用方案7的效果不明显,这与§3对该判别函数的分析一致。方案5的解算结果与方案2相同,确定的相对权比为0.50,此时相对权比并没有发挥作用。

      (4) 本文所采用的两个算例,模型的系数矩阵均为部分元素为含有随机误差的观测数据,另一部分为非随机的固定值;针对此类系数矩阵,在进行模型参数估计时,一般需要采用结构总体最小二乘方法,或选择特定的定权方式,将系数矩阵中的非随机元素固定,这些方法在处理过程中都比较复杂。本文算法是以PEIV模型为基础,可以很好地解决这一问题,与上述两类方法相比更加简单快捷。

    • 本文针对观测向量和系数矩阵权分配不合理、验前随机模型不准确的情况,提出了附有相对权比的PEIV模型的总体最小二乘平差算法,依托PEIV模型的优势, 该算法普遍适用于一般性的系数矩阵和权矩阵;通过在平差准则中加入相对权比,自适应调整观测向量和系数矩阵随机元素对模型参数估计的贡献。通过算例表明,当观测向量和系数矩阵中随机元素的验前单位权方差信息已知时,可采用验前单位权方差法进行参数估计,效果较好;当观测向量和系数矩阵中随机元素的验前单位权方差信息未知时,可以采用判别函数最小化迭代算法进行参数求解,在确定相对权比时选取方案4中的判别函数,算法效果是最好的。相对权比λ的确定是本文算法的关键,如何从理论方面给出严密的确定相对权比的公式是下一步的研究方向。

参考文献 (16)

目录

    /

    返回文章
    返回