丁 娟，葛永慧

（太原理工大學(xué) 礦業(yè)工程學(xué)院測繪科學(xué)與技術(shù)系，太原030024）

遙感圖像在投入使用前，需要對圖像中產(chǎn)生的畸變進(jìn)行幾何糾正。遙感圖像幾何糾正的方法有多項(xiàng)式糾正法、共線方程糾正法、空間投影法等［1］。多項(xiàng)式糾正法原理直觀、計(jì)算簡單，適用于地面平坦地區(qū)；共線方程糾正法引用了地面高程信息，在地形起伏較大的地區(qū)，糾正精度要高于多項(xiàng)式糾正法［2］；空間投影法精確地模擬了衛(wèi)星動態(tài)成像的幾何關(guān)系，它只需少量地面控制點(diǎn)，就能直接建立影像與投影之間、像點(diǎn)與地面點(diǎn)之間嚴(yán)格的數(shù)學(xué)關(guān)系［1］。幾何糾正，即圖像級別中的幾何精糾正，主要方法是多項(xiàng)式糾正法［3］。

多項(xiàng)式糾正法通常用最小二乘（LS）估計(jì)求取轉(zhuǎn)換參數(shù)。但是當(dāng)觀測值中不可避免地含有粗差時，LS估計(jì)也會受到很大的干擾［4］，結(jié)果便會受到歪曲。G．E．P．BOX提出的穩(wěn)健估計(jì)的概念［5］，可以消除或減弱粗差對參數(shù)估計(jì)的影響，然而不同穩(wěn)健估計(jì)方法減弱和消除粗差對參數(shù)估計(jì)影響的能力不同。Li［6］認(rèn)為Fair法的效果要優(yōu)于殘差絕對和最小法（L1法），而Danish法的計(jì)算量相對其他方法較大。Sharmishtha Mitra等［7］的研究表明，由于重尾噪聲分布的存在，Andrews法和Huber法的穩(wěn)健性要優(yōu)于殘差絕對和最小法。Pennacchi［8］用算例證明，迭代100次以后，Cauchcy法比 Welsch法、German－McClure法和Tukey法的效果都要好。焦偉利等［5］提出利用穩(wěn)健理論自動檢測和消除圖像校正中的控制點(diǎn)粗差，實(shí)驗(yàn)結(jié)果表明，穩(wěn)健估計(jì)可以有效消除觀測值中的粗差，使驗(yàn)后單位權(quán)方差變小，同時提高控制點(diǎn)的擬合精度。那么，多項(xiàng)式糾正法在求取轉(zhuǎn)換參數(shù)時，哪些穩(wěn)健估計(jì)方法相對更為有效呢？

筆者采用多項(xiàng)式糾正法建立模型，通過仿真實(shí)驗(yàn)（1 000次）對不同的觀測值數(shù)量、粗差數(shù)量和粗差數(shù)值的13種穩(wěn)健估計(jì)方法消除或減弱粗差的能力進(jìn)行了比較，確定了遙感圖像幾何精糾正中相對更為有效的穩(wěn)健估計(jì)方法。

1 材料與方法

1．1 13種常用的穩(wěn)健估計(jì)方法及其權(quán)函數(shù)

13種常用的穩(wěn)健估計(jì)方法如下。在下列各式中，u代表標(biāo)準(zhǔn)化的殘差指標(biāo)表示權(quán)函數(shù)；a，b和c表示調(diào)和系數(shù)，它們的值均采用有關(guān)文獻(xiàn)的推薦值。

1）Huber法［9］：

式中：c＝1．5．

2）L1法（殘差絕對和最小法）［8］：

式中：k為很小的數(shù)。

3）L1－L2法［8］：

4）Andrews法［10］：

式中：c＝1．339．

5）Hampel法［11］：

式中：a＝1．6；b＝3；c＝6．5．

6）Welsch法［8］：

式中：c＝2．984 6．

7）Tukey法［8］：

式中：c＝4．685 1．

8）Danish法［12－13］：式中：c＝1．5或c＝1．4，在本文中，c＝1．5．

9）Fair法［8］：

式中：c＝1．399 8．

10）German－McClure法［8］：

11）IGG 方案［13－14］：

式中：b＝1．5；c＝2．5；k為很小的數(shù)。

12）IGGⅢ方案［15－16］：

式中：b＝1．0～1．5；c＝2．5～3．0．在本文中，b＝1．5，c＝3．0．

13）Cauchy法［8］：

式中：c＝2．384 9．

1．2 兩種穩(wěn)健估計(jì)方法比較指標(biāo)

1．2．1 殘余真誤差均方誤差與相對增益

文獻(xiàn)［17］提出了殘余真誤差均方誤差和相對增益兩種指標(biāo)，用仿真實(shí)驗(yàn)的方法比較了任意兩種參數(shù)估計(jì)。筆者采用文獻(xiàn)［17］提出的方法比較13種穩(wěn)健估計(jì)方法對于遙感圖像精糾正消除或減弱粗差的能力。

兩種穩(wěn)健估計(jì)方法比較的絕對指標(biāo)——?dú)堄嗾嬲`差均方誤差（＾σf）：

式中：fk為殘余真誤差為觀測值Lk的真值為通過參數(shù)估計(jì)方法得到觀測值Lk的估值，并隨著參數(shù)估計(jì)方法的不同而不同；n為觀測值的數(shù)量。

兩種穩(wěn)健估計(jì)方法比較的相對指標(biāo)——相對增益（RG）：式中和為A 和B 兩種參數(shù)估計(jì)方法對于同一個參數(shù)估計(jì)問題1 000次仿真實(shí)驗(yàn)殘余真誤差均方誤差的平均值。當(dāng)RG＞0時，B方法相對于A方法的相對增益提高；當(dāng)RG＝0（或接近0）時，B方法與A方法參數(shù)估計(jì)結(jié)果等價(jià)；當(dāng)RG＜0時，B方法相對于A方法的相對增益降低。因此，RG從實(shí)質(zhì)上可以說明兩種參數(shù)估計(jì)方法哪個更有效。在本文中，各種穩(wěn)健估計(jì)方法相對于LS法的相對增益簡稱為相對增益。

1．2．2 觀測值中包含粗差的仿真實(shí)驗(yàn)

用觀測值的真值減去包含粗差的隨機(jī)誤差得到S組包含粗差的模擬觀測值：

用各種穩(wěn)健估計(jì)方法和最小二乘法分別計(jì)算S組中每一組模擬觀測值的MSRTE，并取它們的平均值作為觀測值中包含g個粗差ε時此方法的MSRTE。穩(wěn)健估計(jì)方法的終止條件是相鄰兩次殘差之差的絕對值的最大值小于0．1像素。用最小二乘法和不同穩(wěn)健估計(jì)方法得到的MSRTE分別計(jì)算每一種穩(wěn)健估計(jì)方法相對于最小二乘法的增益。

1．3 遙感圖像幾何精糾正模型

遙感圖像幾何精糾正的主要方法是多項(xiàng)式糾正法，其建立了控制點(diǎn)像坐標(biāo)（x，y）與其控制點(diǎn)地面空間坐標(biāo)（u，v）之間的關(guān)系，多項(xiàng)式糾正模型的數(shù)學(xué)表達(dá)式［18］為：

式中：aij、bij為多項(xiàng)式系數(shù)；N 是多項(xiàng)式次數(shù)，取決于圖像變形的程度和地面控制點(diǎn)的數(shù)量，建立N次多項(xiàng)式模型所需控制點(diǎn)數(shù)至少要（N＋1）×（N＋2）／2個。一般情況采用二次多項(xiàng)式，表達(dá)為：

式中，x、y的表達(dá)式具有相同的形式。在本研究的仿真實(shí)驗(yàn)中，以x為例，計(jì)算坐標(biāo)轉(zhuǎn)換參數(shù)a0，a1，a2，a3，a4，a5。仿真實(shí)驗(yàn)中的重合點(diǎn)如圖1所示，重合點(diǎn)的坐標(biāo)真值是模擬值，9～15個重合點(diǎn)的選取方式如下：

9個重合點(diǎn)：P11，P13，P15，P31，P33，P35，P51，P53，P55；

10個重合點(diǎn)：P11，P13，P15，P31，P32，P34，P35，P51，P53，P55；

11個重合點(diǎn)：P11，P12，P14，P15，P31，P33，P35，P51，P52，P54，P55；

12個重合點(diǎn)：P11，P12，P14，P15，P31，P32，P34，P35，P51，P52，P54，P55；

13個重合點(diǎn)：P11，P13，P15，P22，P24，P31，P33，P35，P42，P44，P51，P53，P55；

14個重合點(diǎn)：P11，P13，P15，P22，P23，P24，P31，P35，P42，P43，P44，P51，P53，P55；

15個重合點(diǎn)：P11，P13，P15，P22，P23，P24，P31，P33，P35，P42，P43，P44，P51，P53，P55．

圖1 重合點(diǎn)分布圖

2 結(jié)果與討論

2．1 遙感圖像幾何精糾正仿真實(shí)驗(yàn)

當(dāng)觀測值數(shù)量n＝9，粗差數(shù)量g＝1時，取粗差ε＝5．0σ0進(jìn)行1 000次的仿真實(shí)驗(yàn)，分別計(jì)算LS法和13種穩(wěn)健估計(jì)方法的殘余真誤差均方誤差，并取它們的平均值作為相應(yīng)的殘余真誤差均方誤差（見表1）。對于粗差ε＝10．0σ0作同樣的計(jì)算。相應(yīng)地，計(jì)算出13種穩(wěn)健估計(jì)方法相對于LS法的相對增益（見表2）。

表1 不同穩(wěn)健估計(jì)方法的殘余真誤差均方誤差（n＝9，g＝1） 像素

在后續(xù)的仿真實(shí)驗(yàn)中，分別計(jì)算在n（9～15）個觀測值中含有g(shù)（1～3）個粗差、ε為5．0σ0和10．0σ0時，LS法和13種穩(wěn)健估計(jì)方法的殘余真誤差均方誤差，并分別計(jì)算出相對應(yīng)的13種穩(wěn)健估計(jì)方法相對于LS法的相對增益。

表1中，MLS表示LS法得到的殘余真誤差均方誤差；M1，M2，M3，M4，M5，M6，M7，M8，M9，M10，M11，M12，M13分別表示 Huber法、L1法、L1－L2法、Andrews法、Hampel法、Welsch法、Tukey 法、Danish法、Fair法、German－McClure法、IGG 方案、IGGIII方案、Cauchy法得到的殘余真誤差均方誤差。

表2 不同穩(wěn)健估計(jì)方法的相對增益（n＝9，g＝1） %

表2中，R1，R2，R3，R4，R5，R6，R7，R8，R9，R10，R11，R12，R13分別表示 Huber法、L1法、L1－L2法、Andrews法、Hampel法、Welsch法、Tukey法、Danish法、Fair法、German－McClure法、IGG方案、IGGIII方案、Cauchy法相當(dāng)于LS法的相對增益。

2．2 仿真實(shí)驗(yàn)結(jié)果與討論

圖2是當(dāng)n＝9～15并且g＝1，n＝13～15并且g＝2，n＝15并且g＝3時，13種常用穩(wěn)健估計(jì)方法相對于LS法的相對增益的平均值。

由圖2可知，L1法（R2）和 German－McClure法（R10）是相對更為有效的穩(wěn)健估計(jì)方法。當(dāng)粗差ε＝5．0σ0時，L1法和 German－McClure法的相對增益平均值分別是13%和11%，其他穩(wěn)健估計(jì)方法的相對增益平均值均小于等于10%。當(dāng)粗差ε＝10．0σ0時，L1法和German－McClure法的相對增益平均值分別是39%和38%，其他穩(wěn)健估計(jì)方法的相對增益平均值均小于等于24%。

圖2 不同穩(wěn)健估計(jì)方法相對于LS法的平均相對增益

3 結(jié)論

筆者通過仿真實(shí)驗(yàn)（1 000次），對不同的重合點(diǎn)（觀測值）數(shù)量（9～15個）、粗差數(shù)量（1～3個）和粗差數(shù)值（5．0σ0、10．0σ0）的13種穩(wěn)健估計(jì)方法消除或減弱粗差的能力進(jìn)行了比較。仿真實(shí)驗(yàn)結(jié)果說明，13種常用穩(wěn)健估計(jì)方法都有消除或減弱粗差影響的能力，但是能力大小不盡相同?？傮w上，對于遙感圖像幾何精糾正，L1法和German－McClure法是相對更為有效的穩(wěn)健估計(jì)方法，它們比其他常用穩(wěn)健估計(jì)方法能更有效地消除或減弱粗差的影響。

［1］ 程志梅，陸玲，邱霞明．利用空間投影法實(shí)現(xiàn)遙感圖像的幾何精糾正［J］．科技廣場，2005（10）：103－105．

［2］ 劉云峰，李若．不同 DEM 數(shù)據(jù)對衛(wèi)星遙感影像糾正精度的影響［J］．測繪通報(bào)，2002（7）：26－28．

［3］ 王學(xué)平．遙感圖像幾何校正原理及效果分析［J］．計(jì)算機(jī)應(yīng)用與軟件，2008，25（9）：102－105．

［4］ 歸慶明，張建軍，郭建鋒．壓縮型抗差估計(jì)［J］．測繪工程，2000，29（3）：224－228．

［5］ 焦偉利，何國金，王威，等．穩(wěn)健估計(jì)方法求解遙感圖像幾何校正模型研究［D］．北京：中國科學(xué)院對地觀測與數(shù)字地球科學(xué)中心，2009．

［6］ Huifen Li，Xaingqian Jiang，Zhu Li．Robust estimation in Gaussian filtering for engineering surface characterization［J］．Precision Engineering，2004（28）：186－193．

［7］ Sharmishtha Mitra，Amit Mitra，Debasis Kundu．Genetic algorithm and M－estimator based robust sequential estimation of parameters of nonlinear sinusoidal signals［J］．Communications in Nonlinear Science and Numencal Simulation，2011，16（7）：2796－2809．

［8］ Pennacchi P．Robust Estimate of Excitations in Mechanical Systems Using M－estimators—Theoretical Background and Numerical Applications［J］．Journal of Sound and Vibration，2008，310（4－5）：923－946．

［9］ Baselga S．Global Optimization Solution of Robust Estimation［J］．Journal of Surveying Engineering，2007，133（3）：123－128．

［10］ El－Hawary F，Mbamalu G A N．Fair and Andrews＇s weighting－based IRWLS algorithms for time－delay estimation in underwater target tracking［J］．Ieee Journal of Oceanic Engineering，1993，18（2）：142－150．

［11］ 李浩軍，唐詩華，黃杰．抗差估計(jì)中幾種選權(quán)迭代法常數(shù)選取的探討［J］．測繪科學(xué)，2006，31（6）：70－72．

［12］ Knight L，Wang J L．A Comparison of Outlier Detection Procedures and Robust Estimation Methods in GPS Positioning［J］．Journal of Navigation，2009，62（4）：699－709．

［13］ 王新洲，陶本藻，邱衛(wèi)寧．高等測量平差［M］．北京：測繪出版社，2006．

［14］ 周江文．經(jīng)典誤差理論與抗差估計(jì)［J］．測繪學(xué)報(bào)，1989，18（2）：115－120．

［15］ 常志巧，郝金明，張成軍，等．GPS快速定位中病態(tài)問題的正則化抗差解法［J］．大地測量與地球動力學(xué)，2008，28（3）：83－86．

［16］ 文援蘭，楊元喜，王威．衛(wèi)星精密軌道抗差估計(jì)的研究［J］．空間科學(xué)學(xué)報(bào)，2001，21（4）：341－350．

［17］ Ningning Jia，Yonghui Ge．Remainder reliability and robust estimation：A case study using twelve simulated leveling networks，2013（263－266）：3163－3167．

［18］ 黃世存，章文毅，何國金，等．幾種不同矩陣算法的遙感圖像幾何精糾正效果比較［J］．國土資源遙感，2005（3）：18－22，43．