整体最小二乘法及其在测量数据处理中的应用

整体最小二乘法及其在测量数据处理中的应用丁克良[1]，欧吉坤[2]，陈义[2]1北京建筑工程学院测绘与城市空间信息学院；2中国科学院测量与地球物理研究所3同济大学测量与国土信息工程系DingKe-liang,OUJi-kun，Chenyi(1BeijingInstituteofCivilEngineeringAndArchitecture,SchoolofGeomaticsandurbaninformation,1000442InstituteofGeodesyandGeophysics,ChineseAcademyofSciences,Wuhan,4300773.DepartmentofSurveyingandGeo-Informatics,TongjiUniversity,Shanghai200092,China)摘要在参数求解中，针对参数估计模型的观测向量和系数矩阵都可能存在误差情况，20世纪80年代提出了整体最小二乘方法。本文系统阐述了整体最小二乘基本原理、思想，基于奇异值分解原理详细阐述了整体最小二乘的解算方法。在对普通最小二乘和整体最小二乘深入比较的基础上提出要注意整体最小二乘应用条件问题。关键词：普通最小二乘，整体最小二乘，奇异值分解；测量数据处理Abstract:Forparametersolution,thetotalleastsquare(TLS)methodwasproposedinearly80’softwentycenturyastodealwithsomeestimationmodelsthatboththeobservationanddesignmatrixaresubjecttoerrors.ThebasicconceptofTLS,principleofcomputationalalgorithmsforTLSbasedonthesingularvaluedecompositiontheoryaredescribedsystematicallyinthispaper.TheapplicabilityandlimitationofTLSinsurveyingdataprocessareproposedbasedoncomparisonsandanalysesbetweentheordinaryleastsquaresandthetotalleastsquares.Keywords:ordinaryleastsquares;totalleastsquares;singulardecomposition;surveyingdataprocess1绪论最小二乘法经历了百余年的发展考验，已经成为许多领域数据处理广泛应用的方法。测量数据的处理方法，通常是指按最小二乘法进行测量平差，它是测量数据处理中最基本、最广泛的应用方法[1]，尤其是近几十年来得到了充分的发展和应用。最小二乘平差的基本思想是在最小二乘准则下进行测量数据的调整。测量平差模型均可归结线性方程组=+AXL∆的求解问题。最小二乘准则要求残差的范数平方和极小，它主要是针对观测值中的偶然误差的。然而，实际问题中参数估计中的观测值和系数阵都可能存在误差，针对这种更复杂的情况，20世纪80年代提出了整体最小二乘法。尽管最初的称呼不同，总体最小二乘实际上已有相当长时间了。整体最小二乘的基本思想可以追溯到19世纪，较早的文献诸如Adcock（1878）、Pearson（1901）、Linnik（1961）、Williamson（1968）、York（1966）[2]。当时考399虑的是A和L同时存在误差时，矩阵方程AxL=的近似求解方法，但是，直到1980年Golub和VanLoan才从数值分析的观点首次对这种方法进行了整体分析，并正式称之为整体最小二乘（Golub，VanLoan，1980）[3]。VanHuffel和Vandewalle深入了研究整体最小二乘问题[4,5,6]。此后，整体最小二乘作为一种新的数据处理方法已成为一个研究热点，诸多学者对其进行了广泛研究[6,7,8,9,10,11]。目前、该方法已经成功应用于谱分析、参数估计、自动控制、图像处理、系统辨识、信号处理等相关领域。本文首先详细阐述了整体最小二乘的基本理论、思想、解算方法；结合系数矩阵和观测向量对参数求解的影响，提出了整体最小二乘应用条件问题，对其应用前景进行了讨论。2整体最小二乘基本理论2.1整体最小二乘基本思想对于线性方程组=AxL，普通最小二乘的基本思想是在残差平方和极小的准则约束下求解最佳参数。这里有一个前提，系数矩阵A是没有误差的精确值，但是多数情况系数阵A和观测向量L同时存在误差，若同时考虑二者的误差，此时，线性方程组可表示为[3]()AL+=+AExLE（2.1）其中mn×∈AR，m∈LR，n∈xR，()rankn=A,()ranknm=A；m为观测值个数，n为待估参数个数，AE为系数阵的噪声，LE为观测噪声，误差矩阵[]ALEE属于相互独立的白噪声误差。这一模型称为EIV（Errors-in-Variables）模型[4]。解决这类问题的适宜方法是整体最小二乘法（TotalLeastSquares,TLS）[3,4,5]。对于线性方程组=AxL，整体最小二乘问题就是在以下准则约束下[2,3](1)ˆˆ;ˆˆ[][]minmnFALR×+⎡⎤∈⎣⎦−ALAL（2.2）ˆˆ()R∈LA寻求ˆA、ˆL，任何满足ˆˆ=AxL（2.3）的ˆx均称为线性方程=AxL的整体最小二乘解。ˆˆ[][][]AL=−EEALAL为相应整体最小二乘改正数。式中，FM为Frobenius范数，简称为F范数。2.2整体最小二乘基本算法整体最小二乘的求解是通过奇异值分解来实现的。将线性=AxL改写为10TT⎡⎤⎡⎤−=⎣⎦⎣⎦ALx（2.4）400记增广矩阵[]=CAL，对增广矩阵C进行奇异值分解T=∑CUV其中()121,,,,nndiagσσσσ+∑=（2.5）1210nnσσσσ+≥≥≥≥≥则整体最小二乘解可由增广矩阵右奇异向量的最后一列1n+v得到，即整体最小二乘解为[3]1,1,11,11ˆ,,Tnnnnnvvv++++⎡⎤=−⎣⎦x（2.6）当A为列满秩时，整体最小二乘还有另一种解的形式[3,4,5,6]211n(I)TTtlsnσ−+=−xAAAL（2.7）2.3混合最小二乘法整体最小二乘的基本思想是同时考虑设计矩阵和观测向量的误差，而在许多情况下，设计矩阵的某一列或某几列是常数，如在直线拟合、曲面拟合、GPS非差定位等模型中都存在这种情况。因此，在这种情况下对A的不同列就应区别对待，与此相应的参数可分别采用最小二乘法和整体最小二乘法求解，简称为混合最小二乘[4]将线性方程=AxL表示为[]1122⎡⎤=⎢⎥⎣⎦xAALx（2.10）其中[]12=AAA，11mn×∈AR，22mn×∈AR12TTT⎡⎤=⎣⎦xxx，11n∈xR，22n∈xR，12nnn=+m为观测值个数，n为待估参数个数，1n、2n分别为1A、2A对应的参数个数，1A的元素为常数。和整体最小二乘相比混合最小二乘问题就是在[](21)222ˆˆ[]ˆˆminmnALRF×+∈⎡⎤−⎣⎦ALAL（2.11）准则下，寻求2ˆˆ⎡⎤⎣⎦AL，任何满足1122ˆˆˆˆˆˆ=+=AxAxAxL401的12ˆˆˆ[]TTT=xxx均称为混合最小二乘解。[]222ˆˆAL⎡⎤⎡⎤=−⎣⎦⎣⎦EEALAL为相应的混合最小二乘改正量。混合最小二乘解的求解基本思路是首先采用QR分解法[4]，或者约化的方法将系数矩阵分为常数部分和非常数部分，后者采用整体最小二乘法求解，后者采用普通最小二乘法求解。3整体最小二乘普通最小二乘的比较整体最小二乘和普通最小二乘的区别在与同时考虑观测向量和系数矩阵的误差，参数求解通过增广矩阵奇异值分解来实现，可以说整体最小二乘是普通最小二乘的进一步发展；混合最小二乘是一种更广泛的形式，为便于说明问题，表3.1给出了整体最小二乘、普通最小二乘和混合最小二乘的算法比较。表1普通最小二乘、整体最小二乘和混合最小二乘的比较项目LSTLSLS-TLS数学模型=+AxLV()AL+=+AExLE21122()Axx++=AAEL误差条件A无误差、L含误差A、L均含有误差1A无误差，2A、L均含有误差，约束准则22min−xAxL21minALF⎡⎤⎢⎥−⎣⎦xEE[](1)2222ˆˆ[]ˆˆminmnALRF×+∈⎡⎤−⎣⎦ALAL参数解1ˆ()TTls−=xAAAL211nˆ()TTtlsnσ−+=−xAAIAL221121200ˆ0ˆTTnnσ−+=−⎛⎞⎡⎤⎡⎤⎜⎟⎢⎥⎢⎥⎜⎟⎣⎦⎣⎦⎝⎠xAAALIx4算例分析抛开线性方程的具体意义，给定精确系数矩阵A0和参数X0，直接由方程L0=A0X0模拟观测向量的精确值L0。系数矩阵和观测向量的精确值为402011112213914121525⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦A，00.711.542.092.362.35⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦L，00.41.250.14X−⎡⎤⎢⎥=⎢⎥⎢⎥−⎣⎦构造增广矩阵00[]=0CAL，0()iσA、0()jσC，(1,2,31,,4)ij==分别表示矩阵A0、C0的奇异值。然后分别对系数矩阵、观测向量各元素模拟正态分布随机误差。模拟误差的情况分为仅观测向量受到随机误差的影响、仅系数矩阵受到随机误差的影响、系数矩阵及观测向量都受到误差的影响等三种情况。A和L分别为施加随机误差之后的系数阵和观测向量。为便于和混合最小二乘比较，系数矩阵的第一列元素均为精确值常数1，仅对系数阵的第2、第3列施加随机误差。分别对各种情况采用普通最小二乘、经典整体最小二乘、混合最小二乘计算。表1中，FA、F∆A分别表示系数矩阵及其整体最小二乘改正量的F范数，FC、F∆C为相应增广矩阵的F范数。矩阵范数比值/FF∆AA、/FF∆CC可以从一个侧面反映系数阵误差对参数求解结果的影响。由计算结果可知，对于精确的线性方程L0=A0X0，增广矩阵最小奇异值为零，普通最小二乘、整体最小二乘、混合最小二乘解完全一致。在观测向量施加随机误差之后，前述三种解法的参数解虽然接近，但是已有所差别。受观测向量误差的影响，增广矩阵的最小奇异值由零增大为0.096313，由式（2.7）可知，由于21nnσ+I的影响导致普通最小二乘，整体最小二乘的参数解产生差异。这里，系数矩阵不受误差影响，采用最小二乘法处理比较合理，从计算结果来看，普通最小二乘解更接近参数的真实值。仅当系数阵受到误差影响时，从计算结果看，三种解法差别较大，系数矩阵的误差对参数求解影响较大。仅系数矩阵受误差的影响，这种情况最不适合用普通最小二乘法处理，从计算结果来看，采用整体最小二乘法效果相对较好。第三种情况，系数矩阵、观测向量均受到随机误差的影响。这种情况适于整体最小二乘法处理，注意到系数矩阵的第一列为精确值，不受误差影响，因此采用混合最小二乘法求解比较合理。比较而言，普通最下二乘解的结果偏离真实值相对较远。从/FF∆AA、/FF∆CC范数的比值来看，系数矩阵、观测向量的误差越大，比值越大。而普通最小二乘解偏离真实值越来越远，整体最小二乘解的结果要相对好一些。但误差较小时，矩阵的扰动量小，增广矩阵的最小奇异值较小，三种方法的解算结果接近。在参数的求解中，我们一般采用普通最小二乘法，当系数矩阵误差对参数影响较小时，这种情况比较合理。但是，当系数矩阵误差较大时，矩阵的扰动对参数求解影响较大，就要顾及系数矩阵对参数求解的影响，这种情况采用整体最二乘法可望获得较好的结果。但是，应该看到，整体最小二乘、普通最小二乘算法的区别在于引入增广矩阵最小奇异值，系数矩阵、观测向量误差对增广矩阵最小奇异值的大小都有影响，但是二者对奇异值的大小的影响是不同的。当系数矩阵的扰动对增广矩阵最小奇异值大小贡献较大时，采用整体最小二乘法比较合理，如表中第三种情况；当系数矩阵的扰动对增广矩阵最小奇异403值影响较小，或