Baarda粗差探測和M抗差估計的優(yōu)劣性比較

閆占瑞

（鐵正檢測科技有限公司 山東濟南 250101）

基于均值平移模型的粗差探測是處理觀測值粗差的主要方法之一［1］，它主要采用逐步搜索的方法來實現(xiàn)粗差的定位與定值，通常有向前、向后搜索法。著名的Baarda數(shù)據(jù)探測法就是均值漂移模型的一個典型的實際應用。1968年，巴爾達（Baarda）提出了一套粗差檢測的理論和方法。時至今日，粗差探測仍是測量數(shù)據(jù)研究熱點之一［2］。

基于方差膨脹模型的抗差估計也叫穩(wěn)健估計，是指在測量數(shù)據(jù)存在粗差時，采用一定的估計準則，使參數(shù)估值盡可能接近最優(yōu)的參數(shù)估值［3］。目前，在測量平差研究中，M 估計是使用最廣泛、最簡明的穩(wěn)健估計法［4］，M 抗差估計的抗差性和效率與等價權函數(shù)及其臨界值的合理性和參數(shù)初值的可靠性有關［5］。常用的等價權函數(shù)有L1 法、L1-L2 法、Tukey 法、Danish 法、Fair 法、Cauchy 法、Hampel 法、IGG 方案和IGG3方案等［6］，其中，IGG3 方案是楊元喜院士提出的相關抗差估計解式，近年來得到了廣泛應用［7-11］。

1 Baarda數(shù)據(jù)探測法粗差探測

Baarda 數(shù)據(jù)探測法由荷蘭的Baarda 提出［12］，Baarda 粗差探測的是以標準化殘差為統(tǒng)計量，根據(jù)標準化殘差大小，判斷觀測值是否包含粗差，將粗差剔除后重新平差。在最小二乘平差中，當一個觀測值含有粗差時，會導致多個觀測值的標準化殘差超限。實際計算時，是將絕對值最大且超限的觀測值剔除，然后再進行平差和數(shù)據(jù)探測，直到所有的殘差均不超限，最后用沒有粗差的觀測值進行平差。

設誤差方程式為：

式中：V是n維觀測值殘差向量；A是n×t階系數(shù)矩陣；是t維參數(shù)估值向量；L是n維觀測值向量。設觀測值的權矩陣為P，P是對角陣。

由式（1），參數(shù)的最小二乘解為：

將代入式（1），可求得最小二乘殘差V。又知V的權逆陣為

QV的對角線元素qv1，qv2，…，qvn是觀測值殘差v1，v2，…，vn的權倒數(shù)。當觀測值獨立時：

構造標準化殘差統(tǒng)計量：

式中，σ0為先驗單位權中誤差，vi為第i觀測值的殘差，qvi為殘差協(xié)因數(shù)陣QV主對角線上第i個元素的值。觀測值沒有粗差時，wi是服從標準正態(tài)分布的隨機變量，當|wi|大于給定的限值（如2或3）時，即認定觀測值Li存在粗差。

2 抗差M估計

設有誤差方程式（1）為當觀測值中含有粗差時，參數(shù)的最小二乘（LS）解必然受到歪曲。為求抗差解，一般將中的V2i用其他穩(wěn)健函數(shù)代替，如|V2i|或ρ(Vi)，選擇不同的函數(shù)，其抗差能力也不一樣。

M估計極值函數(shù)：

式中，pi是Li的權。由誤差方差（1）和極值函數(shù)（6）可得參數(shù)的抗差解為：

式中，為等價權矩陣。以表示的對角線元素，為：

式中，w是權因子函數(shù)，常用的權因子函數(shù)有Huber 函數(shù)、IGG1、IGG3 函數(shù)等。w()稱為等價權因子；是標準化殘差，即：

式中，vi是觀測值殘差，mvi是vi的中誤差。mvi由下式計算：

式中：σ0是單位權中誤差，可采用理論值或經驗值；pi是觀測值的權；Ai是矩陣A的第i行，即第i個誤差方程的系數(shù)向量。

抗差估計實質是一個選權迭代的過程，具體步驟如下所示：

Step1：取k=0，參數(shù)的初值為X(k)=[x(1k)，x(2k)，…，x(tk)]T

Step2：計算殘差V(k)=A X∧(k)-L

的權逆陣為：

單位權中誤差公式為：

式中，n0是等價權因子等于零的觀測值個數(shù)。

3 實驗環(huán)節(jié)

3.1 實驗說明

本次實驗主要目的是對處理觀測值中粗差常用的兩種方法優(yōu)劣性進行對比。對于粗差的探測，本次實驗采用的是Baarda數(shù)據(jù)探測法進行粗差探測。對于本次實驗，其中，實驗參數(shù)σ0=1.4''，標準化殘差|wi|的限差值設置為3，當標準化殘差大于3 時，判定觀測值存在粗差，將含有粗差的觀測值剔除。

本次抗差估計采用的是M 抗差估計，關于抗差估計中的權函數(shù)選擇的是由楊元喜院士提出的IGG3 方案，其中，觀測數(shù)據(jù)按照質量高低劃分為3 類：有效信息、可利用信息和有害信息。對于有效信息，采用最小二乘（LS）估計法，采用降權估計處理可利用的信息，采用零權估計處理有害信息。

其中，IGG3權函數(shù)如下所示：

式中，k0、k1為調和系數(shù)，將分成3段，依次對應有效信息、可利用信息和有害信息。在實際計算中，k0的取值范圍為1.0～1.5，k1的取值范圍為2.5～3.0，是標準化殘差，，其中，σ0是單位權方差因子。對于本次實驗，調和系數(shù)k0=1.5，k1=2.5，σ0=1.4''，且迭代收斂條件ε=0.002。

3.2 算例分析

為了比較Baarda 粗差探測和M 抗差估計效果，本文模擬測角網進行了實驗和分析，測角網示意圖如圖1所示。圖中，A、B、C、D這4點為已知點，P1、P2為未知點，為獲得未知點的坐標，獨立等精度觀測了18 個角度，測角中誤差為m=±1.4"，以坐標平差法平差該網，并對誤差方程線性化處理［5］。

圖1 測角網示意圖

3.3 實驗方案實施及結果分析

本實驗主要目的是比較粗差探測和M抗差估計粗差探測的能力。因此，首先，最小二乘平差進行平差，得出最佳估值；然后，按照預設方案，向觀測數(shù)據(jù)中加入粗差，再使用最小二乘平差，并同時采用Baarda數(shù)據(jù)探測法進行粗差探測和IGG3方案進行M抗差估計；最后，分別將各方案結果與正常模式下最佳估值進行比較。

本實驗的粗差加入方案為：（1）在觀測值L2中加入5''的粗差；（2）在觀測值L2、L9中分別加入5''和10''的粗差；（3）在L2、L9、L15中分別加入5''、10''和15''的粗差；（4）在觀測值L13、L14、L15中分別加入5''、10''和15''的粗差。

依據(jù)粗差加入方案進行的實驗結果如表1所示。在本次實驗中，進行方案（4）時發(fā)現(xiàn)在進行粗差探測時探測失敗，粗差探測只探測出了L14、L15觀測值含有粗差，沒有探測出L13，出現(xiàn)了漏判的現(xiàn)象。為尋求原因，又增加兩組實驗進行粗差的探測：（5）在觀測值L8、L14、L15加入粗差；（6）在觀測值L5、L14、L15加入粗差。實驗（5）和實驗（6）的結果如表2所示。

表1 不同方案粗差處理結果與正常模式下的參數(shù)估值的比較

表2 實驗（5）、實驗（6）結果

分析上述實驗結果，可得出以下結論。

（1）隨著加入粗差數(shù)目的增加，由最小二乘平差（LS）得到的參數(shù)的估值偏差也越來越大，結果表明，最小二乘平差（LS）已經不再適用于對未知參數(shù)進行估計。

（2）粗差探測和M 抗差估計都可以有效地減小粗差對觀測數(shù)據(jù)的干擾，獲得較為可靠參數(shù)估計值，但是，粗差數(shù)目的增加使得兩種方法的效果都在下降。

（3）另外，從前3 組實驗結果來看，M 抗差估計比Baarda 粗差探測效果好，得到的結果和正常模式下的參數(shù)估值更加接近，且在第4 組實驗中進行粗差探測時出現(xiàn)了漏判的現(xiàn)象，導致探測失敗。為尋求原因，又增加兩組實驗進行粗差的探測。由表2可知，新增的兩組實驗結果很好分析的原因，粗差探測的效果和含有粗差的觀測值的分布有關，第4 組實驗中加入粗差的3個觀測值是同一個三角網的3個角度觀測值，而新增的實驗粗差加入的位置更加的均勻，所以，粗差探測的效果更好。

4 結語

本文選用Baarda數(shù)據(jù)探測法和M抗差估計并同時選擇IGG3作為權函數(shù)進行粗差探測，并結合實例進行了實驗對比。實驗結果表明：粗差探測的效果和含有粗差的觀測值的分布有關，M 抗差估計的效果與等價權函數(shù)及其臨界值的合理性有關，總體來講，M抗差估計粗差探測效果優(yōu)于Baarda粗差探測。