在直线、平面拟合、空间直角坐标变换、自回归模型求解中,系数矩阵和右端观测量由随机和非随机元素组成,且同一随机元素会在不同的位置出现^[1-4]。常用的变量含误差（errors-in-variables,EIV）模型及相应的整体最小二乘（total least squares,TLS）方法假设系数阵中所有元素含误差。因此,顾及增广误差矩阵元素的随机特性及其相互关系的结构EIV（structured EIV,SEIV）模型和结构加权整体最小二乘（structured weighted TLS,SWTLS）方法应运而生^[5-10]。文献[5]采用极大似然估计方法从含有误差的时间序列数据中确定动态线性系统的参数,被认为是SWTLS问题研究的开端。文献[6]定义了结构整体最小二乘（structured TLS,STLS）这一术语,并将其转换为非线性广义奇异值分解问题求解。文献[7]将结构增广误差矩阵用其中的独立随机元素表示,并命名为约束TLS方法。对于某些列为固定元素而另外一些列为随机元素的混合LS-TLS问题,文献[8]引入正交三角分解计算模型参数。对于误差为非等权相关的情形,文献[9]发展了WLS-WTLS的迭代算法。文献[10]将SEIV模型用部分EIV（partial EIV,PEIV）模型表达,在加权最小二乘准则下推导了参数估计的Gauss-Newton型迭代算法,并进一步给出了参数的一阶近似方差、非线性信赖域区间及偏差。文献[11]给出了PEIV模型的两种迭代算法,其计算效率与系数矩阵中随机元素的数量有关。文献[12]提出了顾及系数阵元素和右端向量相关的广义PEIV模型。文献[13]将系数误差矩阵表示成已知矩阵和独立误差向量的乘积,推导了附有线性和二次约束的STLS问题的迭代算法。文献[14]和文献[15]采用变量投影法将增广系数矩阵表示成仿射结构矩阵与独立随机变量的乘积,然后将STLS问题转换为非线性等式约束优化问题求解。文献[16]研究了含多个右端观测向量的STLS问题。文献[17]提出了结构整体最小范数（structured total least norm,STLN）方法,这一方法可以最小化误差向量的范数,是STLS问题的重要拓展^[17-18]。文献[19-20]研究了STLN问题的快速算法。文献[21]将STLN拓展到系数阵和观测向量具有共同元素的情形。针对自回归（auto-regression, AR）模型这一特定STLS问题,文献[3-4]提出了两种新解法。

在STLS算法推导中,不同文献使用了不同的平差准则,其差异在于是否考虑独立误差的重复次数。如二维坐标转换中,源系统的纵横坐标值在系数矩阵中出现了2次,而目标系统的纵横坐标值在右端向量中只出现了1次。文献[1,7,10,13-15]给出的平差准则未考虑重复次数,而文献[3-4,17,21]均考虑了重复次数。不同的平差准则必定会得到不同的平差结果,其参数估值在统计意义上孰优孰劣,目前尚未给出明确结果。本文从函数模型和数值模拟两方面入手,证明了不考虑独立误差元素重复次数的平差准则能够得到统计意义下更优的解。

1 STLS平差模型及平差准则

EIV模型的函数表达式为^[1]：

式中,$\mathit{y}$和${\mathit{e}}_{\mathit{y}}$分别表示$n$维观测向量及其误差；$\mathit{A}$和${\mathit{E}}_{\mathit{A}}$分别表示$n\times m$维系数矩阵及其误差矩阵；$\mathit{x}$为$m$维参数向量。误差向量的随机模型为：

式中,${\mathit{e}}_{\mathit{A}}=\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{E}}_{\mathit{A}}\right)$,$\mathrm{v}\mathrm{e}\mathrm{c}\left(\bullet \right)$表示矩阵向量化算子,即将$n\times m$维矩阵的每一列从左至右叠加成一个$nm$维列向量；${\sigma }_{\mathrm{0}}^{\mathrm{2}}$是单位权方差；${\mathit{Q}}_{\mathit{A}\mathit{A}}$和${\mathit{Q}}_{\mathit{y}\mathit{y}}$分别是${\mathit{e}}_{\mathit{A}}$和${\mathit{e}}_{\mathit{y}}$的对称正定协因数矩阵；${\mathit{Q}}_{\mathit{y}\mathit{A}}={\mathit{Q}}_{\mathit{A}\mathit{y}}^{\mathrm{T}}$表示两者的协因数矩阵。当向量e的协因数矩阵${\mathit{Q}}_{\mathit{e}}$可逆时,采用WTLS准则求参数的最优估值：

式中,${\mathit{P}}_{\mathit{e}}={\mathit{Q}}_{\mathit{e}}^{-\mathrm{1}}$为$\mathit{e}$的权矩阵。当系数矩阵具有某种结构时,${\mathit{Q}}_{\mathit{e}}$为秩亏矩阵无凯利逆。许多学者将EIV函数模型进行改化,一般是从结构误差矩阵或增广误差矩阵中提取独立误差向量,进一步构造目标函数求解。常用的几种EIV模型修正方法如下：

1） PEIV模型及其平差准则。PEIV模型选取系数阵中独立随机量的真值$\overline{\mathit{a}}$作为待求量,其函数模型为^[10]：

式中,$\mathit{a}$是系数矩阵中随机元素所构成的$t$维列向量；$\overline{\mathit{a}}$和${\mathit{e}}_{\mathit{a}}$分别是相应的真值和误差向量；${\mathit{I}}_n$表示$n$维单位矩阵；$\mathit{h}$是已知的$nm$维常数向量,其元素包含系数阵中的非随机元素；$\mathit{M}$是$nm\times t$维常数矩阵,其形式由系数矩阵中非随机元素的个数及元素间的相关性确定。若$\mathit{a}$和$\mathit{y}$相互独立^[10],即$\mathrm{c}\mathrm{o}\mathrm{v}({\mathit{e}}_{\mathit{a}},{\mathit{e}}_{\mathit{y}})=\mathrm{0}$,且两者的方差矩阵为$\mathit{D}({\mathit{e}}_{\mathit{a}})={\sigma }_{\mathrm{0}}^{\mathrm{2}}{\mathit{\omega }}^{-\mathrm{1}}$,$\mathit{D}({\mathit{e}}_{\mathit{y}})={\sigma }_{\mathrm{0}}^{\mathrm{2}}{\mathit{W}}^{-\mathrm{1}}$,其中$\mathit{\omega }$和$\mathit{W}$分别为${\mathit{e}}_{\mathit{a}}$和${\mathit{e}}_{\mathit{y}}$的权矩阵。相应的平差准则为：

可见PEIV的平差准则中没有考虑独立随机误差的重复次数。如果将${\mathit{e}}_{\mathit{a}}$和${\mathit{e}}_{\mathit{y}}$组合成独立随机误差向量$\mathit{\gamma }={\left[\begin{array}{cc}{\mathit{e}}_{\mathit{a}}^{\mathrm{T}}& {\mathit{e}}_{\mathit{y}}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$,权矩阵${\mathit{P}}_{\gamma }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathit{\omega },\mathit{W})=(\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}({\mathit{\omega }}^{-\mathrm{1}},{\mathit{W}}^{-\mathrm{1}}{))}^{-\mathrm{1}}$,则平差准则式（4）可以写为：

因此PEIV模型的SWTLS解是令独立随机误差的加权平方和最小的参数估值。文献[12]中的GPEIV模型、文献[13]中约束结构SWTLS方法和文献[14-15]中的变量投影法,其本质上都是采用式（5）所示的平差准则。

2） STLN模型及其平差准则。STLN方法是定义一个与参数$\mathit{x}$相关的矩阵$\mathit{X}$,提取系数矩阵中的独立随机误差${\mathit{e}}_{\mathit{a}}$,使下式成立^[17]：

联立式（1a）和式（6）,将观测值残差$\mathit{r}$表达成系数阵独立误差${\mathit{e}}_{\mathit{a}}$和参数$\mathit{x}$的函数：

将式（7）在${\mathit{e}}_{\mathit{a}}$和$\mathit{x}$的近似值处线性化,舍去二次项后得到线性模型,然后采用下列平差准则^[17]：

式中,${\mathit{D}}_a=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(d_{\mathrm{1}},d_{\mathrm{2}},\cdots ,d_t\right)$,$d_i$（$i=\mathrm{1,2},\cdots ,t$）表示${\mathit{e}}_{\mathit{a}}$中的第$i$个元素在误差矩阵${\mathit{E}}_{\mathit{A}}$中重复出现的次数。尽管STLN 方法没有考虑误差的权值,但可以很方便地将观测误差的权纳入平差准则。仍令$\mathit{\gamma }={\left[\begin{array}{cc}{\mathit{e}}_a^{\mathrm{T}}& {\mathit{e}}_y^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}$,定义${\mathit{D}}_s=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\mathit{D}}_a,{\mathit{I}}_n\right)$,由于式（8）中常数项对求极值无影响,则平差准则式（8）等价于：

可见,STLN方法考虑了独立误差${\mathit{e}}_{\mathit{a}}$的重复次数,且以重复次数的平方将独立误差纳入平差模型。

3） 虚拟误差模型及其平差准则。文献[3]以AR模型参数估计为背景,将式（1a）所示的EIV模型在观测值的近似值处线性化。设系数阵真值$\tilde{\mathit{A}}$的近似值为${\mathit{A}}_{\mathrm{0}}$,改正数为$∆\mathit{A}$,参数$\mathit{x}$的近似值为${\mathit{x}}_{\mathrm{0}}$,改正数为$∆\mathit{x}$,则线性化方程为：

式中,$\mathit{v}$表示$\mathit{y}$对应的改正数。然后通过矩阵等价变换得到$∆\mathit{A}{\mathit{x}}_{\mathrm{0}}=\left[\begin{array}{cc}{\mathit{A}}_{\mathrm{10}}& {\mathit{A}}_{\mathrm{20}}\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_{\mathit{B}}\\ \mathit{v}\end{array}\right]$,${\mathit{A}}_{\mathrm{10}}$和${\mathit{A}}_{\mathrm{20}}$表示由${\mathit{x}}_{\mathrm{0}}$按照一定规则构造的近似矩阵；${\mathit{x}}_B$表示设计矩阵中独立观测值对应的改正数向量。将式（10）化为：

组成虚拟观测值误差方程为：

联立式（11）和式（12）得到：

式中,${\mathit{v}}_g$为所有随机观测值改正数向量；${\mathit{B}}_g$为等效设计矩阵；${\mathit{x}}_g$为所有待估参数和虚拟参数的改正数；${\mathit{l}}_g$表示$\mathit{y}-{\mathit{A}}_{\mathrm{0}}{\mathit{x}}_{\mathrm{0}}$的线性变换^[3]。采用如下平差准则：

式中,${\mathit{v}}_g$实际上是$\left({\mathit{E}}_{\mathit{A}}|{\mathit{e}}_{\mathit{y}}\right)$中的独立误差向量$\mathit{\gamma }$；P_g=D_gP由两部分组成,其中${\mathit{D}}_g=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(d_{\mathrm{1}},d_{\mathrm{2}},\cdots ,d_{t+n}\right)$,$d_i$（$i=\mathrm{1,2},\cdots ,t+n$,且t+n为独立误差的个数）表示${\mathit{v}}_g$中的第$i$个元素在$({\mathit{E}}_{\mathit{A}}|{\mathit{e}}_{\mathit{y}})$中的重复次数,$\mathit{P}$为观测值向量的权矩阵。可见,虚拟误差模型解AR模型这一典型SWTLS问题考虑了重复次数。与STLN方法的差异在于,它采用的是重复次数本身,而不是STLN中重复次数的平方。

上述3种典型的STLS方法采用了不同的平差准则,不同的平差准则一定会得到不同的平差结果。为了比较不同准则下的平差结果,下面将上述3种典型的平差准则下的解纳入到一个统一的平差模型中。

2 STLS模型的通用表达式及算法推导

EIV模型（式（1a））又可以表示为^[1]：

式中,$\mathit{B}=\left[\begin{array}{cc}{\mathit{x}}^{\mathrm{T}}\otimes {\mathit{I}}_n& -{\mathit{I}}_n\end{array}\right]$为$n\times \left(nm+n\right)$矩阵,其中$\otimes$为Kronecker积符号,定义为$\mathit{G}\otimes \mathit{H}=\left[g_{ij}\bullet \mathit{H}\right]$,其中$\mathit{G}=\left[g_{ij}\right]$,且$\mathit{H}$为任意矩阵。设$\mathit{\gamma }$为$\left({\mathit{E}}_{\mathit{A}}|{\mathit{e}}_{\mathit{y}}\right)$中的$t$个独立随机误差向量,则${\mathit{e}}_{\mathit{A}}$和${\mathit{e}}_{\mathit{y}}$可以分别表示为：

式中,${\mathit{H}}_{\mathrm{1}}$和${\mathit{H}}_{\mathrm{2}}$分别为$nm\times t$和$n\times t$矩阵,其元素均为常数。联立式（16a）和式（16b）：

将式（17）代入式（15）,并令$\mathit{G}=\mathit{B}\mathit{H}$,则EIV模型可以写为：

独立随机误差的随机模型为：

式中,${\mathit{Q}}_{\mathit{\gamma }}$表示$\mathit{\gamma }$的协因数矩阵；${\mathit{P}}_{\mathit{\gamma }}$表示其权矩阵。令$\mathit{\gamma }={\left[{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}\cdots {\gamma }_t\right]}^{\mathrm{T}}$,且${\gamma }_i$出现的次数为$d_i$（$i=\mathrm{1,2},\cdots ,t$）。定义$\mathit{D}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(k_{\mathrm{1}},k_{\mathrm{2}},\cdots ,k_t\right)$,组成综合权矩阵：

相应的平差准则为：

易知,当$k_i=\mathrm{1}$时,式（21）等价于以PEIV模型为代表的准则式（5）；当$k_i=d_i^{\mathrm{2}}$时,式（21）等价于以STLN模型为代表的准则式（9）；当$k_i=d_i$时,式（21）等价于以虚拟误差模型为代表的准则式（14）。以结构化函数模型（18）为基础,基于平差准则式（21）,构造如下Lagrange目标函数：

式中,$\mathit{\lambda }$为$n$维Lagrange乘子向量。令目标函数式（22）对各待定量$\mathit{x}\mathrm{、}\mathit{\gamma }\mathrm{、}\mathit{\lambda }$的偏导数为0,可得：

由式（23b）得到独立误差向量$\mathit{\gamma }$的估值为：

式中,${\mathit{Q}}_{s\mathit{\gamma }}={\mathit{P}}_{s\mathit{\gamma }}^{-\mathrm{1}}$,表示综合权矩阵的逆；$\widehat{\mathit{x}}\mathrm{、}\widehat{\mathit{\gamma }}\mathrm{、}\widehat{\mathit{\lambda }}$分别为$\mathit{x}\mathrm{、}\mathit{\gamma }\mathrm{、}\mathit{\lambda }$的估值。将式（24）代入式（23c）并移项可得：

将式（25）代入式（23a）可以得到法方程：

则参数的估值为：

由于式（26）法方程矩阵不对称,若在式（26）两端加上${\left(\mathit{A}+{\widehat{\mathit{E}}}_{\mathit{A}}\right)}^{\mathrm{T}}{\left(\mathit{G}{\mathit{Q}}_{s\mathit{\gamma }}{\mathit{G}}^{\mathrm{T}}\right)}^{-\mathrm{1}}{\widehat{\mathit{E}}}_{\mathit{A}}\widehat{\mathit{x}}$,那么$\widehat{\mathit{x}}$可由下式估计：

根据上述推导过程,可以得出通用模型计算SWTLS问题的迭代步骤为：

1） 给定$\mathit{A}$、$\mathit{y}$、${\mathit{Q}}_{\mathit{\gamma }}$,根据$\mathit{\gamma }$的重复次数和不同平差准则的计入方式确定对角阵$\mathit{D}$,根据$({\mathit{E}}_{\mathit{A}}|{\mathit{e}}_{\mathit{y}})$的结构确定${\mathit{H}}_{\mathrm{1}}$和${\mathit{H}}_{\mathrm{2}}$。采用式（20）计算综合权阵${\mathit{P}}_{s\mathit{\gamma }}$和对应的协因数矩阵${\mathit{Q}}_{s\mathit{\gamma }}={\mathit{P}}_{s\mathit{\gamma }}^{-\mathrm{1}}$。计算${\mathit{Q}}_{\mathit{y}\mathit{y}}={\mathit{H}}_{\mathrm{2}}{\mathit{Q}}_{\mathit{\gamma }}{\mathit{H}}_{\mathrm{2}}^{\mathrm{T}}$和参数的初值${\widehat{\mathit{x}}}^{\left(\mathrm{0}\right)}={\widehat{\mathit{x}}}_{\mathrm{L}\mathrm{S}}={\left({\mathit{A}}^{\mathrm{T}}{\mathit{Q}}_{\mathit{y}\mathit{y}}^{-\mathrm{1}}\mathit{A}\right)}^{-\mathrm{1}}{\mathit{A}}^{\mathrm{T}}{\mathit{Q}}_{\mathit{y}\mathit{y}}^{-\mathrm{1}}\mathit{y}$；

2） 根据初值${\widehat{\mathit{x}}}^{\left(\mathrm{0}\right)}$计算矩阵$\mathit{B}$和$\mathit{G}$。采用式（25）和式（24）分别计算$\widehat{\mathit{\lambda }}$和$\widehat{\mathit{\gamma }}$,然后由式（16a）计算${\widehat{\mathit{e}}}_{\mathit{A}}$的值,从而有${\widehat{\mathit{E}}}_{\mathit{A}}=\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{c}\left({\widehat{\mathit{e}}}_{\mathit{A}}\right)$,$\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{c}\left(\bullet \right)$是$\mathrm{v}\mathrm{e}\mathrm{c}\left(\bullet \right)$的逆运算,表示将$nm\mathrm{维}$向量恢复成$n\times m$维矩阵；

3） 将$\mathit{G}$、${\mathit{Q}}_{s\mathit{\gamma }}$、${\widehat{\mathit{E}}}_{\mathit{A}}$、${\widehat{\mathit{x}}}^{\left(\mathrm{0}\right)}$分别代入式（27）或式（28）,计算参数估值$\widehat{\mathit{x}}$；

4） 如果最后两次的估值$\widehat{\mathit{x}}$足够接近于给定的阈值,终止迭代。否则,转向步骤2）。

3 STLS模型平差准则的优化选取

尽管由式（18）和式（19）组成的SWTLS一般模型能通过选取不同的$\mathit{D}$矩阵获得不同准则下的迭代最优解,但无法获得参数的方差或均方误差（mean square error,MSE）等精度评定指标,只能在模拟实验中检验参数精度。由文献[22-23]可知,WTLS并不是一种新的平差方法,仅仅是LS框架下的另外一种平差模型,且经典平差理论中的Gauss-Helmert（GH）模型（附有参数的条件平差模型）是EIV模型的一个特例。因此WTLS问题的解可由非线性GH模型导出。若给定$\mathit{x}$的近似值${\mathit{x}}_{\mathrm{0}}$和$\mathit{\gamma }$的近似值${\mathit{\gamma }}_{\mathrm{0}}$,式（18）可线性化为：

式中,$\mathrm{\Delta }\mathit{x}=\mathit{x}-{\mathit{x}}_{\mathrm{0}}$；$\mathrm{\Delta }\mathit{\gamma }=\mathit{\gamma }-{\mathit{\gamma }}_{\mathrm{0}}$；${\mathit{G}}_{\mathrm{0}}$和${\mathit{E}}_{\mathit{A}}^{\mathrm{0}}$为近似值${\mathit{x}}_{\mathrm{0}}$和${\mathit{\gamma }}_{\mathrm{0}}$处的$\mathit{G}$和${\mathit{E}}_{\mathit{A}}$的值。且有：

将式（30）和式（31）代入式（29）并作适当变换可以得到：

可见SEIV模型的线性化形式（式（32））是一个标准GH模型,该模型不再体现$\mathit{\gamma }$的重复次数,实际上重复次数已经由投影矩阵$\mathit{H}$表达了。重复次数体现了EIV模型的结构性,而投影矩阵$\mathit{H}$正是描述上述结构性的量。根据经典LS理论,其平差准则应为：

因此,从SWTLS问题的线性化形式证明了其平差准则应满足式（33）。由线性化模型式（32）得到的最小二乘解是在给定${\mathit{x}}_{\mathrm{0}}$和${\mathit{\gamma }}_{\mathrm{0}}$情形下的一个近似解,为了得到更严密的解,将式（32）进一步表达为：

式中,$\mathit{l}=\mathit{y}-\mathit{A}{\mathit{x}}_{\mathrm{0}}$。为了与§2中符号一致,记${\mathit{P}}_{s\mathit{\gamma }}=I_{\mathrm{t}}{\mathit{P}}_{\mathit{\gamma }}={\mathit{P}}_{\mathit{\gamma }}$（$I_{\mathrm{t}}$为t维单位矩阵）。根据式（34）及平差准则式（33）,构造如下Lagrange乘子函数：

分别求$\phi \left(∆\mathit{x},∆\mathit{\gamma }\right)$对$∆\mathit{x}$、$∆\mathit{\gamma }$和$\mathit{\lambda }$的偏导数,并令其值为0,可得：

由式（37）可得：

将式（39）代入式（38）得：

由式（40）可得Lagrange乘子向量为：

式（40）两边同乘以${\left(\mathit{A}+{\mathit{E}}_{\mathit{A}}^{\mathrm{0}}\right)}^{\mathrm{T}}$并顾及式（36）,则有：

由GH模型推导得到的参数估计（式（43））与§2中通用模型得到的解（式（28））具有完全相同的形式。GH模型的算法流程如下：

1） 给定$\mathit{A}$、$\mathit{y}$、${\mathit{Q}}_{\gamma }$,根据$({\mathit{E}}_{\mathit{A}}|{\mathit{e}}_{\mathit{y}})$的结构确定${\mathit{H}}_{\mathrm{1}}$、${\mathit{H}}_{\mathrm{2}}$。计算${\mathit{Q}}_{\mathit{y}\mathit{y}}={\mathit{H}}_{\mathrm{2}}{\mathit{Q}}_{\mathit{\gamma }}{\mathit{H}}_{\mathrm{2}}^{\mathrm{T}}$和参数的初值${\widehat{\mathit{x}}}^{\left(\mathrm{0}\right)}={\widehat{\mathit{x}}}_{\mathrm{L}\mathrm{S}}={\left({\mathit{A}}^{\mathrm{T}}{\mathit{Q}}_{\mathit{y}\mathit{y}}^{-\mathrm{1}}\mathit{A}\right)}^{-\mathrm{1}}{\mathit{A}}^{\mathrm{T}}{\mathit{Q}}_{\mathit{y}\mathit{y}}^{-\mathrm{1}}\mathit{y}$,使用上标（i）$(i=\mathrm{0,1},\mathrm{2},\cdots ,N,N\mathrm{为迭}\mathrm{代次}\mathrm{数})$对不同迭代次数下的各参数进行区分,则$\mathit{\gamma }$和${\mathit{E}}_{\mathit{A}}$的初值分别为${\mathit{\gamma }}^{\left(\mathrm{0}\right)}=\mathrm{0}$和${\mathit{E}}_{\mathit{A}}^{\left(\mathrm{0}\right)}=\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{c}\left({\mathit{H}}_{\mathrm{1}}{\mathit{\gamma }}^{\left(\mathrm{0}\right)}\right)=\mathrm{0}$；

2） 将${\widehat{\mathit{x}}}^{\left(i\right)}$代入式（31）计算${\mathit{B}}^{\left(i\right)}$和${\mathit{G}}^{\left(i\right)}$,由式（41）计算${\mathit{\lambda }}^{\left(i\right)}$,由式（39）计算$∆{\mathit{\gamma }}^{\left(i\right)}$并更新${\mathit{\gamma }}^{\left(i\right)}={\mathit{\gamma }}^{\left(i-\mathrm{1}\right)}+∆{\mathit{\gamma }}^{\left(i\right)}$,由式（30）计算${\mathit{E}}_{\mathit{A}}^{\left(i\right)}$；

3） 由式（42）计算$\mathrm{\Delta }{\mathit{x}}^{\left(i-\mathrm{1}\right)}$；

4） 如果$\mathrm{\Delta }\mathit{x}$小于给定的阈值,则终止迭代；否则,更新${\mathit{x}}^{\left(i\right)}=\mathrm{\Delta }{\mathit{x}}^{\left(i-\mathrm{1}\right)}+{\mathit{x}}^{\left(i-\mathrm{1}\right)}$,转向步骤2）。

注意,若不限定式（35）中${\mathit{P}}_{s\mathit{\gamma }}={\mathit{P}}_{\mathit{\gamma }}$,那么采用线性化的GH模型,在准则式（21）下,即可选取不同类型$\mathit{D}$矩阵的条件下,可以得到和§2中通用模型完全一致的解。采用GH模型进行算法推导有3个方面的意义：（1）证明了SEIV模型的平差准则中不应该再考虑重复矩阵。（2）证明了由通用模型（式（18））和线性化GH模型作为条件,在相同的平差准则下得到的参数估值是一致的。（3）将SEIV模型线性化后,便于采用协方差传播定律计算参数估值的近似精度。

4 STLS模型参数估值的近似方差

根据§2、§3的算法得到SWTLS解以后,可计算独立误差向量的估值$\widehat{\mathit{\gamma }}$,单位权方差计算式为：

SWTLS解的近似协因数矩阵可以根据线性化GH模型的解（式（42））求解。虚拟观测值$\mathit{l}$可以写成：

式中,$\mathit{\eta }$为$(\mathit{A}|\mathit{y})$中的$t$个独立随机观测向量。根据协因数传播律,向量$\mathit{l}$的协因数矩阵为：

实际计算中近似值${\mathit{x}}_{\mathrm{0}}$可以用SWTLS估值$\widehat{\mathit{x}}$代替,相应的残差矩阵近似值${\mathit{E}}_{\mathit{A}}^{\mathrm{0}}$用对应的估值${\tilde{\mathit{E}}}_{\mathit{A}}$代替。如果忽略式（42）中$\tilde{\mathit{A}}=\mathit{A}+{\mathit{E}}_{\mathit{A}}^{\mathrm{0}}$这一项中$\mathit{A}$的随机性,那么可由协因数传播律得到SWTLS解的近似协因数矩阵为：

从而得到SWTLS解的近似方差为：

5 数值实验

首先采用二维仿射变换实例验证算法的正确性和可行性。设第$i$个公共点$P_i$在源坐标系和目标坐标系中的坐标分别为$\left(x_i,y_i\right)$和$\left(X_i,Y_i\right)$,仿射变换模型为$X_i=a_{\mathrm{0}}+a_{\mathrm{1}}{x}_i+a_{\mathrm{2}}{y}_i$,$Y_i=b_{\mathrm{0}}+b_{\mathrm{1}}{x}_i+b_{\mathrm{2}}{y}_i$,其中$a_{\mathrm{0}}$和$b_{\mathrm{0}}$为两个坐标系的平移量,$a_{\mathrm{1}}$、$a_{\mathrm{2}}$、$b_{\mathrm{1}}$、$b_{\mathrm{2}}$分别表示两个坐标系之间的旋转和尺度变换参数。假设有3个或以上的公共点,仿射变换模型可以用EIV模型（式（1））表示,且有：

给定12个公共点在源坐标系中的坐标分别为（-1,1）、（1.2,-3.0）、（-2.6,3.0）、（3.0,1.5）、（-4.8,-1.0）、（5.2,0.2）、（6.0,5.5）、（-7.2,2.2）、（7.8,-2.0）、（8.5,2.2）、（-9.5,-5.0）、（10.0,-0.8）,参数的真值为$\tilde{\mathit{x}}$=[10 4 -2 -10 1 3]^T,根据仿射变换模型计算出12个点在目标系中的真实坐标分别为（4,-8）、（20.8,-17.8）、（-6.4,-3.6）、（19.0,-2.5）、（-7.2,-17.8）、（30.4,-4.2）、（23.0,12.5）、（-23.2,-10.6）、（45.2,-8.2）、（39.6,5.1）、（-18.0,-34.5）、（51.6,-2.4）。对所有坐标值添加方差为${\mathit{D}}_{\mathrm{48}}={\sigma }_{\mathrm{0}}^{\mathrm{2}}{\mathit{I}}_{\mathrm{48}}={\mathit{I}}_{\mathrm{48}}$（${\mathit{I}}_{\mathrm{48}}$为$\mathrm{48}\times \mathrm{48}$的单位矩阵）的随机误差。首先,不考虑系数矩阵中源坐标重复2次而目标坐标重复1次的差异,取重复矩阵$\mathit{D}={\mathit{D}}_{\mathrm{1}}={\mathit{I}}_{\mathrm{48}}$时,分别采用本文方法、PEIV模型^[10]和变量投影法^[14]进行求解,得到的参数估计结果如表1所示。

表  1  不同平差方法的SWTLS解（$\mathit{D}={\mathit{D}}_{\mathrm{1}}$）

Table  1.  SWTLS Solutions with Different Adjustment Methods （$\mathit{D}={\mathit{D}}_{\mathrm{1}}$）

平差方法	${\widehat{a}}_0$	${\widehat{a}}_1$	${\widehat{a}}_2$	${\widehat{b}}_0$	${\widehat{b}}_1$	${\widehat{b}}_2$
PEIV模型[10]	10.412 055	4.065 602	-2.136 739	-9.641 747	1.116 170	2.912 190
变量投影法[14]	10.412 055	4.065 602	-2.136 739	-9.641 747	1.116 170	2.912 190
本文方法	10.412 055	4.065 602	-2.136 739	-9.641 747	1.116 170	2.912 190

下载: 导出CSV 

| 显示表格

当取$\mathit{D}={\mathit{D}}_{\mathrm{2}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathrm{4},\cdots ,\mathrm{4,1},\cdots ,\mathrm{1})$（4和1的个数都为24）时,采用本文方法和STLN方法^[17]进行求解,得到的参数估计结果如表2所示。

表  2  不同平差方法的SWTLS解（$\mathit{D}={\mathit{D}}_{\mathrm{2}}$）

Table  2.  SWTLS Solutions with Different Adjustment Methods （$\mathit{D}={\mathit{D}}_{\mathrm{2}}$）

平差方法	${\widehat{a}}_0$	${\widehat{a}}_1$	${\widehat{a}}_2$	${\widehat{b}}_0$	${\widehat{b}}_1$	${\widehat{b}}_2$
STLN方法[17]	10.418 681	4.056 913	-2.101 082	-9.639 440	1.118 336	2.855 364
本文方法	10.418 681	4.056 913	-2.101 082	-9.639 440	1.118 336	2.855 364

下载: 导出CSV 

| 显示表格

从表1可以看出,由于PEIV模型^[10]、变量投影法^[14]均未考虑源系统中坐标的重复次数,两者的平差准则等价于本文提出的通用模型中重复矩阵取单位矩阵时的准则（式（21））,尽管三者对结构误差的处理采用了不同的函数模型,但在相同的准则下都得到了完全一致的平差结果,说明本文提出的通用平差算法是可行有效的。由表2可以看出,STLN方法^[17]顾及了系数阵中随机观测值的重复次数,并且是将观测重复数的平方纳入平差模型,和本文通用模型在同样的平差准则下也得到了一致的平差结果,进一步证明了本文算法的可行性和有效性。

为了验证不同平差准则下SWTLS解的统计性质,模拟系数阵元素和右端项含有公共元素的SEIV模型,系数阵和右端项的真值及误差的结构分别为：

待估参数的真值为$\tilde{\mathit{x}}={[\begin{array}{ccc}\mathrm{1}& \mathrm{5}& \mathrm{2}\end{array}]}^{\mathrm{T}}$,独立误差向量$\mathit{\gamma }=[{\gamma }_{\mathrm{1}}{\gamma }_{\mathrm{2}}\cdots {\gamma }_{\mathrm{25}}{]}^{\mathrm{T}}$是零均值独立同分布的随机变量,其方差为${\mathit{D}}_{\mathit{\gamma }}={\sigma }_{\mathrm{0}}^{\mathrm{2}}{\mathit{I}}_{\mathrm{25}}$,其中${\mathit{I}}_{\mathrm{25}}$表示25×25阶单位矩阵,${\gamma }_i$（$i=\mathrm{1,2},\cdots ,\mathrm{25}$）的重复次数$d_i\mathrm{可以}$由$\left({\mathit{E}}_{\mathit{A}}|{\mathit{e}}_{\mathit{y}}\right)$的结构获得。分别采用3种平差准则计算参数的估值,准则式（21）中对应的$\mathit{D}$矩阵分别取${\mathit{D}}_{\mathrm{1}}={\mathit{I}}_{\mathrm{25}}$,${\mathit{D}}_{\mathrm{2}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(d_{\mathrm{1}},d_{\mathrm{2}},\cdots ,d_{\mathrm{25}})=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathrm{1,3},\mathrm{1,3},\cdots ,\mathrm{1,3},\mathrm{1})$,D₃=diag（$d_{\mathrm{1}}^{\mathrm{2}}$,$d_{\mathrm{2}}^{\mathrm{2}}$,$\cdots$,$d_{\mathrm{25}}^{\mathrm{2}}$）=diag（1,9,1,9,$\cdots$,1,9,1）。根据式（44）和式（48）计算单位权方差和参数估值的方差。由于参数真值已知,可以计算均方误差$\mathrm{M}\mathrm{S}\mathrm{E}(\widehat{\mathit{x}})={\left(\widehat{\mathit{x}}-\tilde{\mathit{x}}\right)}^{\mathrm{T}}\left(\widehat{\mathit{x}}-\tilde{\mathit{x}}\right)$,以此来检核3种平差准则下估值的精度,并与近似方差对照。分别取误差水平${\sigma }_{\mathrm{0}}^{\mathrm{2}}=\mathrm{0.25}$和${\sigma }_{\mathrm{0}}^{\mathrm{2}}=\mathrm{1}$进行计算。为了体现结果的统计性,将不同误差水平的实验各进行10 000次,将每次运算的结果取平均值,结果分别见表3和表4。

表  3  不同平差准则下的SWTLS结果（${\sigma }_{\mathrm{0}}^{\mathrm{2}}=\mathrm{0.25}$）

Table  3.  SWTLS Results Under Different Adjustment Principles （${\sigma }_{\mathrm{0}}^{\mathrm{2}}=\mathrm{0.25}$）

统计项	统计量	$\mathit{D}={\mathit{D}}_1$	$\mathit{D}={\mathit{D}}_2$	$\mathit{D}={\mathit{D}}_3$
估计参数值	${\widehat{x}}_1$	0.999 998	1.000 281	0.999 975
	${\widehat{x}}_2$	5.002 389	5.007 359	5.009 953
	${\widehat{x}}_3$	1.999 596	1.999 320	1.998 190
单位权方差	${\widehat{\sigma }}_0^2$	0.239 737	0.537 349	1.463 591
参数方差估值	${\widehat{\sigma }}_{x_1}^2$	9.43×10-6	1.72×10-5	5.76×10-5
	${\widehat{\sigma }}_{x_2}^2$	5.23×10-4	9.14×10-4	3.19×10-3
	${\widehat{\sigma }}_{x_3}^2$	8.28×10-6	1.60×10-5	5.05×10-5
	${\widehat{\sigma }}_{\mathrm{\Sigma }}^2$	5.40×10-4	9.47×10-4	3.30×10-3
MSE	${\widehat{m}}_{x_1}^2$	9.76×10-6	1.10×10-5	1.63×10-5
	${\widehat{m}}_{x_2}^2$	5.44×10-4	6.21×10-4	9.40×10-4
	${\widehat{m}}_{x_3}^2$	8.50×10-6	9.39×10-6	1.42×10-5
	${\widehat{m}}_{\mathrm{\Sigma }}^2$	5.63×10-4	6.41×10-4	9.71×10-4

下载: 导出CSV 

| 显示表格

表  4  不同平差准则下的SWTLS结果（${\sigma }_{\mathrm{0}}^{\mathrm{2}}=\mathrm{1}$）

Table  4.  SWTLS Results Under Different Adjustment Principles （${\sigma }_{\mathrm{0}}^{\mathrm{2}}=\mathrm{1}$）

统计项	统计量	$\mathit{D}={\mathit{D}}_1$	$\mathit{D}={\mathit{D}}_2$	$\mathit{D}={\mathit{D}}_3$
估计参数值	${\widehat{x}}_1$	1.000 220	1.003 770	1.007 659
	${\widehat{x}}_2$	4.997 838	5.022 963	5.049 989
	${\widehat{x}}_3$	2.000 725	2.004 000	2.007 714
单位权方差	${\widehat{\sigma }}_0^2$	0.960 341	2.146 473	5.839 121
参数方差估值	${\widehat{\sigma }}_{x_1}^2$	3.78×10-5	6.87×10-5	2.30×10-4
	${\widehat{\sigma }}_{x_2}^2$	2.09×10-3	3.65×10-3	1.27×10-2
	${\widehat{\sigma }}_{x_3}^2$	3.32×10-5	6.41×10-5	2.02×10-4
	${\widehat{\sigma }}_{\mathrm{\Sigma }}^2$	2.17×10-3	3.78×10-3	1.32×10-2
MSE	${\widehat{m}}_{x_1}^2$	4.02×10-5	4.50×10-5	6.67×10-5
	${\widehat{m}}_{x_2}^2$	2.21×10-3	2.53×10-3	3.85×10-3
	${\widehat{m}}_{x_3}^2$	3.53×10-5	3.88×10-5	5.82×10-5
	${\widehat{m}}_{\mathrm{\Sigma }}^2$	2.29×10-3	2.62×10-3	3.97×10-3

下载: 导出CSV 

| 显示表格

从表3可以看出,不考虑重复次数的平差准则（采用${\mathit{D}}_{\mathrm{1}}$）获得的参数估值更接近于真值。参数各分量的MSE最小,从数值上验证了不考虑系数阵误差重复系数的平差准则是最优的。单位权方差的估值略小于真值,是由于式（43）没有考虑非线性模型线性化引入的偏差项。同理,各参数分量的方差均小于对应的MSE,原因是SEIV模型的非线性特性会引入估计偏差,方差加上偏差的平方和等于MSE。但是,对比方差和MSE的差异可知,参数3个分量的方差估值和MSE之间的差异（$\left|({\widehat{\sigma }}_{x_i}^{\mathrm{2}}-{\widehat{m}}_{x_i}^{\mathrm{2}})/{\widehat{m}}_{x_i}^{\mathrm{2}}\right|$）分别为3.38%、3.86%和2.59%,说明当误差较小时,参数方差估值和MSE的差异很小,可以作为精度评定的指标。因此,本文不再采用非线性最小二乘偏差修正或Monte Carlo模拟的思想求解偏差值^[10,24]。

从上述理论分析可知,将SEIV模型线性化后,其平差准则不应顾及误差重复次数。式（44）和式（48）只适用于$\mathit{D}$为单位阵的情形。从后面两列可以看出,当平差准则中计入误差重复的次数时,单位权方差估值与真值以及参数分量的方差估值与MSE均有较大的差异。当$\mathit{D}={\mathit{D}}_{\mathrm{3}}$时,单位权方差估值约为真值的6倍,方差估值约为对应MSE的3.5倍,这是由于单位权方差公式中重复计算残差平方和所致,由此进一步证明不能采用这两种平差准则用线性近似方法求单位权中误差及方差。