总体最小二乘求解线性模型的一种新算法*

2014-04-14冷顺绿杨根新杨新云

地矿测绘 2014年1期

冷顺绿，杨根新，杨新云

(1．云南省地图院，云南昆明 650034;2．云南国土资源职业学院，云南昆明 650217;3．云南省地矿测绘院，云南昆明 650218)

0 引言

在测量数据处理中，常常采用最小二乘法来进行线性模型的参数估计。这是在不考虑系数矩阵含有误差而只考虑观测向量含有误差的情况下得到的参数估值，即理论上的最优值。但当线性模型中系数矩阵也含有误差时，应采用总体最小二乘法来求取参数估值。Golub 等人提出运用奇异值分解法(SVD)来求取总体最小二乘参数估值［1］，但SVD 的复杂性限制了总体最小二乘在测量数据处理中的应用［3，4］。鉴于此，文献［3，4］推导了总体最小二乘的迭代算法。本文在此基础上，根据线性间接平差模型，运用拉格朗日原理推导了一种线性模型的总体最小二乘算法。并以文献［5］中实例数据加以验证，其结果与文献［4］中的迭代算法计算的结果完全一致，说明该算法正确合理。

1 算法推导

总体最小二乘线性模型为:

式中:L 为m ×1 的观测向量;V 为m ×1 的观测值改正向量;A为m × n 的系数矩阵;EA为系数改正矩阵;^x 为n ×1 的参数估值向量。要求取最佳参数估值，即根据总体最小二乘原理要满足如下条件:

式中:VA= vec(EA)，为mn ×1 的向量。其中，vec 为矩阵拉值运算，vec(EA)即将矩阵EA从左至右逐列拉直成mn × 1 的向量。构造拉格朗日目标函数如下:

式中:K 为m ×1 的拉格朗日乘数向量;EA^x = (^xT⊗I)TVA，其中⊗为矩阵的Kronecker 积，I 为m 阶单位矩阵。根据拉格朗日目标函数求极值原理，分别对V、VA、^x、K 求导并令其等于0，化简可得:

根据式(4)前两式再联立式(1)可得:

式中:F = (^xT⊗I)。

联立式(5)与式(4)的第3 式，可得:

式中:EA= reshape(- FTK)，其中reshape 为重构矩阵运算，将一列元素重新构成m × n 的矩阵，即为vec 的反运算。根据式(6)即可得到参数^x 的表达式:

由上述推导可知:

则单位权中误差计算公式为:

式中:m 为观测方程的个数;n 为必要观测数。由上述推导，总体最小二乘的参数^x 估值可以采用迭代法，其具体步骤如下:

1)对参数赋初值为^x(0);

2)按式(10)计算F、K 值和新的参数估值:

4)输出参数估值，按式(9)求得单位权中误差。

2 实例分析

为检验本文迭代算法的正确性，采用文献［5］中表3．18 的数据，见表1。拟合一曲线函数y = ax3+bx2+cx+d，曲线函数的参数真值为1、-2、-3、4。分别采用最小二乘法、文献［5］TLS 法、文献［4］TLS 法和本文所述方法来进行计算。拟合结果，如表2 所示。

表1 曲线拟合数据Tab．1 Curve fitting data

表2 不同方法拟合结果Tab．2 The fitting results of different methods

从表2 的拟合结果可以看出，由于总体最小二乘迭代算法考虑了系数矩阵的误差，所求得的参数值与真值较为接近，拟合精度比采用最小二乘法精度高。采用本文方法拟合的结果与采用文献［4］方法拟合的结果相同，得到的参数值与文献［5］计算的结果较为接近，且拟合精度较高，说明本文所提出的总体最小二乘迭代算法正确合理。