部分变量误差模型的整体抗差最小二乘估计

赵　俊,归庆明

1. 信息工程大学地理空间信息学院，河南 郑州 450001; 2. 信息工程大学理学院，河南 郑州 450001



部分变量误差模型的整体抗差最小二乘估计

赵俊1,归庆明2

1. 信息工程大学地理空间信息学院，河南 郑州 450001; 2. 信息工程大学理学院，河南 郑州 450001

Foundationsupport：TheNationalNaturalScienceFoundationofChina(Nos. 41174005; 41474009)

部分变量误差模型(partialEIVmodel)的加权整体最小二乘(weightedtotalleast-squares,WTLS)估计不具备抵御粗差的能力。鉴于粗差可能同时出现在观测值和系数矩阵中，本文在提出部分变量误差模型WTLS估计的两步迭代解法的基础上，运用抗差M估计的等价权方法，发展了一种整体抗差最小二乘(TRLS)估计方法，并采用一致最大功效统计量确定降权因子。针对WTLS估计两步迭代解法的特点，设计了两个不同的降权方案：第1个方案是在估计系数矩阵元素时，不对观测值降权，仅对系数矩阵降权；第2个方案是在估计系数矩阵元素时，既对系数矩阵降权，同时也对观测值降权。通过对模拟2D仿射变换和线性拟合实例进行计算和分析，结果表明第1方案优于第2方案，并且优于基于残差和验后单位权方差的抗差估计和现有的变量误差模型抗差估计。

部分变量误差模型; 加权整体最小二乘估计；粗差；整体抗差最小二乘估计；一致最大功效统计量

最近几年，作为变量误差模型(errors-in-variables,EIV)的参数估计方法，整体最小二乘(totalleast-squares,TLS)或加权整体最小二乘(weightedTLS,WTLS)方法在大地测量领域获得了广泛关注[1-2]，但大多数学者的工作主要集中于WTLS估计的算法研究[3-8]。WTLS估计是在一般的变量误差模型中引入观测误差和系数矩阵误差的方差-协方差矩阵的参数估计方法[3]。事实上，变量误差模型中并不是所有的系数矩阵元素都是随机的。例如，在坐标变换中，一些元素不含误差，可以认为是固定的，而另外一些元素则是由观测值所组成的，可以认为是随机的。顾及这一事实，通常考虑的变量误差模型被推广到更一般的部分变量误差模型[9]。部分变量误差模型的WTLS估计方法，其优势在于不用拉格朗日乘子且法方程未知参数的个数大大减少，因而获得了更多关注[10-11]。但是，WTLS估计与经典的LS估计一样，对粗差很敏感，不具有抵御粗差影响的能力[12-17]。

如何有效地处理观测数据中的粗差一直是测量领域关注和研究的热点[18-19]，抗差估计是一种处理粗差的有效方法，并且针对Gauss-Markov模型的抗差估计成果较多[20-26]。虽然针对变量误差模型，文献[14—17]提出了为数不多的几个抗差估计，但这些估计均采用残差和验后单位权方差确定降权因子，易受粗差的影响，进而影响其抗差功效。由于部分变量误差模型的系数矩阵也可能含有粗差，所以经典的抗差估计方法并不能直接应用于部分变量误差模型，更为有效的方法值得研究。为此，本文将在提出部分变量误差模型WTLS估计的两步迭代解法的基础上，运用抗差LS估计的等价权方法，提出一种整体抗差最小二乘(TRLS)估计方法，并构建一致最大功效统计量确定降权因子。针对WTLS估计两步迭代解法的特点，设计两个不同的降权方案：第1方案是在估计系数矩阵元素时，不对观测值降权，仅对系数矩阵降权；第2方案是在估计系数矩阵元素时，既对系数矩阵降权，同时也对观测值降权。

1　WTLS估计的两步迭代解法

考虑如下部分变量误差模型[9]

(1)

(2)

下面提出部分变量误差模型(1)的WTLS估计的两步迭代解法。首先，将任意给定的X的初值X(0)代入式(1)，得到

(3)

(4)

再令

(5)

(6)

对应的残差向量为

(7)

(8)

(9)

L=A*X+Δ

(10)

再次应用LS原理[19]，可得未知参数向量X的估计为

(11)

对应的残差向量为

(12)

最后，可得单位权方差σ2的估计为

(13)

值得指出的是，上述WTLS估计的两步迭代解法本质上与现有算法[9]等价。鉴于上述公式是在LS框架下推导得到的，所以经典LS估计和抗差LS估计的相关理论和方法[21-24]都可以直接应用于部分变量误差模型。

2　整体抗差最小二乘(TRLS)方法

(14)

和

(15)

根据基于敏感度分析的抗差估计理论[26]，本文将按下列方式确定式(14)、(15)中的等价权阵

(16)

式中，ωi为降权因子，可通过Huber权函数[20]进行确定

(17)

(18)

(19)

式中，m为给定的常数。

众所周知，残差和单位权方差的验后估计易受粗差的影响[26]；因此，按式(19)确定降权因子，会造成抗差功效的损失。基于敏感度分析的抗差估计理论表明，基于一致最大功效统计量的抗差估计优于基于残差或标准化残差的抗差估计[26]。令人遗憾的是，通过现有的WTLS估计方法均无法构造此一致最大功效统计量，为此本文利用第1节提出的WTLS估计的两步迭代解法，构造如下统计量，以将基于一致最大功效统计量的抗差估计拓展到部分变量误差模型

(20)

(21)

Fig. 3 is the chip photo of the 245 GHz 2nd subharmonic receiver with on-chip antenna with each circuit block indicated in the figure.

(22)

(23)

对应的降权因子分别为

(24)

(25)

(26)

(27)

注意到式(14)中既包含系数矩阵的权，也包含观测值的权，为此提出如下两个降权方案，以便作比较分析。

方案1，在式(14)中仅对Pe进行降权。

综上所述，计算TRLS估计的迭代算法设计如下：

3　算例及其分析

算例1，2D仿射变换的数学模型为

(28)

假定变换参数的参考值为

(29)

而无误差的观测数据如表1所示。

表1　无误差的原始观测数据

根据MonteCarlo策略[15]，将协方差矩阵为

(30)QS=0.005diag([12315427218 36])

QT=0.005diag([13611843654 52])

的随机误差分别加在13个点的无误差的原始观测值上，得到模拟观测值，并且在第4个点的旧坐标分量xs上加入大小为2 m的粗差。临界值ε和m分别选为0.000 01和1.5，最大迭代次数为200，权函数取Huber权函数。为消除随机因素的影响，随机模拟运行了200次。单个粗差情形下，由各方案计算得到的变换参数估值与参考值之间的欧氏距离如图1所示，图2则给出了单位权方差的估值。

与此同时，对上述200次模拟结果进行了统计分析，其结果如表2、表3和表4所示。从表中不难看出，当系数矩阵元素受到粗差污染时，WTLS估计扭曲得十分厉害， 其抗差性很差； 而

基于方案1的TRLS估计则较为有效地削弱了粗差对未知参数估值的影响，获得了可靠的参数估值；但方案2并没有达到预期的抗差效果，其稳定性还不如WTLS估计，这说明在估计系数矩阵时，无须进行相应的降权。如果将TRLS估计抵御粗差的成功率定义为

(31)则两种方案的成功率分别为100%和29.5%，这表明方案1的TRLS估计更具稳定性和可靠性。

表2　单个粗差情形下200次模拟的变换参数估值与参考值之间的平均距离

表3　单个粗差情形下单个变换参数估值与参考值之间的绝对距离

表4　单个粗差情形下200次模拟的单位权方差估值的统计结果

下面考虑多个粗差的情形。为此，继续在第4个点的旧坐标分量ys上加入大小为2 m的粗差，同时在第2个点的新坐标两分量xt和yt上也都加入大小为2 m的粗差。模拟计算200次，其变换参数估值与参考值之间的欧氏距离以及单位权方差估计分别如图3和图4所示，而变换参数估值与参考值之间的平均距离如表5所示。不难看出，方案2的TRLS估计基本上不具有抵御粗差的能力，而方案1的TRLS估计则克服了粗差的影响，并且成功率达到100%。

表5　多个粗差情形下200次模拟的变换参数与参考值之间欧氏距离的平均距离

综上可知，不论是单个粗差还是多个粗差，方案2的TRLS估计均不具有较强的抗差能力，而方案1的TRLS估计则能够有效削弱单个或多个粗差对未知参数估值的不良影响，这说明在处理部分变量误差模型中的粗差时，应采用方案1的TRLS估计方法。为了显示本文算法的优越性，在相同模拟条件下，仍然应用前面的相关数据，将本文算法与基于残差和验后单位权方差确定降权因子(式(19))的抗差估计进行比较，其计算结果如表6和表7所示。相关结果表明本文算法优于现有算法。

表6　单个粗差情形下200次的模拟变换参数估值与参考值之间距离

表7　多个粗差情形下200次模拟的变换参数估值与参考值之间的距离

算例2，考虑如表8所示的实测数据[29]，该数据是关于20个孩子的血液样本中血清中卡那霉素的水平。其中一项数据是通过足跟方式获取，而另外一项则是用导管的方式获取。文献[29]指出表8中的第2个观测值含有粗差，据此本文将去掉该观测值所得到的WTLS估值作为未知参数的参考值。在计算过程中，迭代终止的条件为相邻迭代估值之间的距离小于0.000 01，或者是迭代次数超过200次，权函数取Huber权函数。分别采用WTLS估计、方案1和方案2的TRLS估计进行解算，并与文献[15]中的结果进行对比，如表9所示。表10则给出了由各种方法得到的单位权方差估值。图5表示的是不同参数估值对应的拟合结果。

图1　各方案计算得到的变换参数估值与参考值之间的欧氏距离(单个粗差情形)Fig.1　Euclidean distances between estimated transformation parameters by different schemes and reference values with single outlier

图2　单位权方差的估值(单个粗差情形)Fig.2　Estimate of variance of unit weight with single outlier

图3　各方案计算得到的变换参数估值与参考值之间的欧氏距离(多个粗差情形)Fig.3　Euclidean distance between estimated transformation parameters by different schemes and reference values with multiple outliers

图4　200次模拟的单位权方差估值(多个粗差情形)Fig.4　Estimate of variance of unit weight for 200 simulations with multiple outliers

图5　线性拟合图Fig.5　Linear fitting

表8　血液样本数据

表9　各种方法的计算结果

表10　单位权方差的估值

从表9可以看出，方案1和方案2的TRLS方法均具有一定的抵御粗差的功效，但就未知参数估值与参考值之间的欧氏距离来看，方案1明显优于方案2。与文献[15]中的结果进行比较发现，本文的TRLS估值更稳健。从表10可以看出，本文算法的计算精度高于文献[15]，这也有力地验证了本文方法的有效性和优良性。

4　结　论

(1) 部分变量误差模型的WTLS估计同LS估计一样，对粗差很敏感，即抗差性能很差。

(2) 提出了部分变量误差模型WTLS估计的一种两步迭代解法，该解法的计算公式均是在LS框架下推导得到的，具有易于理解和可操作性强等特点，为构建部分变量误差模型的整体抗差估计奠定了基础。

(3) 利用部分变量误差模型的两步迭代解法，构建了与抗差LS估计相类似的TRLS解式，保持了传统的数据处理方式，并构造了一致最大功效统计量以确定降权因子，其中未知的单位权方差可采用具有高崩溃污染率的LMS方法进行估计。

(4) 结合部分变量误差模型两步迭代解法的特点，设计了两个降权方案：第1个方案是在估计系数矩阵元素时，不对观测值进行降权，仅对系数矩阵降权，而第2个方案则是在估计系数矩阵元素时对观测值和系数矩阵同时降权。

(5) 2D相似变换模拟算例的计算结果表明，单个粗差情形下两种不同降权方案的TRLS估计均具有一定的抗差能力，但方案1的抗差功效更优。多个粗差情形下，方案2的TRLS估计基本上不具有抗差性。因此，处理部分变量误差模型中的粗差时，应采用方案1的TRLS估计。线性拟合实测数据的计算结果表明本文方法优于现有变量误差模型的抗差估计方法。

(6) 需要注意的是，统计量(式(20)、式(21))并非严格意义上的一致最大功效统计量，因此，部分变量误差模型的整体抗差估计及其算法值得进一步研究。

[1]FANGX.WeightedTotalLeastSquaresSolutionsforApplicationinGeodesy[D].Hanover:LeibnizUniversityHanover, 2011.

[2]SNOWKB.TopicsinTotalLeast-squaresAdjustmentwithintheErrors-in-variablesModel:SingularCofactorMatricesandPriorInformation[D].Ohio:OhioStateUniversity, 2012.

[3]SCHAFFRINB,WIESERA.OnWeightedTotalLeast-squaresAdjustmentforLinearRegression[J].JournalofGeodesy, 2008, 82(7): 415-421.

[4]SHENYunzhong,LIBofeng,CHENYi.AnIterativeSolutionofWeightedTotalLeast-squaresAdjustment[J].JournalofGeodesy, 2011, 85(4): 229-238.

[5]MAHBOUBV.OnWeightedTotalLeast-squaresforGeodeticTransformations[J].JournalofGeodesy, 2012, 86(5): 359-367.

[6]AMIRI-SIMKOOEIAR,JAZAERIS.WeightedTotalLeastSquaresFormulatedbyStandardLeastSquaresTheory[J].JournalofGeodeticScience, 2012, 2(2): 113-124.

[7]FANGXing.WeightedTotalLeastSquares:NecessaryandSufficientConditions,FixedandRandomParameters[J].JournalofGeodesy, 2013, 87(8): 733-749.

[8]JAZAERIS,AMIRI-SIMKOOEIAR,SHARIFIMA.IterativeAlgorithmforWeightedTotalLeastSquaresAdjustment[J].SurveyReview, 2014, 46(334): 19-27.

[9]XUPeiliang,LIUJingnan,SHIChuang.TotalLeastSquaresAdjustmentinPartialErrors-in-variablesModels:AlgorithmandStatisticalAnalysis[J].JournalofGeodesy, 2012, 86(8): 661-675.

[10]SHIYun,XUPeiliang,LIUJingnan,etal.AlternativeFormulaeforParameterEstimationinPartialErrors-in-variablesModels[J].JournalofGeodesy, 2015, 89(1): 13-16.

[11]ZENGWenxian,LIUJingnan,YAOYibin.OnPartialErrors-in-variablesModelswithInequalityConstraintsofParametersandVariables[J].JournalofGeodesy, 2015, 89(2): 111-119.

[12]SCHAFFRINB,UZUNS.Errors-in-variablesforMobileMappingAlgorithmsinthePresenceofOutliers[J].ArchivesofPhotogrammetry,CartographyandRemoteSensing, 2011, 22: 377-387.

[13]SCHAFFRINB,UZUNS.OntheReliabilityofErrors-in-variablesModels[J].ActaetCommentationesUniversitatisTartuensisDeMathematica, 2012, 16(1): 69-81.

[14]陈义, 陆珏. 以三维坐标转换为例解算稳健总体最小二乘方法[J]. 测绘学报, 2012, 41(5): 715-722.CHENGYi,LUJue.Performing3DSimilarityTransformationbyRobustTotalLeastSquares[J].ActaGeodaeticaetCartographicaSinica, 2012, 41(5): 715-722.

[15]MAHBOUBV,AMIRI-SIMKOOEIAR,SHARIFIMA.IterativelyReweightedTotalLeastSquares:ARobustEstimationinErrors-in-variablesModels[J].SurveyReview, 2013, 45(329): 92-99.

[16]LUJ,CHENY,LIBF,FANGX.RobustTotalLeastSquareswithReweightingIterationforThree-dimensionalSimilarityTransformation[J].SurveyReview, 2014, 46(334): 28-36.

[17]TAOYQ,GAOJX,YAOYF.TLSAlgorithmforGPSHeightFittingBasedonRobustEstimation[J].SurveyReview, 2014, 46(336): 184-188.

[18]ROUSSEEUWPJ,LEROYAM.RobustRegressionandOutlierDetection[M].NewYork:JohnWiley&Sons, 1987.

[19]KOCHKR.ParameterEstimationandHypothesisTestinginLinearModels[M]. 2nded.NewYork:Springer, 1999.

[20]HUBERPJ.RobustStatistics[M].NewYork:JohnWiley&Sons, 1981.

[21]周江文. 经典误差理论与抗差估计[J]. 测绘学报, 1989, 18(2): 115-120.

ZHOUJiangwen.ClassicalTheoryofErrorsandRobustEstimation[J].ActaGeodaeticaetCartographicaSinica, 1989, 18(2): 115-120.

[22]YANGY,SONGL,XUT.RobustEstimatorforCorrelatedObservationsBasedonBifactorEquivalentWeights[J].JournalofGeodesy, 2002, 76(6-7): 353-358.

[23]YANGY.RobustEstimationforDependentObservations[J].ManuscriptaGeodaetica, 1994, 19(1): 10-17.

[24]YANGY.RobustEstimationofGeodeticDatumTransformation[J].JournalofGeodesy, 1999, 73(5): 268-274.

[25]XUPeiliang.Sign-constrainedRobustLeastSquares,SubjectiveBreakdownPointandtheEffectofWeightsofObservationsonRobustness[J].JournalofGeodesy, 2005, 79(1-3): 146-159.

[26]GUOJianfeng,OUJikun,WANGHaitao.RobustEstimationforCorrelatedObservations:TwoLocalSensitivity-BasedDownweightingStrategies[J].JournalofGeodesy, 2010, 84(4): 243-250.

[27]黄维彬. 近代平差理论及其应用[M]. 北京: 解放军出版社, 1992.HUANGWeibin.TheoryandApplicationsofModernAdjustment[M].Beijing:PLAPublishingHouse, 1992.

[28]郭建锋, 赵俊. 粗差探测与识别统计检验量的比较分析[J]. 测绘学报, 2012, 41(1): 14-18.GUOJianfeng,ZHAOJun.ComparativeAnalysisofStatisticalTestsUsedforDetectionandIdentificationofOutliers[J].ActaGeodaeticaetCartographicaSinica, 2012, 41(1): 14-18.

[29]KELLYG.TheInfluenceFunctionintheErrorsinVariablesProblem[J].TheAnnalsofStatistics, 1984, 12(1): 87-100.

(责任编辑：丛树平)

收稿日期： 2015-07-16

修回日期： 2015-12-26

第一作者简介： 赵俊(1987—)，男，博士生，研究方向为测量数据处理。

Firstauthor：ZHAOJun(1987—),male,PhDcandidate,majorsinsurveyingdataprocessing.

E-mail：zhaojun4368@163.com

LUOZhicai

TotalRobustifiedLeastSquaresEstimationinPartialErrors-in-variablesModel

ZHAOJun1,GUIQingming2

1.InstituteofGeospatialInformation,InformationEngineeringUniversity,Zhengzhou450001,China; 2.InstituteofScience,InformationEngineeringUniversity,Zhengzhou450001,China

Theweightedtotalleast-squares(WTLS)estimateforthepartialerrors-in-variables(EIV)modelisverysusceptibletooutliers.BecausetheobservationsandcoefficientmatrixinthepartialEIVmodelmaybecontaminatedwithoutlierssimultaneously,atotalrobustifiedleastsquares(TRLS)estimationforthepartialEIVmodelisproposedbycombiningatwo-stepiteratedalgorithmoftheWTLSestimatewiththeequivalentweightmethodofrobustM-estimation.Andtheuniformlymostpowerfulteststatisticsareconstructedtodeterminethedown-weightingfactors.Forthecharacteristicsofthetwo-stepiteratedmethod,twodifferentdown-weightingschemesarepresented.Inthefirstschemedown-weightingisonlyimplementedforthecoefficientmatrixandnotfortheobservationswhensomeelementsofthecoefficientmatrixareestimated,andthesecondschemeiscontrary.Asimulatedtwo-dimensionalaffinetransformationandalinearfittingwithrealdataareanalyzed.TheresultsshowthattheTRLSwiththefirstschemeissuperiortoonewiththesecondscheme,anditoutperformstheexistingrobustmethodswithresidualandposteriorestimateofvarianceofunitweightandexistingrobustmethodsforthegeneralEIVmodel.

partialEIVmodel；weightedtotalleast-squares；outlier；TRLS；uniformlymostpowerfulteststatistics

2015-04-27

2016-03-01

吴怿昊(1987—)，男，博士生，研究方向为物理大地测量学。Firstauthor：WUYihao(1987—),male,PhDcandidate,majorsinphysicalgeodesy.

E-mail：whuwyh@126.com

罗志才

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

ZHAOJun,GUIQingming.TotalRobustifiedLeastSquaresEstimationinPartialErrors-in-variablesModel[J].ActaGeodaeticaetCartographicaSinica,2016,45(5):552-559.DOI:10.11947/j.AGCS.2016.20150374.

P207

A

1001-1595(2016)05-0552-08

国家自然科学基金(41174005；41474009)

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