-
总体最小二乘平差理论及应用研究是目前测绘数据处理领域的热点研究问题，总体最小二乘准则最早由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}} $)作为判别函数，算法效果显著。
-
$$ \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列向量；B是nm×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]T，hi是n×1的列向量；B=[B1T B2T … BmT]T，Bi是n×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 }} $T⊗In)B；Sh=($\mathit{\boldsymbol{\hat \beta }} $T⊗In)h；Nh、NB、NBh和NhB都是m×m的矩阵，矩阵元素分别为：Nh(i, j)=hiTP1hj，NB(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)。uh、uB分别为：
$$ {{\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模型问题中一般需要求解的未知量包括观测向量y的n个改正数，系数矩阵A的n×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的值。在观测值x和y的坐标上分别加上均值为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)分别为原始坐标系和目标坐标系中控制点的坐标值；sX、sY为沿X轴和Y轴的缩放比例，α为旋转参数；c1、c2为沿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，方差阵为DS和DT的随机误差，其中，
表 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\|$的大小(β为a1、b1、c1、a2、b2、c2组成的参数列向量)，计算结果如表 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中的判别函数，算法效果是最好的。相对权比λ的确定是本文算法的关键，如何从理论方面给出严密的确定相对权比的公式是下一步的研究方向。
Total Least Squares Adjustment of Partial Errors-in-Variables Model with Weight Scaling Factor
-
摘要: 针对观测向量和系数矩阵权分配不合理、验前随机模型不准确的情况，以部分误差变量（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}} $作为判别函数是判别函数最小化迭代算法中效果最好的。Abstract: As a prior stochastic model contains inaccurate information, the weight matrices of observation and coefficient matrix are unreasonable. To address this problem, we investigate the total least squares adjustment of partial errors-in-variables(PEIV)model with a weight scaling factor that adaptively adjusts the contribution of the observation and coefficient matrix to parameter estimation. A prior unite weight variance and minimum discriminate function method are deduced, so the proposed method is valid for a structured coefficient matrix. Some conclusions are drawn from simulations of straight line fitting and coordinate transformation. When the prior unit weight variances of observation and coefficient matrix are known and accurate, the prior unit weight variance method is very effective; the minimum discriminate function method with the ${{\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}} $ as its discriminate function to determine weight scaling factor yielded the best performace.
-
图 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 -
点击查看大图
图(3) / 表(4)
计量
- 文章访问数: 1817
- HTML全文浏览量: 187
- PDF下载量: 296
- 被引次数: 0
