Baarda粗差探測和M抗差估計的優劣性比較

閆占瑞

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

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

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

1 Baarda數據探測法粗差探測

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

設誤差方程式為：

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

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

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

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

構造標準化殘差統計量：

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

2 抗差M估計

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

M估計極值函數：

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

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

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

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

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

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

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

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

的權逆陣為：

單位權中誤差公式為：

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

3 實驗環節

3.1 實驗說明

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

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

其中，IGG3權函數如下所示：

式中，k0、k1為調和系數，將分成3段，依次對應有效信息、可利用信息和有害信息。在實際計算中，k0的取值范圍為1.0～1.5，k1的取值范圍為2.5～3.0，是標準化殘差，，其中，σ0是單位權方差因子。對于本次實驗，調和系數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抗差估計粗差探測的能力。因此，首先，最小二乘平差進行平差，得出最佳估值；然后，按照預設方案，向觀測數據中加入粗差，再使用最小二乘平差，并同時采用Baarda數據探測法進行粗差探測和IGG3方案進行M抗差估計；最后，分別將各方案結果與正常模式下最佳估值進行比較。

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

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

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

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

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

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

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

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

4 結語

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