平差的实际问题中常存在设计矩阵由含有误差的观测值构成，传统最小二乘方法在求解此类模型的过程中，因忽略设计矩阵中实际存在的误差而具有理论的不完备性。考虑设计矩阵误差，将所有观测值改正数平方和最小作为平差准则的求解方法最初由Adcock于1877年提出^[1]。1980年，Golub和Van Loan提出了奇异值分解(SVD)算法，并将这类设计矩阵带有误差的问题命名为整体最小二乘方法(total least square, TLS)^[2]。使用SVD解法虽能获得稳定数值解，但计算复杂度高^[3]，如果设计矩阵中含有误差的元素重复出现，使用SVD解法，同一观测值在不同位置改正数将不一样。另外，使用SVD解法估计参数，难以给出以观测值表达的未知参数平差值表达式，从而无法利用协方差传播率进行精度评定^[4]。近些年，文献[4]提出一种顾及设计矩阵随机误差的最小二乘组合新解法，其将设计矩阵中所有含误差的元素作为虚拟观测值引入虚拟观测方程求解模型，解决了TLS方法的精度评定问题。加权整体最小二乘亦有较大进展^[5-10]。此外，文献[11]研究了三维坐标转换的通用整体最小二乘解法；文献[12]对部分变量误差模型提出了整体抗差最小二乘估计方法；文献[13]对Partial EIV模型提出了一种新解法；文献[14]提出了一种基于中位数法的抗差整体最小二乘估计方法。

AR模型是一类基本的时间序列模型，在工程实践与导航定位等方面有着广泛应用^[15-22]。对于t阶AR模型，任意时刻观测值z_m为自身最近t阶滞后项的线性组合^[23]，如下

式(1)用矩阵形式表达可写为

式中，观测向量L = [z_t+1 z_t+2 … z_t+n]^T，令任意时刻观测值z_i的改正数为v_i；设计矩阵；未知参数向量X =[φ₁ φ₂ … φ_t]^T。

由此可见，含误差的观测值规律地出现在设计矩阵不同位置，且观测向量中元素亦重复出现在设计矩阵中。使用传统最小二乘方法估计AR模型参数，因忽略设计矩阵中存在含误差的元素而理论缺乏严密性。AR模型中，观测向量误差与设计矩阵误差同源，若使用SVD方法估计AR模型参数，则多次出现的同一观测值的改正数不再唯一。文献[24]提出一种非线性的AR模型整体最小二乘迭代算法，但难以应用协方差传播定律给出精度评定公式。

前述方法对求解AR模型存在一定缺陷，本文提出一种顾及设计矩阵误差，适用于AR模型的整体最小二乘新解法，该方法通过引入未在观测向量中出现且含误差的观测值作为虚拟观测值，将设计矩阵对应的改正数与未知参数初值乘积改写为由新设参数改正数及实际观测值改正数表达的线性组合，使得观测向量与设计矩阵中带误差的元素个数相同。由此将整体最小二乘问题转化为经典最小二乘问题，最终使用间接平差方法进行参数求解。本文方法有效地克服了AR模型中同一参数在不同位置改正数不一致的情况，并能依据协方差传播定律进行精度评定。

1 顾及设计矩阵误差的AR模型新解法 1.1 算法推导

同时考虑设计矩阵B与观测向量L误差，式(2)可写为

忽略二阶小项Δ Bx，式(3)可改写为

式中，v = [v_t+1 v_t+2 … v_t+n]^T为观测向量L改正数；X₀= φ₁₀ φ₂₀ … φ_t0]^T，由待求参数初值组成，x为未知参数改正数；B₀由设计矩阵B中各元素初值构成；，由设计矩阵各元素改正数构成。

1.1.1 矩阵等价变换

矩阵Δ B中改正数规律地重复出现，为将不同的改正数单独取出，同时考虑未在设计矩阵中出现的改正数v_t+n，将Δ BX₀作如式(5)所示等价变换

1.1.2 引入虚拟观测值

在AR模型中，观测向量L中元素在设计矩阵B中重复出现，两者相同元素改正数应一致。考虑到B中的[z₁ z₂ … z_t]^T未在L中出现，为使观测向量L与设计矩阵B中带误差的元素个数相同，引入仅在设计矩阵B中出现且带误差的元素作为虚拟观测值L_B，即L_B=[z₁ z₂ … z_t]^T。令x_B=[v₁ v₂ … v_t]^T，由此Δ BX₀可由x_B与观测值改正数v组成的列向量线性表示，即式(6)

式中，；，将式(6)结果代入式(4)可得

1.1.3 函数模型构建

将式(7)中右边B₂₀v项移至等式左边，提取公因式v得(E－B₂₀)v。由于矩阵B₂₀中非零元素集中出现在主对角线下方，主对角线及其上方元素全为0，故行列式|E－B₂₀|恒为1，即矩阵E－B₂₀可逆，从而可得式(8)

虚拟观测值误差方程为

如果令

则有

式中，v_g为z₁、z₂、…、z_t、z_t+1、…、z_t+n改正数构成的列向量；B_g为等效设计矩阵；x_g为所有实际待估参数的和虚拟参数的改正数。算法的估计准则为v_g^TPv_g=min，由此整体最小二乘问题转化为经典最小二乘问题，从而可按经典平差中的间接平差方式求解。

1.2 解算流程

本文方法解算流程为：

(1) 确定观测值及虚拟观测值初值，选用普通最小二乘所得的AR系数作为未知参数初值，即X₀=(B₀^TB₀)^－1B₀^TL。

(2) 按照式(10)所定义的各矩阵，将AR系数初值或更新值代入，构建等效设计矩阵B_g，矩阵l_g，按间接平差原理计算x_g= B_g^TPB_g^-1B_g^TPl_g，同时更新未知参数X。

(3) 根据式(10)计算v_g，更新实际观测值L估值及虚拟观测值L_B估值。

(4) 重复步骤(2)、(3)直到未知参数改正数x小于一定阈值停止迭代，得到最终的未知参数结果。

1.3 精度评定

在经典最小二乘中，权由观测值精度确定，而观测值出现的次数，则会在平差过程中，通过对系数矩阵的作用，进一步影响法方程进而影响平差结果。但是在整体最小二乘中，系数矩阵中的观测值具有双重身份，除了影响系数矩阵本身外，还作为观测值参与平差，本文认为观测值这部分影响与观测值在系数矩阵中出现的次数也有关系。本文所定义的权阵P由两部分构成，即由观测值精度信息构建的权阵P ′，以及一个与次数有关的因子阵K，其中P=P′K。为简化问题，本文先行假定AR模型各期观测值等精度，即P ′为单位阵，P=K，若观测值不等精度，则在按照精度信息得到的权阵P ′的基础上再乘以一个次数有关的因子阵K。

由于在AR模型中，设计矩阵中元素在不同位置中规律地重复出现，观测值亦重复出现于设计矩阵中，在使用整体最小二乘方法进行求解未知参数时，实际上将同一观测值当作多次相同观测建模，为此，在各期观测值等精度的条件下，模型各观测值权可由观测值出现的次数确定，考虑引入虚拟观测值，该算法权阵P定义如下

使用上述方法经过若干次迭代可以获得未知参数X平差值，单位权中误差σ计算如式(12)所示

式中，模型自由度f=(n+u)－(t+u)=n－t。

在本文算法中，参数X₀=[φ₁₀ φ₂₀ … φ_t0]^T的初值为最小二乘所得到的AR系数。需要注意的是，若初值精度不够，则会导致迭代次数多，收敛速度慢，初值精度过差会导致迭代不能收敛。

2 算例分析 2.1 模拟数据

为验证算法的可行性与精度，依据式(13)模拟一个平稳的AR(3)模型，其AR系数分别取为φ₁= 0.8、φ₂=-0.5、φ₃=-0.3, 模拟数据如表 1所示。

式中，e_t为高斯白噪声，表 1模拟数据所附加噪声e_t~N(0, 1)。

对于表 1模拟数据，分别采用本文方法(迭代阈值为0.000 000 1)、SVD方法以及直接最小二乘3种方法求解参数及单位权中误差，列于表 2。

从表 2可知，本文方法所得到的参数估值较SVD方法和传统最小二乘方法更接近于真值，由于考虑了设计矩阵误差及同一观测值出现在不同位置改正数应一致等事实，本文方法的单位权中误差明显小于最小二乘方法和SVD方法，并且更加接近虚拟观测噪声的中误差。当对模拟数据不附加高斯白噪声时，本文方法所得参数估值将与真值一致。

2.2 实测沉降数据

表 2所示结果验证了算法的可行性与精度。为更好地验证本文算法可行性，采用文献[25]中例3-7实测沉降观测数据，该组实测沉降数据可用AR(3)模型拟合，故对实测沉降数据建立AR(3)模型，φ₁、φ₂、φ₃为待求AR系数。在实测沉降数据中，本文算法所选用的参数初值X₀= [φ₁₀ φ₂₀ φ₃₀]^T为最小二乘方法所获得的AR系数，在获取参数初值之后，结合实测沉降观测值代入式(10)，按照各矩阵定义构造模型，迭代求解以获得最终结果。

本文算法、普通最小二乘(LS)以及SVD 3种方法所得参数估值如表 3所示。实测观测数据及3种方法所得平差值如表 4所示。图 1所示为本文算法、普通最小二乘(LS)以及SVD 3种方法平差结果与原始观测值残差绝对值对比。

表 3中本文方法所求结果与SVD法较为接近，但本文方法所得单位权中误差小于SVD法，直接最小二乘结果与本文结果相差较远，查阅文献[23]知，本文所引用实测沉降数据并非平稳的自回归时间序列，这可能是和直接最小二乘结果估值结果相差较大的原因。由表 4及图 1可知，本文方法所得观测值的平差值更契合原始观测序列变化趋势。新方法看起来误差相对平稳，波动小一些。第14期的结果，其他方法较差，新方法较好，说明新方法由于理论严密，所以对于抑制较大误差效果更好。

2.3 预报精度分析

AR模型能利用前若干期观测值预报后若干期观测值，为检验本文算法的外符合精度，以表 4所列前30期实测沉降观测值为基础，分别利用本文算法、普通最小二乘(LS)以及SVD 3种方法求解模型系数并预报第31期至36期沉降数据。依据前30期实测沉降数据所求参数结果及相应方法预报值如表 5所示。图 2所示为本文算法、普通最小二乘(LS)以及SVD 3种方法预报结果与实测沉降观测值残差绝对值对比。

由表 5及图 2可知，在第31期至36期数据中，由SVD方法所求参数的预测效果最差，本文算法整体上预测效果最优，在第34期预报中3种方法预测效果均差，可能是观测噪声较大所致。向后预报期数越小，本文算法所求系数的预报结果与普通最小二乘方法相比，更接近与实际观测值，而在较远期预报中，二者效果相当。总体而言，在本文算法、普通最小二乘(LS)以及SVD 3种方法中，本文算法的外符合精度最优。

3 结语

在AR模型中，观测向量与设计矩阵的误差同源，观测值规律地重复出现在设计矩阵中，且设计矩阵中自身元素亦重复出现。传统最小二乘方法难以解决系数矩阵及观测值向量皆带误差的数据处理问题。本文提出了一种考虑设计矩阵误差的AR模型整体最小二乘新解法，引入仅在设计矩阵中出现且含误差的元素作为虚拟观测值，使观测向量与设计矩阵中带误差的元素个数相同。然后巧妙地对观测方程进行等价变换，实现了最小二乘框架下求解整体最小二乘问题，有效地克服了传统SVD方法的理论缺陷且能应用协方差传播定律给出未知参数的精度评定公式。最后通过对模拟数据及实测数据的验证，发现本文方法具有比SVD方法及经典最小二乘方法更高的精度及更优的外符合精度。