一种加权整体最小二乘估计的高效算法

王建民，倪福泽，赵建军; WANG Jianmin; NI Fuze; ZHAO Jianjun

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

一种加权整体最小二乘估计的高效算法 PDF

- ORCID：
王建民 ^1,2
- ORCID：
倪福泽 ³
- ORCID：
赵建军 ²

1. 太原理工大学矿业工程学院，山西太原 030024； 2. 成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川成都 610059； 3. 中煤（西安）航测遥感研究院有限公司，陕西西安 710100

中图分类号： P207

最近更新：2021-05-12

DOI：10.11908/j.issn.0253-374x.20253

摘要

加权整体最小二乘法（WTLS）是估计errors-in-variables （EIV）模型参数严密的方法，当面临大数据集时，其计算效率有限。针对EIV模型中设计矩阵呈现出的结构性特征，在最小二乘准则的约束条件下，通过仅给设计矩阵的随机列赋予权重，推证了适用于EIV模型参数估计的部分加权整体最小二乘法（PWTLS）。PWTLS无需借助拉格朗日辅助法，能够精确估计EIV模型参数；另外，该算法缩减了矩阵的维数，同时在迭代过程中避免了估计设计矩阵的随机误差，从而减小了矩阵运算量，提升了计算效率。最后以真实数据和模拟数据为例与其他7种同类算法进行对比，结果表明， PWTLS取得了与同类算法相同的精度，但计算效率显著提高，验证了算法的可行性。

关键词

整体最小二乘; 变量误差模型; 计算效率; 坐标转换

高斯‒马尔可夫模型在许多工程实践中得到成功应用，通常认为模型中的设计矩阵是无误差的。然而，在许多情景下，设计矩阵中的元素并不能保证是常数，其本身也是含有随机误差的观测值，设计矩阵和观测向量中含有随机误差的模型通常被定义为变量误差模型（errors-in-variables， EIV）^［

1-2］。

最初，因为求解EIV模型参数的方法不成熟，通常忽略设计矩阵中的随机误差，而采用加权最小二乘法（weighted least-squares，WLS）近似估计EIV参数，这种方法是不严密的。1980年，Golub等^［

3］提出采用整体最小二乘法（total least-squares ，TLS）估计EIV模型的未知参数。当设计矩阵和观测向量均含有随机误差时，TLS就是最小二乘（least-squares，LS）的扩展且更具一般性^{［参考文献 4

百度学术}4］。TLS将设计矩阵中的所有元素按等权对待，但实际上设计矩阵中的元素并非都含有随机误差，且精度也不尽相同，因此TLS被进一步扩展为加权整体最小二乘法（weighted TLS，WTLS），WTLS分别给观测向量和设计矩阵附加了不同的权重，比TLS更具一般性。

一般来说，主要有两类方法实现TLS和WTLS的求解，一类方法是基于奇异值分解（singular value decomposition，SVD）的直接解算方法^［

4–8］，Akyilmaz^{［参考文献 6

百度学术}6］在SVD分解中给设计矩阵和观测向量赋予了相应的权矩阵，但引入的权阵并没有真正起到调配误差的作用，其本质上是对每个参数和观测值进行了缩放^{［参考文献 9

百度学术}9］。另一类常用的方法是基于拉格朗日辅助法构建极值函数分别对未知量进行求偏导数，许多文献据此提出了WTLS的迭代算法^{［参考文献 2

百度学术}2，5，8，10-12］，建立在统计框架下的这些迭代算法估计EIV模型参数是严密的^{［参考文献 13

百度学术}13］。因为WTLS同时顾及了设计矩阵和观测向量的权重，在工程实践中应用迭代法估计EIV模型参数更为普遍。

从现有的研究成果来看，虽然国内外学者就EIV模型提出了许多卓有成效的解算方法，但是与传统的WLS相比，算法比较复杂，不具有操作上的优势^［

14］，而且大多数迭代算法仅考虑了参数估计精度和迭代收敛速率，当面临大数据集时计算效率有限。

为了提高WTLS的效率，一些学者从不同侧面改进了WTLS算法。Shen等^［

10］提出基于Newton-Gauss的非线性WLS平差方法，Jazaeri等^{［参考文献 15

百度学术}15］导出了基于WLS准则下的迭代算法，这些方法都明显提升了迭代收敛速率。但仅通过增强迭代收敛速率，计算效率提升有限。事实上，绝大多数EIV模型的设计矩阵具有结构性特征，仅有部分元素是随机的，其通常的解法是构造特定的权矩阵，使系数矩阵中固定元素位置的改正数为零，但这样会导致部分矩阵维数增大，从而加大了矩阵的运算量，不利于计算效率的提高。为此，Xu等^{［参考文献 9

百度学术}9］提出了partial EIV（PEIV）模型，并给出了相应的算法。PEIV模型一个显著的优点是未知数的总量明显减小，然而多次迭代收敛拖慢了其计算效率^{［参考文献 9

百度学术}9，16-17］。Zhao^{［参考文献 17

百度学术}17］应用拉格朗日辅助法对参量求导增强了PEIV算法的效率，王乐洋等^{［参考文献 16

百度学术}16］应用间接平差的原理进行改进，使算法更加简洁、紧凑。这两种方法改进的结果是迭代快速收敛，效率明显提高，但不足之处是需要根据原设计矩阵的结构解析出一个固定矩阵和固定向量，而且在迭代过程仍需重构设计矩阵，不利于效率的提升。

WTLS在大地测量等众多科学研究和工程领域中得到应用^［

18-19］， WTLS算法基本成熟，渐成理论体系。算法之间的差异主要体现在计算效率和算法原理的复杂程度上。在大数据的背景下，计算效率是决定算法能有效应用的一个重要因素。例如应用WTLS处理点云数据，数据量会显著增加，算法性能就会显得尤为重要。

本文针对EIV模型设计矩阵呈现出的结构性特征，将设计矩阵分成随机和固定部分，仅给随机部分赋予权重，并且不涉及拉格朗日辅助法对参量求偏导，在WLS平差约束准则下，推证了一种估计EIV模型参数的迭代算法，该算法能够显著降低矩阵的运算量，极大提升计算效率。

1 基于WLS准则的WTLS高效算法原理

1.1　EIV模型及WTLS算法

EIV模型的函数模型和随机模型通常表示为^［

8-9］

L = (A - E) θ + e

（1）

[\begin{matrix} e \\ v e c (E) \end{matrix}] ~ ([\begin{matrix} 0 \\ 0 \end{matrix}], σ^{2} [\begin{matrix} Q_{e} & 0 \\ 0 & Q_{E} \end{matrix}])

（2）

式（1）、（2）中：E为（n×t）阶设计矩阵A的随机误差；e为n维观测向量L的随机误差；（n×n）阶矩阵Q_e和（nt×nt）阶矩阵Q_E分别为e和vec（E）的协因数矩阵；vec （⋅）表示矩阵的拉直变换；θ为待估计的t个未知参数；σ²为单位权方差。

对于EIV模型，可直接应用最小二乘准则 $e^{T} {Q_{e}}^{- 1} e = m i n$ ，其参数估值为^［

15］

\hat{θ} = {[{(A - \hat{E})}^{T} Q_{e}^{- 1} (A - \hat{E})]}^{- 1} {(A - \hat{E})}^{T} Q_{e}^{- 1} L

（3）

式中： $\hat{E}$ 和 $\hat{θ}$ 分别表示 E和θ的估计值。

仅对随机误差e加以约束，难以估计E，一般的WTLS的平差准则为^［

10-11，15］

e^{T} {Q_{e}}^{- 1} e + v e c {(E)}^{T} {Q_{E}}^{- 1} v e c (E) = m i n

（4）

在式（4）的约束条件下，利用拉格朗日辅助法可以精确估计未知参数θ^［

8］，即

\begin{array}{l} \hat{θ} = {[{(A - \hat{E})}^{T} Q^{- 1} (A - \hat{E})]}^{- 1} \times \\ {(A - \hat{E})}^{T} Q^{- 1} (L - \hat{E} \hat{θ}) (5) \end{array}

v e c (\hat{E}) = - Q_{E} X Q^{- 1} (L - A \hat{θ})

（6）

式（6）中： $Q = Q_{e} + {\hat{X}}^{T} Q_{E} \hat{X}$ ， $\hat{X} = (\hat{θ} \otimes I_{n})$ ， I_n为单位阵， $\otimes$ 表示两个矩阵的直积运算。

式（6）是WTLS常用的迭代算式，一般的WTLS方法在迭代过程中需要更新矩阵 $\hat{E}$ ，不利于算法效率的提升。接下来，根据设计矩阵A呈现出的结构性特征，推导一种基于WLS准则的EIV模型参数估计高效算法。

1.2　顾及设计矩阵结构的WTLS算法

理论上EIV模型中的设计矩阵的元素均可包含随机误差，但在实际应用中，绝大多数EIV模型的设计矩阵都会存在部分固定列并且不含有随机误差^［

9，17，20］。对于含有固定列的设计矩阵，通常是构建一个特定的协因数矩阵Q_E，在迭代的过程中使得固定元素的改正数为零，但这样会导致Q_E的维数被固定列“扩大”，制约了计算效率。

设A中有t₁列固定列和t₂列随机列，将A中的固定列和随机列进行矩阵分块，对应的参数θ也分成两部分，即

A = [A_{1}, A_{2}]

（7）

θ = [θ_{1}, θ_{2}]^{T}

（8）

式（7）、（8）中：A₁为由固定列组成的（n×t₁）阶矩阵；A₂为由随机列组成的（n×t₂）阶矩阵；与A₁和A₂对应的参数分别为θ₁和θ₂。

根据式（7）和式（8），EIV模型（式（1）和式（2））重新改写为

L = (A_{1}, A_{2} - E_{2}) (\begin{array}{l} θ_{1} \\ θ_{2} \end{array}) + e

（9）

[\begin{matrix} e \\ v e c (E_{2}) \end{matrix}] ~ ([\begin{matrix} 0 \\ 0 \end{matrix}], σ^{2} [\begin{matrix} Q_{e} & 0 \\ 0 & Q_{E_{2}} \end{matrix}])

（10）

式（9）、（10）中： E₂为A₂的随机误差矩阵，其阶数为（n×t₂）； $Q_{E_{2}}$ 为vec（E₂）的协因数矩阵，其阶数为（nt₂×nt₂）。

将式（9）改写为以下误差方程：

L - A_{1} {\hat{θ}}_{1} = {\hat{X}}_{2}^{T} v e c ({\hat{A}}_{2}) + \hat{e}

（11）

式中： ${\hat{X}}_{2}^{T} = ({\hat{θ}}_{2}^{T} \otimes I_{n})$ 。

顾及式（11）， $v e c (A_{2})$ 表示为

v e c (A_{2}) = v e c ({\hat{A}}_{2}) + v e c ({\hat{E}}_{2})

（12）

将式（11）和式（12）组合形成如下的非线性方程，即

H = B v e c ({\hat{A}}_{2}) + \hat{Ω}

（13）

式中： $H = [\begin{array}{l} L - A_{1} {\hat{θ}}_{1} \\ v e c (A_{2}) \end{array}]$ ； $B = [\begin{array}{l} {\hat{X}}_{2}^{T} \\ I_{n t_{2}} \end{array}]$ ； $\hat{Ω} = [\begin{array}{l} \hat{e} \\ v e c ({\hat{E}}_{2}) \end{array}]$ 。

在 ${\hat{Ω}}^{T} P_{H} \hat{Ω} = m i n$ 的约束下，应用WLS原理，式（13）中的 $v e c ({\hat{A}}_{2})$ 的估值为

\begin{array}{l} v e c ({\hat{A}}_{2}) = [B^{T} P_{H} {B]}^{- 1} B^{T} P_{H} H = \\ (Q_{E_{2}}^{- 1} + {\hat{X}}_{2} Q_{e}^{- 1} {\hat{X}}_{2}^{T})^{- 1} \times \\ [{\hat{X}}_{2} Q_{e}^{- 1} (L - A_{1} {\hat{θ}}_{1}) + Q_{E_{2}}^{- 1} v e c (A_{2})] \end{array}

（14）

式中： $P_{H} = [\begin{matrix} Q_{e}^{- 1} & 0 \\ 0 & Q_{_{E_{2}}}^{- 1} \end{matrix}]$ 。

为进一步化简式（14），需利用式（15）的矩阵求逆引理^［

21］，即

(J^{- 1} - K D^{- 1} {Z)}^{- 1} = J + J K {(D - Z J K)}^{- 1} Z J

（15）

令 $J^{- 1} = Q_{E_{2}}^{- 1}, K = - {\hat{X}}_{2}, D = Q_{e}$ 且 $Z = {\hat{X}}_{2}^{T}$ ，式（14）化简为

\begin{array}{l} v e c ({\hat{A}}_{2}) = [Q_{E_{2}}^{- 1} - Q_{E_{2}}^{- 1} {\hat{X}}_{2} (Q_{e} + {\hat{X}}_{2}^{T} Q_{E_{2}}^{- 1} {\hat{X}}_{2})^{- 1} \times \\ {\hat{X}}_{2}^{T} Q_{E_{2}}^{- 1}] [{\hat{X}}_{2} Q_{e}^{- 1} (L - A_{1} {\hat{θ}}_{1}) + Q_{E_{2}}^{- 1} v e c (A_{2})] = \\ v e c (A_{2}) + Q_{E_{2}}^{- 1} {\hat{X}}_{2} (Q_{e} + {\hat{X}}_{2}^{T} Q_{E_{2}}^{- 1} {\hat{X}}_{2})^{- 1} \times \\ (L - A_{1} {\hat{θ}}_{1} - A_{2} {\hat{θ}}_{2}) = \\ v e c (A_{2}) + Q_{E_{2}}^{- 1} {\hat{X}}_{2} Q_{2}^{- 1} (L - A \hat{θ}) \end{array}

（16）

式中： $Q_{2} = (Q_{e} + {\hat{X}}_{2}^{T} Q_{E_{2}} {\hat{X}}_{2})$ 。

由式（16）可知，设计矩阵A₂的随机误差的估值为

v e c ({\hat{E}}_{2}) = - Q_{E_{2}} {\hat{X}}_{2} Q_{2}^{- 1} (L - A \hat{θ})

（17）

由式（7）可知，观测向量L的残差可表示为

\begin{array}{l} \hat{e} = L - \hat{A} \hat{θ} = \\ L - [A_{1}, A_{2} - {\hat{E}}_{2}] [{\hat{θ}}_{1}, {\hat{θ}}_{2}]^{T} = \\ L - A \hat{θ} + {\hat{X}}_{2}^{T} v e c ({\hat{E}}_{2}) = \\ L - A \hat{θ} - {\hat{X}}_{2}^{T} Q_{E_{2}} {\hat{X}}_{2} Q_{2}^{- 1} (L - A \hat{θ}) = \\ (I - {\hat{X}}_{2}^{T} Q_{E_{2}} {\hat{X}}_{2} Q_{2}^{- 1}) (L - A \hat{θ}) = \\ Q_{e} Q_{2}^{- 1} (L - A \hat{θ}) \end{array}

（18）

由式（18）可知

L = \hat{A} \hat{θ} + Q_{e} Q_{2}^{- 1} (L - A \hat{θ})

（19）

将式（19）代入式（3），得

\begin{array}{l} \hat{θ} = ({\hat{A}}^{T} Q_{e}^{- 1} {\hat{A})}^{- 1} {\hat{A}}^{T} Q_{e}^{- 1} [\hat{A} \hat{θ} + Q_{e} Q_{2}^{- 1} (L - A \hat{θ})] = \\ \hat{θ} + ({\hat{A}}^{T} Q_{e}^{- 1} {\hat{A})}^{- 1} {\hat{A}}^{T} Q_{2}^{- 1} (L - A \hat{θ}) \end{array}

（20）

则

{\hat{A}}^{T} λ = {(A - \hat{E})}^{T} λ = 0

（21）

式中： $λ = Q_{2}^{- 1} (L - A \hat{θ})$ ，λ为n×1的向量。

由式（21）可知

A^{T} λ = {\hat{E}}^{T} λ

（22）

根据A的矩阵分块（式（7）），式（22）可分解为两部分，即

A_{1}^{T} λ = {\hat{E}}_{1}^{T} λ = O λ

（23）

A_{2}^{T} λ = {\hat{E}}_{2}^{T} λ

（24）

式中： ${\hat{E}}_{1}^{T} = O$ ，O 为（t₁×n）零矩阵。

式（25）和式（26）是两个克罗内克乘积法则，用于式（24）。

v e c (S_{n \times 1}) = v e c (S_{n \times 1}^{T})

（25）

v e c (A B C) = (C^{T} \otimes A) v e c (B)

（26）

式（24）可表示为

A_{2}^{T} λ = {\hat{E}}_{2}^{T} λ = I_{t_{2}} {\hat{E}}_{2}^{T} λ = v e c (λ^{T} {\hat{E}}_{2} I_{t_{2}})

（27）

则

A_{2}^{T} λ = A_{2}^{T} Q_{2}^{- 1} (L - A \hat{θ}) = (I_{t_{2}} \otimes λ^{T}) v e c ({\hat{E}}_{2})

（28）

将式（17）代入式（28），得

A_{2}^{T} Q_{2}^{- 1} (L - A \hat{θ}) = - (I_{t_{2}} \otimes λ^{T}) Q_{E_{2}} {\hat{X}}_{2} Q_{2}^{- 1} (L - A \hat{θ}) = G (L - A \hat{θ})

（29）

式中： $G = - (I_{t_{2}} \otimes λ^{T}) Q_{E_{2}} {\hat{X}}_{2} Q_{2}^{- 1}$ ，G为（t₂×n）的矩阵。

组合式（23）和式（29），得

[\begin{array}{l} A_{1}^{T} \\ A_{2}^{T} \end{array}] Q_{2}^{- 1} (L - A \hat{θ}) = [\begin{matrix} O \\ G \end{matrix}] (L - A \hat{θ})

（30）

式（30）可化简为

(A^{T} Q_{2}^{- 1} A - U A) \hat{θ} = (A^{T} Q_{2}^{- 1} - U) L

（31）

式中： $U = {[\begin{matrix} O & G \end{matrix}]}^{T}$ 。

式（31）中的 $(A^{T} Q_{2}^{- 1} A - U A)$ 是列满秩矩阵，则未知参数 $\hat{θ}$ 最终的估计结果为

\hat{θ} = (A^{T} Q_{2}^{- 1} {A - U A)}^{- 1} (A^{T} Q_{2}^{- 1} - U) L

（32）

由式（17）和式（18）可估计 $σ^{2}$ ，即

{\hat{σ}}^{2} = \frac{{\hat{e}}^{T} Q_{e}^{- 1} \hat{e} + v e c ({\hat{E}}_{2})^{T} Q_{E_{2}}^{- 1} v e c ({\hat{E}}_{2})}{n - t} = \frac{λ^{T} (A \hat{θ} - L)}{n - t}

（33）

上述推导是基于WLS最小二乘准则，并且式（32）在形式上与间接平差的法方程类似，但方程两侧都存在未知参数 $\hat{θ}$ ，仍需迭代计算。为了加以区别和简写，上述算法称为部分加权整体最小二乘法（partially weighted TLS，PWTLS），PWTLS的算法流程步骤如下：

（1）初始化初值，即

{\hat{θ}}_{(0)} = {(A^{T} Q_{L}^{- 1} A)}^{- 1} A^{T} Q_{L}^{- 1} L

（2）开始迭代，即

{\hat{θ}}_{(i)} = {[\begin{matrix} {\hat{θ}}_{1 (i)} & {\hat{θ}}_{2 (i)} \end{matrix}]}^{T}

X_{2 (i)}^{T} = {\hat{θ}}_{2 (i)}^{T} \otimes I_{n}

{Q_{2}}_{(i)} = (Q_{e} + {\hat{X}}_{2 (i)}^{T} Q_{E_{2}} {\hat{X}}_{2 (i)})

λ_{(i)} = Q_{2 (i)}^{- 1} (L - A {\hat{θ}}_{(i)})

G_{(i)} = - (I_{t_{2}} \otimes λ_{(i)}^{T}) Q_{E_{2}} {\hat{X}}_{2 (i)} Q_{2 (i)}^{- 1}

U_{(i)} = {[\begin{matrix} O & G_{(i)} \end{matrix}]}^{T}

{\hat{θ}}_{(i + 1)} = {[A^{T} Q_{2 (i)}^{- 1} A - U_{(i)} A]}^{- 1} (A^{T} Q_{2 (i)}^{- 1} - U_{(i)}) L

（3）重复步骤（2）直到 $|{\hat{θ}}_{(i + 1)} - {\hat{θ}}_{(i)}| < ε$ ，迭代结束（下标i表示迭代次数，ε为设定的收敛阈值）。

由上述算法可知，在迭代过程中，并不涉及对设计矩阵A拆分成A₁和A₂，仅需把随机列A₂对应的参数θ₂分解出来参与迭代。只要EIV模型的设计矩阵存在固定列，根据列与参数的对应关系，很容易将A组合成固定和随机部分。通常情况下，习惯上总是将EIV模型的固定列排列在一起，而且绝大多数EIV模型至少有一列是固定列，因而PWTLS估计EIV参数具有普适性。另外，在实际应用中，设计矩阵的随机列可能会含有少量固定元素，例如三维七参数坐标转换模型中，有3列固定列组成A₁，在随机列中有一个固定零元素，对于随机列中的个别固定元素赋予值为零的协因数后，仍可组成随机列A₂。

从算法流程可以看出，PWTLS算法不同于一般的WTLS算法，在每次迭代的过程中并不需要估计随机误差E；通过给设计矩阵的随机列赋予权重，使得相关矩阵运算的时间复杂度由O（nt）降低到 O（nt₂），占用的内存空间由O（nt²）减小到O（nt₂²），例如原参数矩阵 $\hat{X}$ 的维数（nt × nt）降为 ${\hat{X}}_{2}$ （nt₂× nt₂），基于这两点，使PWTLS算法的效率显著提高。

2 算例与分析

为了验证PWTLS的性能，本文选择了7种具有代表性的WTLS算法^［

10-11，15-17，20，22］与本文算法进行对比分析，这7种算法从不同角度对WTLS进行效率增强，具体原理见参考文献，有的算法在相关文献中也与其他算法进行过对比研究。

2.1　真实算例

直线拟合简单而有效，经常用于证明WTLS算法的有效性。设有直线方程为 $y = a x + b$ ，其中x 和y是直线的坐标分量，a和b是直线的斜率和截距。实验数据来自于文献［

23］，见表1，其中（x_i，y_i）和（w_xi，w_yi）表示坐标和相应的权值。

表1 坐标观测值和相应的权重

Tab. 1 Coordinate observations and weights

编号	x_i	W_xi	y_i	W_yi
1	0	1 000.0	5.9	1.0
2	0.9	1 000.0	5.4	1.8
3	1.8	500.0	4.4	4.0
4	2.6	800.0	4.6	8.0
5	3.3	200.0	3.5	20.0
6	4.4	80.0	3.7	20.0
7	5.2	60.0	2.8	70.0
8	6.1	20.0	2.8	70.0
9	6.5	1.8	2.4	100.0
10	7.4	1.0	1.5	500.0

根据直线方程和坐标点数量，设计矩阵A，参数向量θ及观测向量L具体形式为

A = [\begin{matrix} 1 & x_{1} \\ 1 & x_{2} \\ ⋮ & ⋮ \\ 1 & x_{10} \end{matrix}], L = [\begin{matrix} y_{1} \\ y_{2} \\ ⋮ \\ y_{10} \end{matrix}], θ = [\begin{matrix} a \\ b \end{matrix}]

（34）

协因数矩阵Q_e 和Q_E 分别为

\begin{array}{l} Q_{e} = d i a g (W_{y_{1}}^{- 1}, W_{y_{2}}^{- 1}, \dots, W_{y_{10}}^{- 1}), \\ Q_{E} = d i a g (\underset{10}{\underset{⏟}{0,0, \dots, 0,}} W_{x_{1}}^{- 1}, W_{x_{2}}^{- 1}, \dots, W_{x_{10}}^{- 1}) \end{array}

（35）

在式（35）中，设计矩阵A可分成A₁和A₂两部分，即

A_{1} = {[\underset{10}{\underset{⏟}{1,1, . . ., 1}}]}^{T}, A_{2} = [x_{1}, x_{2}, . . ., x_{10}]^{T}

（36）

A₂的协因数矩阵 $Q_{E_{2}}$ 的结构取自Q_E的一部分，即

Q_{E_{2}} = d i a g (W_{x_{1}}^{- 1}, W_{x_{2}}^{- 1}, \dots, W_{x_{10}}^{- 1})

（37）

设迭代收敛阈值 ε=10^-10 ，分别应用8种算法估计参数和单位权中误差，实验结果见表2。从表2中可明显看出，8种严密算法估计的参数相差在10^-11量级且中误差相同，表明PWTLS能够获得和另外7种严密算法相同的解。文献［

22］对PEIV模型进行改进，其目的也是提升效率，从本例的结果来看，收敛速度高达454次，其余算法的收敛速度都较快。文献［16］同样是对PEIV模型算法进行优化，取得了理想的效果。比较而言，文献［20］的算法收敛最快，收敛速度仅是影响计算效率一方面的因素。本例中直线拟合的设计矩阵结构简单而且仅有10个点，主要用于PWTLS参数估计精度的验证及收敛速度的初步分析，接下来通过大数据量应用仿真模拟实验来对比分析各个算法的效率。

表2 8种不同算法的参数和中误差的估计结果

Tab. 2 Estimated parameters and standard deviations obtained using eight algorithms

算法	参数a	参数b	中误差σ	迭代次数
PWTLS	5.479 910 224 033	-0.4805 334 074 462	1.217 91	8
文献[11]	5.479 910 224 033	-0.4805 334 074 462	1.217 91	8
文献[17]	5.479 910 224 033	-0.4805 334 074 462	1.217 91	8
文献[16]	5.479 910 22 4032	-0.4805 334 074 461	1.217 91	7
文献[20]	5.479 910 224 032	-0.4805 334 074 462	1.217 91	5
文献[15]	5.479 910 224 032	-0.4805 334 074 461	1.217 91	7
文献[10]	5.479 910 224 030	-0.4805 334 074 456	1.217 91	13
文献[22]	5.479 910 226 342	-0.4805 334 079 855	1.217 91	454

2.2　仿真算例

平面相似四参数转换模型是常用的坐标基准转换模型之一，通常由式（38）表示^［

24］

{[\begin{matrix} x \\ y \end{matrix}]}_{T} = [\begin{matrix} ξ \\ η \end{matrix}] + k [\begin{matrix} c o s a & - s i n a \\ s i n a & c o s a \end{matrix}] {[\begin{matrix} x \\ y \end{matrix}]}_{S}

（38）

式中：ξ 和η 为平移参数；k为尺度参数；a为旋转角；下标S和T分别表示原坐标系和目标坐标系。

假设有d个重合点用于求解转换参数，式（38）经过变换可由EIV模型描述，设计矩阵A，参数θ及观测向量L表示为

\underset{2 d, 4}{A} = {[\begin{matrix} 1 & 0 & x_{1} & y_{1} \\ 0 & 1 & y_{1} & - x_{1} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & 0 & x_{d} & y_{d} \\ 0 & 1 & y_{d} & - x_{d} \end{matrix}]}_{S}, θ = [\begin{matrix} ξ \\ η \\ u \\ w \end{matrix}], \underset{2 d, 1}{L} = {[\begin{matrix} x_{1} \\ y_{1} \\ ⋮ \\ x_{d} \\ y_{d} \end{matrix}]}_{T}

（39）

式中： $u = k c o s a$ ； $w = k s i n a$ 为辅助参数。

参数ξ 和η对应的系数是1和0的常量，设计矩阵A的第1列和第2列是固定列，组成A₁；其他两列组成A₂形成随机列，A₁和A₂表示为

A_{1} = [\begin{matrix} 1 & 0 \\ 0 & 1 \\ ⋮ & ⋮ \\ 1 & 0 \\ 0 & 1 \end{matrix}],

A_{2} = [\begin{matrix} x_{1} & y_{1} \\ y_{1} & - x_{1} \\ ⋮ & ⋮ \\ x_{d} & y_{d} \\ y_{d} & - x_{d} \end{matrix}]

（40）

对应于A₁和A₂，相应的参数可分成两部分，即

θ_{1} = {[ξ, η]}^{T}, θ_{2} = {[u, w]}^{T}

（41）

原坐标系中的坐标由随机函数产生d个坐标点，目标坐标由给定的参数真值ξ= -27.366， η= -71.185， u=1.000 001 092 和w=6.400 15×10^-7通过转换模型反算得到。原坐标和目标坐标中的随机误差由随机函数（式（42）、（43））随机生成，并将随机误差附加到相应的坐标分量中。

Δ_{S} = {[\begin{matrix} Δ_{x_{1}} & Δ_{y_{1}} & \dots & Δ_{x_{d}} & Δ_{y_{d}} \end{matrix}]}^{T} = m v n r n d (0,0.05,2 d)

（42）

Δ_{T} = {[\begin{matrix} Δ_{x_{1}} & Δ_{y_{1}} & \dots & Δ_{x_{d}} & Δ_{y_{d}} \end{matrix}]}^{T} = m v n r n d (0,0.05,2 d)

（43）

式中：Δ_S和Δ_T分别表示原始坐标和目标坐标中的真误差，mvnrnd （0， σ， 2d）生成期望为0、标准差为 σ = 0.05 m、维数为2d×1的随机向量。

由随机误差和标准差定义观测量向量L的协因数阵Q_e

\underset{2 d, 2 d}{Q_{e}} = d i a g {(\begin{matrix} \frac{Δ_{x_{1}}^{2}}{σ^{2}} & \frac{Δ_{y_{1}}^{2}}{σ^{2}} & \dots & \frac{Δ_{x_{d}}^{2}}{σ^{2}} & \frac{Δ_{y_{d}}^{2}}{σ^{2}} \end{matrix})}_{T}

（44）

根据设计矩阵A₂的结构，vec（A₂）的协因数阵 $Q_{E_{2}}$ 定义为

Q_{E_{2}} = [\begin{matrix} q_{x_{1}}^{1,1} & 0 & 0 & 0 & 0 & 0 & q_{x_{1}}^{1,2 d + 2} & 0 & 0 & 0 \\ 0 & q_{y_{1}}^{2,2} & 0 & 0 & 0 & - q_{y_{1}}^{2,2 d + 1} & 0 & 0 & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & q_{x_{d}}^{2 d - 1,2 d - 1} & 0 & 0 & 0 & 0 & 0 & q_{x_{d}}^{2 d - 1,4 d} \\ 0 & 0 & 0 & 0 & q_{y_{d}}^{2 d, 2 d} & 0 & 0 & 0 & - q_{y_{d}}^{2 d, 4 d - 1} & 0 \\ 0 & - q_{y_{1}}^{2 d + 1,2} & 0 & 0 & 0 & q_{y_{1}}^{2 d + 1,2 d + 1} & 0 & 0 & 0 & 0 \\ q_{x_{1}}^{2 d + 2,1} & 0 & 0 & 0 & 0 & 0 & q_{x_{1}}^{2 d + 2,2 d + 2} & 0 & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & 0 & 0 & 0 & - q_{y_{d}}^{4 d - 1,2 d} & 0 & 0 & 0 & q_{y_{d}}^{4 d - 1,4 d - 1} & 0 \\ 0 & 0 & 0 & q_{x_{d}}^{4 d, 2 d - 1} & 0 & 0 & 0 & 0 & 0 & q_{x_{d}}^{4 d, 4 d} \end{matrix}]

（45）

式中： $q_{x_{i}} = \frac{Δ_{x_{i}}^{2}}{σ^{2}}; q_{y_{i}} = \frac{Δ_{y_{i}}^{2}}{σ^{2}} (i = 1,2, . . ., d)$ ；上标表示元素在矩阵中的位置，其中 $Δ_{x_{i}}$ 和 $Δ_{y_{i}}$ 来自于向量Δ_S。

基于上述的模型和数据，随机模拟1 000次实验，固定重合点数d=200，应用均方根误差（RMSE）来衡量参数估计的精度，参数ξ、η、u 和w的RMSE由式（46）、（47）估计

{\hat{σ}}_{ξ} = \sqrt[]{\frac{1}{1 000} \sum_{i = 1}^{1 000} (ξ - {\hat{ξ}}_{i})^{2}}

，

{\hat{σ}}_{η} = \sqrt[]{\frac{1}{1 000} \sum_{i = 1}^{1 000} (η - {\hat{η}}_{i})^{2}}

（46）

{\hat{σ}}_{u} = \sqrt[]{\frac{1}{1 000} \sum_{i = 1}^{1 000} (u - {\hat{u}}_{i})^{2}}

，

{\hat{σ}}_{w} = \sqrt[]{\frac{1}{1 000} \sum_{i = 1}^{1 000} (w - {\hat{w}}_{i})^{2}}

（47）

式（46）、（47）中： ${\hat{ξ}}_{i}$ 、 ${\hat{η}}_{i}$ 、 ${\hat{u}}_{i}$ 和 ${\hat{w}}_{i}$ 为经第i次模拟的对应参数的估计值。

在相同数据和计算平台的实验条件下（收敛阈值ε =10^-10），分别应用7种具有较高收敛速度的算法进行模拟，参数的RMSE的统计结果和1 000次模拟的累计耗时见表3，结果表明，严密算法中的RMSE差异为10×10^-5mm，进一步证明了PWTLS与其他严密算法的参数估计精度相同。

表3 1 000次模拟的RMSE和累计耗时

Tab. 3 Results of parameter estimation and details of the time consumed for different algorithms

算法	${\hat{σ}}_{ξ}$ /mm	${\hat{σ}}_{η}$ /mm	${\hat{σ}}_{u}$ /10^-6	${\hat{σ}}_{w}$ /10^-6	平均迭代次数	总耗时/s
PWTLS	6.600 843 647 7	6.769 967 561 0	7.825 463 890 4	8.310 172 114 1	3.63	7.61
文献[11]	6.600 843 647 7	6.769 967 561 0	7.825 463 890 4	8.310 172 114 1	3.63	19.86
文献[17]	6.600 857 165 2	6.769 954 304 2	7.825 463 137 6	8.310 192 209 9	3.19	20.71
文献[16]	6.600 835 421 5	6.769 953 801 2	7.825 473 558 9	8.310 158 454 4	3.08	22.97
文献[20]	6.600 843 647 6	6.769 967 560 9	7.825 463 889 6	8.310 172 114 3	3.01	24.14
文献[15]	6.600 843 647 1	6.769 967 560 9	7.825 463 888 3	8.310 172 114 5	3.08	38.95
文献[10]	6.600 827 049 8	6.769 950 721 3	7.825 406 436 6	8.310 176 011 5	3.73	39.28

从计算的总耗时可以看出，PWTLS是所有严密算法中计算效率最高的，与其他算法计算效率相比，最少提高了161%，最多提高了近416%。其主要是因为PWTLS使用 $Q_{E_{2}}$ 代替了Q_E，维数由800×800降到400×400；另外，在每次迭代过程中，PWTLS算法并不需要估计随机误差矩阵E。

为了突出PWTLS计算效率的优越性，从另一个角度执行了10次模拟实验，要求重合点数据以间隔为100个点，逐渐从d=100增加到1 000。计算耗时统计结果如图1所示，结果表明，PWTLS的耗时曲线变化相对平缓，意味着随着数据量的增大，PWTLS 的优势越来越明显。

图1 7种算法的计算效率比较

Fig. 1 Computational efficiency of seven algorithms

3 结语

WTLS已经在众多科学和工程技术中得到应用，当面临处理大数据集时，其计算效率有限，降低算法的复杂度和提高计算效率对WTLS算法将显得越来越重要。

本文针对EIV模型的设计矩阵呈现出的结构性特征，提出了基于WLS准则约束条件下的PWTLS算法，无需引入辅助量，迭代方程与经典的间接平差法方程类似，便于理解，易于实现，估计的参数是严密的。PWTLS通过给设计矩阵中的随机列赋予权重，使得设计矩阵的协因数阵维数降低，从而在迭代过程中缩减了所有参与运算矩阵的维数，极大地减少了计算量，促进了效率的提升；另外，PWTLS算法避免了重复估计设计矩阵随机列的误差，这两方面的因素使矩阵的运算量降低，从而显著提升了计算效率。最后，通过算例验证了所提出的PWTLS算法与同类算法参数估计精度相当，但计算效率显著提高。

作者贡献声明

王建民：提出了论文中PWTLS算法原理，完成了论文的整体写作。

倪福泽：编写了计算程序，完成了论文中的实验计算和分析，绘制了图表。

赵建军：提供研究经费的资助，修改了中英文摘要。

参考文献

刘经南，曾文宪，徐培亮. 整体最小二乘估计的研究进展［J］. 武汉大学学报（信息科学版）， 2013， 38（5）： 505. [百度学术]

LIU Jingnan， ZENG Wenxian， XU Peiliang. Overview of total least squares methods［J］. Geomatics and Information Science of Whuhan University， 2013， 38（5）： 505. [百度学术]

FANG X. Weighted total least squares solutions for applications in geodesy［D］. Hannover：Gottfried Wilhelm Leibniz Universität Hannover， 2011. [百度学术]

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. [百度学术]

MARKOVSKY I， VAN HUFFEL S . Overview of total least-squares methods［J］. Signal Processing， 2007， 87（10）： 2283. [百度学术]

TONG X， JIN Y， LI L. An improved weighted total least squares method with applications in linear fitting and coordinate transformation［J］. Journal of Surveying Engineering， 2011， 137（4）： 120. [百度学术]

AKYILMAZ O. Total least squares solution of coordinate transformation［J］. Survey Review， 2007， 39（303）： 68. [百度学术]

FELUS Y A. Application of total least squares for spatial point process analysis［J］. Journal of Surveying Engineering， 2004， 130（3）： 126. [百度学术]

SCHAFFRIN B， WIESER A. On weighted total least-squares adjustment for linear regression［J］. Journal of Geodesy， 2008， 82（7）： 415. [百度学术]

XU P， LIU J， SHI C. Total least squares adjustment in partial errors-in-variables models： algorithm and statistical analysis［J］. Journal of Geodesy， 2012， 86（8）： 661. [百度学术]

SHEN Y， LI B， CHEN Y. An iterative solution of weighted total least-squares adjustment［J］. Journal of Geodesy， 2011， 85（4）： 229. [百度学术]

MAHBOUB V. On weighted total least-squares for geodetic transformations［J］. Journal of Geodesy， 2012， 86（5）： 59. [百度学术]

MAHBOUB V， SHARIFI M A. On weighted total least-squares with linear and quadratic constraints［J］. Journal of Geodesy， 2013， 87（3）： 279. [百度学术]

XU P， LIU J， ZENG W， et al. Effects of errors-in-variables on weighted least squares estimation［J］. Journal of Geodesy， 2014， 88（7）： 705. [百度学术]

姚宜斌，孔建. 顾及设计矩阵随机误差的最小二乘组合新解法［J］. 武汉大学学报（信息科学版）， 2014， 39（9）： 1028. [百度学术]

YAO Yibin， KONG Jian. A new combined LS method considering random errors of design matrix［J］. Geomatics and Information Science of Whuhan University， 2014， 39（9）： 1028. [百度学术]

JAZAERI S， AMIRI-SIMKOOEI A R， SHARIFI M A. Iterative algorithm for weighted total least squares adjustment［J］. Survey Review， 2014， 46（334）： 19. [百度学术]

王乐洋，余航，陈晓勇.Partial EIV模型的解法［J］．测绘学报，2016 ，45（1）：22. [百度学术]

WANG Leyang， YU Hang， CHEN Xiaoyong. An algorithm for partial EIV model［J］. Acta Geodaetica et Cartographica Sinica， 2016， 45（1）： 22. [百度学术]

ZHAO J. Efficient weighted total least-squares solution for partial errors-in-variables model［J］. Survey Review， 2017， 49（356）： 346. [百度学术]

王乐洋.地壳应变参数反演的总体最小二乘方法［J］.大地测量与地球动力学， 2013， 33（3）： 106. [百度学术]

WANG Leyang. Inversion of crustal strain parameters based on total least squares［J］. Journal of Geodesy and Geodynamics，2013，33（3）：106. [百度学术]

方兴，曾文宪，刘经南，等.三维坐标转换的通用整体最小二乘算法［J］. 测绘学报， 2014， 43（11）： 1139. [百度学术]

FANG Xing ， ZENG Wenxian ， LIU Jingnan ， et al.A general total least squares algorithm for three- dimensional coordinate transformations［J］. Acta Geodaetica et Cartographica Sinica， 2014， 43（11）： 1139. [百度学术]

ZHOU Y， KOU X， LI J， et al. Comparison of structured and weighted total least-squares adjustment methods for linearly structured errors-in-variables models［J］. Journal of Surveying Engineering， 2017， 143（1）： 04016019. [百度学术]

HENDERSON H V， SEARLE S R. On deriving the inverse of a sum of matrices［J］. SIAM Review， 1981， 23（1）： 53. [百度学术]

SHI Y， XU P， LIU J， et al. Alternative formulae for parameter estimation in partial errors-in-variables models［J］. Journal of Geodesy， 2015， 89（1）： 13. [百度学术]

NERI F， SAITTA G， CHIOFALO S. An accurate and straightforward approach to line regression analysis of error-affected experimental data［J］. Journal of Physics E： Scientific Instruments， 1989， 22（4）： 215. [百度学术]

WANG J. Locating and estimating multiple gross errors during coordinate transformation［J］. Survey Review， 2015， 47（345）： 458. [百度学术]