吳光明,魯鐵定,2,3,鄧小淵,4,吳建江,5

(1.東華理工大學(xué) 測繪工程學(xué)院,江西 南昌 330013;2.流域生態(tài)與地理環(huán)境監(jiān)測國家測繪地理信息局重點(diǎn)實(shí)驗(yàn)室,江西 南昌 330013;3.江西省數(shù)字國土重點(diǎn)實(shí)驗(yàn)室,江西 南昌 330013;4.浙江省地理信息中心，浙江 杭州 310000;5.浙江省地球物理技術(shù)應(yīng)用研究所,浙江 杭州 310014)

病態(tài)問題在測量數(shù)據(jù)處理過程例如在控制網(wǎng)平差[1]、GPS快速定位[2-3]、航空重力向下延拓[4-6]等領(lǐng)域時(shí)常碰見。處理病態(tài)問題常見方法有Tikhonov正則化法[7-8]、嶺估計(jì)[9-10]等，然而嶺估計(jì)和正則化法改變了方程的等量關(guān)系，使得估計(jì)結(jié)果有偏并且?guī)X參數(shù)和正則化參數(shù)確定困難[11-12]。

1 譜修正迭代法

對于經(jīng)典平差模型

L+Δ=AX.

(1)

式中：L為n×1觀測值向量，A為m×n系數(shù)矩陣，X為n×1待定參數(shù)，Δ為觀測向量的噪聲，Δ～N(0,σ02I)。其最小二乘估計(jì)及估計(jì)的協(xié)方差為

(2)

(3)

當(dāng)法矩陣AΤA出現(xiàn)病態(tài)(條件數(shù)一般大于103)，則法矩陣求逆將會表現(xiàn)得不穩(wěn)定，導(dǎo)致求解出的參數(shù)估值不可靠[12]。

(4)

或者

(5)

(6)

2 靶向譜修正迭代法

譜修正迭代中的(ATA+αI)-1與嶺估計(jì)相似，嶺估計(jì)式為

(7)

(8)

整理后

(9)

式中：R是靶向矩陣，其構(gòu)造方法是基于較小特征值對應(yīng)的特征向量而構(gòu)成，

(10)

式中：R1是正則化矩陣，Gi是法矩陣ATA的特征向量。小特征值的判定方法可用特征值標(biāo)準(zhǔn)差分量之和占標(biāo)準(zhǔn)差比重達(dá)到95%以上[19]，即

(11)

式中：Λi是法矩陣ATA的特征值。將R1作為靶向矩陣并代入式(8)則整理得到本文靶向譜修正迭代式

(12)

譜修正迭代結(jié)果是無偏估計(jì)，根據(jù)文獻(xiàn)[13]證明方法，本文靶向修正迭代結(jié)果也是無偏的。只要迭代矩陣α(ATA+αR)-1·R的譜半徑小于1，迭代就能收斂[20]。在對α(ATA+αR)-1·R進(jìn)行譜分解，得到

α(ATA+αR)-1·R=

(13)

譜修正參數(shù)α選擇有多種方法[16-18],沒有統(tǒng)一結(jié)論。因此本文根據(jù)文獻(xiàn)[18]的方法，選擇一個(gè)較長區(qū)間、一定步長的參數(shù)，討論不同參數(shù)條件下，兩種方法迭代結(jié)果的比較。

3 算例及分析

3.1 算例1

采用文獻(xiàn)[4]中的模擬病態(tài)問題算例，法矩陣條件數(shù)是4.1847×105,嚴(yán)重病態(tài)，其中未知參數(shù)的真值為X=[111111111]T。為比較兩種方法的估計(jì)結(jié)果，取迭代初值為X=[0.80.80.80.80.80.80.80.80.8]T或其他初值均可，不能選擇最小二乘估值作為迭代初值，若選擇最小二乘估值將無法迭代。根據(jù)文獻(xiàn)[18]的方法，取譜修正參數(shù)α=1，10，20，30，分別用最小二乘估計(jì)、譜修正迭代法、靶向譜修正迭代法對這個(gè)問題進(jìn)行解算，結(jié)果見表1。

表1 兩種算法的解算結(jié)果

圖1 兩種方法結(jié)果

從圖1可以看出，兩種方法的差值范數(shù)基本相等，結(jié)果基本一致；見圖2(a)圖是α=1～50，由于豎軸刻度較大，兩種方法迭代次數(shù)也近乎相等；由圖2(b)看出，在α>11時(shí)，本文方法迭代次數(shù)少于譜修正迭代法。根據(jù)這個(gè)算例，靶向譜修正迭代法與譜修正迭代法相比，在參數(shù)估計(jì)和差值范數(shù)上沒有較大變化；但在參數(shù)逐漸增大時(shí)，迭代次數(shù)降低，計(jì)算效率提高，體現(xiàn)出本文方法的優(yōu)勢。

圖2 兩種方法迭代次數(shù)

3.2 算例2

表2 兩種算法的解算結(jié)果

從圖3可以看出，圖3(a)是α=1～17，兩種方法的差值范數(shù)基本相等；由圖3(b)看出，在α>17時(shí)，本文方法差值范數(shù)小于譜修正迭代法。見圖4，1<α<17時(shí)，兩種方法迭代次數(shù)基本相等；而α>17時(shí)，本文方法迭代次數(shù)遠(yuǎn)小于譜修正迭代法。這個(gè)算例表明，參數(shù)逐漸增大，兩種方法在參數(shù)估計(jì)、差值范數(shù)、迭代次數(shù)上先是沒有較大變化，但參數(shù)繼續(xù)增大時(shí)，本文方法解算的差值范數(shù)和迭代次數(shù)均降低，進(jìn)一步驗(yàn)證靶向譜修正迭代法高效、快速解算。

圖3 兩種方法解算結(jié)果

圖4 兩種方法迭代次數(shù)

4 結(jié)束語

譜修正迭代法是修正法矩陣所有譜且結(jié)果是無偏估計(jì)，而一般病態(tài)問題是法矩陣的幾個(gè)譜奇異，因此存在譜多余修正問題。針對該問題，本文提出靶向譜修正迭代法，即將譜修正迭代矩陣中的單位陣變換成靶向矩陣，迭代過程中只修正法矩陣奇異的譜。通過模擬算例對比兩種方法解算的估計(jì)結(jié)果、偏差范數(shù)和迭代次數(shù)，得到以下結(jié)論：

1)譜修正參數(shù)相同條件下，兩種方法迭代矩陣的譜半徑相等，且估計(jì)結(jié)果均是無偏的；

2)伴隨著參數(shù)增大，本文方法的解算的參數(shù)估值將優(yōu)于譜修正迭代，并且偏差范數(shù)更低、迭代次數(shù)相比也更少。

因此，在探討處理病態(tài)問題時(shí)，相比譜修正迭代法，本文方法在估計(jì)結(jié)果更優(yōu)且快速迭代，計(jì)算效率高。本文遞增選取參數(shù)，參數(shù)選擇未深入分析，因而參數(shù)如何選取有待進(jìn)一步研究。