閆占瑞

（鐵正檢測(cè)科技有限公司 山東濟(jì)南 250101）

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

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

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

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

設(shè)誤差方程式為：

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

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

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

QV的對(duì)角線元素qv1，qv2，…，qvn是觀測(cè)值殘差v1，v2，…，vn的權(quán)倒數(shù)。當(dāng)觀測(cè)值獨(dú)立時(shí)：

構(gòu)造標(biāo)準(zhǔn)化殘差統(tǒng)計(jì)量：

式中，σ0為先驗(yàn)單位權(quán)中誤差，vi為第i觀測(cè)值的殘差，qvi為殘差協(xié)因數(shù)陣QV主對(duì)角線上第i個(gè)元素的值。觀測(cè)值沒(méi)有粗差時(shí)，wi是服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量，當(dāng)|wi|大于給定的限值（如2或3）時(shí)，即認(rèn)定觀測(cè)值Li存在粗差。

2 抗差M估計(jì)

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

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

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

式中，為等價(jià)權(quán)矩陣。以表示的對(duì)角線元素，為：

式中，w是權(quán)因子函數(shù)，常用的權(quán)因子函數(shù)有Huber 函數(shù)、IGG1、IGG3 函數(shù)等。w()稱為等價(jià)權(quán)因子；是標(biāo)準(zhǔn)化殘差，即：

式中，vi是觀測(cè)值殘差，mvi是vi的中誤差。mvi由下式計(jì)算：

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

抗差估計(jì)實(shí)質(zhì)是一個(gè)選權(quán)迭代的過(guò)程，具體步驟如下所示：

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

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

的權(quán)逆陣為：

單位權(quán)中誤差公式為：

式中，n0是等價(jià)權(quán)因子等于零的觀測(cè)值個(gè)數(shù)。

3 實(shí)驗(yàn)環(huán)節(jié)

3.1 實(shí)驗(yàn)說(shuō)明

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

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

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

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

3.2 算例分析

為了比較Baarda 粗差探測(cè)和M 抗差估計(jì)效果，本文模擬測(cè)角網(wǎng)進(jìn)行了實(shí)驗(yàn)和分析，測(cè)角網(wǎng)示意圖如圖1所示。圖中，A、B、C、D這4點(diǎn)為已知點(diǎn)，P1、P2為未知點(diǎn)，為獲得未知點(diǎn)的坐標(biāo)，獨(dú)立等精度觀測(cè)了18 個(gè)角度，測(cè)角中誤差為m=±1.4"，以坐標(biāo)平差法平差該網(wǎng)，并對(duì)誤差方程線性化處理［5］。

圖1 測(cè)角網(wǎng)示意圖

3.3 實(shí)驗(yàn)方案實(shí)施及結(jié)果分析

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

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

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

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

表2 實(shí)驗(yàn)（5）、實(shí)驗(yàn)（6）結(jié)果

分析上述實(shí)驗(yàn)結(jié)果，可得出以下結(jié)論。

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

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

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

4 結(jié)語(yǔ)

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