留言板

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

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

带有界不确定性的加权混合估计方法

宋迎春 宋采薇 左廷英

宋迎春, 宋采薇, 左廷英. 带有界不确定性的加权混合估计方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
引用本文: 宋迎春, 宋采薇, 左廷英. 带有界不确定性的加权混合估计方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
SONG Yingchun, SONG Caiwei, ZUO Tingying. A Weighted Mixed Estimation Method with Bounded Uncertainty[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
Citation: SONG Yingchun, SONG Caiwei, ZUO Tingying. A Weighted Mixed Estimation Method with Bounded Uncertainty[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300

带有界不确定性的加权混合估计方法

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

国家自然科学基金 41674009

国家自然科学基金 41574006

国家自然科学基金 41674012

详细信息
    作者简介:

    宋迎春,博士,教授,主要研究方向为测量数据处理的理论与方法。csusyc@csu.edu.cn

    通讯作者: 左廷英,博士,教授。dengmin208@tom.com
  • 中图分类号: P207

A Weighted Mixed Estimation Method with Bounded Uncertainty

Funds: 

The National Natural Science Foundation of China 41674009

The National Natural Science Foundation of China 41574006

The National Natural Science Foundation of China 41674012

More Information
    Author Bio:

    SONG Yingchun, PhD, professor, specializes in theory and method of measuring data processing. E-mail:csusyc@csu.edu.cn

    Corresponding author: ZUO Tingying, PhD, professor. E-mail: dengmin208@tom.com
图(1) / 表(3)
计量
  • 文章访问数:  1001
  • HTML全文浏览量:  278
  • PDF下载量:  63
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-05-23
  • 刊出日期:  2020-07-30

带有界不确定性的加权混合估计方法

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

    国家自然科学基金 41674009

    国家自然科学基金 41574006

    国家自然科学基金 41674012

    作者简介:

    宋迎春,博士,教授,主要研究方向为测量数据处理的理论与方法。csusyc@csu.edu.cn

    通讯作者: 左廷英,博士,教授。dengmin208@tom.com
  • 中图分类号: P207

摘要: 大地测量中各种异方差多源观测模型进行融合都需要进行混合估计。由于附加信息和样本信息在估计过程中作用是不均等的,需要建立新的加权平差准则,平衡先验约束和观测信息对参数估计的影响。把多源观测数据看成是观测信息和一些随机约束信息,首先利用椭球近似描述有界不确定信息,建立了基于外接椭球特征矩阵迹最小的平差准则,然后提出了一个新的观测信息融合方法,并给出了一种优化权的计算方法,使得加权混合估计方法能有效应用于大地测量数据处理,最后,通过算例验证了算法的有效性,说明了集员估计解与带权混合估计的关系。

English Abstract

宋迎春, 宋采薇, 左廷英. 带有界不确定性的加权混合估计方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
引用本文: 宋迎春, 宋采薇, 左廷英. 带有界不确定性的加权混合估计方法[J]. 武汉大学学报 ● 信息科学版, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
SONG Yingchun, SONG Caiwei, ZUO Tingying. A Weighted Mixed Estimation Method with Bounded Uncertainty[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
Citation: SONG Yingchun, SONG Caiwei, ZUO Tingying. A Weighted Mixed Estimation Method with Bounded Uncertainty[J]. Geomatics and Information Science of Wuhan University, 2020, 45(7): 949-955. doi: 10.13203/j.whugis20180300
  • 随着测量理论与技术的发展,测量平差的对象已从过去的单一同类观测扩展为同类不同精度或异构多源观测,通过多源数据融合可以获取更全面、更有效的观测信息,它比单一信息源更精确、更可靠[1-3]。然而,多源观测有不同的函数模型和随机模型,多源数据融合中常出现一些复杂的先验信息,如参数间某些随机先验信息可以构成随机模型约束[4],又如噪声在某一范围内变化,参数值有一个已知的取值范围,它们可以构成参数的可行域利用集员估计方法进行计算[5-7]。Fogel[8]首先提出基于有界噪声的参数集员估计方法,其集合的Chebyshev中心可作为参数真实值的一个点估计。Fogel和Huang[9]又进一步改善了其算法,将椭球引入到参数可行集的近似描述中来,提出了椭球集员估计算法。近年来,椭球集员估计算法得到了迅速的发展[10-12],利用这些算法可以更好地将多源观测中的各类先验信息建模,获取最优估计值。大地测量中,许多学者提出了联合平差方法处理多源观测信息,可是,此类方法大多需要知道各类观测量噪声的准确先验方差,从而确定各类融合数据的权[13-15]。在无法准确知道各类观测量的先验方差时,有学者提出了迭代求权值方法[14-15, 26],但是迭代重加权最小二乘法是参数的一个非线性估计,其协方差阵计算困难,无法对其进行精度评估。在数学上,许多学者把多源观测数据看成是在新观测得到的样本信息上加上一些先验的随机约束信息[16-18],Durbin[19]、Theil和Goldberger[20]、Theil[21]早期研究了这类问题,提出了混合估计。由于附加信息和样本信息在估计过程中作用是不均等的,Schaffrin和Toutenburg[22]在混合估计的基础上,提出了加权混合估计。这些算法大多注重于算法的效率和估计的优良性,不能直接用于大地测量数据处理。本文针对大地测量中各种异构多源观测信息,利用椭球近似描述有界噪声,提出了椭球集员估计算法。以两个椭球交集的外接椭球特征矩阵的迹最小的平差准则,建立了一个新的观测信息融合方法,利用最优化权平衡先验约束和观测信息对参数估计的影响,使得加权混合估计方法能有效应用于大测量数据处理。并通过算例验证了算法的有效性,说明了集员估计解与带权混合估计的关系。

    • e为有界噪声,其有界性可以用范数形式来描述,如$\left\| \mathrm{e} \right\|\le \gamma $,也可用区间来描述,如e ∈ [eleu],这些描述方法都可以统一用椭球集合来表示:

      $$ \mathit{\boldsymbol{E}}\left( \mathit{\boldsymbol{e}} \right) = \left\{ {\mathit{\boldsymbol{e}}:{\rm{}}{\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{e}} \le 1} \right\} $$ (1)

      式中,P为正定矩阵,它是椭球的特征矩阵,用来刻画椭球的形状特征,类似于e的协方差阵描述e的特征。有许多学者研究了矩阵P的构造[10, 12, 23],在几何上,椭球的扁平程度以及椭球的体积是由椭球的特征矩阵来确定的,由于本文研究的不确定性是一种有界约束(椭球约束),这个有界性是通过矩阵P来刻画的,它相当于e的协方差阵来刻画e的特征一样[25-26]。对于平差模型:

      $$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}} $$ (2)

      式中,X = (x1x2x m) Tm维未知向量;Ln维观测向量;An × m的设计矩阵,秩(A) = men维有界随机误差向量,满足$\left\| \mathit{\boldsymbol{e}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 = {\mathit{\boldsymbol{e}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{e}} < 1$,其中P e-1为正定对角矩阵。显然,

      $$ \begin{array}{*{20}{c}} {\left\| \mathit{\boldsymbol{e}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 = \left\| {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 = }\\ {{{(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}})}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}) < 1} \end{array} $$ (3)

      在模型式(2)中,取X L = (AT A)-1 AT L,由AX - L = (AX - P A L)+(P A - I) L,有:

      $$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 = }\\ {\left\| {\mathit{\boldsymbol{AX}} - {\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 + \left\| {({\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}} - \mathit{\boldsymbol{I}})\mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 + }\\ {2{{[({\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}} - \mathit{\boldsymbol{I}})\mathit{\boldsymbol{L}}]}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}(\mathit{\boldsymbol{AX}} - {\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{L}})} \end{array} $$ (4)

      式中,P A = A(AT A)-1 AT。又因P e-1为对角矩阵,可使得AP e-1 = P e-1 A,而(P A - I) A = 0,故

      $$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{L}}^{\rm{T}}}({\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}} - \mathit{\boldsymbol{I}})\mathit{\boldsymbol{P}}_e^{ - 1}(\mathit{\boldsymbol{AX}} - {\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{L}}) = \\ {\mathit{\boldsymbol{L}}^{\rm{T}}}({\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}} - \mathit{\boldsymbol{I}})\mathit{\boldsymbol{AP}}_e^{ - 1}(X - {({\mathit{\boldsymbol{A}}^\mathit{\boldsymbol{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^\mathit{\boldsymbol{T}}}\mathit{\boldsymbol{L}}) = 0 \end{array} $$

      由式(4)有:

      $$ \begin{array}{*{20}{c}} {\left\| {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 = \left\| {\mathit{\boldsymbol{AX}} - {\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 + }\\ {\left\| {({\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}} - \mathit{\boldsymbol{I}})\mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2} \end{array} $$ (5)

      故有:

      $$ \begin{array}{*{20}{c}} {{{(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_L})}^{\rm{T}}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}}(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_L}) = }\\ {\left\| {\mathit{\boldsymbol{AX}} - {\mathit{\boldsymbol{P}}_\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 < \left\| {\mathit{\boldsymbol{AX}} - \mathit{\boldsymbol{L}}} \right\|_{\mathit{\boldsymbol{P}}_e^{ - 1}}^2 < 1} \end{array} $$ (6)

      这说明对任何的XE (e)= { X:(L - AX) T P e-1(L - AX) < 1 },都会有:

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{X}} \in \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{X}}_L}, {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}}) = }\\ {\{ \mathit{\boldsymbol{X}}:{{(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_L})}^{\rm{T}}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}}(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_L}) < 1\} } \end{array} $$

      $$ \mathit{\boldsymbol{E}}(\mathit{\boldsymbol{e}}) \subseteq \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{X}}_L}, {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}}) $$ (7)

      即椭球E (X LAT P e-1 A)包含了椭球E (e)。此结论也可以由下面的推理得到:若L = AX是相容方程组,取X 0使得L = AX 0,当L = AX是不相容时,取X 0 = X LS = (AT A)-1 AT L,这时,LAX 0,利用式(1)有:

      $$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{e}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{e}} = {{(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}})}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}) = }\\ {{{(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}})}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{AX}}) \approx }\\ {{{(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_0})}^{\rm{T}}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}^{ - 1}}\mathit{\boldsymbol{A}}(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_0})} \end{array} $$

      e的有界不确定性也可以近似地表示为:E (e)= { X:(X - X 0) T AT P -1 A(X - X 0) ≤ 1 }

      为了后面研究的方便,本文使用不确定性椭球E (X LAT P e-1 A)代替有界误差椭球集合E (e)。式(2)中的X是没有任何约束的,但在一些实际问题中,往往要求X满足某种约束条件,如X是非负的,或X满足某个线性等式约束。Theil和Goldberger[20]提出了混合回归模型,它是将线性约束HX = c随机化得到的随机化约束:

      $$ \mathit{\boldsymbol{h}} = \mathit{\boldsymbol{HX}} + \mathit{\boldsymbol{w}} $$ (8)

      式中,hp维随机向量;Hp × m矩阵;wp维有界随机误差向量,满足w T P w-1w < 1。这样就形成了带有随机约束的有界误差混合平差模型:

      $$ \left\{ {\begin{array}{*{20}{l}} {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{e}}}\\ {\mathit{\boldsymbol{ }}{\rm{s}}{\rm{.t}}{\rm{.}}\mathit{\boldsymbol{ h}} = \mathit{\boldsymbol{HX}} + \mathit{\boldsymbol{w}}} \end{array}} \right. $$ (9)

      文献[16-24]给出了许多的混合估计算法,相关的算法已在实例2中列出,这些算法大多注重于算法的效率和估计的优良性。类似前面的推导,对于模型式(4),同样可以定义如下:

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{E}}({\mathit{\boldsymbol{X}}_\mathit{\boldsymbol{h}}}, {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}}) = }\\ {\{ \mathit{\boldsymbol{X}}:{{(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_\mathit{\boldsymbol{h}}})}^T}({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}})(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{X}}_\mathit{\boldsymbol{h}}}) \le 1\} } \end{array} $$ (10)

      其中,X h = (H T P w-1 H)-1 H T P w-1 h。下面将利用有界不确定性来推导带权混合估计中的权值的一个新确定方法。

    • E (z 1P 1)和E (z 2P 2)为两个椭球,它们分别定义为:

      $$ \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_i}, {\mathit{\boldsymbol{P}}_i}) = \{ \mathit{\boldsymbol{X}}:{(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{z}}_i})^{\rm{T}}}\mathit{\boldsymbol{P}}_i^{ - 1}(\mathit{\boldsymbol{X}} - {\mathit{\boldsymbol{z}}_i}) \le 1] $$

      其中,i = 1,2。它们的交定义为:

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_1}, {\mathit{\boldsymbol{P}}_1}) \cap \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_2}, {\mathit{\boldsymbol{P}}_2}) = }\\ {\{ \mathit{\boldsymbol{X}}:\mathit{\boldsymbol{X}} \in \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_1}, {\mathit{\boldsymbol{P}}_1}), \mathit{\boldsymbol{X}} \in \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_2}, {\mathit{\boldsymbol{P}}_2})\} } \end{array} $$ (11)

      显然,两个椭球的交集不一定是一个椭球,它可以用一个外包椭球来近似[11, 23],如图 1所示,设其外包椭球族为E (zP),有:

      图  1  两个椭球交的最小外包椭球

      Figure 1.  Minimum Circumscribed Ellipsoid with Two Ellipsoid Intersections

      $$ {\mathit{\boldsymbol{E}}(\mathit{\boldsymbol{z}}, \mathit{\boldsymbol{P}}) \supseteq \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_1}, {\mathit{\boldsymbol{P}}_1}) \cap \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{z}}_2}, {\mathit{\boldsymbol{P}}_2})} $$ (12)
      $$ {{\mathit{\boldsymbol{P}}^{ - 1}} = \frac{1}{{1 - \rho }}[a\mathit{\boldsymbol{P}}_1^{ - 1} + (1 - a)\mathit{\boldsymbol{P}}_2^{ - 1}]} $$ (13)
      $$ {\mathit{\boldsymbol{z}} = [a\mathit{\boldsymbol{P}}_1^{ - 1} + (1 - a)\mathit{\boldsymbol{P}}_2^{ - 1}]{\mathit{\boldsymbol{z}}_0}} $$ (14)
      $$ {{\mathit{\boldsymbol{z}}_0} = [a\mathit{\boldsymbol{P}}_1^{ - 1}{\mathit{\boldsymbol{z}}_1} + (1 - a)\mathit{\boldsymbol{P}}_2^{ - 1}{\mathit{\boldsymbol{z}}_2}]} $$ (15)
      $$ \begin{array}{*{20}{c}} {\rho = a\mathit{\boldsymbol{z}}_1^{\rm{T}}\mathit{\boldsymbol{P}}_1^{ - 1}{\mathit{\boldsymbol{z}}_1} + (1 - a)\mathit{\boldsymbol{z}}_2^{\rm{T}}\mathit{\boldsymbol{P}}_2^{ - 1}{\mathit{\boldsymbol{z}}_2} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{z}}_0^{\rm{T}}[a\mathit{\boldsymbol{P}}_1^{ - 1} + (1 - a)\mathit{\boldsymbol{P}}_2^{ - 1}]{\mathit{\boldsymbol{z}}_0}} \end{array} $$ (16)

      其中,0 < a < 1。上面的交运算式(12)至式(16)的推导可以参考文献[23],通过优化系数a,可以得到最小迹或最小体积外包椭球[11, 23],见图 1

      现假设式(9)的可行解集$\mathit{\boldsymbol{E}} = \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{X}}_L}, {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}}) \cap \mathit{\boldsymbol{E}}({\mathit{\boldsymbol{X}}_\mathit{\boldsymbol{h}}}, {\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}})$的外包椭球$\mathit{\boldsymbol{E}}({{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}, {\mathit{\boldsymbol{P}}_{{\rm{ME}}}})$为:

      $$ \{ \mathit{\boldsymbol{X}}:{(\mathit{\boldsymbol{X}} - {{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}})^{\rm{T}}}\mathit{\boldsymbol{P}}_{{\rm{ME}}}^{ - 1}(\mathit{\boldsymbol{X}} - {{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}) \le 1] $$ (17)

      由式(12)和式(13)有:

      $$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{P}}_{{\rm{ME}}}^{ - 1} = \frac{1}{{1 - \rho }}[a{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}} + (1 - a){\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}}]}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (0 < a < 1)} \end{array} $$ (18)

      Q = aAT P e-1 A + (1 - a) H T P w-1 H,由文献[23]给出的证明可知0 ≤ ρ < 1,则有:

      $$ \begin{array}{l} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}} = {\mathit{\boldsymbol{Q}}^{ - 1}}[a{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_L} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (1 - a){\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}}{\mathit{\boldsymbol{X}}_h}] = \\ {\mathit{\boldsymbol{Q}}^{ - 1}}[a{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}L + (1 - a){\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{h}}] \end{array} $$ (19)

      由文献[16]可知,${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$是式(19)的一个解,从式(19)可以看出,${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$是岭估计算法的一种推广,形式上就是Schaffrin和Toutenburg [22]提出的加权混合估计。如果能确定a,不仅可以得到带有椭球不确定性的参数估计方法,还可以从不确定性的角度得到一个加权混合估计的权确定方法。

    • 通过一个简单计算式(18)可以化为:

      $$ {\mathit{\boldsymbol{P}}_{{\rm{ME}}}} = (1 - \rho )[a{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}} + (1 - a){\mathit{\boldsymbol{H}}^T}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}}] $$ (20)

      由式(16)有:

      $$ \begin{array}{l} \rho = a{\mathit{\boldsymbol{L}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{L}} + (1 - a){\mathit{\boldsymbol{h}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{h}} - [a{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{L}} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {(1 - a){\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{h}}{]^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{ - 1}}[a{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{L}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} (1 - a){\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{h}}]} \end{array} \end{array} $$ (21)

      不同的a可以得到不同的外包椭球,为了保证椭球交的外包椭球的最小性,可以通过优化系数a,得到最小迹外包椭球:

      $$ a = {\rm{arg}}{\kern 1pt} {\kern 1pt} {\rm{min}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{tr}}{\kern 1pt} {\kern 1pt} {\kern 1pt} ({\mathit{\boldsymbol{P}}_{{\rm{ME}}}}) $$ (22)

      许多的文献给出了式(21)中a的计算方法,但通常较为复杂[10-12]。本文的方法是直接搜索,因为0 < a < 1,0 ≤ ρ < 1,可知a必须使得0 ≤ ρ < 1成立的条件下,mintr (P ME)取最小值。直接利用搜索算法(a从0开始到1止,通过增量Δa,逐步搜索得到使tr (P ME)达到最小的a)可以求出a的值,从而求出参数估计值${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$和P ME

    • 算例1:为了便于画出椭圆进行分析说明,设计如式(9)所示二维的平差模型,其中:$\mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {3.485{\kern 1pt} {\kern 1pt} {\kern 1pt} 4}&{0.964{\kern 1pt} {\kern 1pt} 0}\\ {8.542{\kern 1pt} {\kern 1pt} {\kern 1pt} 1}&{0.799{\kern 1pt} {\kern 1pt} 0}\\ { - 3.123{\kern 1pt} {\kern 1pt} {\kern 1pt} 7}&{0.385{\kern 1pt} {\kern 1pt} {\kern 1pt} 5} \end{array}} \right]$;$\mathit{\boldsymbol{L}} = {[21.312{\kern 1pt} {\kern 1pt} 9, 40.165{\kern 1pt} {\kern 1pt} {\kern 1pt} 9, - 8.903{\kern 1pt} {\kern 1pt} {\kern 1pt} 2]^{\rm{T}}}$;$\mathit{\boldsymbol{h}} = {\left[ {\begin{array}{*{20}{c}} { - 10.219{\kern 1pt} {\kern 1pt} 6}&{ - 27.988{\kern 1pt} {\kern 1pt} {\kern 1pt} 5} \end{array}} \right]^{\rm{T}}}$,参数的真值为${\mathit{\boldsymbol{X}}_{{\rm{ true }}}} = {\left[ {\begin{array}{*{20}{c}} 4&7 \end{array}} \right]^{\rm{T}}}$;误差向量为:$\mathit{\boldsymbol{e}} = \mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{X}}_{{\rm{ true }}}} = {\left[ {\begin{array}{*{20}{l}} {0.622{\kern 1pt} {\kern 1pt} 9}&{0.404{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}&{0.892{\kern 1pt} {\kern 1pt} 9} \end{array}} \right]^{\rm{T}}}$,$\mathit{\boldsymbol{w}} = \mathit{\boldsymbol{h}} - \mathit{\boldsymbol{H}}{\mathit{\boldsymbol{X}}_{{\rm{ true }}}} = {\left[ {\begin{array}{*{20}{c}} {0.236{\kern 1pt} {\kern 1pt} {\kern 1pt} 2}&{ - 0.524{\kern 1pt} {\kern 1pt} {\kern 1pt} 1} \end{array}} \right]^{\rm{T}}}$, 满足eT P e- 1 e < 1和w T P w-1 w < 1,此处有:

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{P}}_e^{ - 1} = {\rm{diag}} (0.5998, 1.7121, 0.2241)}\\ {\mathit{\boldsymbol{P}}_w^{ - 1} = {\rm{diag}} (2.2317, 2.4205)} \end{array} $$

      利用式(19)至式(21)进行计算得到权值$a = 0.529{\kern 1pt} {\kern 1pt} 0, {\mathit{\boldsymbol{\hat X}}_M} = {[4.0485\quad 7.2166]^{\rm{T}}}$,

      $$ {\mathit{\boldsymbol{P}}_M} = \left[ {\begin{array}{*{20}{l}} {0.004{\kern 1pt} {\kern 1pt} 3}&{0.005{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}\\ {0.005{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}&{0.314{\kern 1pt} {\kern 1pt} {\kern 1pt} 5} \end{array}} \right]。 $$

      图 1中,E (X LAT P e-1 A)为有界噪声e的误差椭球,E (X hH T P w-1 H)为有界噪声w的误差椭球,E (X MP M-1)为E (X LAT P e-1 A)和E (X hH T P w-1 H)交的最小外包椭球。

      算例2:设有二次曲面f (xy) = -2.735x + 1.543y+3.648x 2-2.741y 2-2.681xy,现对f (xy)模拟了两类观测数据,一类为带有协方差为σ e2 P e-1的观测噪声,另一类为带有协方差为σ w2 P w-1的观测噪声,见表 1表 2,权矩阵P eP w为对角矩阵,单位权方差σ e2 = 0.81,σ w2 = 1.00,P e-1= diag (1.616 8,1.988 6,1.562 1,1.703 6,1.321 2,0.544 2, 0.481 7), P w-1 = diag (2.437 3, 1.209 8, 1.980 2,0.858 2,2.884 2,1.729 7,2.899 9,1.541 8,2.592 7)。采用如下的6种方案进行计算,并分析各种算法的效率,设参数估计的偏差平方和$f(\mathit{\boldsymbol{\hat X}}) = {(\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{X}}_{{\rm{ true }}}})^{\rm{T}}}(\mathit{\boldsymbol{\hat X}} - {\mathit{\boldsymbol{X}}_{{\rm{ true }}}})$,其计算结果见表 3

      表 1  带有界噪声(eT P e-1e < 1)的观测向量和设计矩阵

      Table 1.  Observation Vector and Design Matrix with Bounded Noise(eT P e-1e < 1)

      L e A
      10.180 0 −0.222 6 2.070 7 −0.114 7 4.287 9 0.013 2 −0.237 6
      −3.594 0 −0.281 9 0.325 0 −0.897 1 0.105 6 0.804 8 −0.291 5
      280.066 9 0.390 2 9.545 2 0.993 3 91.111 2 0.986 6 9.481 2
      112.733 0 0.193 9 −5.056 0 0.388 8 25.563 5 0.151 1 −1.965 5
      30.451 3 0.236 4 3.363 9 0.230 0 11.315 8 0.052 9 0.773 8
      205.951 4 −0.127 8 −7.065 4 0.234 3 49.919 4 0.054 9 −1.655 5
      22.035 8 0.530 1 −1.897 6 0.659 5 3.600 7 0.434 9 −1.251 4

      表 2  带有界噪声(w T P w-1w < 1)的观测向量和设计矩阵

      Table 2.  Observation Vector and Design Matrix with Bounded Noise(w T P w-1w < 1)

      h w H
      13.553 4 −0.431 7 −1.481 5 0.449 8 2.194 9 0.202 3 −0.666 4
      9.244 0 −0.018 4 1.951 2 −0.233 1 3.807 0 0.054 3 −0.454 8
      16.727 8 0.117 9 −1.599 2 0.797 7 2.557 4 0.636 4 −1.275 7
      98.788 1 −0.019 4 5.350 5 −0.863 1 28.627 7 0.745 0 −4.618 1
      341.244 3 0.309 8 −9.486 2 −0.470 0 89.987 3 0.220 9 4.458 5
      41.485 8 0.155 4 3.855 9 0.249 1 14.867 7 0.062 1 0.960 6
      75.261 7 0.098 7 5.068 8 0.360 6 25.692 4 0.130 0 1.827 6
      264.017 4 0.055 0 9.071 2 0.473 8 82.285 9 0.224 5 4.297 9
      29.351 7 0.149 1 −2.253 6 0.866 0 5.078 5 0.749 9 −1.951 5

      表 3  几种不同方案的参数估计结果

      Table 3.  Parameter Estimation Results of Several Different Schemes

      ${\mathit{\boldsymbol{X}}_{{\rm{ true }}}}$ ${{\mathit{\boldsymbol{\hat X}}}_L}$ ${{\mathit{\boldsymbol{\hat X}}}_h}$ ${{\mathit{\boldsymbol{\hat X}}}_M}$ ${{\mathit{\boldsymbol{\hat X}}}_{{\rm{MS}}}}$ $\mathit{\boldsymbol{\hat X}}_{{\rm{MS}}}^{\left( k \right)}$ ${{\mathit{\boldsymbol{\hat X}}}_E}$
      −2.735 −2.736 1 −2.713 3 −2.742 7 −2.743 1 −2.742 6 −2.742 3
      1.543 2.233 1 0.621 5 1.682 3 1.702 1 1.677 4 1.652 3
      3.648 3.641 4 3.621 4 3.650 2 3.650 1 3.650 2 3.650 2
      −2.741 −2.342 9 −0.246 7 −2.784 4 −2.790 7 −2.782 2 −2.767 4
      −2.681 −2.687 8 −2.247 7 −2.666 0 −2.665 2 −2.666 2 −2.666 6
      $f(\mathit{\boldsymbol{\hat X}})$ 0.634 8 7.259 6 0.021 6 0.028 1 0.020 0 0.012 9

      方案1:采用观测数据表 1进行最小二乘计算。估计解为:

      $$ {{{\mathit{\boldsymbol{\hat X}}}_L} = {{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}})}^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{L}}} $$ (23a)
      $$ {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_L}}} = {{({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}})}^{ - 1}}} $$ (23b)

      方案2:采用观测数据表 2进行最小二乘计算。估计解为:

      $$ {{{\mathit{\boldsymbol{\hat X}}}_h} = {{({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}})}^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{h}}} $$ (24a)
      $$ {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_h}}} = {{({\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}})}^{ - 1}}} $$ (24b)

      方案3:采用混合加权算法(单位权方差σe2σw2已知)[16]。估计解为:

      $$ \begin{array}{l} {{\mathit{\boldsymbol{\hat X}}}_M} = {\left( {\frac{1}{{\sigma _e^2}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{A}} + \frac{1}{{\sigma _w^2}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{H}}} \right)^{ - 1}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {{\kern 1pt} \frac{1}{{\sigma _e^2}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}\mathit{\boldsymbol{L}} + \frac{1}{{\sigma _w^2}}{\mathit{\boldsymbol{H}}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}\mathit{\boldsymbol{h}}} \right) \end{array} $$ (25a)
      $$ {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_M}}} = {\left( {\frac{1}{{\sigma _\varDelta ^2}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_\varDelta }\mathit{\boldsymbol{A}} + \frac{1}{{\sigma _e^2}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_e}\mathit{\boldsymbol{H}}} \right)^{ - 1}} $$ (25b)

      方案4:采用混合加权算法(单位权方差σe2σw2未知)[16]。估计解为:

      $$ \begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\hat X}}}_{{\rm{MS}}}} = {{\left( {\frac{1}{{\mathit{\boldsymbol{s}}_1^2}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_\varDelta }\mathit{\boldsymbol{A}} + \frac{1}{{\mathit{\boldsymbol{s}}_2^2}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_e}\mathit{\boldsymbol{H}}} \right)}^{ - 1}} \times }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\frac{1}{{\mathit{\boldsymbol{s}}_1^2}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_\varDelta }\mathit{\boldsymbol{L}} + \frac{1}{{\mathit{\boldsymbol{s}}_2^2}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_e}\mathit{\boldsymbol{h}}} \right)} \end{array} $$ (26)

      其中,s12s2 2分别为σe2σw2的估计值,计算式分别为:

      $$ \begin{array}{l} \mathit{\boldsymbol{s}}_1^2 = \frac{{{{(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_{\rm{L}}})}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_{\rm{L}}})}}{{n - m}}\\ \mathit{\boldsymbol{s}}_2^2 = \frac{{{{(h - \mathit{\boldsymbol{H}}{{\mathit{\boldsymbol{\hat X}}}_h})}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}(\mathit{\boldsymbol{h}} - \mathit{\boldsymbol{H}}{{\mathit{\boldsymbol{\hat X}}}_h})}}{{p - m}} \end{array} $$

      方案5:采用迭代重加权算法,首先计算出σe2σw2的初始估计值:

      $$ \mathit{\boldsymbol{s}}_1^2\left( 0 \right) = \frac{{{{(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_L})}^{\rm{T}}}\mathit{\boldsymbol{P}}_e^{ - 1}(\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{\hat X}}}_L})}}{{n - m}} $$ (27)
      $$ \mathit{\boldsymbol{s}}_2^2\left( 0 \right) = \frac{{{{({\rm{h}} - \mathit{\boldsymbol{H}}{{\mathit{\boldsymbol{\hat X}}}_h})}^{\rm{T}}}\mathit{\boldsymbol{P}}_w^{ - 1}(\mathit{\boldsymbol{h}} - \mathit{\boldsymbol{H}}{{\mathit{\boldsymbol{\hat X}}}_h})}}{{p - m}} $$ (28)

      然后,利用下面的式子进行迭代计算:

      $$ \begin{array}{l} \mathit{\boldsymbol{\hat X}}_{{\rm{MS}}}^{(k + 1)} = {\left( {\frac{1}{{\mathit{\boldsymbol{s}}_1^2(k)}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_\varDelta }\mathit{\boldsymbol{A}} + \frac{1}{{\mathit{\boldsymbol{s}}_2^2(k)}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_e}\mathit{\boldsymbol{H}}} \right)^{ - 1}} \times \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left( {\frac{1}{{\mathit{\boldsymbol{s}}_1^2(k)}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_\varDelta }\mathit{\boldsymbol{L}} + \frac{1}{{\mathit{\boldsymbol{s}}_2^2(k)}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_e}\mathit{\boldsymbol{h}}} \right) \end{array} $$ (29)
      $$ {\mathit{\boldsymbol{s}}_1^2(k + 1) = \frac{{{{[\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A\hat X}}_M^{(k)}]}^{\rm{T}}}{\mathit{\boldsymbol{P}}_\Delta }[\mathit{\boldsymbol{L}} - \mathit{\boldsymbol{A\hat X}}_M^{(k)}]}}{{n - m}}} $$ (30)
      $$ {\mathit{\boldsymbol{s}}_2^2(k + 1) = \frac{{{{[\mathit{\boldsymbol{h}} - \mathit{\boldsymbol{H\hat X}}_{{\rm{MS}}}^{(k)}]}^{\rm{T}}}{\mathit{\boldsymbol{P}}_e}[\mathit{\boldsymbol{h}} - \mathit{\boldsymbol{H\hat X}}_{{\rm{MS}}}^{(k)}]}}{{p - m}}} $$ (31)

      直到$\left\| {\mathit{\boldsymbol{s}}_1^2(k + 1) - \mathit{\boldsymbol{s}}_1^2(k)} \right\| < \varepsilon $和$\left\| {\mathit{\boldsymbol{s}}_2^2(k + 1) - \mathit{\boldsymbol{s}}_2^2(k)} \right\| < \varepsilon $同时成立时终止迭代,ε = 1.0 × 10- 7为预先给定的正数。

      方案6:采用本文有界不确定性加权混合估计算法。把e看作为m维有界随机误差噪声,满足eT P e-1 e < 1,把w看作为p维有界随机误差噪声,满足w T P w-1 w < 1,利用式(19)至式(21)进行计算。算法步聚为:

      1)由于观测向量L中,观测数n较小,${{\mathit{\boldsymbol{\hat X}}}_L}$和${{\mathit{\boldsymbol{\hat X}}}_h}$的精度都不是最好的,但把两个观测进行融合得到的混合估计${{\mathit{\boldsymbol{\hat X}}}_M}$精度有了较大的提高。由于在测绘工程实际中,单位权方差一般无法准确地知道,算例中提供的多种混合估计方案的精度都有了明显的提高。

      2)普通混合估计${{\mathit{\boldsymbol{\hat X}}}_M}$是线性无偏估计,它们的均方误差阵可以由前面推导的协方差公式计算出来,为了比较它们的优良性,本算例中分别计算它们的迹为$ {\rm{tr}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_L}}}) = 2.512{\kern 1pt} {\kern 1pt} 2$,$ {\rm{tr}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_h}}}) = 17.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 4$和$ {\rm{tr}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_M}}}) = 0.365{\kern 1pt} {\kern 1pt} {\kern 1pt} 6$。

      另外,从表 3中的偏差平方和也可以看出${{\mathit{\boldsymbol{\hat X}}}_M}$要优于${{\mathit{\boldsymbol{\hat X}}}_L}$和${{\mathit{\boldsymbol{\hat X}}}_h}$。一般情形下都会有:

      $$ \begin{array}{*{20}{l}} { {\rm{tr }}[ {\rm{MSE}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_L}}})] \le {\rm{tr }} [ {\rm{MSE}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_M}}})]}\\ { {\rm{tr }} [ {\rm{MSE}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_{\rm{e}}}}})] \le {\rm{tr }} [ {\rm{MSE}} ({\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_{{{\mathit{\boldsymbol{\hat X}}}_M}}})]} \end{array} $$

      即${{\mathit{\boldsymbol{\hat X}}}_L}$和${{\mathit{\boldsymbol{\hat X}}}_h}$的均方误差大于普通混合估计${{\mathit{\boldsymbol{\hat X}}}_M}$的均方误差[16]

      3)${{\mathit{\boldsymbol{\hat X}}}_{{\rm{MS}}}}$和重迭代加权估计$\mathit{\boldsymbol{\hat X}}_{{\rm{MS}}}^{\left( {k + 1} \right)}$都是利用两步估计计算得到的,它们解决了σ e2σ w2未知时的参数估计问题,但它们的协方差阵求解比较困难。无法进行均方误差的比较,其精度评估无法进行,文献[16]给出了${{\mathit{\boldsymbol{\hat X}}}_{{\rm{MS}}}}$协方差的一个计算公式,但非常复杂。表 3给出的偏差平方和只能作为一个参考,在实际应用中由于不知道真值,其实无法比较其算法的效果。

      4)采用迭代重加权法解决了σ e2σ w2未知时的参数估计问题[24],但是,迭代的每一步估计$\mathit{\boldsymbol{\hat X}}_{{\rm{MS}}}^{\left( k \right)}$是一个非线性估计,它的协方差阵计算非常复杂,而且迭代过程未必收敛,有时只能以两步估计的值为估计结果。因此,迭代重加权最小二乘法是参数的一个粗略估计,无法对其精度进行评估。

      5)利用本文算法中的式(21)可以计算出混合估计的权因子为a = 0.446 0,计算得到的椭球中心${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$ =(-2.742 3,1.652 3,3.650 2,-2.767 4,-2.666 6)T,由于a可以看作是一个非随机变量,则${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$也是一个线性估计。椭球的特征矩阵P ME为:

      $$ {\mathit{\boldsymbol{P}}_{{\rm{ME}}}} = \left[ {\begin{array}{*{20}{r}} {0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}&{ - 0.002{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}&{0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 1}&{ - 0.003{\kern 1pt} {\kern 1pt} {\kern 1pt} 6}&{ - 0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 6}\\ { - 0.002{\kern 1pt} {\kern 1pt} 5}&{0.121{\kern 1pt} {\kern 1pt} {\kern 1pt} 8}&{0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 2}&{ - 0.062{\kern 1pt} {\kern 1pt} 6}&{ - 0.003{\kern 1pt} {\kern 1pt} {\kern 1pt} 0}\\ {0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 1}&{0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 2}&{0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 1}&{ - 0.002{\kern 1pt} {\kern 1pt} 4}&{ - 0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 7}\\ { - 0.003{\kern 1pt} {\kern 1pt} {\kern 1pt} 6}&{ - 0.062{\kern 1pt} {\kern 1pt} 6}&{ - 0.002{\kern 1pt} {\kern 1pt} 4}&{0.275{\kern 1pt} {\kern 1pt} {\kern 1pt} 9}&{0.021{\kern 1pt} {\kern 1pt} {\kern 1pt} 4}\\ { - 0.001{\kern 1pt} {\kern 1pt} {\kern 1pt} 6}&{ - 0.003{\kern 1pt} {\kern 1pt} {\kern 1pt} 0}&{ - 0.000{\kern 1pt} {\kern 1pt} {\kern 1pt} 7}&{0.021{\kern 1pt} {\kern 1pt} {\kern 1pt} 4}&{0.011{\kern 1pt} {\kern 1pt} {\kern 1pt} 3} \end{array}} \right] $$

      特征矩阵的最小迹为:tr(P ME) = 0.410 6,根据Fogel的观点[8],椭球集合的中心${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$可作为参数真实值X的一个点估计,PME可以作为${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$的协方差。tr(P ME)可以看作是${{\mathit{\boldsymbol{\hat X}}}_{{\rm{ME}}}}$的均方差与前面的结果比较,它的精度非常接近${{\mathit{\boldsymbol{\hat X}}}_M}$的精度,但${{\mathit{\boldsymbol{\hat X}}}_M}$的计算需要已知σ e2σ w2,这在测量实际中很难达到。

    • 联合平差的原理和方法与线性模型中的混合估计理论是一致,异质观测、前期观测都可以看成是一种随机约束。大地测量中各种异方差多源观测模型进行融合时,先验信息的融合都需要进行混合估计,由于附加信息和样本信息在估计过程中作用是不均等的,需要建立新的加权平差准则,平衡先验约束和观测信息对参数估计的影响。大地测量中各种异构多源观测信息,可以用椭球近似来描述。椭球集员估计算法可以用来进行观测信息的融合,通过椭球特征迹的最小化可以求得混合估计的权来平衡先验约束和观测信息对参数估计的影响,使得加权混合估计方法能有效应用于大测量数据处理。

参考文献 (26)

目录

    /

    返回文章
    返回