变量误差(errors-in-variables，EIV)模型及其总体最小二乘(total least squares，TLS)估计研究已逾百年，取得了一系列重要成果，并在科学和工程中得到了广泛应用^[1]。测绘领域对TLS估计的应用集中于线性回归、基准变换和摄影测量空间后方交会等问题^[2-3]，在这些应用中，EIV模型的系数矩阵存在结构特征，即系数矩阵存在固定量和随机量以及不同位置的随机量之间存在相关性。对于系数矩阵存在固定列的情况，可以采用基于矩阵QR分解的混合LS-TLS（least squares-total least squares）估计求解^[1]。而更一般的情况，可以采用各种结构总体最小二乘(structural total least squares, STLS)估计^[4-7]求解。Xu等^[8]提出了Partial EIV (PEIV)模型，通过构造虚拟观测值推导出一种STLS问题的算法，王乐洋等^[9-10]对上述算法进行了改进，提高了算法的计算效率，且适用于相关观测PEIV模型的参数估计。也可以通过构造秩亏的观测权阵，利用加权总体最小二乘(weighted total least squares，WTLS)估计^[11-12]求解STLS问题。

自回归(autoregressive, AR)模型是一类重要的时间序列模型，它在科学研究和工程实践中有广泛的应用，如形变分析与预报、导航定位、极移预测、陀螺漂移分析等^[13-14]。AR模型的求解实际上是一个STLS问题。文献[13]推导了一种非线性的AR模型TLS迭代算法，但难以采用协方差传播定律进行精度评定；文献[14]通过引入虚拟观测值，将TLS问题转化为经典的最小二乘问题，采用间接平差方法进行参数求解。AR模型观测值可能来源于不同类型的观测或平差，甚至来自不同期观测，因而具有不同的精度。已有算法^{[10, 13-14]}认为，系数矩阵和观测向量中的随机量拥有相同的验前单位权方差并不合适，需要确定合理的观测权矩阵以确保参数的最优估计，因此需要进行方差分量估计，对随机模型进行修正。

对于高斯-马尔可夫(Gauss-Markov，GM)模型，文献[15]提出了最小二乘方差分量估计(least squares variance component estimation, LS-VCE)，文献[16]将LS-VCE应用到GPS的数据处理中，说明了这种方法的有效性。对于EIV模型，文献[17]给出了EIV模型的LS-VCE算法；文献[18]基于LS-VCE原理提出了一种严密的TLS估计算法，实现了参数和方差分量同步估计，但这种算法并不适用于包含固定列的系数矩阵，因而作者进一步给出了两种替代算法；文献[19]推导了EIV模型的MINQUE(minimum norm quadratic unbiased estimator)算法；文献[20]研究了PEIV模型的Helmert方差分量估计算法；文献[21]以PEIV模型为基础，推导了附有相对权比的TLS平差算法，并给出了两种相对权比的确定方法。上述各种算法并未考虑系数矩阵和观测向量的相关性，故不适用于协方差分量的估计，不能直接应用于AR模型的方差分量估计。文献[22]进一步分析了相关观测PEIV模型的LS-VCE算法；文献[23]基于非负最小二乘理论给出了PEIV模型的非负LS-VCE算法。文献[20-23]均是将系数矩阵和观测向量中的随机量作为两类数据在PEIV模型下进行方差分量估计，需要进行简单的多类观测数据扩展，否则也不能直接应用于AR模型的方差分量估计。

为此，本文提出了一种STLS问题的LS-VCE算法，并将其应用于AR模型的方差分量估计。该方法通过变量投影法^[7]将系数矩阵和观测向量构成的增广矩阵展开成仿射矩阵形式，将EIV模型在参数近似值和增广矩阵中随机量近似值处展开成非线性GH(Gauss-Helmert)模型，由此将TLS问题转化为LS问题，并采用附有系统参数的条件平差法进行迭代求解。同时根据LS-VCE原理对STLS问题的随机模型进行分析，计算随机量的方差和协方差分量估值，在迭代过程中更新观测权阵，达到修正参数估值的效果，并通过实测算例验证了本文算法的有效性。

1 STLS法 1.1 变量投影法

EIV模型考虑了系数矩阵和观测向量误差，其函数模型如下：

$ \mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}} - \left( {A - {\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}} \right) \cdot \mathit{\boldsymbol{x}} = 0 $

(1)

式中，$\mathit{\boldsymbol{y}} \in {{\bf{R}}^{m \times 1}}$表示观测向量；${\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}$表示观测向量的误差向量；$\mathit{\boldsymbol{A}} \in {{\bf{R}}^{m \times n}}$表示系数矩阵；${{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}}$表示系数矩阵的误差矩阵；$\mathit{\boldsymbol{x}} \in {{\bf{R}}^{n \times 1}}$表示参数向量。对于STLS问题，其系数矩阵和观测向量构成的增广矩阵$\mathit{\boldsymbol{C}} = \left[ {\mathit{\boldsymbol{A}}\;\;\mathit{\boldsymbol{y}}} \right]$具有某种结构，其误差矩阵${\rm{\Delta }}\mathit{\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}}}&{{\mathit{\boldsymbol{e}}_\mathit{\boldsymbol{y}}}} \end{array}} \right]$应具有相同结构。假设C仅存在某种仿射结构，即不考虑随机量之间存在非线性关系的情况，对于随机量之间存在非线性关系时，可以先将其进行线性化然后再求解。将C中独立随机量组成的向量设为$\mathit{\boldsymbol{p}} \in {{\bf{R}}^{{n_p}}}$，${{n_p}}$为独立随机量个数。故C的仿射函数形式如下^[7]：

$ C = \mathit{\boldsymbol{S}}\left( \mathit{\boldsymbol{p}} \right) = {\mathit{\boldsymbol{S}}_0} + \sum\limits_{j = 1}^{{n_p}} {{\mathit{\boldsymbol{S}}_j}} {p_j}\;\;\;\;\;\left( {{\mathit{\boldsymbol{S}}_0}, {\mathit{\boldsymbol{S}}_j} \in {{\rm{R}}^{m \times \left( {n + 1} \right)}}} \right) $

(2)

式中，${\mathit{\boldsymbol{S}}_0}$为固定量仿射结构矩阵，${\mathit{\boldsymbol{S}}_0}$中的非零元素为C中的固定量；${\mathit{\boldsymbol{S}}_j} = \left( {j = 1, 2 \cdots {n_p}} \right)$为随机量${p_j} = \left( {j = 1, 2 \cdots {n_p}} \right)$对应的仿射结构矩阵，它描述了${p_j}$在C中的分布。故可将STLS问题的数学模型变形为：

$ \left\{ {\begin{array}{*{20}{l}} {f\left( {{\rm{\Delta }}\mathit{\boldsymbol{p}}, \mathit{\boldsymbol{x}}} \right) = \mathit{\boldsymbol{S}}\left( {\mathit{\boldsymbol{p}} - {\rm{\Delta }}\mathit{\boldsymbol{p}}} \right){\mathit{\boldsymbol{x}}_{{\rm{ext}}}} = 0}\\ {{\rm{\Delta }}\mathit{\boldsymbol{p}} \sim \left( {0, {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}} \right)} \end{array}} \right. $

(3)

式中，${\mathit{\boldsymbol{x}}_{{\rm{ext}}}} = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{x}}\\ { - 1} \end{array}} \right]$为x的增广向量；${{\rm{\Delta }}\mathit{\boldsymbol{p}}}$为p的误差向量；Q_P为P的协方差阵。对式(3)中的函数模型进一步变形得：

$ \mathit{\boldsymbol{S}}\left( {\mathit{\boldsymbol{p}} - {\rm{\Delta }}\mathit{\boldsymbol{p}}} \right){\mathit{\boldsymbol{x}}_{{\rm{ext}}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{x}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}} = 0 $

(4)

其中，

$ \mathit{\boldsymbol{G}}\left( \mathit{\boldsymbol{x}} \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{S}}_1}{\mathit{\boldsymbol{x}}_{{\rm{ext}}}}}&{{\mathit{\boldsymbol{S}}_2}{\mathit{\boldsymbol{x}}_{{\rm{ext}}}}}& \ldots &{{\mathit{\boldsymbol{S}}_{{n_\mathit{\boldsymbol{p}}}}}{\mathit{\boldsymbol{x}}_{{\rm{ext}}}}} \end{array}} \right] = \left( {\mathit{\boldsymbol{x}}_{{\rm{ext}}}^{\rm{T}} \otimes {\mathit{\boldsymbol{I}}_m}} \right)\mathit{\boldsymbol{\psi }} $

(5)

式中，$\mathit{\boldsymbol{\psi }} = \left[ {\begin{array}{*{20}{c}} {{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_1}} \right)}&{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_2}} \right)}& \ldots &{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_{{n_\mathit{\boldsymbol{p}}}}}} \right)}\end{array}} \right]$；$ \otimes $为克罗内克乘积算子；I_m表示m × m的单位矩阵；vec(·)为矩阵向量化算子。

1.2 基于变量投影的STLS算法

EIV模型是一种非线性模型，故将式(4)在参数近似值${{{\hat x}^{i - 1}}}$和增广矩阵C中随机量误差近似值${\rm{\Delta }}{{\mathit{\boldsymbol{\hat p}}}^{i - 1}}$处进行一阶泰勒展开，得：

$ f\left( {{\rm{\Delta }}\mathit{\boldsymbol{p}}, \mathit{\boldsymbol{x}}} \right) \approx \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}} - {{\hat A}^{i - 1}}\mathit{\boldsymbol{x}} + \hat y_c^{i - 1} = 0 $

(6)

式中，${{\hat A}^{i - 1}} = \mathit{\boldsymbol{A}} - \hat E_\mathit{\boldsymbol{A}}^{i - 1}$ ($\hat E_\mathit{\boldsymbol{A}}^{i - 1}$的构造方法见式(2))；$\hat y_c^{i - 1} = \mathit{\boldsymbol{y}} - \hat E_\mathit{\boldsymbol{A}}^{i - 1}{{\hat x}^{i - 1}}$。式(6)可视为经典平差中附加系统参数的条件平差模型，可采用经典平差原理进行求解。令${\rm{\Delta }}\hat y_c^{i - 1} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}}$，故${\rm{\Delta }}\hat y_c^{i - 1} \sim \left( {0, \;\mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{i - 1}}} \right)} \right)$，将式(6)变为经典平差中的参数平差模型：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{\Delta }}\hat y_c^{i - 1} = {{\hat A}^{i - 1}}\mathit{\boldsymbol{x}} - \hat y_c^{i - 1}}\\ {{\rm{\Delta }}\hat y_c^{i - 1} \sim \left( {0, \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{i - 1}}} \right)} \right)} \end{array}} \right. $

(7)

根据参数平差法得：

$ {{\hat x}^i} = {[{{\hat A}^{i - 1}}^{\rm{T}}{({\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}})^{ - 1}}{{\hat A}^{i - 1}}]^{ - 1}}{{\hat A}^{i - 1}}^{\rm{T}}{({\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}})^{ - 1}}\hat y_c^{i - 1} $

(8)

其中，

$ {\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}} = {\mathit{\boldsymbol{Q}}_{{\rm{\Delta }}\hat y_c^{i - 1}}} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^{i - 1}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{i - 1}}} \right) $

(9)

同时，可得$\hat y_c^{i - 1}$的估值$\hat y_c^{i}$及其误差的估值$\;{\mathit{\boldsymbol{V}}_{\hat y_c^i}}$为：

$ \left\{ \begin{array}{l} \hat y_c^i = {\mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}}\hat y_c^{i - 1} = {{\hat A}^{i - 1}}{{\hat x}^i} = \left( {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{E}}_A^{i - 1}} \right){{\hat x}^i}\\ {\mathit{\boldsymbol{V}}_{\hat y_c^i}} = \mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}^ \bot \hat y_c^{i - 1} = \hat y_c^{i - 1} - \hat y_c^i = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{{\rm A}}}{{\hat x}^i} = \hat e_\mathit{\boldsymbol{y}}^{i - 1} - \mathit{\boldsymbol{E}}_\mathit{\boldsymbol{A}}^{i - 1}{{\hat x}^i} \end{array} \right. $

(10)

式中，${\mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}} = {{\hat A}^{i - 1}}{({{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}{{\hat A}^{i - 1}})^{ - 1}}{{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}$；$\mathit{\boldsymbol{P}}_{{{\hat A}^{i - 1}}}^ \bot = {\mathit{\boldsymbol{I}}_m} - {{\hat A}^{i - 1}}{({{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}{{\hat A}^{i - 1}})^{ - 1}}{{\hat A}^{i - 1}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{\hat y_c^{i - 1}}^{ - 1}$，均为正交投影矩阵。

已知x估值${{\hat x}^i}$后，则式(6)变为经典平差中的条件平差模型：

$ \left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{G}}\left( {{{\hat x}^i}} \right){\rm{\Delta }}\mathit{\boldsymbol{p}} + \left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\hat x}^i}} \right) = 0}\\ {{\rm{\Delta }}\mathit{\boldsymbol{p}} \sim \left( {0, {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}} \right)} \end{array}} \right. $

(11)

根据条件平差法得：

$ {\rm{\Delta }}{{\mathit{\boldsymbol{\hat p}}}^i} = - {\mathit{\boldsymbol{Q}}_p}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^i}} \right)\mathit{\boldsymbol{Q}}_{\hat y_c^i}^{ - 1}\left( {\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}{{\hat x}^i}} \right) $

(12)

式中，${\mathit{\boldsymbol{Q}}_{\hat y_c^i}} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^i}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^i}} \right)$。已知${\rm{\Delta }}\mathit{\boldsymbol{p}}$的估值${\rm{\Delta }}{{\mathit{\boldsymbol{\hat p}}}^i}$后，由式(2)得：

${\rm{\Delta }}{\hat C^i} = \left[ {\begin{array}{*{20}{c}} {\hat E_\mathit{\boldsymbol{A}}^i}&{\hat e_\mathit{\boldsymbol{y}}^i} \end{array}} \right] = \mathit{\boldsymbol{S}}\left( {{\rm{\Delta }}{{\hat p}^i}} \right) = \sum\limits_{j = 1}^{{n_\mathit{\boldsymbol{p}}}} {{\mathit{\boldsymbol{S}}_j}} {\rm{\Delta }}\hat p_j^i$

(13)

则${{\hat A}^i} = \mathit{\boldsymbol{A}} - \hat E_\mathit{\boldsymbol{A}}^i$。

根据最小二乘原理，上述STLS估计的单位权方差可以表示为：

${{\hat \sigma }^2} = {\rm{\Delta }}{{\hat p}^i}^{\rm{T}}\mathit{\boldsymbol{P}}{\rm{\Delta }}{{\hat p}^i}/r$

(14)

式中，P为权矩阵；r为自由度，即r = m - n。

对于参数估值的精度，由参数平差模型式(7)和参数估值迭代表达式(8)，根据协因数传播定律，可得：

${\mathit{\boldsymbol{Q}}_{\hat x\hat x}} = {[{{\hat A}^i}^{\rm{T}}{(\mathit{\boldsymbol{Q}}_{\hat y}^i)^{ - 1}}{{\hat A}^i}]^{ - 1}}$

(15)

上述推导过程仅考虑到了一阶误差的影响，对于更高阶误差所引起的参数估计以及精度评定偏差参见文献[24-26]。

2 STLS问题的最小二乘方差分量估计

LS-VCE由Teunissen等提出^[15-16]，该方法采用最小二乘准则进行方差分量估计，在线性GM模型假设下，可以得到方差分量的线性最优无偏估计。各种方差分量估计的目的均是对随机模型进行修正，区别在于所采用的估计原则以及对统计信息的依赖程度等方面的差异，并没有哪一种方差分量估计方法的估计特性较其他方法绝对占优。LS-VCE的优点在于，将方差分量估计问题纳入到最小二乘平差理论体系中，有利于方差分量估计的精度评定以及结合最小二乘理论进行其他方面的拓展研究。

式(3)中的随机模型可以分解为更一般的分量形式：

${\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Q}}_0} + \sum\limits_{k = 1}^d {{\sigma _k}{\mathit{\boldsymbol{Q}}_k}} $

(16)

式中，${{\sigma _k}}(k=1, 2…d)$表示随机量p的方差或协方差分量；d为方差分量的个数；Q₀为随机模型中的已知部分；Q_k$(k=1, 2…d)$为${{\sigma _k}}$对应的协因数矩阵且为线性无关的对称矩阵，并要求式(16)中的线性组合Q_p为正定矩阵。

在实际平差中，系数矩阵和观测向量中的随机量可能具有不同的方差分量，如果考虑到随机量之间的相关性，也可能具有不同的协方差分量，此时采用唯一的单位权方差进行平差并不合理。Teunissen提出的LS-VCE适用于GM模型，并可以进一步推广到EIV模型。将STLS问题转换为LS问题，在式(8)和式(12)的迭代计算过程中，Q_p为不断修正的协方差阵，可以通过方差分量估计确定未知的方差和协方差分量，以修正下一轮迭代中所给定的权，从而实现优化参数估值的目的。假设STLS迭代结束后得到的参数估值和随机量误差估值分别为${\hat x}$和${\rm{\Delta }}\hat p$，将它们代入式(7)得：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{\Delta }}{{\hat y}_c} = \hat A\mathit{\boldsymbol{x}} - {{\hat y}_c}}\\ {{\rm{\Delta }}{{\hat y}_c} \sim \left( {0, \sigma _0^2\mathit{\boldsymbol{G}}\left( {\hat x} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right)} \right)} \end{array}} \right. $

(17)

式中，$\hat A = \mathit{\boldsymbol{A}} - {{\hat E}_A}$；${\rm{\Delta }}{{\hat y}_c} = \mathit{\boldsymbol{G}}\left( {\hat x} \right){\rm{\Delta }}\mathit{\boldsymbol{p}}$；${{\hat y}_c} = \mathit{\boldsymbol{y}} - {{\hat E}_A}\hat x$。

式（17）可视为GM模型，其中${{\hat y}_c}$、${\rm{\Delta }}{{\hat y}_c}$、${\hat A}$分别为GM模型的观测向量、观测向量的误差向量和系数矩阵，因而可以直接将LS-VCE应用于式(17)中的平差问题。

将式(16)代入式(9)使${{\hat y}_c}$的协方差阵变形为：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_{_{{{\hat y}_c}}}} = {\mathit{\boldsymbol{Q}}_{{\rm{\Delta }}{{\hat y}_c}}} = \mathit{\boldsymbol{G}}\left( {\hat x} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right) = }\\ {\mathit{\boldsymbol{G}}\left( {\hat x} \right){\mathit{\boldsymbol{Q}}_0}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right) + \mathit{\boldsymbol{G}}\left( {\hat x} \right)\left( {\sum\limits_{k = 1}^d {{\sigma _k}{\mathit{\boldsymbol{Q}}_k}} } \right){\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {\hat x} \right) = }\\ {\mathit{\boldsymbol{Q}}_0^\mathit{\boldsymbol{G}} + \sum\limits_{k = 1}^d {{\sigma _k}} \mathit{\boldsymbol{Q}}_k^\mathit{\boldsymbol{G}}} \end{array} $

(18)

设矩阵${{\hat B}^{\rm{T}}}$为${\hat A}$的零空间，即${{\hat B}^{\rm{T}}}\hat A = 0$，${{\hat A}^{\rm{T}}}\hat B = 0$，同时引入向量t。由文献[15-16]可知，将${{\hat B}^{\rm{T}}}$左乘式(17)中的函数模型可得${{\hat B}^{\rm{T}}}$、t和${{\hat y}_c}$的关系如下：

$ \mathit{\boldsymbol{t}} = {{\hat B}^{\rm{T}}}{{\hat y}_c}, \;{\rm{E}}\left( \mathit{\boldsymbol{t}} \right) = 0 $

(19)

式中，E(·)表示数学期望。由协因数传播定律得：

$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{t}}} = {\rm{E}}\left( {\mathit{\boldsymbol{t}}{\mathit{\boldsymbol{t}}^{\rm{T}}}} \right) = {{\hat B}^{\rm{T}}}{\mathit{\boldsymbol{Q}}_{_{{{\hat y}_c}}}}\hat B $

(20)

式（20）即为方差分量估计的函数模型。为了消除式(20)中的重复观测方程，将vh(·)算子(由对称矩阵所对应的上三角矩阵或者下三角矩阵按列排列构成列向量)应用于上式，得：

$ {\rm{E}}\left( {{\mathit{\boldsymbol{y}}_{{\rm{vh}}}}} \right) = {\mathit{\boldsymbol{A}}_{{\rm{vh}}}}\mathit{\boldsymbol{\sigma }} $

(21)

式中，y_vh和A_vh为对式(20)进行vh(·)操作后的向量和矩阵；$\mathit{\boldsymbol{\sigma }}$为未知的方差和协方差分量组成的向量。由最小二乘准则并经过一系列推导和简化^[15-16]得到方差和协方差分量估值及其协因数阵如下：

$ \left\{ {\begin{array}{*{20}{l}} {\hat \sigma = {\mathit{\boldsymbol{N}}^{ - 1}}\mathit{\boldsymbol{l}}}\\ {{\mathit{\boldsymbol{Q}}_{\hat \sigma }} = {\mathit{\boldsymbol{N}}^{ - 1}}} \end{array}} \right. $

(22)

其中，N、l分别为计算方差分量的法方程的系数矩阵和常数向量，其元素分别为：

$ \left\{ \begin{array}{l} {n_{{j_1}{j_2}}} = \frac{1}{2}{\rm{tr}}\left( {\mathit{\boldsymbol{Q}}_{{j_1}}^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot \mathit{\boldsymbol{Q}}_{{j_2}}^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot } \right)\\ {l_{{j_1}}} = \frac{1}{2}\mathit{\boldsymbol{V}}_{{{\hat y}_c}}^{\rm{T}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{Q}}_{{j_1}}^{\bf{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}{\mathit{\boldsymbol{V}}_{{{\hat y}_c}}} - \frac{1}{2}{\rm{tr}}\left( {\mathit{\boldsymbol{Q}}_0^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot \mathit{\boldsymbol{Q}}_{{j_1}}^\mathit{\boldsymbol{G}}\mathit{\boldsymbol{Q}}_{{{\mathit{\boldsymbol{\hat y}}}_c}}^{ - 1}\mathit{\boldsymbol{P}}_{\hat A}^ \bot } \right) \end{array} \right. $

(23)

式中，$\mathit{\boldsymbol{P}}_{\hat A}^ \bot = {\mathit{\boldsymbol{I}}_m} - \hat A{({{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\hat A)^{ - 1}}{{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}$；$\mathit{\boldsymbol{P}}_{\hat A}^ \bot = {\mathit{\boldsymbol{I}}_m} - \hat A{({{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\hat A)^{ - 1}}{{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\;\;\;{\mathit{\boldsymbol{V}}_{{{\hat y}_c}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}\hat x$；${j_1}, {j_2} = 1, 2 \ldots d$。由以上推导可知，STLS问题的最小二乘方差分量估计法是由内外两个迭代循环过程构成，其中内循环用于计算方差和协方差分量的估值，外循环用于计算参数和随机误差的估值。具体的迭代步骤为：

1) 数据准备：A、y、Q_p、ψ，内循环收敛条件ε = 1 × 10^-10，外循环收敛条件τ = 1 × 10^-8。

2) 给定方差和协方差分量的初值$\mathit{\boldsymbol{\hat \sigma }} = {{\mathit{\boldsymbol{\hat \sigma }}}^0} = \left[ {\begin{array}{*{20}{c}} {\sigma _1^0}&{\sigma _2^0}& \ldots &{\sigma _t^0} \end{array}} \right]$。

3) 给定参数初值$\hat x = {{\hat x}^0} = {({\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{y}}$。

4) 设置外循环迭代次数s₁ = 0，并进入外循环。

5) 用式(16)更新${\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Q}}_0} + \sum\limits_{k = 1}^d {\hat \sigma _k^{{s_1}}} {\mathit{\boldsymbol{Q}}_k}$。

6) 设置内循环迭代次数s₂ = 0，并进入内循环：

(a) 用式(5)更新$\mathit{\boldsymbol{G}}\left( {{{\hat x}^{{s_2}}}} \right)$，并用式(9)更新${\mathit{\boldsymbol{Q}}_{\hat y_c^{{S_2}}}} = \mathit{\boldsymbol{G}}\left( {{{\hat x}^{{s_2}}}} \right){\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}}{\mathit{\boldsymbol{G}}^{\rm{T}}}\left( {{{\hat x}^{{s_2}}}} \right)$。

(b) 用式(12)计算得${\rm{\Delta }}{{\hat p}^{{s_2}}}$。

(c) 用式(13)计算得$\hat E_\mathit{\boldsymbol{A}}^{{s_2}}$，则有${{\hat A}^{{s_2}}} = A - \hat E_\mathit{\boldsymbol{A}}^{{s_2}}$。

(d) 用式(8)计算得${{\hat x}^{{s_2} + 1}}$。

(e) 重复内循环，直到$\;\left\| {{{\hat x}^{{s_2} + 1}} - {{\hat x}^{{s_2}}}} \right\| \le \varepsilon $。

7) 用式(18)更新${\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}}$。

8) 计算$\mathit{\boldsymbol{P}}_{\hat A}^ \bot = {\mathit{\boldsymbol{I}}_m} - \hat A{({{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}\hat A)^{ - 1}}{{\hat A}^{\rm{T}}}\mathit{\boldsymbol{Q}}_{{{\hat y}_c}}^{ - 1}$，${\mathit{\boldsymbol{V}}_{{{\hat y}_c}}} = \mathit{\boldsymbol{y}} - \mathit{\boldsymbol{A}}\hat x$。

9) 用式(23)计算得N和l。

10) 用式(22)计算得${{\hat \sigma }^{{s_1} + 1}} = {\mathit{\boldsymbol{N}}^{ - 1}}\mathit{\boldsymbol{l}}$，${\mathit{\boldsymbol{Q}}_{{{\hat \sigma }^{{S_1} + 1}}}} = {\mathit{\boldsymbol{N}}^{ - 1}}$。

11) 重复外循环，直到$\left\| {{{\mathit{\boldsymbol{\hat \sigma }}}^{{s_1} + 1}} - {{\mathit{\boldsymbol{\hat \sigma }}}^{{s_1}}}} \right\| \le \tau $。

3 算例分析

现以某一建筑物某个时段的36期沉降观测^[27]为例，高程观测值如表 1所示。文献[27]通过对参数的显著性进行检验，确定采用3阶AR模型对沉降数据进行分析和预测。为了便于比较，本文同样采用3阶AR模型，其对应的EIV模型如下：

$ \left\{ \begin{array}{l} H_4^{\rm{*}} = H_1^{\rm{*}}{\xi _1} + H_2^{\rm{*}}{\xi _2} + H_3^{\rm{*}}{\xi _3}\\ H_5^{\rm{*}} = H_2^{\rm{*}}{\xi _1} + H_3^{\rm{*}}{\xi _2} + H_4^{\rm{*}}{\xi _3}\\ H_6^{\rm{*}} = H_3^{\rm{*}}{\xi _1} + H_4^{\rm{*}}{\xi _2} + H_5^{\rm{*}}{\xi _3}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \vdots \\ H_{36}^{\rm{*}} = H_{33}^{\rm{*}}{\xi _1} + H_{34}^{\rm{*}}{\xi _2} + H_{35}^{\rm{*}}{\xi _3} \end{array} \right. $

(24)

式中，${\xi _1}$、${\xi _2}$、${\xi _3}$表示AR模型的模型参数；$H_i^{\rm{*}}\left( {i = 1, 2 \ldots 36} \right)$表示36期沉降观测的高程；*表示随机量。

利用本文算法进行解算时，还需要给出ψ。首先由系数矩阵和观测向量构成的增广矩阵C可表示为：

$ C = \left[ {\begin{array}{*{20}{c}} {H_1^{\rm{*}}}&{H_2^{\rm{*}}}&{H_3^{\rm{*}}}&{H_4^{\rm{*}}}\\ {H_2^{\rm{*}}}&{H_3^{\rm{*}}}&{H_4^{\rm{*}}}&{H_5^{\rm{*}}}\\ {H_3^{\rm{*}}}&{H_4^{\rm{*}}}&{H_5^{\rm{*}}}&{H_6^{\rm{*}}}\\ \vdots & \vdots & \vdots & \vdots \\ {H_{33}^{\rm{*}}}&{H_{34}^{\rm{*}}}&{H_{35}^{\rm{*}}}&{H_{36}^{\rm{*}}} \end{array}} \right] $

(25)

然后根据式(2)构造C的仿射结构矩阵S₁~S₃₀。不难发现，依次保留C各主对角线元素并置为1，即为S₁~S₃₀，进而可得$\mathit{\boldsymbol{\psi }} = \left[ {\begin{array}{*{20}{c}} {{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_1}} \right)}&{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_2}} \right)}& \ldots &{{\rm{vec}}\left( {{\mathit{\boldsymbol{S}}_{30}}} \right)} \end{array}} \right]$。

假设上述36组沉降数据来源于3个阶段的观测，其中第1阶段观测了12组数据，即${\mathit{\boldsymbol{H}}_f} = {\left[ {\begin{array}{*{20}{c}} {{H_1}}&{{H_2}}& \ldots &{{H_{12}}} \end{array}} \right]^{\rm{T}}}$；第2阶段也观测了12组数据，即${\mathit{\boldsymbol{H}}_s} = {\left[ {\begin{array}{*{20}{c}} {{H_{13}}}&{{H_{14}}}& \ldots &{{H_{24}}} \end{array}} \right]^{\rm{T}}}$；第3阶段同样观测了12组数据，即${\mathit{\boldsymbol{H}}_t} = {\left[ {\begin{array}{*{20}{c}} {{H_{25}}}&{{H_{26}}}& \ldots &{{H_{36}}} \end{array}} \right]^{\rm{T}}}$。参照式(16)，假定函数模型式(24)对应的随机模型如下：

$ {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{p}}} = {\mathit{\boldsymbol{Q}}_0} + \sigma _f^2{\mathit{\boldsymbol{Q}}_f} + \sigma _s^2{\mathit{\boldsymbol{Q}}_s} + \sigma _t^2{\mathit{\boldsymbol{Q}}_t} $

(26)

式中，$\sigma _f^2$、$\sigma _s^2$、$\sigma _t^2$分别为第1阶段、第2阶段、第3阶段沉降观测的方差，其对应的协因数矩阵分别为${\mathit{\boldsymbol{Q}}_f}$、${\mathit{\boldsymbol{Q}}_s}$、${\mathit{\boldsymbol{Q}}_t}$；${\mathit{\boldsymbol{Q}}_0}$表示随机模型的已知部分，其中，

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{Q}}_0} = {0_{36}}, \;{\mathit{\boldsymbol{Q}}_f} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{I}}_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}} \end{array}} \right], }\\ {{\mathit{\boldsymbol{Q}}_s} = \left[ {\begin{array}{*{20}{c}} {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{\mathit{\boldsymbol{I}}_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}} \end{array}} \right], {\mathit{\boldsymbol{Q}}_t} = \left[ {\begin{array}{*{20}{c}} {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{0_{12}}}\\ {{0_{12}}}&{{0_{12}}}&{{\mathit{\boldsymbol{I}}_{12}}} \end{array}} \right]} \end{array} $

(27)

式中，I₁₂表示12×12的单位矩阵；0₁₂表示12×12的零矩阵；0₃₆表示36×36的零矩阵。

文献[17-19]的算法并不适用于AR模型的方差分量估计，在此对文献[17]中的算法进行改进，称为Amiri-Simkooei改进算法，即考虑系数矩阵和观测向量之间的相关性，建立132×132的观测权阵Q，同时进一步修改文献[17]中${\mathit{\boldsymbol{Q}}_{\tilde y}}$的表达式为${\mathit{\boldsymbol{Q}}_{\tilde y}} = {{\hat B}^{\rm{T}}}\mathit{\boldsymbol{Q}}\hat B$(其中，$\hat B = \left[ {{{\hat x}^T} \otimes {\mathit{\boldsymbol{I}}_m} - {\mathit{\boldsymbol{I}}_m}} \right]$])^[17]，由于AR模型的系数矩阵和观测向量误差同源，观测权阵Q将变成非对角秩亏矩阵。

假设方差分量的迭代初值分别为$\sigma _f^2\left( 0 \right) = 1$，$\sigma _s^2\left( 0 \right) = 1$，$\sigma _t^2\left( 0 \right) = 1$，分别采用最小二乘估计(LS估计)、PEIV模型解法^[10]、Amiri-Simkooei改进算法^[17]、本文算法进行求解，参数估计的迭代收敛阈值设为ε = 1 × 10^-10，方差分量估计的迭代收敛阈值设为τ = 1 × 10^-8。

表 2列出了各种算法得到的参数估值及其精度（标准差），表 3列出了Amiri-Simkooei改进算法和本文算法得到的方差分量估值及其精度，表 4列出了LS估计、PEIV模型解法得到的单位权方差估值及其精度，图 1给出了本文算法的方差分量收敛曲线。

由表 2可知，LS估计和PEIV模型解法所得的参数估值差异明显，在假定36期观测等精度的前提下，PEIV模型解法给出的模型估值更加可靠。在引入3个方差分量σ_f²、σ_s²、σ_t²后，Amiri-Simkooei改进算法和本文算法所得的参数估值产生了进一步的改变，这说明观测权阵对参数估计具有显著影响，合理的观测权阵是参数估计可靠性的前提。同时，本文算法和Amiri-Simkooei改进算法的计算结果一致，说明了本文算法的有效性以及WTLS估计^[12]和本文给出的STLS估计具有等价性，但本文算法显著减少了待估参数个数，也无需重新构造秩亏的观测权阵，因而在一定程度上可以提高算法的计算效率。通过表 3和表 4可知，本例引入的3个方差分量σ_f²、σ_s²、σ_t²估值存在显著差异，并且与PEIV模型解法得到的单位权方差估值也存在差异，因而可以得出结论：文献[27]对沉降数据进行处理时做出的等精度观测假设并不合适，由于对沉降观测没有进一步的了解，因此无法给出合理的随机模型对表 1中的沉降数据进行可靠的分析。

需要注意的是，本文给出的AR模型方差分量估计方法最终归结为一个非凸优化问题，所以并不能保证得到其全局最优解，其全局收敛性与模型的非线性程度、观测方程的不吻合程度等相关，并且在算法设计上采用了高斯法，因而收敛性也与初值相关^[28]。

4 结语

AR模型是一类重要的时间序列模型，它在测量数据处理中有广泛应用，AR模型的求解本质上是一个STLS问题。针对AR模型平差时观测权阵分配不合理、随机模型不准确的情况，本文首先基于变量投影法给出了一种STLS估计算法，然后将其与LS-VCE相结合，推导出了STLS问题的一种方差分量估计算法。该算法可以直接应用于AR模型的方差分量估值，并通过实测算例验证了它的有效性。结论如下：

1）本文给出的STLS估计算法可以认为是文献[7]的一种替代算法，该算法在一定程度上简化了推导过程，提高了计算效率，并给出了参数的一阶近似精度。

2）本文将STLS问题转化为附加参数的条件平差问题进行求解，即将TLS问题纳入到最小二乘平差理论体系中，因而能够直接和现有的方差分量估计算法相结合，如最小二乘方差分量估计。本文算法适用于多类数据的方差分量估计，并且也可以进行协方差分量的估计。

3）本文只讨论了EIV模型系数矩阵状态良好、观测量只存在随机误差时的方差分量估计问题，EIV模型下方差分量的稳健估计算法和正则化估计算法是下一步的研究重点。