穩(wěn)健估計(jì)的一種改進(jìn)迭代算法

方 興，黃李雄，曾文憲，吳 云

武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079

當(dāng)觀測(cè)值不含粗差、觀測(cè)誤差服從零均值的分布時(shí)，經(jīng)典最小二乘估計(jì)具有最優(yōu)無(wú)偏性。然而，在實(shí)際的數(shù)據(jù)采集中經(jīng)常會(huì)出現(xiàn)粗差，由于最小二乘不具有抗差性，對(duì)粗差觀測(cè)值十分敏感，少量的粗差就可能對(duì)最小二乘估計(jì)結(jié)果產(chǎn)生破壞性影響[1-6]。因此，在粗差干擾出現(xiàn)的情況下，選擇適當(dāng)?shù)墓烙?jì)方法，盡量降低或者去除粗差對(duì)參數(shù)估值的影響，得出正常模式下的最優(yōu)或接近最優(yōu)的參數(shù)估值，即為測(cè)量數(shù)據(jù)處理領(lǐng)域穩(wěn)健估計(jì)的基本思想。對(duì)于現(xiàn)代空間測(cè)量技術(shù)如GNSS、攝影測(cè)量與遙感等，利用穩(wěn)健估計(jì)處理自動(dòng)化采集、海量的測(cè)量數(shù)據(jù)中混入的粗差具有重要的理論和實(shí)際意義。

1953年，文獻(xiàn)[7]首次提出了穩(wěn)健估計(jì)(又稱為抗差估計(jì))這一專業(yè)術(shù)語(yǔ)。文獻(xiàn)[8]提出了污染分布模式。文獻(xiàn)[9]闡述了定位參數(shù)的抗差估計(jì)。文獻(xiàn)[10]提出了影響函數(shù)和崩潰點(diǎn)的概念。這些都為穩(wěn)健估計(jì)理論奠定了基礎(chǔ)。1980年，文獻(xiàn)[11]將穩(wěn)健估計(jì)理論引入測(cè)量界，提出了著名的“丹麥法”。此外，國(guó)內(nèi)學(xué)者在穩(wěn)健估計(jì)方面也取得重大突破，如文獻(xiàn)[12]提出等價(jià)權(quán)的概念；文獻(xiàn)[13]提出自適應(yīng)的穩(wěn)健估計(jì)算法等。目前，穩(wěn)健估計(jì)可以分為3大類：M估計(jì)、L估計(jì)、R估計(jì)。其中，M估計(jì)是經(jīng)典極大似然估計(jì)的推廣，是測(cè)量數(shù)據(jù)處理中穩(wěn)健估計(jì)最重要的估計(jì)方法。

穩(wěn)健估計(jì)的選權(quán)迭代算法屬于M估計(jì)，在測(cè)量中應(yīng)用廣泛?；驹頌椋菏紫炔捎米钚《朔ǖ玫接^測(cè)值的估計(jì)殘差；根據(jù)殘差確定各觀測(cè)值新的權(quán)因子，進(jìn)行下一輪平差計(jì)算；反復(fù)迭代，直至前后兩次估計(jì)值的變化值小于限差為止。算法每次迭代求解法方程均需進(jìn)行矩陣的求逆運(yùn)算，矩陣求逆的復(fù)雜度是其維數(shù)的三次方。因此，平差模型的待估參數(shù)多，即法方程正交矩陣維數(shù)高，算法的計(jì)算量和計(jì)算時(shí)間急劇上升，算法效率難以滿足現(xiàn)代測(cè)量中大規(guī)模測(cè)量方案、海量數(shù)據(jù)及平差模型實(shí)時(shí)解算等要求。

本文基于矩陣逆的運(yùn)算法則，提出了一種有效改進(jìn)現(xiàn)有選權(quán)迭代算法的方法。迭代計(jì)算中，僅需計(jì)算由前一次殘差確定的新權(quán)陣下的模型未知參數(shù)估計(jì)值的變化量，避免了法方程正交矩陣的求逆，顯著提高了算法的效率，并能保持原穩(wěn)健權(quán)函數(shù)的抗差性。本文以水準(zhǔn)網(wǎng)為例，通過(guò)增加估計(jì)參數(shù)的數(shù)量和粗差數(shù)量，驗(yàn)證了當(dāng)模型的估計(jì)量增加時(shí)，改進(jìn)的選權(quán)迭代算法的計(jì)算效率要遠(yuǎn)高于經(jīng)典選權(quán)迭代算法。

1 經(jīng)典選權(quán)迭代算法

設(shè)平差的線性函數(shù)模型和隨機(jī)模型為

(1)

式中，l為觀測(cè)值向量；B為系數(shù)矩陣并且列滿秩；x為未知參數(shù)向量；Δ為誤差向量；D為對(duì)角方差-協(xié)方差矩陣；σ0為單位權(quán)中誤差；Q為對(duì)角協(xié)因數(shù)矩陣；P為對(duì)角權(quán)陣。

獨(dú)立不等權(quán)觀測(cè)情況下的M估計(jì)準(zhǔn)則為

(2)

(3)

(4)

(5)

將式(3)代入式(5)，得到M估計(jì)下的法方程

(6)

解式(6)得

(7)

根據(jù)殘差v確定各觀測(cè)值新的權(quán)因子，再進(jìn)行下一輪平差，反復(fù)迭代直至前后兩次估計(jì)值的變化小于限差，計(jì)算結(jié)束。算法的核心是確定權(quán)函數(shù)ρ(或φ)，即由每次迭代的殘差估計(jì)值確定新的權(quán)陣，以保證估值的抗差性和效率[13]。以下介紹幾種常見(jiàn)的權(quán)函數(shù)。

(1) Huber法[14-17]

(2) Danish法[18-20]

(3) IGGⅢ方案[21-23]

c1=1.5,c2=3.0

以上權(quán)函數(shù)中，w(u)為權(quán)函數(shù)，c1和c2為相應(yīng)的調(diào)和系數(shù)。

2 改進(jìn)選權(quán)迭代算法

(8)

式中

(9)

將式(9)代入式(8)得

(10)

根據(jù)矩陣求逆引理[24-25]

(A+UCV)-1=A-1-A-1U(C-1+VA-1U)-1VA-1

(11)

將式(10)展開(kāi)

(12)

(13)

式(7)中，令wold、wnew分別表示權(quán)陣更新前后的權(quán)陣，則

wnew=BT(P+dP)l=BTPl+BTdPl=wold+dw

(14)

式中，dw為

(15)

令

(16)

由式(12)—式(16)得

[wold+bi·dpili]=

diλ-diγ-1β-diγ-1αλ=

[λ(γ-α)-β]γ-1di

(17)

則式(13)可表示為

(18)

改進(jìn)選權(quán)迭代算法步驟總結(jié)如下：

(1) 估計(jì)最小二乘初值：

(2) 選權(quán)迭代計(jì)算：

我的姐姐們整天跟著父親下地，不是薅草，就是栽苗挖苕，總有干不完的活兒。母親雖然不允許我走出院子，但她有時(shí)候竟忘了我的存在。她喂了更多的雞和豬，每天要打幾背簍的草，還要切細(xì)剁碎，哄著那些牲口吃光。我沉悶極了，孤零零地坐在院子邊，等著有人經(jīng)過(guò)時(shí)和我說(shuō)上幾句話。

3 實(shí)例分析

為了從數(shù)值上說(shuō)明改進(jìn)選權(quán)迭代算法在計(jì)算效率上的優(yōu)勢(shì)，筆者設(shè)計(jì)了一個(gè)條狀、結(jié)構(gòu)規(guī)整的水準(zhǔn)網(wǎng)實(shí)例，如圖1所示。點(diǎn)A、B、C、D為已知控制點(diǎn)，其余為未知點(diǎn)，未知點(diǎn)個(gè)數(shù)為t(t為偶數(shù))，觀測(cè)值個(gè)數(shù)為n,且n=1.5t+2。

水準(zhǔn)網(wǎng)設(shè)計(jì)試驗(yàn)數(shù)據(jù)數(shù)據(jù)說(shuō)明如下：

控制點(diǎn)高程：HA=1 m,HB=0 m,HC=1.5 m,HD=0.5 m

高差觀測(cè)值：

h1=h2=h7=h8=…=(0.5+Δ)m

h4=h5=h10=h11=…=(-0.5+Δ)m

h3=h6=h9=h12=…=(1+Δ)m

式中，Δ為符合正態(tài)分布(μ=0，σ=0.005)的隨機(jī)

誤差，設(shè)計(jì)粗差的量級(jí)約為1 m。

圖1 模擬的水準(zhǔn)網(wǎng)Fig.1 Simulated leveling network

當(dāng)水準(zhǔn)網(wǎng)包含不同數(shù)量級(jí)的未知點(diǎn)t=100、500、1000、5000、10 000時(shí)，且分別包含1到5個(gè)粗差的情況下，筆者采用Matlab編程對(duì)改進(jìn)選權(quán)迭代算法與經(jīng)典算法的計(jì)算效率進(jìn)行了比較。程序采用Matlab自帶的tic/toc函數(shù)，記錄了算法改進(jìn)前后程序運(yùn)行所需時(shí)間，試驗(yàn)結(jié)果見(jiàn)表1。

表1 運(yùn)行時(shí)間對(duì)比

試驗(yàn)結(jié)果表明：

(1) 改進(jìn)算法計(jì)算效率要高于經(jīng)典選權(quán)迭代算法(圖2、圖3)，尤其是在平差模型中待求參數(shù)個(gè)數(shù)t較大或者多個(gè)粗差的情況下，計(jì)算效率顯著提升。當(dāng)未知數(shù)個(gè)數(shù)t=10 000且含5個(gè)粗差時(shí)，改進(jìn)算法所需計(jì)算時(shí)間僅為經(jīng)典算法的1/30。

(2) 改進(jìn)算法的計(jì)算效率取決于采用最小二乘估計(jì)初值的計(jì)算時(shí)間，迭代計(jì)算無(wú)需求逆，當(dāng)粗差個(gè)數(shù)增加時(shí)，計(jì)算時(shí)間僅有微小增長(zhǎng)。例如，當(dāng)未知數(shù)個(gè)數(shù)t=10 000，粗差從1個(gè)增長(zhǎng)到5個(gè)時(shí)，經(jīng)典算法計(jì)算時(shí)間急劇增長(zhǎng)，由829 s增加4倍，達(dá)到3285 s；而改進(jìn)算法從100 s僅增加到109 s，改進(jìn)算法對(duì)于處理多個(gè)粗差情況非常高效。

4 結(jié) 論

經(jīng)典選權(quán)迭代算法屬于M估計(jì)，是當(dāng)前運(yùn)用最為廣泛的穩(wěn)健估計(jì)方法之一。本文采用矩陣逆運(yùn)算法則，提出了一種改進(jìn)的選權(quán)迭代算法。改進(jìn)算法在迭代過(guò)程中，僅需計(jì)算新權(quán)陣下的最小二乘估計(jì)改正值，避免了經(jīng)典算法中每步對(duì)平差模型法方程的求逆運(yùn)算，顯著提高了算法的計(jì)算效率。改進(jìn)選權(quán)迭代算法適用于現(xiàn)代空間測(cè)量技術(shù)下對(duì)自動(dòng)觀測(cè)的海量數(shù)據(jù)的實(shí)時(shí)或者準(zhǔn)實(shí)時(shí)的處理需求，如GNSS、遙感等觀測(cè)數(shù)據(jù)，當(dāng)這些數(shù)據(jù)含有一定數(shù)量的粗差時(shí)，改進(jìn)算法可以大大降低經(jīng)典選權(quán)迭代算法的計(jì)算時(shí)間。本文算法基于獨(dú)立不等精度觀測(cè)值的情況，筆者下一步的研究目標(biāo)是針對(duì)海量數(shù)據(jù)，討論如何提高相關(guān)觀測(cè)值[5，17，26]的選權(quán)迭代算法的計(jì)算效率。

圖2 改進(jìn)前后運(yùn)行時(shí)間比較(未知點(diǎn)數(shù)為10 000)Fig.2 Computational cost comparison when the unknown parameters is 10 000

圖3 改進(jìn)前后運(yùn)行時(shí)間比較(誤差數(shù)為5)Fig.3 Computational cost comparison when there are 5 errors