一种基于丹麦法的改进型双步M 估计

林国标 刘立龙 蔡成辉 黎峻宇 黄良珂

1 广西矿冶与环境科学实验中心，桂林市建干路12号，541004

2 桂林理工大学测绘地理信息学院，桂林市建干路12号，541004

3 广西空间信息与测绘重点实验室，桂林市建干路12号，541004

抗差估计的抗差性及有效性主要取决于初值的准确性与权函数的合理性［1］。目前，一般选取最小二乘估计的结果作为抗差估计的初值［2］，但最小二乘估计对粗差具有均衡性和不敏感性，致使粗差观测值的残差并非最大，导致选权迭代时出现错误的判断。同时，合理的权函数应该分3段：正常段、可疑段、淘汰段，其临界值和均方差因子应根据实际情况选取，具有可变性［1，3］。在传统的抗差估计方法中，Andrews法、Tukey 法不含正常段，Huber 法、丹麦法没有淘汰段，IGG法、Hampel法虽然具有正常段、可疑段及淘汰段，但其临界值往往取某一固定的先验值。这些方法仅采用一步抗差估计［4］，具有一定的局限性。

针对以上问题，本文提出一种基于丹麦法的改进型双步M 估计方法。采用抗差性较强的L1－范估计的平差值作为抗差估计的初值，在丹麦法等价权函数的基础上引入淘汰域，将改进的丹麦法权函数作为抗差权，选取可变动的临界值，通过双步M 估计，提高抗差估计的精度。

1 改进型双步M 估计

1.1 初值选取

由于LS估计不具抗差性，LS估计的残差不能有效反映粗差的位置。而L1－范估计的平差值具有较好的抗差性，能较为准确地体现真实残差［5］，本文将其作为双步M 估计的初值。

首先，列误差方程：

解得残差：

然后，以L1－范权函数重新定权，再进行迭代计算，收敛时得残差初值：

式中，Pi为观测权，为等价权，V0为初始残差，l为常数项，c为一个很小的正常数。

1.2 等价权函数选择

丹麦法的权函数如下［3］：

丹麦法权函数缺乏淘汰段，导致抗差效率降低。加入淘汰段，对丹麦法进行改进：

1.3 第一步——综合抗差阶段

该阶段有两个关键问题需要解决：一是确定临界值k0、k1，二是选取等价权。k0、k1一般分别取某一固定值［3］，k0＝1～2.5，k1＝3～6。临界值在k0、k1基础上再乘以一个可变因子（r为多余度），可解决临界值取值不当的问题，实现系数阵空间与观测空间的同时抗差［6］。r的平均值为（n－m）／n（n、m分别为观测值个数与参数个数），故可变因子为。此时临界值会随具体问题而变化，增强了抗差估计的灵活性。本文取k0＝2d，k1＝4d。

选取改进后的丹麦法等价权作为抗差权，单位权中误差σ0采用第一次或第二次抗差估计所求得的观测残差来计算，能够获得比较准确的均方差因子［7］：

1.4 第二步——优化阶段

将第一步所得的估计值作为该阶段的初值，利用截尾LS法重新平差计算。取权函数［8］：

临界值取k＝3。根据此时的v与σ2，将大于3倍均方差的观测值删去，然后利用剩余的观测值进行加权LS 平差，得到最终的参数解。精度评定采用：

式中，n、m、t分别为总观测数、必要观测数和删去的粗差观测数；为最终权。

2 算例分析

本文结合一个测边网算例和一个水准网算例［9］，对比本文提出的改进型双步M 估计与丹麦法、Huber法、L1法、IGGⅠ法及两步抗差估计5种方法的平差结果。将无粗差时的LS法平差值当作真值，传统抗差估计方法和两步抗差估计的临界值取先验值，均方差因子通过初次平差的残差由中位数法计算，即0.674 5。

算例1 如图1，测边网中A、B、C、D为已知点，P1、P2、P3、P4为未知点，观测13条边长，测距精度σS＝3mm＋1×10－6S。起算数据及观测边长见表1。在第3、12条边上分别附上－30dm、50 dm 的粗差，各种抗差估计方法的计算结果见表2。

图1 测边网图Fig.1 Trilateration network

算例中，测边网的污染率为15.4%，属于较严重污染。由表2可知：

1）在参数估计方面，几种传统抗差估计方法的估计结果与真值之差较大：丹麦法、Huber法、L1法和IGGⅠ法最大的差值分别为11.2cm、30.4cm、20.3cm、43.1cm；两步抗差估计的结果与真值之差较小，最大为－3.2cm；本文方法与真值之差较小，最大差值为－2.5cm，整体比两步抗差估计小一些。

2）在最终权方面，P3、P12为粗差权，其余为正常观测值的最终权，传统抗差估计方法的P3、P12未能同时降到最小或0，如丹麦法、Huber法、L1法和IGGⅠ法的P3分别为0.03、0.11、0.03、0.13，L1法正常观测值的最终权受粗差影响波动很大；两步抗差估计的P3、P12均降为0，而P6、P8、P13分别降为0.42、0.40、0.73；本文方法的P3、P12均降为0，其余权值均为观测权。

若在上述粗差基础上再加上一个粗差，即在第3、10、12条边上分别附上－30dm、－20dm、50dm 的3个粗差，再一次进行平差计算，则各种抗差估计方法的计算结果见表3。此时测边网的污染率为23.1%。由表3可知：

表1 测边网观测数据及起算数据Tab.1 Observed value and known data

表2 各种抗差估计方法的计算结果Tab.2 Results by various robust estimation methods

表3 各种抗差估计方法的计算结果Tab.3 Results by various robust estimation methods

1）在参数估计方面，传统抗差估计方法和两步抗差估计的结果与真值之差比较大，如丹麦法、Huber法、L1法、IGGⅠ法和两步抗差估计的最大差值分别为109.0、175.6、115.1、121.8、121.8 cm；本文方法的结果与真值之差较小，最大仅为－2.7cm。

2）在最终权方面，P3、P10、P12为粗差权，其余为正常观测值的最终权。传统抗差估计方法和两步抗差估计的P3、P10仍为观测权，分别为1.86、14.2，L1法的粗差权虽然最后有所降低，但是正常观测权之比波动很大；本文方法的P3、P10、P12均降为0，其余权仍为观测权。

综上所述，当测边网含有2个粗差时，传统抗差估计方法的抗差效果较差，两步抗差估计和本文方法的抗差效果较好，而本文方法的整体抗差效果优于两步抗差估计；当测边网含有3个粗差时，传统抗差估计方法和两步抗差估计基本丧失了抗差估计的能力，而本文方法的抗差效果依然比较显著。这是由于传统抗差估计方法和两步抗差估计往往采用对粗差具有均衡性的LS法的平差结果作为抗差初值，不能准确反映粗差的位置。而前者临界值是固定的，不能有效适应实际问题，因此它们抵御多维粗差的效果并不显著。

算例2 如图2，水准网中A、B为已知点，HA＝8.016m，HB＝9.016m，点P1、P2、P3为待定点。进行7条水准路线观测，见表4。在第1、4路线上分别附加－20mm、40mm 的粗差，各种方法的计算结果见表5。

图2 水准网图Fig.2 Leveling network

表4 水准网观测高差值及路线长Tab.4 Observed value and length

表5 各种抗差估计方法的计算结果Tab.5 Results by various robust estimation methods

算例中，多余观测数为4，污染率为28.6%，属于比较严重的污染，因此本算例只进行包含2个粗差的比较。由表5可知：

1）在参数估计方面，传统抗差估计方法和两步抗差估计的估计结果与真值之差较大，丹麦法、Huber法、L1法、IGGⅠ法和两步抗差估计的最大差值分别为－8.7 mm、－8.0 mm、－6.7 mm、－8.3mm、－8.8mm；本文方法的结果与真值之差较小，最大仅为－0.3mm。

2）在最终权方面，P1、P4为粗差权，其余为正常观测值的最终权，传统抗差估计和两步抗差估计的P1、P4未能同时降为0，如丹麦法、Huber法、IGGⅠ法及两步抗差估计的P1观测权均为2.73，L1法的P1为0.21，L1法的P3、P7有所降低；本文方法的P1、P4均降为0，其余最终权维持不变。

和算例1同理，因为传统抗差估计和两步抗差估计常常采用不能准确反映粗差位置的LS法平差结果作为抗差初值，前者临界值是固定的，所以往往不能抵抗多维粗差的干扰。而本方法以L1－范估计结果作为初值，第一步的抗差权采用改进的丹麦法等价权，临界值可变；第二步采用截尾LS法，取得了较好的估计结果。

3 结 语

基于丹麦法的改进型双步M 估计是一种以L1－范估计的平差值作为初值、以改进后的丹麦等价权作为综合抗差阶段权因子的抗差估计方法，比较有效地解决了抗差估计初值与临界值选取不当的问题，提高了抗差估计的精度。算例表明，无论是测边网平差，还是水准网平差，该方法都能有效抵御多维粗差的干扰，获得优于传统抗差估计和两步抗差估计的平差结果，具有较强的抗差性。

［1］杨元喜，吴富梅.临界值可变的抗差估计等价权函数［J］.测绘科学技术学报，2006（5）：320－324（Yang Yuanxi，Wu Fumei.Modified Equivalent Weight Function with Variable Criterion for Robust Estimation［J］.Journal of Geomatics Science and Technology，2006，（5）：320－324）

［2］邱卫宁.具有稳健初值的选权迭代法［J］.武汉大学学报：信息科学版，2003（4）：452－454（Qiu Weining.Method for Selecting Weight Iteration with Robust Initial Value［J］.Geomatics and Information Science of Wuhan University，2003，（4）：452－454）

［3］李浩军，唐诗华，黄杰.经典选权迭代法研究与两步抗差估计的提出［J］.海洋测绘，2007（1）：17－20（Li Haojun，Tang Shihua，Huang Jie.The Discussion about Some Kinds of Selecting Weight Iteration Method［J］.Hydrographic Surveying and Charting，2007（1）：17－20）

［4］靳奉祥.抗差估计理论与方法研究［J］.山东科技大学学报：自然科学版，2003（4）：1－6（Jin Fengxiang.Study on Robust Estimation Theory and Method［J］.Journal of Shandong University of Science and Technology：Natural Science，2003（4）：1－6）

［5］王振杰，欧吉坤，曲国庆.抗差估计的初值选择［J］.地壳形变与地震，2001（3）：32－35（Wang Zhenjie，Ou Jikun，Qu Guoqing.Choice of Initial Value for Robust Estimation［J］.Crustal Deformation and Earthquake，2001，（3）：32－35）

［6］欧吉坤.一种三步抗差方案的设计［J］.测绘学报，1996（3）：173－179（Ou Jikun.Design of a New Scheme of Robust Estimation by Three Steps［J］.Acta Geodaetica et Cartographica Sinica，1996（3）：173－179）

［7］高晓，戴吾蛟.抗差Helmert方差分量估计在GPS／BDS组合定位中的应用［J］.大地测量与地球动力学，2014，34（1）：173－176（Gao Xiao，Dai Wujiao.Application of Robusthelmert Variance Component Estimation to Position in Combination of GPS and BDS［J］.Journal of Geodesy and Geodynamics，2014，34（1）：173－176）

［8］王福丽，成英燕，韦铖，等.利用抗差多项式拟合法探测修复GNSS周跳［J］.大地测量与地球动力学，2013，33（3）：129－132（Wang Fuli，Cheng Yingyan，Wei Cheng，et al.GNSS Cycle Slip Detection and Correction Using Robust Polynomial Fitting［J］.Journal of Geodesy and Geodynamics，2013，33（3）：129－132）

［9］武汉大学测绘学院测量平差学科组.误差理论与测量平差基础［M］.武汉：武汉大学出版社，2009（Research Group of Surveying Adjustment in School of Geodesy and Geomatics of Wuhan University.Error Theory and Fundation of Surveying Adjustment［M］.Wuhan：Wuhan University Press，2009）