-
变量误差(errors‐in‐variables,EIV)模型[1]是一种同时顾及系数矩阵和观测向量误差的回归模型。对此模型可以采用总体最小二乘(total least squares,TLS)估计进行求解,TLS估计本质上是对最小二乘(least squares,LS)估计的扩展,是LS估计的去正则化处理[2]。在系数矩阵存在误差的情况下,TLS估计优于LS估计,并且具有弱一致性和渐进无偏性[2-3]。
目前,测绘领域对TLS及其各种扩展算法的应用集中于直线拟合、基准变换、自回归分析等问题[4-5],在这些应用中,EIV模型存在结构特征,即系数矩阵和观测向量存在固定量和随机量及不同位置的随机量之间存在线性相关性。可以采用基于广义奇异值分解或者矩阵QR分解的混合最小二乘-总体最小二乘(least squares-total least squares,LS-TLS)算法[3, 6-7]进行求解,这些算法只针对系数矩阵存在固定列的情况;也可以采用观测权阵生成法则[8],将随机量纳入平差的随机模型,通过构造奇异的观测权阵并利用加权总体最小二乘(weighted total least squares,WTLS)估计进行求解,这种算法增加了待估量的个数;甚至可以采用各种约束总体最小二乘(constrained total least squares,CTLS)算法进行求解[9-10],不过这将增加求解的复杂程度。
EIV模型存在结构特征本质上是结构总体最小二乘(structured total least squares,STLS)问题。目前,已存在诸多算法可以解决这一问题,如文献[11-12]引入伪观测值将随机量纳入平差的函数模型;文献[13]将STLS问题转化为非线性非约束优化问题,然后采用牛顿法进行求解;文献[14]将参数和随机量视为两组待估量,进而将STLS问题转换为两个LS问题;文献[10]分析了线性及二次约束下的STLS问题。上述各种算法仅考虑了系数矩阵的结构特征。
文献[15]利用变量投影法生成结构矩阵,将STLS问题转化为非线性双重最小化问题,再采用牛顿法进行求解,但是这种算法推导过程相对复杂,迭代效率较低,并且不能给出参数的近似精度估计。对于这一问题,首先,本文顾及系数矩阵和观测向量中随机量的观测权阵,针对系数矩阵和观测向量存在结构特征的情况,利用变量投影法[15]将其展开成仿射函数形式;然后,将STLS问题转化为非线性约束优化问题,并通过拉格朗日法推导出了一种STLS估计算法,进一步对其特征进行了分析;最后,给出了单位权方差和参数协因数阵的一阶近似估计公式。通过两个实例,验证了本文算法的有效性,并分析了本文算法与已有处理STLS问题的相关算法之间的关系。
-
EIV模型如下:
$$ \mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}=(\mathit{\boldsymbol{A}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}})\mathit{\boldsymbol{x}}, m>n=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\left(\mathit{\boldsymbol{A}}\right) $$ (1) 式中,$ \mathit{\boldsymbol{y}}\in {\mathit{\boldsymbol{R}}}^{m\times 1} $为观测向量;ey为其误差向量;$ \mathit{\boldsymbol{A}}\in {\mathit{\boldsymbol{R}}}^{m\times n} $为系数矩阵;EA为其误差矩阵;$ \mathit{\boldsymbol{x}}\in {\mathit{\boldsymbol{R}}}^{m\times 1} $为参数向量。
STLS问题的函数模型如式(1)所示,该模型要求系数矩阵的增广矩阵$ \mathit{\boldsymbol{C}}=\left[\mathit{\boldsymbol{A}}\;\;\mathit{\boldsymbol{y}}\right] $与其误差矩阵$ \mathrm{\Delta }\mathit{\boldsymbol{C}}=\left[{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\right] $具有相同的结构。在此假设$ \mathit{\boldsymbol{C}} $仅存在某种仿射结构,由$ \mathit{\boldsymbol{C}} $中独立随机量组成的向量为$ \mathit{\boldsymbol{p}}\in {\mathit{\boldsymbol{R}}}^{{n}_{\mathit{p}}} $,$ {n}_{\mathit{p}} $为独立随机量个数。故根据变量投影法可以将$ \mathit{\boldsymbol{C}} $表示为如下仿射函数[15]:
$$ \mathit{\boldsymbol{C}}=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right)={\mathit{\boldsymbol{S}}}_{0}+\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{S}}}_{i}{\mathit{\boldsymbol{p}}}_{i}, {\mathit{\boldsymbol{S}}}_{0}, {\mathit{\boldsymbol{S}}}_{i}\in {\mathit{\boldsymbol{R}}}^{m\times (n+1)} $$ (2) 式中,$ {\mathit{\boldsymbol{S}}}_{0} $为固定量仿射结构矩阵,$ {\mathit{\boldsymbol{S}}}_{0} $中非零元素为$ \mathit{\boldsymbol{C}} $中固定量;$ {\mathit{\boldsymbol{S}}}_{i} $为随机量$ {\mathit{\boldsymbol{p}}}_{i} $对应的仿射结构矩阵,一般为稀疏矩阵。STLS问题的数学模型为:
$$ \left\{\begin{array}{l}\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{p}}-\mathrm{\Delta }\mathit{\boldsymbol{p}}){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=0\\ \mathrm{\Delta }\mathit{\boldsymbol{p}}~N(0, {\sigma }_{0}^{2}\mathit{\boldsymbol{Q}})\end{array}\right. $$ (3) 式中,$ {\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=\left[\begin{array}{c}\mathit{\boldsymbol{x}}\\ -1\end{array}\right] $为$ \mathit{\boldsymbol{x}} $的增广向量;$ \mathit{\boldsymbol{Q}} $为$ \mathit{\boldsymbol{p}} $的正定协因数阵;$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $为$ \mathit{\boldsymbol{p}} $的观测误差。STLS问题的目标函数为:
$$ \mathrm{m}\mathrm{i}\mathrm{n}{‖\mathrm{\Delta }\mathit{\boldsymbol{p}}‖}_{{\mathit{\boldsymbol{Q}}}^{-1}}^{2} $$ (4) 对STLS问题的函数模型,进一步变形为[15]:
$$ \begin{array}{l}\mathit{\boldsymbol{S}}(\mathit{\boldsymbol{p}}-\mathrm{\Delta }\mathit{\boldsymbol{p}}){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=0\iff \\ \;\;\;\;\;\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{S}}}_{i}\mathrm{\Delta }{\mathit{\boldsymbol{p}}}_{i}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\iff \\ \;\;\;\;\;\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{S}}}_{i}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\mathrm{\Delta }{\mathit{\boldsymbol{p}}}_{i}=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\iff \\ \;\;\;\;\;\;\;\;\mathit{\boldsymbol{G}}\left(\mathit{\boldsymbol{x}}\right)\mathrm{\Delta }\mathit{\boldsymbol{p}}=\mathit{\boldsymbol{r}}\left(\mathit{\boldsymbol{x}}\right)\end{array} $$ (5) 式中,$ \mathit{\boldsymbol{r}}\left(\mathit{\boldsymbol{x}}\right)=\mathit{\boldsymbol{S}}\left(\mathit{\boldsymbol{p}}\right){\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}=-(\mathit{\boldsymbol{y}}-\mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}}) $;结构矩阵$ \mathit{\boldsymbol{G}}\left(\mathit{\boldsymbol{x}}\right)=\left[\begin{array}{cccc}{\mathit{\boldsymbol{S}}}_{1}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}& {\mathit{\boldsymbol{S}}}_{2}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}& \cdots & {\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}{\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}\end{array}\right]=\mathit{\boldsymbol{B}}\mathit{\boldsymbol{\psi }} $,其中,$ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{cccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)& \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}\right)\end{array}\right] $,$ \mathit{\boldsymbol{B}}={\mathit{\boldsymbol{x}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}\otimes {\mathit{\boldsymbol{I}}}_{m} $,$ \mathrm{v}\mathrm{e}\mathrm{c}(\cdot ) $为矩阵向量化算子,$ \otimes $为克罗内克乘积算子。
此时,STLS问题可以表示为:$ \underset{\mathit{\boldsymbol{x}}, \mathrm{\Delta }\mathit{\boldsymbol{p}}}{\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathrm{\Delta }\mathit{\boldsymbol{p}}‖}_{{\mathit{\boldsymbol{Q}}}^{-1}}^{2} $。其中,
$$ \mathit{\boldsymbol{G}}\left(\mathit{\boldsymbol{x}}\right)\mathrm{\Delta }\mathit{\boldsymbol{p}}=\mathit{\boldsymbol{r}}\left(\mathit{\boldsymbol{x}}\right), \mathrm{ }\mathrm{ }\mathrm{\Delta }\mathit{\boldsymbol{p}}~N(0, {\sigma }_{0}^{2}\mathit{\boldsymbol{Q}}) $$ (6) 式(6)本质上是一个非线性约束优化问题,可以从目标函数出发采用拉格朗日乘数法进行求解。
-
对于STLS问题,由拉格朗日乘数法构造其目标函数为:
$$ \mathit{\Phi }(\mathit{\boldsymbol{x}}, \mathrm{\Delta }\mathit{\boldsymbol{p}}, \mathit{\boldsymbol{\lambda }})=\mathrm{\Delta }{\mathit{\boldsymbol{p}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}\mathrm{\Delta }\mathit{\boldsymbol{p}}-2{\mathit{\boldsymbol{\lambda }}}^{\mathrm{T}}\left[\mathit{\boldsymbol{G}}\right(\mathit{\boldsymbol{x}})\mathrm{\Delta }\mathit{\boldsymbol{p}}-\mathit{\boldsymbol{r}}(\mathit{\boldsymbol{x}}\left)\right] $$ (7) 将式(7)分别对$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $、$ \mathit{\boldsymbol{\lambda }} $、$ \mathit{\boldsymbol{x}} $求偏导并令其等于零得:
$$ \frac{\partial \mathit{\Phi } }{\partial \mathrm{\Delta }\mathit{\boldsymbol{p}}}=2{\mathit{\boldsymbol{Q}}}^{-1}\mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}-2\mathit{\boldsymbol{G}}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)}^{\mathrm{T}}\widehat{\mathit{\boldsymbol{\lambda }}}=0 $$ (8) $$ \frac{\partial \mathit{\Phi } }{\partial \mathit{\boldsymbol{\lambda }}}=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}-\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=0 $$ (9) $$ \frac{\partial \mathit{\Phi } }{\partial \mathit{\boldsymbol{x}}}=(\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}{)}^{\mathrm{T}}\widehat{\mathit{\boldsymbol{\lambda }}}=0 $$ (10) 式中,$ \left(\widehat{\cdot }\right) $表示相应参数的估值,并根据式(2)重构系数矩阵误差矩阵的估值为:
$$ {\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}=\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}^{\mathit{\boldsymbol{A}}}}{\mathit{\boldsymbol{S}}}_{i}^{\mathit{\boldsymbol{A}}}\left(\mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}_{i}^{\mathit{\boldsymbol{A}}}\right) $$ (11) 式中,$ \mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}_{i}^{\mathit{\boldsymbol{A}}} $和$ {\mathit{\boldsymbol{S}}}_{i}^{\mathit{\boldsymbol{A}}} $分别为系数矩阵中独立随机量估值和其对应的仿射结构矩阵;$ {n}_{\mathit{\boldsymbol{p}}}^{\mathit{\boldsymbol{A}}} $为系数矩阵中独立随机量的个数。
假设$ \widehat{\mathit{\boldsymbol{x}}} $已知,式(9)可视为经典平差中的条件平差模型,由条件平差法可得:
$$ \mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}=\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{\Gamma }}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)}^{-1}\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $$ (12) 其中,$ \mathit{\boldsymbol{\Gamma }}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $。
将式(12)代入式(8),并将方程两侧同乘$ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}} $得:
$$ 2\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)-2\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{\lambda }}}=0 $$ (13) $ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $具有凯利逆,故可得:
$$ \widehat{\mathit{\boldsymbol{\lambda }}}={\mathit{\pmb{\Gamma }}}^{-1}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{r}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $$ (14) 将式(14)代入式(10),变形,得:
$$ (\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}{)}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-\mathit{\boldsymbol{A}}\widehat{\mathit{\boldsymbol{x}}}+{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}})=0 $$ (15) 进一步可得:
$$ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $$ (16) 式中,$ \widehat{\mathit{\boldsymbol{A}}}=\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}} $。式(16)具有经典平差中参数平差模型最小二乘解的形式。
算法的迭代过程如下:
1)初始值:假设$ \mathit{\boldsymbol{A}} $无误差,通过LS估计求得参数的初值$ {\widehat{\mathit{\boldsymbol{x}}}}^{0} $。
2)由$ {\widehat{\mathit{\boldsymbol{x}}}}^{t} $通过式(12)可得$ \mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}^{t+1} $,进而通过式(11)重构$ {\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1} $,即$ {\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1}=\sum\limits _{i=1}^{{n}_{\mathit{\boldsymbol{p}}}^{\mathit{\boldsymbol{A}}}}{\mathit{\boldsymbol{S}}}_{i}^{\mathit{\boldsymbol{A}}}(\mathrm{\Delta }{\widehat{p}}_{i}^{\mathit{\boldsymbol{A}}}{)}^{t+1} $,进一步得$ {\widehat{\mathit{\boldsymbol{A}}}}^{t+1}=\mathit{\boldsymbol{A}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1} $,然后得到$ {\mathit{\boldsymbol{\Gamma }}}^{-1}\left({\widehat{\mathit{\boldsymbol{x}}}}^{t}\right) $,最终计算出新的参数估值为:
$$ {\widehat{\mathit{\boldsymbol{x}}}}^{t+1}=\left[\right({\widehat{\mathit{\boldsymbol{A}}}}^{t+1}{)}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\left({\widehat{\mathit{\boldsymbol{x}}}}^{t}\right){\widehat{\mathit{\boldsymbol{A}}}}^{t+1}{]}^{-1}·\\ \left({\widehat{\mathit{\boldsymbol{A}}}}^{t+1}{)}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\right({\widehat{\mathit{\boldsymbol{x}}}}^{t}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}^{t+1}{\widehat{\mathit{\boldsymbol{x}}}}^{t}) $$ (17) 3)重复步骤2),直到满足收敛条件$ ‖{\widehat{\mathit{\boldsymbol{x}}}}^{t+1}-{\widehat{\mathit{\boldsymbol{x}}}}^{t}‖<\epsilon $,其中,$ \epsilon $为收敛阈值。
对于STLS问题,假设只有$ \mathit{\boldsymbol{y}} $存在误差,可将EIV模型转化为高斯-马尔可夫(Gauss-Markov,GM)模型,即:
$$ \mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}}=\mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}} $$ (18) 当$ \mathit{\boldsymbol{A}} $存在误差$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $,将其代入式(18)有:
$$ (\mathit{\boldsymbol{A}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}})\mathit{\boldsymbol{x}}=(\mathit{\boldsymbol{y}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{x}})-({\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\mathit{\boldsymbol{x}})={\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}-{\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}} $$ (19) 顾及GM模型的假设,式(19)左侧的系数矩阵$ \mathit{\boldsymbol{A}}-{\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $无误差。同时,式(16)中ey受到$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $的增益,故$ {\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}} $的协因数阵为$ {\mathit{\boldsymbol{Q}}}_{{\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}}}=\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}=\mathit{\boldsymbol{\Gamma }}\left(\mathit{\boldsymbol{x}}\right) $。进而由参数平差法得:
$$ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}){\widehat{\mathit{\boldsymbol{y}}}}_{\mathit{\boldsymbol{c}}} $$ (20) 将式(12)和式(20)相结合用于求解STLS问题。由上述分析可得,TLS问题的求解可以转化为两个LS问题的求解。这一思想在文献[14]中得到了体现,该文献将STLS问题中的待估量分为$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $和$ \mathit{\boldsymbol{x}} $两类,利用式(12)求解$ \mathrm{\Delta }\mathit{\boldsymbol{p}} $,通过直接改正系数矩阵误差求解$ \mathit{\boldsymbol{x}} $,即:
$$ \widehat{\mathit{\boldsymbol{x}}}=({\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}{\widehat{\mathit{\boldsymbol{A}}})}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}\mathit{\boldsymbol{y}} $$ (21) 当$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\to 0 $,式(20)趋近于式(21),因而,式(21)可以视为式(20)的简化形式。在建立求解$ \mathit{\boldsymbol{x}} $的GM模型时,式(21)并没有顾及$ \mathit{\boldsymbol{A}} $的误差$ {\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}} $改正后对GM模型的影响。这一影响不断被$ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}} $吸收,影响迭代效率和迭代稳定性,同时,也无法直接获得TLS估计下$ \mathit{\boldsymbol{y}} $的改正数。
-
根据TLS原理,上述STLS估计算法的单位权方差可以表示为:
$$ {\sigma }^{2}=\frac{\mathrm{\Delta }{\widehat{\mathit{\boldsymbol{p}}}}^{\mathrm{T}}{\mathit{\boldsymbol{Q}}}^{-1}\mathrm{\Delta }\widehat{\mathit{\boldsymbol{p}}}}{d} $$ (22) 式中,$ d=m-n $为自由度,其数值与观测方程个数和参数个数有关。本文得到的参数$ \mathit{\boldsymbol{x}} $的估计表达式具有LS解的形式,因而,可以采用协因数传播律求解参数的协因数阵。由§1.2的分析可知,$ {\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}} $对应的误差向量$ {\mathit{\boldsymbol{e}}}_{{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{c}}}} $的协因数阵为$ \mathit{\boldsymbol{\Gamma }}\left(\mathit{\boldsymbol{x}}\right) $,所以根据协因数传播律得:
$$ {\mathit{\boldsymbol{Q}}}_{\widehat{\mathit{\boldsymbol{x}}}\widehat{\mathit{\boldsymbol{x}}}}=[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\pmb{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1} $$ (23) 式(23)即为STLS问题参数协因数阵的一阶近似估计公式。
-
选取的直线拟合算例[16]给出了10组坐标$ ({u}_{i}, {v}_{i}) $用于拟合直线的截距$ a $和斜率$ b $,其中,i=1,2…10。其对应的EIV模型为:
$$ \left[\begin{array}{c}{v}_{1}\\ {v}_{2}\\ \vdots\\ {v}_{10}\end{array}\right]-\left[\begin{array}{c}{e}_{{v}_{1}}\\ {e}_{{v}_{2}}\\ \vdots \\ {e}_{{v}_{10}}\end{array}\right]=\left(\left[\begin{array}{cc}\begin{array}{c}1\\ 1\\ \vdots\\ 1\end{array}& \begin{array}{c}{u}_{1}\\ {u}_{2}\\ \vdots\\ {u}_{10}\end{array}\end{array}\right]-\left[\begin{array}{cc}\begin{array}{c}0\\ 0\\ \vdots\\ 0\end{array}& \begin{array}{c}{e}_{{u}_{1}}\\ {e}_{{u}_{2}}\\ \vdots\\ {e}_{{u}_{10}}\end{array}\end{array}\right]\right)\cdot \left[\begin{array}{c}a\\ b\end{array}\right] $$ (24) 式中,$ ({e}_{{u}_{i}}, {e}_{{v}_{i}}) $为10组坐标对应的误差,其中,i=1,2…10。由式(24)可知,截距$ a $对应的系数矩阵元素的误差理论上为零,因而本质上直线拟合问题是STLS问题。10组坐标及其对应的权如表 1所示[16]。
表 1 观测点和对应的权
Table 1. Observed Points and Corresponding Weights
点号 坐标 坐标的权 u v Wu Wv 1 0.0 5.9 1 000.0 1.0 2 0.9 5.4 1 000.0 1.8 3 1.8 4.4 500.0 4.0 4 2.6 4.6 800.0 8.0 5 3.3 3.5 200.0 20.0 6 4.4 3.7 80.0 20.0 7 5.2 2.8 60.0 70.0 8 6.1 2.8 20.0 70.0 9 6.5 2.4 1.8 100.0 10 7.4 1.5 1.0 500.0 将表 1中的数据分别采用3种WTLS估计算法(文献[17]算法、文献[18]算法和文献[19]算法)和2种STLS估计算法(文献[15]算法和本文算法)进行求解。各种迭代算法收敛阈值设为$ 1\times {10}^{-10} $,相应的参数估值$ \widehat{a} $和$ \widehat{b} $、验后单位权方差估值$ {\widehat{\sigma }}_{0}^{2} $、参数估计精度$ {\widehat{\sigma }}_{\widehat{a}}^{2} $和$ {\widehat{\sigma }}_{\widehat{b}}^{2} $以及迭代次数如表 2所示。
表 2 不同算法所得直线拟合结果
Table 2. Line Fitting Results from Different Methods
参数 文献[17]算法 文献[18]算法 文献[19]算法 文献[15]算法 本文算法 $ \widehat{a} $ 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 $ \widehat{b} $ -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 $ \widehat{\sigma }{}_{0}^{2} $ 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 $ \widehat{\sigma }{}_{\widehat{a}}^{2} $ 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0 $ \widehat{\sigma }{}_{\widehat{b}}^{2} $ 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5 迭代次数 13 8 7 9 7 由表 2可知,各种算法均能得到一致的估计结果,即本文算法与现有的WTLS估计算法结果一致,说明本文算法是有效的。它们的收敛次数有差别:文献[17]算法需要13次迭代收敛;文献[15]算法优于文献[17]算法,需要9次迭代收敛;文献[18]算法的迭代收敛减少到8次;文献[19]算法和本文算法迭代收敛次数最少,均为7次;本文算法相对文献[15]算法迭代效率有提高。
通过构造秩亏的观测权阵[8],各种WTLS估计算法[17-19]可以解决STLS问题,但是WTLS估计需要求解系数矩阵的增广矩阵中的所有元素,因此,其待估量个数和观测权阵维数大于STLS估计。参数估计表达式的差异是决定迭代效率的重要因素,上述各种WTLS估计算法中,文献[19]算法迭代效率最高,并且参数估计表达式具有LS解的形式,本文给出的参数估计表达式结构与之相同,这也是本文算法在上述几种加权或者结构总体最小估计算法中迭代效率最高的原因。
-
自回归分析是统计上处理时间序列的一种常用方法。文献[20]给出了线性自回归模型的一个测量实例。该实例为对某一建筑物进行了36期沉降观测(如表 3所示[20]),并且认为不同期观测独立等精度。
表 3 高程数据
Table 3. Elevation Data
点号 高程/mm 点号 高程/mm 点号 高程/mm 点号 高程/mm 1 26.33 10 25.46 19 28.09 28 26.57 2 26.27 11 26.12 20 26.78 29 28.36 3 26.43 12 27.28 21 28.66 30 27.94 4 25.56 13 26.67 22 26.75 31 26.81 5 26.82 14 27.95 23 27.24 32 28.50 6 26.56 15 26.74 24 28.02 33 27.68 7 25.93 16 27.53 25 26.81 34 26.57 8 26.43 17 25.31 26 28.50 35 28.36 9 26.52 18 26.90 27 27.68 36 27.94 为便于比较,参照文献[20],采用三阶线性自回归模型进行建模。对应的EIV模型如下:
$$ \left\{\begin{array}{l}{H}_{4}^{\mathrm{*}}={H}_{1}^{\mathrm{*}}{\xi }_{1}+{H}_{2}^{\mathrm{*}}{\xi }_{2}+{H}_{3}^{\mathrm{*}}{\xi }_{3}\\ {H}_{5}^{\mathrm{*}}={H}_{2}^{\mathrm{*}}{\xi }_{1}+{H}_{3}^{\mathrm{*}}{\xi }_{2}+{H}_{4}^{\mathrm{*}}{\xi }_{3}\\ ⋮\\ {H}_{36}^{\mathrm{*}}={H}_{33}^{\mathrm{*}}{\xi }_{1}+{H}_{34}^{\mathrm{*}}{\xi }_{2}+{H}_{35}^{\mathrm{*}}{\xi }_{3}\end{array}\right. $$ (25) 式中,$ {\xi }_{1} $、$ {\xi }_{2} $、$ {\xi }_{3} $表示自回归模型的参数;$ {H}_{i}^{\mathrm{*}} $表示第i期沉降观测值,其中,i=1,2…36;“*”表示随机量。
利用本文算法进行解算时,还需要给出$ \mathit{\boldsymbol{\psi }} $,系数矩阵的增广矩阵C可表示为:
$$ \mathit{\boldsymbol{C}}=\left[\begin{array}{cccc}{H}_{1}^{\mathrm{*}}& {H}_{2}^{\mathrm{*}}& {H}_{3}^{\mathrm{*}}& {H}_{4}^{\mathrm{*}}\\ {H}_{2}^{\mathrm{*}}& {H}_{3}^{\mathrm{*}}& {H}_{4}^{\mathrm{*}}& {H}_{5}^{\mathrm{*}}\\ ⋮& ⋮& ⋮& ⋮\\ {H}_{33}^{\mathrm{*}}& {H}_{34}^{\mathrm{*}}& {H}_{35}^{\mathrm{*}}& {H}_{36}^{\mathrm{*}}\end{array}\right] $$ (26) 根据式(2)构造$ \mathit{\boldsymbol{C}} $的仿射结构矩阵$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{36} $,$ \mathit{\boldsymbol{C}} $具有Hankel矩阵结构。不难发现,依次保留$ \mathit{\boldsymbol{C}} $各主对角线元素并置为1,即为$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{36} $,进而可得$ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{ccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{36}\right)\end{array}\right] $。
分别采用1种WTLS估计算法(文献[19]算法)和4种STLS估计算法(文献[10]算法、文献[21]算法、文献[15]算法、本文算法)对表 3中数据进行求解,设置收敛阈值为$ 1\times {10}^{-10} $,相应的参数估值($ {\widehat{\xi }}_{1} $、$ {\widehat{\xi }}_{2} $和$ {\widehat{\xi }}_{3} $)、验后单位权方差估值$ {\widehat{\sigma }}_{0}^{2} $、参数估计精度($ {\widehat{\sigma }}_{{\widehat{\xi }}_{1}}^{2} $、$ {\widehat{\sigma }}_{{\widehat{\xi }}_{2}}^{2} $和$ {\widehat{\sigma }}_{{\widehat{\xi }}_{3}}^{2} $)及迭代次数如表 4所示。
表 4 线性自回归模型不同算法所得结果
Table 4. Results from Different Methods of Linear Auto-Regression Model
参数 文献[19]算法 文献[10]算法 文献[21]算法 文献[15]算法 本文算法 $ {\widehat{\xi }}_{1} $ 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 $ {\widehat{\xi }}_{2} $ 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 $ {\widehat{\xi }}_{3} $ -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 $ \widehat{\sigma }{}_{0}^{2} $ 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 $ \widehat{\sigma }{}_{{\widehat{\xi }}_{1}}^{2} $ 0.012 658 067 0 0.012 658 067 0 0.012 658 067 0 ‒ 0.012 658 067 0 $ \widehat{\sigma }{}_{{\widehat{\xi }}_{2}}^{2} $ 0.009 481 774 9 0.009 481 774 9 0.009 481 774 9 ‒ 0.009 481 774 9 $ \widehat{\sigma }{}_{{\widehat{\xi }}_{3}}^{2} $ 0.009 074 553 9 0.009 074 553 9 0.009 074 553 9 ‒ 0.009 074 553 9 迭代次数 43 43 43 60 43 算例1中的直线拟合问题仅是系数矩阵存在固定列,为进一步说明本文算法的有效性,选取了线性自回归问题,该问题中系数矩阵的增广矩阵为Hankel矩阵。由表 4可知,各种算法估计结果具有一致性,验证了本文算法的有效性。在迭代效率方面,文献[15]算法需要60次迭代收敛;其余各种算法优于文献[15]算法,均需要43次迭代收敛,可见,本文给出的STLS估计算法迭代效率同样优于文献[15]算法。
通过提取独立随机量建立结构矩阵,STLS问题可以利用各种STLS估计算法[10, 15, 21]进行求解。文献[10]算法、文献[21]算法以及本文算法参数估计表达式和结构矩阵提取方式如表 5所示。
表 5 参数估计表达式比较
Table 5. Comparison of Parameter Estimation Expression
算法 参数估计表达式 结构矩阵提取方法 文献[10]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \widehat{\mathit{\boldsymbol{B}}}={\widehat{\mathit{\boldsymbol{x}}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}\otimes {\mathit{\boldsymbol{I}}}_{m} $,$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{l}}l}=\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\mathit{\boldsymbol{S}}}^{\mathrm{T}} $;$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵;$ \widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{S}} $为结构矩阵 $ \mathit{\boldsymbol{S}}\cdot \mathit{\boldsymbol{e}}=\left[\begin{array}{cc}{\mathit{\boldsymbol{C}}}_{\mathit{\boldsymbol{A}}}& 0\\ 0& {I}_{m}\end{array}\right]\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $ 文献[21]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,结构矩阵$ {\widehat{\mathit{\boldsymbol{C}}}}_{2}=\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{C}}}_{1} $ $ {\mathit{\boldsymbol{C}}}_{1}\cdot \mathit{\boldsymbol{e}}={\mathit{\boldsymbol{C}}}_{1}\cdot \left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $ 本文算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \mathit{\boldsymbol{\Gamma }}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $,其中,结构矩阵$ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{\psi }} $ $ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{cccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)& \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}\right)\end{array}\right] $通过变量投影法建立$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}} $ 文献[19]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\right) $;$ {\mathit{\boldsymbol{Q}}}_{1} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵 ‒ 由表 5可知,文献[10]算法、文献[21]算法和本文算法参数估计表达式的结构与文献[19]算法均是相同的,所以它们的迭代效率相同,这种参数估计表达式更早见于文献[22-23]。但是3种方法的结构矩阵提取方法存在差别,其中,文献[10]算法和文献[21]算法结构矩阵提取时,并不考虑系数矩阵和观测向量包含相同随机量的情况,因而,在处理自回归模型时需要利用文献[8]中给出的方法生成秩亏的观测权阵。理论上,文献[21]算法给出的矩阵$ {\mathit{\boldsymbol{C}}}_{1} $形式并不唯一,其中一种形式即为文献[10]算法中的矩阵$ \mathit{\boldsymbol{S}} $。当然各种算法公式推导过程等方面也存在差异,详情参考文献[10, 21]。
-
文献[15]引入变量投影法构造结构矩阵,这种方法原理简单,规律性强,独立于求解问题。但是文献[15]的算法推导过程相对复杂、计算效率较低,并且不能给出参数的近似精度估计。为此,本文给出了一种将系数矩阵的增广矩阵中随机量纳入随机模型的新思路,通过变量投影法对函数模型进行重构,将系数矩阵的增广矩阵表示为独立随机量的仿射函数形式,采用非线性拉格朗日乘数法得到了一种STLS估计算法,并给出了单位权方差和参数的一阶近似协因数阵估计公式。
1)两个算例表明,本文算法与已有能够解决系数矩阵和观测向量存在结构特征的加权或者结构总体最小二乘算法估计结果一致,说明了本文算法的有效性,并且相比文献[15]的算法提高了迭代效率。
2)与通过构造秩亏观测权阵的各种WTLS估计算法相比,本文算法仅需要建立系数矩阵的增广矩阵中独立随机量的正定观测权阵,降低了待估量个数。与其他STLS估计算法相比,本文引入了变量投影法这一结构矩阵构造方法,顾及了EIV模型中系数矩阵和观测向量的整体结构,即考虑了系数矩阵和观测向量包含相同随机量的情况,因而,在处理诸如自回归模型等问题时,无需再建立秩亏的观测权阵。
3)通过对各种算法的比较可知,STLS估计和WTLS估计分别从从函数模型和随机模型的角度确保TLS估计的统计最优性,WTLS估计和STLS估计具有等价性。
4)本文给出的STLS算法具有一般性,适用于系数矩阵和观测向量包含随机量、固定量以及固定列等情况,当系数矩阵全为固定值,本文算法转变成加权LS估计算法,当系数矩阵和观测向量仅包含随机量,本文算法等价于WTLS估计算法。
-
摘要: 测绘领域诸多实际应用中系数矩阵和观测向量具有结构特征,即系数矩阵和观测向量中包含固定量(甚至固定列) 和随机量,并且不同位置的随机量线性相关。针对这个问题,从变量误差(errors‐in‐variables,EIV)函数模型出发,首先,将系数矩阵和观测向量构成的增广矩阵表示为仿射函数形式,并采用变量投影法对函数模型进行重构;然后,利用拉格朗日法推导出了一种结构总体最小二乘(structured total least squares,STLS) 估计算法。算例分析结果表明,该算法与已有能够解决系数矩阵和观测向量存在结构特征的加权或结构总体最小二乘算法估计结果一致,说明了该算法的有效性,同时阐明了该算法与已有相关算法的关系。Abstract:
Objectives To solve the problem that the coefficient matrix and the observation vector have structured characteristics for many practical applications in the field of surveying and mapping, that is, the coefficient matrix and the observation vector contain fixed quantities (or even fixed columns) and random quantities, and the random quantities at different positions are linearly related. Methods Starting from the errors-in-variables (EIV) function model, the augmented matrix composed of the coefficient matrix and the observation vector is expressed as an affine function form, and the function model is reconstructed by the variable projection method. Then, a structured total least squares (STLS) estimation algorithm is derived by the Lagrange method. Results The example results show that the proposed algorithm is consistent with the existing weighted or structured total least squares estimation algorithms that can solve structured problem in the coefficient matrix and the observation vector. Compared with the weighted total least squares (WTLS) estimation algorithms, the proposed algorithm only needs to establish a positive definite weight matrix of independent random variables and reduces the number of estimates. Compared with other STLS estimation algorithms, the proposed algorithm takes into account the overall structure of the coefficient matrix and the observation vector. Conclusions This shows the effectiveness of the proposed algorithm. The STLS estimation and the WTLS estimation ensure the statistical optimality from the perspective of the function model and the stochastic model, respectively. -
表 1 观测点和对应的权
Table 1. Observed Points and Corresponding Weights
点号 坐标 坐标的权 u v Wu Wv 1 0.0 5.9 1 000.0 1.0 2 0.9 5.4 1 000.0 1.8 3 1.8 4.4 500.0 4.0 4 2.6 4.6 800.0 8.0 5 3.3 3.5 200.0 20.0 6 4.4 3.7 80.0 20.0 7 5.2 2.8 60.0 70.0 8 6.1 2.8 20.0 70.0 9 6.5 2.4 1.8 100.0 10 7.4 1.5 1.0 500.0 表 2 不同算法所得直线拟合结果
Table 2. Line Fitting Results from Different Methods
参数 文献[17]算法 文献[18]算法 文献[19]算法 文献[15]算法 本文算法 $ \widehat{a} $ 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 5.479 910 224 0 $ \widehat{b} $ -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 -0.480 533 407 4 $ \widehat{\sigma }{}_{0}^{2} $ 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 1.483 294 149 3 $ \widehat{\sigma }{}_{\widehat{a}}^{2} $ 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0 0.129 058 064 0 $ \widehat{\sigma }{}_{\widehat{b}}^{2} $ 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5 0.004 987 222 5 迭代次数 13 8 7 9 7 表 3 高程数据
Table 3. Elevation Data
点号 高程/mm 点号 高程/mm 点号 高程/mm 点号 高程/mm 1 26.33 10 25.46 19 28.09 28 26.57 2 26.27 11 26.12 20 26.78 29 28.36 3 26.43 12 27.28 21 28.66 30 27.94 4 25.56 13 26.67 22 26.75 31 26.81 5 26.82 14 27.95 23 27.24 32 28.50 6 26.56 15 26.74 24 28.02 33 27.68 7 25.93 16 27.53 25 26.81 34 26.57 8 26.43 17 25.31 26 28.50 35 28.36 9 26.52 18 26.90 27 27.68 36 27.94 表 4 线性自回归模型不同算法所得结果
Table 4. Results from Different Methods of Linear Auto-Regression Model
参数 文献[19]算法 文献[10]算法 文献[21]算法 文献[15]算法 本文算法 $ {\widehat{\xi }}_{1} $ 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 1.179 081 343 2 $ {\widehat{\xi }}_{2} $ 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 0.041 899 550 4 $ {\widehat{\xi }}_{3} $ -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 -0.214 448 099 2 $ \widehat{\sigma }{}_{0}^{2} $ 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 0.413 991 225 1 $ \widehat{\sigma }{}_{{\widehat{\xi }}_{1}}^{2} $ 0.012 658 067 0 0.012 658 067 0 0.012 658 067 0 ‒ 0.012 658 067 0 $ \widehat{\sigma }{}_{{\widehat{\xi }}_{2}}^{2} $ 0.009 481 774 9 0.009 481 774 9 0.009 481 774 9 ‒ 0.009 481 774 9 $ \widehat{\sigma }{}_{{\widehat{\xi }}_{3}}^{2} $ 0.009 074 553 9 0.009 074 553 9 0.009 074 553 9 ‒ 0.009 074 553 9 迭代次数 43 43 43 60 43 表 5 参数估计表达式比较
Table 5. Comparison of Parameter Estimation Expression
算法 参数估计表达式 结构矩阵提取方法 文献[10]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{ll}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \widehat{\mathit{\boldsymbol{B}}}={\widehat{\mathit{\boldsymbol{x}}}}_{\mathrm{e}\mathrm{x}\mathrm{t}}^{\mathrm{T}}\otimes {\mathit{\boldsymbol{I}}}_{m} $,$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{l}}l}=\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\mathit{\boldsymbol{S}}}^{\mathrm{T}} $;$ {\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵;$ \widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{S}} $为结构矩阵 $ \mathit{\boldsymbol{S}}\cdot \mathit{\boldsymbol{e}}=\left[\begin{array}{cc}{\mathit{\boldsymbol{C}}}_{\mathit{\boldsymbol{A}}}& 0\\ 0& {I}_{m}\end{array}\right]\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $ 文献[21]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left({\widehat{\mathit{\boldsymbol{C}}}}_{2}{\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{e}}}{\widehat{\mathit{\boldsymbol{C}}}}_{2}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,结构矩阵$ {\widehat{\mathit{\boldsymbol{C}}}}_{2}=\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{C}}}_{1} $ $ {\mathit{\boldsymbol{C}}}_{1}\cdot \mathit{\boldsymbol{e}}={\mathit{\boldsymbol{C}}}_{1}\cdot \left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{a}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]=\left[\begin{array}{c}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{A}}}\\ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right] $ 本文算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}{\left(\widehat{\mathit{\boldsymbol{x}}}\right)\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}{\mathit{\boldsymbol{\Gamma }}}^{-1}\right(\widehat{\mathit{\boldsymbol{x}}}\left)\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ \mathit{\boldsymbol{\Gamma }}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{G}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{x}}}\right) $,其中,结构矩阵$ \mathit{\boldsymbol{G}}\left(\widehat{\mathit{\boldsymbol{x}}}\right)=\widehat{\mathit{\boldsymbol{B}}}\mathit{\boldsymbol{\psi }} $ $ \mathit{\boldsymbol{\psi }}=\left[\begin{array}{cccc}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{1}\right)& \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{2}\right)& \cdots & \mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}}\right)\end{array}\right] $通过变量投影法建立$ {\mathit{\boldsymbol{S}}}_{1}~{\mathit{\boldsymbol{S}}}_{{n}_{\mathit{\boldsymbol{p}}}} $ 文献[19]算法 $ \widehat{\mathit{\boldsymbol{x}}}=\left[{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\right(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}{\widehat{\mathit{\boldsymbol{A}}}]}^{-1}{\widehat{\mathit{\boldsymbol{A}}}}^{\mathrm{T}}\left(\widehat{\mathit{\boldsymbol{B}}}{\mathit{\boldsymbol{Q}}}_{1}{\widehat{\mathit{\boldsymbol{B}}}}^{\mathrm{T}}{)}^{-1}\right(\mathit{\boldsymbol{y}}-{\widehat{\mathit{\boldsymbol{E}}}}_{\mathit{\boldsymbol{A}}}\widehat{\mathit{\boldsymbol{x}}}) $式中,$ {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{\boldsymbol{E}}}_{\mathit{\boldsymbol{A}}}\right) $;$ {\mathit{\boldsymbol{Q}}}_{1} $为$ {\left[\begin{array}{cc}{\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{{\rm A}}}}& {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{y}}}\end{array}\right]}^{\mathrm{T}} $的协因数阵 ‒ -
[1] Fuller W A. Measurement Error Models[M]. New York: Wiley, 1987 [2] Golub G H, Van Loan C F. An Analysis of the Total Least Squares Problem[J]. SIAM Journal on Numerical Analysis, 1980, 17(6): 883-893 doi: 10.1137/0717073 [3] Van Huffel S, van Dewalle J. The Total Least Squares Problem: Computational Aspects and Analysis[M]. Philadelphia: SIAM, 1991 [4] 刘经南, 曾文宪, 徐培亮. 整体最小二乘估计的研究进展[J]. 武汉大学学报·信息科学版, 2013, 38(5): 505-512 http://ch.whu.edu.cn/article/id/2641 Liu Jingnan, Zeng Wenxian, Xu Peiliang. Overview of Total Least Squares Methods[J]. Geomatics and Information Science of Wuhan University, 2013, 38(5): 505-512 http://ch.whu.edu.cn/article/id/2641 [5] 王乐洋, 许才军. 总体最小二乘研究进展[J]. 武汉大学学报·信息科学版, 2013, 38(7): 850-856 http://ch.whu.edu.cn/article/id/2703 Wang Leyang, Xu Caijun. Progress in Total Least Squares[J]. Geomatics and Information Science of Wuhan University, 2013, 38(7): 850-856 http://ch.whu.edu.cn/article/id/2703 [6] Dunne B E, Williamson G A. QR-Based TLS and Mixed LS-TLS Algorithms with Applications to Adaptive ⅡR Filtering[J]. IEEE Transactions on Signal Processing, 2003, 51(2): 386-394 doi: 10.1109/TSP.2002.806980 [7] Yan Shijian, Fan Jinyan. The Solution Set of the Mixed LS-TLS Problem[J]. International Journal of Computer Mathematics, 2011, 77(4): 545-561 [8] Mahboub V. On Weighted Total Least-Squares for Geodetic Transformations[J]. Journal of Geodesy, 2012, 86(5): 359-367 doi: 10.1007/s00190-011-0524-5 [9] Schaffrin B. A Note on Constrained Total Least-Squares Estimation[J]. Linear Algebra and Its Applications, 2006, 417(1): 245-258 doi: 10.1016/j.laa.2006.03.044 [10] Fang Xing. A Structured and Constrained Total Least-Squares Solution with Cross-Covariances[J]. Studia Geophysica et Geodaetica, 2014, 58(1), DOI: 10.1007/s11200-012-0671-z [11] Xu Peiliang, Liu Jiangnan, Shi Chuang. Total Least Squares Adjustment in Partial Errors-in-Variables Models: Algorithm and Statistical Analysis[J]. Journal of Geodesy, 2012, 86(8): 661-675 doi: 10.1007/s00190-012-0552-9 [12] Shi Yun, Xu Peiliang, Liu Jingnan, et al. Alternative Formulae for Parameter Estimation in Partial Errors-in-Variables Models[J]. Journal of Geodesy, 2015, 89(1): 13-16 doi: 10.1007/s00190-014-0756-2 [13] Rosen J B, Park H, Glick J. Total Least Norm Formulation and Solution for Structured Problems[J]. SIAM Journal on Matrix Analysis and Applications, 1996, 17(1): 110-126 doi: 10.1137/S0895479893258802 [14] Yalamov P Y, Yuan J Y. A Successive Least Squares Method for Structured Total Least Squares[J]. Journal of Computational Mathematics, 2003, 21(4): 463-472 [15] Markovsky I, Van Huffel S. High-performance Numerical Algorithms and Software for Structured Total Least Squares[J]. Journal of Computational and Applied Mathematics, 2005, 180(2): 311-331 doi: 10.1016/j.cam.2004.11.003 [16] Neri F, Saitta G, Chiofalo S. An Accurate and Straight Forward Approach to Line Regression Analysis of Error-Affected Experimental Data[J]. Journal of Physics E: Scientific Instruments, 1989, 22(4): 215-217 doi: 10.1088/0022-3735/22/4/002 [17] Shen Yunzhong, Li Bofeng, Chen Yi. An Iterative Solution of Weighted Total Least-Squares Adjustment[J]. Journal of Geodesy, 2011, 85(4): 229-238 doi: 10.1007/s00190-010-0431-1 [18] Snow K. Topics in Total Least-Squares Adjustment within the Errors-in-Variables Model: Singular Cofactor Matrices and Prior Information[D]. Columbus: The Ohio State University, 2012 [19] Amiri-Simkooei A, Jazaeri S. Weighted Total Least Squares Formulated by Standard Least Squares Theory[J]. Journal of Geodetic Science, 2012, 2(2): 113-124 doi: 10.2478/v10156-011-0036-5 [20] 王新洲, 陶本藻, 邱卫宁, 等. 高等测量平差[M]. 北京: 测绘出版社, 2006 Wang Xinzhou, Tao Benzao, Qiu Weining, et al. Higher Surveying Adjustment[M]. Beijing: Surveying and Mapping Press, 2006 [21] 王乐洋, 许光煜, 温贵森. 一种相关观测的Partial EIV模型求解方法[J]. 测绘学报, 2017, 46(8): 978-987 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201708008.htm Wang Leyang, Xu Guangyu, Wen Guisen. A Method for Partial EIV Model with Correlated Observations[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(8): 978-987 https://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201708008.htm [22] Fang Xing. Weighted Total Least Squares Solution for Application in Geodesy[D]. Hanover: Leibniz University Hanover, 2011 [23] Fang Xing. Weighted Total Least Squares: Necessary and Sufficient Conditions, Fixed and Random Parameters[J]. Journal of Geodesy, 2013, 87(8): 733-749 doi: 10.1007/s00190-013-0643-2 -