总体最小二乘(total least squares，TLS)方法是近30多年来发展起来的一种能同时顾及观测值误差和模型系数矩阵误差的数学方法，其理论及应用研究是目前国内外研究的热点问题^[1]。文献^[2]首次提出总体最小二乘概念，随着变量误差模型(errors-in-variables，EIV)由单变量模型扩展到多元变量模型^[3]，总体最小二乘也相应扩展到多元总体最小二乘(multivariate TLS，MTLS)^{[4, 5]}。某些模型的系数矩阵会因未知参数写成向量形式含有重复出现的随机元素和大量的零元素；多元总体最小二乘将参数向量和观测向量改写成矩阵形式，避免了这种现象，进而削弱了系数矩阵中行列之间的相关性，而且使得系数矩阵的构造更加简单。

在国外相关研究中，文献^{[5, 6, 7]}认为观测矩阵和系数矩阵中的随机元素是等精度的，给出了基于奇异值分解法和拉格朗日乘数法的多元总体最小二乘算法；文献^[8]讨论了基于拉格朗日乘数法的多元加权总体最小二乘算法，将随机元素的方差阵分解成两个矩阵取直积的形式，但是当观测数据含有两个以上的方差分量或是观测数据含有多类观测值时，随机变量的方差协方差阵就很难写成或是不能写成两个矩阵取直积的形式；文献^[9]推导了基于奇异值分解法的多元加权总体最小二乘算法，但是这种算法只使用了行尺度权阵，其权阵不具有一般性。以上文献的所有算法都不适用于系数矩阵与观测矩阵中随机元素具有相关性的情况。文献^[10]给出了基于拉格朗日乘数法的多元加权总体最小二乘算法，其权阵具有一般性，但没有具体给出系数矩阵含有常数元素的解法以及没有解决算法的精度评定问题，这种算法的解算效率也还需要验证。在国内相关研究中，文献^[11]将基于奇异值分解法的多元总体最小二乘算法应用在三维坐标转换中；文献^[12]给出的多元总体最小二乘算法受到正交Procrustes算法制约仅适用于观测值误差对称独立的情况，随机元素相关性问题也没有解决；文献^[13]只给出了观测值等精度同方差情况下的多元总体最小二乘算法。对于系数矩阵含有常量的情况，文献^{[10, 14, 15, 16]}将系数矩阵拆分成分别包含随机元素和常数元素的系数矩阵；文献^{[11, 12, 13]}采用消元法，分两步计算随机元素对应的未知参数和常数元素对应的未知参数；文献^{[8, 17]}使常数元素的方差和协方差为零，然后进行加权解算；文献^[18]首次将EIV模型改写成更具一般性的PEIV(partial EIV)模型，解决了系数矩阵任意位置含有常数元素的平差问题。但以上方法都没有把系数矩阵和观测向量(或矩阵)中的所有元素纳入到一个协因数阵中进行处理。

基于以上问题，本文推导了顾及系数矩阵含有常数元素情况下的多元总体最小二乘的牛顿解法，并讨论了多元总体最小二乘的牛顿法解和拉格朗日乘数法解之间的关系以及给出了多元总体最小二乘的精度评定公式，最后通过模拟数据和实测数据验证了本文算法的可行性和有效性。

多元变量EIV模型^{[6, 7, 8]}为

点击浏览公式

点击浏览公式

顾及系数矩阵A中某列含有常数元素，在文献^{[8, 17, 18]}的基础上，系数矩阵和相应的协因数阵定义为

点击浏览公式

点击浏览公式

式中，A_s是系数矩阵A的随机元素部分；A_d是系数矩阵A的常数元素部分；E_s为随机元素的误差矩阵；Q_i(i=1,2,…，t)∈R^n×n(t+m)。

此时，协因数阵Q为奇异阵，但在计算时并没有直接求取Q的逆，因此并不影响解算过程。

本文在文献^[10]推导牛顿法思路的基础上，进一步推导了一般情形下多元总体最小二乘问题的牛顿解法。根据总体最小二乘问题平差准则V^TPV=min，由求条件极值的拉格朗日乘数法构造如下目标函数

点击浏览公式

式中，λ∈R^nm×1为拉格朗日乘数向量；$\otimes $为Kronecker-Zehfuss积。 对式(5)中的元素分别求导得

点击浏览公式

点击浏览公式

式中，矩阵$\sum $=($\overset{\Lambda }{\mathop \Xi }\,$^T$\otimes $I_n-I_m$\otimes $I_n)，这与文献^[10]的表达一致，符号“＾ ”表示最优估值。

由式(6)和式(7)，得到改正数最优估值为^[10]

点击浏览公式

进而得到牛顿法求解无约束问题的目标函数

点击浏览公式

式中

式(9)对vec(Ξ)求一阶导数，得

点击浏览公式

经推导可得式(10)第1部分

点击浏览公式

由逆矩阵求导原理，得式(10)第2部分

点击浏览公式

式中

点击浏览公式

点击浏览公式

式中，e_i=[0…010…0]^T为m维列向量，除了第i个元素为1外，其余的元素全部为0；i=[(j-1)/t]+1；[·]为取整运算；${\bar{Q}}$_j∈R^nm×n(t+m)。

由式(11)和式(12)得

点击浏览公式

点击浏览公式

并令

点击浏览公式

式中，i,j=1,2,…,tm；${\bar{Q}}$∈R^nmtm×nmtm；M∈R^tm×tm；sqrt(·)为开方运算。

则

点击浏览公式

依据牛顿法原理得到未知参数的改正数

点击浏览公式

以及

点击浏览公式

式中，vec($\overset{\wedge }{\mathop \Xi }\,$)为未知参数最佳估值；vec(Ξ⁰)为未知参数近似值。

多元总体最小二乘问题的牛顿解法的迭代步骤为(k为迭代次数)： 第1步，根据给定的A、Y、Q，利用式(14)和式(17)构造${\bar{Q}}$和${\bar{\bar{Q}}}$，计算最小二乘解

第2步，由上一步获得的数值，分别计算

由式(15)和式(18)，进一步计算G^k+1、H^k+1和

第3步，根据式(20)更新未知参数估值。

第4步，重复第2和3步，直到当

拉格朗日乘数法的多元总体最小二乘解为^[10]

点击浏览公式

点击浏览公式

式中，拉格朗日乘数法未知参数改正数为

点击浏览公式

令

点击浏览公式

利用矩阵反演公式^{[10, 19]}，得

点击浏览公式

又因为

点击浏览公式

联合式(23)、式(25)和式(26)，得

点击浏览公式

式中，I_tm∈R^tm×tm为单位矩阵；$R=D_{1}^{-1}{{(D_{1}^{-1}+D_{2}^{-1})}^{-1}}$。

由式(27)可知，牛顿法未知参数改正数与拉格朗日乘数法未知参数改正数的区别在于R矩阵。拉格朗日乘数法把平差准则的条件极值转化成式(5)的自由极值，对目标函数式(5)中的变量分别求一阶偏导数，并令其为零来求取目标函数的极小值；牛顿法判断极小值不仅需要目标函数对变量的一阶导数为零，还需要二阶导数大于零，即 Hessian 矩阵正定，这使得牛顿法对全局的判断能力更强，寻找极小值的速度也更快^[20]；vec($\Delta \overset{\wedge }{\mathop{\Xi }}\,$)_L是一阶导数意义下的未知参数改正数，vec($\Delta \overset{\wedge }{\mathop{\Xi }}\,$)_N是二阶导数意义下的未知参数改正数，目标函数的二阶导数部分组成R矩阵。一般情形下，当目标函数的二阶导数存在时，R矩阵也是存在的，Rvec($\Delta \overset{\wedge }{\mathop{\Xi }}\,$)_L可以理解为拉格朗日乘数法未知参数改正数的修正量，这也正是牛顿法比拉格朗日乘数法迭代效率高的原因。当R或Rvec($\Delta \overset{\wedge }{\mathop{\Xi }}\,$)_L为零时，牛顿法就等价于拉格朗日乘数法。 2 多元总体最小二乘解法精度评定 2.1 单位权方差估值公式

根据文献^{[6, 7]}，多元总体最小二乘解法的单位权方差估值公式为

点击浏览公式

式中，m(n-t)为自由度，mn为观测矩阵Y元素的个数，mt为参数矩阵Ξ元素的个数。 当系数矩阵含有常数元素时， $\mathop {{\rm{ }}V}\limits^ \wedge {\mkern 1mu} = \left[{\matrix{ {vec(\mathop {{E_{As}}}\limits^ \wedge {\mkern 1mu} )} \cr {vec(\mathop {{\rm{ }}E}\limits^ \wedge {{\mkern 1mu} _Y})} \cr } } \right]P = {\left[{\matrix{ {{Q_{{A_S}}}_{{A_S}}} & {{Q_{{A_S}}}_Y} \cr {{Q_{Y{A_S}}}} & {{Q_{YY}}} \cr } } \right]^{ - 1}}$；式(28)也与文献^[6]的式(3.7)和文献^[7]的式(37)相一致，但本式的表达具有一般性。

多元总体最小二乘法中，数据多以矩阵形式存在。为计算协因数阵，把矩阵写成向量形式，得到

点击浏览公式

式(29)将向量$\overset{\wedge }{\mathop V}\,$、$\overset{\wedge }{\mathop L}\,$和$\xi $近似写成L的函数，根据协因数传播律得到各个向量间的一阶近似协因数阵^[21]。本文求出的多元总体最小二乘问题的未知参数的协因数阵也与文献^{[19, 22, 23]}中给出的一般总体最小二乘问题的未知参数的协因数阵在表达形式上相一致。所有协因数阵如表 1所示，其中， ${{\overset{\wedge }{\mathop{X}}\,}_{3}}={{{\bar{A}}}^{T}}\overset{\wedge }{\mathop{X}}\,_{2}^{-1}\bar{A};{{\overset{\wedge }{\mathop{X}}\,}_{4}}=\overset{\wedge }{\mathop{X}}\,_{3}^{-1}{{{\bar{A}}}^{T}}\overset{\wedge }{\mathop{X}}\,_{2}^{-1};{{\overset{\wedge }{\mathop{X}}\,}_{5}}=Q\overset{\wedge }{\mathop{{{\sum }^{T}}}}\,\overset{\wedge }{\mathop{X}}\,_{2}^{-1}\overset{\wedge }{\mathop{\sum }}\,Q$

根据文献^[17]中给出的源坐标系的真值数据，本文通过给定的二维仿射变换参数真值计算得到目标坐标系的真值数据。在得到两组坐标系的真值数据后，依据文献^[17]，本文在模拟数据中源坐标系所加误差的协因数阵为

模拟数据中目标坐标系所加误差的协因数阵为

共模拟500次，分别使用本文提出的牛顿法(N法)和拉格朗日-牛顿法(L-N法)，以及Fang(2011)法^[10]、Schaffrin(2009)法^[8]和加权最小二乘法(WLS)计算二维仿射变换参数，迭代法的终止条件均为10^-12，将获得的二维仿射变换参数的平均值，以及二维仿射变换参数平均值与其真值的差值2范数和平均迭代次数列于表 2中。

参数为矩阵时，文献^[10]并没有给出在系数矩阵含有常数列时的具体公式算法。依据其提出的在参数为向量时系数矩阵含有常数列的加权总体最小二乘算法，本文编程实现了这种将系数矩阵分成两部分的多元加权总体最小二乘算法，并简称为Fang(2011)法。拉格朗日-牛顿法(L-N法)是指把Fang(2011)法(拉格朗日乘数法)获得的第一次多元加权总体最小二乘参数解作为本文牛顿法迭代的初始值；在选取迭代初始值时，如果使用多元最小二乘参数解作为初始值计算式(20)，则为牛顿法(N法)，如果使用多元加权总体最小二乘参数解作为初始值计算式(20)，则为拉格朗日-牛顿法(L-N法)；在能够获得多元加权总体最小二乘参数解的情况下，建议使用收敛速度更快的L-N法。

由表 2可知，N法、L-N法、Fang(2011)法和Schaffrin(2009)法这4种多元加权总体最小二乘法获得的参数结果一致，并且参数估值的差值范数要小于WLS法。从平均迭代次数来看，Schaffrin(2009)法的迭代次数最多，Fang(2011)法次之，L-N法的迭代次数要稍小于N法，为最小。本文提出的多元总体最小二乘牛顿解法的迭代效率要高于Fang(2011)法和Schaffrin(2009)法这两种解法。由文献^[20]可知，因为牛顿法是一种用H(Hessian)矩阵导出的最速下降法，牛顿步径也即vec($\Delta \overset{\wedge }{\mathop{\Xi }}\,$)_N是好的搜索方向，而且牛顿法也具有二次收敛速度，所以，牛顿法的迭代次数最少。本文提出的牛顿法比Fang(2011)法(拉格朗日乘数法)有着更高的收敛速度，它们的区别主要是由Rvec($\Delta \overset{\wedge }{\mathop{\Xi }}\,$)_L矩阵决定的，而且当${{\Xi }^{0}}$靠近$\overset{\wedge }{\mathop{\Xi }}\,$时，牛顿法有着更高的收敛速度，这也解释了为何L-N法比N法有着更少的迭代次数。

以上算例没有考虑系数矩阵元素与观测矩阵元素的相关性。为了验证本文的牛顿法能够解决更一般的多元总体最小二乘问题，在文献^[19]的基础上，使用Matlab 7.11.0 软件，利用函数randn(100,13*4)^T×randn(100,13*4)生成协因数阵Q_r，利用函数mvnrnd生成均值均为0，方差-协方差阵为0.001 Q_r的随机误差矩阵，并将其施加到真值数据中以作为计算的模拟数据。

由于模拟数据使得两套坐标系具有了相关性，此时Schaffrin(2009)法不再适用，选用本文提出的牛顿法和拉格朗日-牛顿法，以及Fang(2011)法模拟计算10 000次，计算时取Q=Q_r，迭代法的终止条件均为10^-12，并由获得的10 000次的参数估值绘制出的参数估值频率直方图及其相应的拟合出的正态分布曲线图如图 1所示。

图 1 所估参数的直方图及其正态分布拟合曲线 Fig. 1 Histogram and fitting curve of normal distribution of estimated parameters

图选项

计算结果表明，牛顿法的平均迭代次数为6.972 6，拉格朗日-牛顿法的平均迭代次数为6.828 7，Fang(2011)法的平均迭代次数为12.188 2。牛顿法和拉格朗日-牛顿法使用的协因数阵Q为奇异矩阵，但这并不影响未知参数的解算。Fang(2011)法的迭代次数最多，拉格朗日-牛顿法迭代次数要少于牛顿法，这也与上文的结论相符合。参数估值的频率直方图近似于正态分布，相应的正态分布拟合曲线大约也以参数的真值为对称轴。以上的结论证明了多元总体最小二乘牛顿解法的可行性和有效性。

选取文献^[19]中二维坐标转换数据作为计算数据，共有5组在源坐标系和目标坐标系下的公共点。 源坐标系和目标坐标系权阵为^[19]

使用本文提出的牛顿法和拉格朗日-牛顿法，以及Fang(2011)法计算二维仿射变换参数，迭代法的终止条件均为10^-12。3种方法的未知参数估值列于表 3中。利用表 3中的未知参数估值结果和表 1中的公式计算未知参数估值协因数阵，并将其列于表 4中。

表 3为3种方法的共同结果，其精度一致，且$\overset{\wedge }{\mathop{\sigma _{0}^{2}}}\,$为0.147 294 126 638。结果表明，牛顿法迭代了2次，拉格朗日-牛顿法迭代了1次，Fang(2011)法迭代了3次，这也说明了本文算法的合理性。

本文推导了多元总体最小二乘问题的牛顿解法，并比较了基于牛顿法的多元加权总体最小二乘解和基于拉格朗日乘数法的多元加权总体最小二乘解之间的关系，给出了多元总体最小二乘算法精度评定的16种协因数阵的近似计算公式。通过算例验证，并与现有的解法作了比较，得到以下几点结论： (1) 牛顿法有着更少的迭代次数，更高的解算效率。 (2) 在观测矩阵和系数矩阵元素对应的协因数阵中，本文取常数元素的协因数为0。这虽然导致协因数阵Q奇异，但并不影响参数的估计结果；本文的算法也能解决观测矩阵与系数矩阵元素具有相关性的问题。 (3) 牛顿法的解算效率会受到迭代初始值的影响，本文采用的拉格朗日-牛顿法因其迭代初始值更接近参数最优估值，而使得解算效率比起牛顿法有所提高。 本文给出的多元总体最小二乘解法精度评定的公式只是近似公式，如何推导更为严密的总体最小二乘精度评定公式还需进一步研究。