稳健估计

Posted 2020-08-15 王小会

测量中粗差不可避免，但是常用的最小二乘估计却不具备抗干扰的能力，对含粗差的观测量相当敏感，因此如果拿最小二乘原理去处理韩有粗差的数据，会造成严重的后果。

在现代广义平差测量中，常常使用稳健估计来探测和处理粗差，其中最经常用的就是就是M估计，这里p函数我选择Huber函数，来实现。

具体原理：

 

程序中我使用了我前两篇博客中的矩阵库和间接平差类。

 

 

using CMath;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Adjust
{
    /// <summary>
    /// 稳健估计，M估计
    /// </summary>
    public static class Huber
    {
        /// <summary>
        /// 得到粗差剔除后的改正值
        /// </summary>
        /// <param name="fun">观测方程</param>
        /// <param name="e">迭代精度</param>
        /// /// <param name="size">中误差调整系数</param>
        /// <returns>得到粗差剔除后的改正值</returns>
        public static Matrix GetHuber(LinearEquation fun,double size=1, double e=0.1)
        {
            //构造一个参数矩阵，默认为0
            double[,] xite = new double[fun.B.Cols, 1];
            for (int i = 0; i < xite.GetLength(0); i++)
            {
                xite[i,0] = 0.0;
            }
            Matrix Xite = new Matrix(xite);
            Matrix X = fun.GetNW();
            //构造一个精度矩阵
            double[,] eite=new double[fun.B.Cols,1];
             for (int i = 0; i < eite.GetLength(0); i++)
            {
                xite[i,0] = e;
            }
            Matrix Eite = new Matrix(xite);
            //开始迭代
            while (Matrix.MAbs(X - Xite) > Eite)
            {
                X = Xite;
                bool iserror = false;
                //包含含有粗差的索引
                List<int> index = new List<int>();
                //检测是否有粗差
                //可能含有多个粗差
                try
                {
                    for (int i = 0; i < fun.GetV().Rows; i++)
                    {
                        if (Math.Abs(fun.GetV()[i, 0]) > size * fun.GetSigam())
                        {
                            iserror = false;
                            index.Add(i);
                        }
                    }
                }
                catch
                {
                    throw new Exception("误差方程构造错误");
                }
                //未含有粗差
                if (iserror) break;
                //含有粗差
                else
                {
                    //调整观测权
                    //需要调整多个
                    for (int i = 0; i < index.Count; i++)
                    {    
                        //这里使用Huber来改正权
                        fun.P[index[i], index[i]] = size*fun.GetSigam() / Math.Abs(fun.GetV()[index[i],0]);
                    }
                }
                //重新解算
                Xite = fun.GetNW();
            }
            return Xite;
        }
    }
   
}