顾及设计矩阵误差的AR模型新解法

姚宜斌，熊朝晖，张 豹，张 良，孔 建

1.武汉大学测绘学院，湖北 武汉 430079；2.武汉大学地球空间环境与大地测量教育部重点实验室，湖北 武汉 430079；3.地球空间信息技术协同创新中心，湖北 武汉 430079；4.武汉大学中国南极测绘研究中心，湖北 武汉 430079

顾及设计矩阵误差的AR模型新解法

姚宜斌1,2,3，熊朝晖1，张 豹1，张 良1，孔 建4

1.武汉大学测绘学院，湖北 武汉 430079；2.武汉大学地球空间环境与大地测量教育部重点实验室，湖北 武汉 430079；3.地球空间信息技术协同创新中心，湖北 武汉 430079；4.武汉大学中国南极测绘研究中心，湖北 武汉 430079

在自回归模型求解中，设计矩阵和观测值均存在误差，传统的最小二乘法不能很好地解决这一问题。本文提出了一种顾及设计矩阵误差的AR模型新解法，通过引入虚拟观测值，使观测向量与设计矩阵不仅同源而且带误差的元素个数相同，然后通过对观测方程进行等价变换巧妙实现了在最小二乘框架下求解自回归问题。利用模拟数据及实测数据分别对新算法进行了内符合精度检验，并利用实测数据对新算法进行外符合精度检验，结果表明新算法得到的结果显著优于奇异值分解(singular value decomposition,SVD)解法及传统最小二乘解法,验证了算法的精度和有效性。

AR模型；设计矩阵误差；整体最小二乘；虚拟观测值；奇异值分解

平差的实际问题中常存在设计矩阵由含有误差的观测值构成，传统最小二乘方法在求解此类模型的过程中，因忽略设计矩阵中实际存在的误差而具有理论的不完备性。考虑设计矩阵误差，将所有观测值改正数平方和最小作为平差准则的求解方法最初由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模型，任意时刻观测值zm为自身最近t阶滞后项的线性组合[23]，如下

(1)

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

L=BX

(2)

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

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

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

1.1 算法推导

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

L+v=(B0+ΔB)(X0+x)=B0X0+

B0x+ΔBX0+ΔBx

(3)

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

v=B0X0+B0x+ΔBX0-L

(4)

1.1.1 矩阵等价变换

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

(5)

1.1.2 引入虚拟观测值

(6)

v=B0X0+B0x+B10xB+B20v-L

(7)

1.1.3 函数模型构建

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

v=(E-B20)-1B0x+(E-B20)-1B10xB-(E-B20)-1(L-B0X0)

(8)

虚拟观测值误差方程为

vB=xB

(9)

如果令

则有

vg=Bgxg-lg

(10)

1.2 解算流程

本文方法解算流程为：

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

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

1.3 精度评定

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

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

t+1,t,t-1,…,2,1)

(11)

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

(12)

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

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

2 算例分析

2.1 模拟数据

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

表1 AR(3)模型模拟数据Tab.1 Simulation date of AR(3) model

zt=φ1zt-1+φ2zt-2+φ3zt-3+et

(13)

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

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

表23种不同方法所求参数值及单位权中误差
Tab.2Theparametervalueandtheerrorofunitweightofthreemethods

参数SVDLS本文方法真值φ10．83630．77090．80430．8φ2-0．5461-0．4667-0．5084-0．5φ3-0．2725-0．3246-0．3068-0．3σ01．44551．41600．9042

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

2.2 实测沉降数据

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

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

表33种方法所求沉降数据参数估值
Tab.3Theparametervalueofobservationdataofsettlementofthreemethods

参数SVDLS本文方法φ1-0．37480．0411-0．4863φ20．28980．32780．4264φ31．09010．63511．0670σ00．89540．76730．7320

表4 沉降数据实测值及3种方法平差值Tab.4 Adjustment of observation of three method and the original data mm

图1 3种方法平差结果与原始观测值比较Fig.1 Comparative results between original data and adjustment of observation of three value

2.3 预报精度分析

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

表5 3种方法预报结果Tab.5 Forecast result of three method

图2 3种方法预测结果与实测沉降观测值比较Fig.2 Comparative results between original data and forecast result of observation of three value

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

3 结 语

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

[1] ADCOCK R J.Note on the Method of Least Squares[J].Analyst，1877，4(6)：183-184.

[2] 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-893.

[3] VAN HUFFEL S，ZHA Hongyuan.An Efficient Total Least Squares Algorithm Based on a Rank-revealing Two-sided Orthogonal Decomposition[J].Numerical Algorithms，1993，4(1)：101-133.

[4] 姚宜斌，孔建.顾及设计矩阵随机误差的最小二乘组合新解法[J].武汉大学学报(信息科学版)，2014，39(9)：1028-1032.

YAO Yibin，KONG Jian.A New Combined LS Method Considering Random Errors of Design Matrix[J].Geomatics and Information Science of Wuhan University，2014，39(9)：1028-1032.

[5] 曾文宪，方兴，刘经南，等.附有不等式约束的加权整体最小二乘算法[J].测绘学报，2014，43(10)：1013-1018.DOI：10.13485/j.cnki.11-2089.2014.0173.

ZENG Wenxian，FANG Xing，LIU Jingnan，et al.Weighted Total Least Squares Algorithm with Inequality Constraints[J].Acta Geodaetica et Cartographica Sinica，2014，43(10)：1013-1018.DOI：10.13485/j.cnki.11-2089.2014.0173.

[6] SCHAFFRIN B，WIESER A.On Weighted Total Least-squares Adjustment for Linear Regression[J].Journal of Geodesy，2008，82(7)：415-421.

[7] SCHUERMANS M，MARKOVSKY I，VAN HUFFEL S.An Adapted Version of the Element-wise Weighted Total Least Squares Method for Applications in Chemometrics[J].Chemometrics and Intelligent Laboratory Systems，2007，85(1)：40-46.

[8] VAN HUFFEL S，VANDEWALLE J.Analysis and Properties of the Generalized Total Least Squares Problem AX≈B When Some or All Columns in A are Subject to Error[J].SIAM Journal on Matrix Analysis and Applications，1989，10(3)：294-315.

[9] XU Peiliang，LIU Jingnan，SHI Chuang.Total Least Squares Adjustment in Partial Errors-in-variables Models：Algorithm and Statistical Analysis[J].Journal of Geodesy，2012，86(8)：661-675.

[10] FANG X.Weighted Total Least Squares Solutions for Applications in Geodesy[D].Hannover，Germany：Leibniz University，2011.

[11] 方兴，曾文宪，刘经南，等.三维坐标转换的通用整体最小二乘算法[J].测绘学报，2014，43(11)：1139-1143.DOI：10.13485/j.cnki.11-2089.2014.0193.

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-1143.DOI：10.13485/j.cnki.11-2089.2014.0193.

[12] 赵俊，归庆明.部分变量误差模型的整体抗差最小二乘估计[J].测绘学报，2016，45(5)：552-559.DOI：10.11947/j.AGCS.2016.20150374.

ZHAO Jun，GUI Qingming.Total Robustified Least Squares Estimation in Partial Errors-in-variables Model[J].Acta Geodaetica et Cartographica Sinica，2016，45(5)：552-559.DOI：10.11947/j.AGCS.2016.20150374.

[13] 王乐洋，余航，陈晓勇.Partial EIV模型的解法[J].测绘学报，2016，45(1)：22-29.DOI：10.11947/j.AGCS.2016.20140560.

WANG Leyang，YU Hang，CHEN Xiaoyong.An Algorithm for Partial EIV Model[J].Acta Geodaetica et Cartographica Sinica，2016，45(1)：22-29.DOI：10.11947/j.AGCS.2016.20140560.

[14] 陶叶青，高井祥，姚一飞.基于中位数法的抗差总体最小二乘估计[J].测绘学报，2016，45(3)：297-301.DOI：10.11947/j.AGCS.2016.20150234.

TAO Yeqing，GAO Jingxiang，YAO Yifei.Solution for Robust Total Least Squares Estimation Based on Median Method[J].Acta Geodaetica et Cartographica Sinica，2016，45(3)：297-301.DOI：10.11947/j.AGCS.2016.20150234.

[15] 吴富梅，杨元喜.基于高阶AR模型的陀螺随机漂移模型[J].测绘学报，2007，36(4)：389-394.

WU Fumei，YANG Yuanxi.Gyroscope Random Drift Model Based on the Higher-order AR Model[J].Acta Geodaetica et Cartographica Sinica，2007，36(4)：389-394.

[16] 潘国荣，刘大杰.顾及邻近点变形因素项的动态模型辨识及预测[J].测绘学报，2001，30(1)：32-35.

PAN Guorong，LIU Dajie.Dynamic Modeling Identification and Predication in Consideration of the Adjacent Point Deformation[J].Acta Geodaetica et Cartographica Sinica，2001，30(1)：32-35.

[17] 杨元喜，崔先强.动态定位有色噪声影响函数——以一阶AR模型为例[J].测绘学报，2003，32(1)：6-10.

YANG Yuanxi，CUI Xianqiang.Influence Functions of Colored Noises on Kinematic Positioning：Taking the AR Model of First Class as an Example[J].Acta Geodaetica et Cartographica Sinica，2003，32(1)：6-10.

[18] 叶志伟，尹晖，张守建.AR模型谱在超导重力数据信号检测中的分析研究[J].武汉大学学报(信息科学版)，2007，32(6)：536-539.

YE Zhiwei，YIN Hui，ZHANG Shoujian.Using AR Model Spectrum Algorithms to Detect Superconducting Gravimetric Signals[J].Geomatics and Information Science of Wuhan University，2007，32(6)：536-539.

[19] 张昊，王琪洁，朱建军，等.对钱德勒参数进行时变修正的CLS+AR模型在极移预测中的应用[J].武汉大学学报(信息科学版)，2012，37(3)：286-289.

ZHANG Hao，WANG Qijie，ZHU Jianjun，et al.Application of CLS+AR Model Polar Motion to Prediction Based on Time-varying Parameters Correction of Chandler Wobble[J].Geomatics and Information Science of Wuhan University，2012，37(3)：286-289.

[20] 王乐洋，许才军，鲁铁定.边长变化反演应变参数的总体最小二乘方法[J].武汉大学学报(信息科学版)，2010，35(2)：181-184.

WANG Leyang，XU Caijun，LU Tieding.Inversion of Strain Parameter Using Distance Changes Based on Total Least Squares[J].Geomatics and Information Science of Wuhan University，2010，35(2)：181-184.

[21] 魏二虎，殷志祥，李广文，等.虚拟观测值法在三维坐标转换中的应用研究[J].武汉大学学报(信息科学版)，2014，39(2)：152-156.

WEI Erhu，YIN Zhixiang，LI Guangwen，et al.On 3D Coordinate Transformations with Virtual Observation Method[J].Geomatics and Information Science of Wuhan University，2014，39(2)：152-156.

[22] 姚宜斌，黄书华，孔建，等.空间直线拟合的整体最小二乘算法[J].武汉大学学报(信息科学版)，2014，39(5)：571-574.

YAO Yibin，HUANG Shuhua，KONG Jian，et al.Total Least Squares Algorithm for Fitting Spatial Straight Lines[J].Geomatics and Information Science of Wuhan University，2014，39(5)：571-574.

[23] CRYER J D，CHAN K S.时间序列分析及应用[M].潘红宇，译.北京：机械工业出版社，2011.

CRYER J D，CHAN K S.Time Series Analysis with Applications in R[M].PAN Hongyu，tran.Beijing：China Machine Press，2011.

[24] 姚宜斌，黄书华，陈家君.求解自回归模型参数的整体最小二乘新方法[J].武汉大学学报(信息科学版)，2014，39(12)：1463-1466.

YAO Yibin，HUANG Shuhua，CHEN Jiajun.A New Method of TLS to Solving the Autoregressive Model Parameter[J].Geomatics and Information Science of Wuhan University，2014，39(12)：1463-1466.

[25] 王新洲，陶本藻，邱卫宁，等.高等测量平差[M].北京：测绘出版社，2013.

WANG Xinzhou，TAO Benzao，QIU Weining，et al.Advanced Surveying Adjustment[M].Beijing：Mapping Publishing Company，2013.

A New Method to Solving AR Model Parameters Considering Random Errors of Design Matrix

YAO Yibin1,2,3，XIONG Zhaohui1，ZHANG Bao1，ZHANG Liang1，KONG Jian4

1.School of Geodesy and Geomatics，Wuhan University，Wuhan 430079，China; 2.Key Laboratory of Geospace Environment and Geodesy，Ministry of Education,Wuhan University，Wuhan 430079，China; 3.Collaborative Innovation Center for Geospatial Technology，Wuhan 430079，China; 4.Chinese Antarctic Center of Surveying and Mapping，Wuhan 430079，China

The ordinary least square method could not solve the problem that the error exist both in design matrix and observation vector while compute parameter values of AR model.In this article, a new method is proposed which consider the random errors of design matrix.The source of design matrix and observation vector is same and the amount of parameters contain error can be equal by introducing virtual observation.Then, this problem could be solved under the framework of normal least square by equivalence transformation of observation equation.The result of this new method is superior to SVD method and normal least square method by simulation date and observation data which verify the feasibility and effectiveness of this method.

AR model；design matrix error；TLS；virtual observations；SVD method

The General Program of National Natural Science Foundation of China(Nos.41274022；41574028)；Natural Science Foundation for Distinguished Young Scholars of Hubei Province of China(No.2015CFA036)

YAO Yibin(1976—),male,professor,majors in geodetic data processing,GNSS space environment science,etc.

XIONG Zhaohui

姚宜斌，熊朝晖，张豹，等.顾及设计矩阵误差的AR模型新解法[J].测绘学报，2017,46(11)：1795-1801.

10.11947/j.AGCS.2017.20170004.

YAO Yibin，XIONG Zhaohui，ZHANG Bao，et al.A New Method to Solving AR Model Parameters Considering Random Errors of Design Matrix[J].Acta Geodaetica et Cartographica Sinica,2017,46(11):1795-1801.DOI:10.11947/j.AGCS.2017.20170004.

P228

A

1001-1595(2017)11-1795-07

国家自然科学基金(41274022;41574028)；湖北省杰出青年科学基金(2015CFA036)

(责任编辑：宋启凡)

2017-01-03

修回日期：2017-08-18

姚宜斌(1976—)，男，教授，主要从事测量数据处理理论与方法、GNSS空间环境学研究。

E-mail：ybyao@sgg.whu.edu.cn

熊朝晖

E-mail：cehui_xiong@whu.edu.cn