穩健估計

測量中粗差不可避免，但是常用的最小二乘估計卻不具備抗干擾的能力，對含粗差的觀測量相當敏感，因此如果拿最小二乘原理去處理韓有粗差的數據，會造成嚴重的后果。

在現代廣義平差測量中，常常使用穩健估計來探測和處理粗差，其中最經常用的就是就是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;
        }
    }
   
}