楊秋偉,陳華,周聰,李翠紅

(1.紹興文理學院 土木工程系 浙江 紹興 312000;2. 寧波工程學院 建筑與交通工程學院,浙江 寧波 315211;3. 浙江省土木工程工業(yè)化建造工程技術研究中心(寧波工程學院)，浙江 寧波 315211)

0 引 言

在測量平差問題中,最小二乘估計是求解線性方程組最常用的方法之一,例如,考慮一個常見的線性估計模型:

y=A·x,

(1)

式中：向量y為已知的觀測向量(m×1維);A為系數(shù)矩陣(m×n維列滿秩矩陣);向量x為未知的參數(shù)向量(n×1維).利用最小二乘法求解方程(1)的過程如下:

對方程(1)兩邊乘以AT,可得法方程為

z=B·x,

(2)

式中,B=ATA,z=ATy.由方程(2)可得x的最小二乘估計為:

xlse=B-1·z.

(3)

由方程(3)所得的最小二乘解滿足最優(yōu)線性無偏性,因此最小二乘估計又稱為無偏估計.然而,當線性方程的系數(shù)矩陣A存在復共線性時,即A中某些列向量相關度很高時,由方程(3)所得的解可能誤差很大甚至完全失真,這種情況稱為病態(tài)最小二乘問題[1-5].為了得到準確的x估值,必須研究其他的穩(wěn)健算法.

1 奇異值截斷方法

如前所述,針對病態(tài)最小二乘問題,學者們已提出各種穩(wěn)健估計算法.其中,奇異值截斷方法[6-11]是最具代表性的一類方法.這類方法普適性強,既可適用于方程數(shù)目大于或等于未知量數(shù)目的情況(系數(shù)矩陣A列滿秩),也可適用于方程數(shù)目小于未知量數(shù)目的情況(系數(shù)矩陣A行滿秩),還可適用于系數(shù)矩陣為虧損矩陣的情況.其理論依據(jù)是,在方程組系數(shù)矩陣的所有非零奇異值中,較大的奇異值對數(shù)據(jù)噪聲不敏感,而較小奇異值(接近于0的奇異值)體現(xiàn)了系數(shù)矩陣的復共線性,且對數(shù)據(jù)噪聲很敏感.因此,在方程求解時忽略較小的奇異值而只保留較大的奇異值,相當于一定程度上消除了系數(shù)矩陣的復共線性,增強了抵抗數(shù)據(jù)噪聲干擾的能力,從而可以獲得比較穩(wěn)定的計算結果.

利用奇異值截斷方法求解線性方程組的過程簡述如下.不失一般性,考慮一個任意的線性方程組如下：

y=C·x.

(4)

與方程(1)不同,方程(4)中的系數(shù)矩陣C對矩陣維數(shù)和是否滿秩沒有特殊要求.接下來首先對方程(4)中的系數(shù)矩陣C做奇異值分解,即

C=UΛVT,

(5)

U=[u1,u2,…],V=[v1,v2,…],

(6)

(7)

式中,σ1,σ2,…,σp為矩陣C的p個非零奇異值,且滿足σ1≥σ2≥…≥σp.方程(5)代入方程(4)中可得x的解為

(8)

由方程(8)可見,越小的奇異值倒數(shù)越大,其微小變化將引起解的很大變化.因此,為了提高解的穩(wěn)定性,在方程(8)中可以舍去較小的奇異值,而只保留較大的奇異值.假設只保留前t個奇異值,t又稱為截斷閾值,則截斷奇異值后所得的解為

(9)

方程(9)即為奇異值截斷方法計算x的公式.顯然,決定該方法計算精度的關鍵是確定合適的截斷閾值t.不少學者就此開展研究工作,代表性的方法有F檢驗法[10],廣義交互確認法(GCV)[12],L曲線法[13]等.但這些方法均需要花費額外的計算量,操作步驟上比較復雜,且容易出現(xiàn)誤判.

2 奇異值修正方法

鑒于奇異值截斷閾值選擇過程比較復雜,增加的額外計算量較多,有些學者提出了奇異值修正的方法[9,12-16],其核心思想是對較小的奇異值進行適量的修正(增大),從而使得結果更為穩(wěn)定.這些奇異值修正方法所提的修正公式基本都是分段函數(shù),比如文獻[9]中所提修正公式為

(10)

文獻[14-16]中所用的修正公式為

(11)

顯然,這些修正公式中仍然需要確定一個修正閾值q,來對不同的奇異值選擇相應的修正函數(shù). 確定修正閾值q的計算量仍然與確定截斷閾值t的計算量相當,在應用上仍然不太方便,尤其是對于大規(guī)模的線性方程組,計算過程很復雜.因此,本文提出一種新的奇異值修正公式,對所有的非零奇異值都采用同樣的修正公式,形式如下

(12)

(13)

顯然,采用本文所提的奇異值修正公式,并不需要確定任何閾值,基本沒有增加額外的計算量,有利于工程應用.由后繼的數(shù)值算例可見,本文方法可以獲得比奇異值截斷法精度更高的計算結果.

3 數(shù)值算例

3.1 病態(tài)平差模型

以文獻[17]所用的病態(tài)平差模型y=A·x為例,驗證所提的奇異值修正法,該方程的系數(shù)矩陣A和觀測向量y的真值為

(14)

未知參數(shù)向量x的真值為

yc=[-10.178 1, 10.640 4, 1.372 4, 12.05, 13.607 7,-0.112 2, 20.942 2, 1.554, 9.328 6, 12.112 3]T.

采用本文所提奇異值修正法,即利用方程(12)和(13),所得的估值xnew的元素值列于表1中.另外,為了比較本文方法、最小二乘法和截斷奇異值方法,表1中同時給出了最小二乘估計xlse的元素值和取最優(yōu)截斷閾值t=3時的x估值xt的元素值.為了評估這三種解與真值之間的接近程度,表1中還給出了各個估值與真值之差的2-范數(shù)eΔx.

表1 本文方法與最小二乘估計和截斷奇異值方法的比較

由表1可見,最小二乘估計xlse與真值差別很大,已完全失真.而采用本文所提奇異值修正方法,所得結果(xnew)的元素值要比奇異值截斷法(xt)的元素值更準確,且本文方法計算公式十分簡單快捷,并不需要任何額外的運算,方便工程應用.

3.2 Hilbert病態(tài)矩陣

Hilbert病態(tài)矩陣[18]經(jīng)常用于測試各種穩(wěn)健估計算法的可靠性,該矩陣由如下公式產(chǎn)生:

(15)

式中,i、j為該矩陣中元素hij的行號和列號.隨著矩陣維數(shù)的增大,Hilbert矩陣的病態(tài)性也會越來越嚴重.以Hilbert矩陣為系數(shù)矩陣,建立線性方程如下:

Hm×n·x=y.

(16)

不失一般性,假設該方程中x的真值取為全部元素等于1的列向量,而觀測向量y取為Hm×n·(1,1,…,1)T.現(xiàn)利用本文方法來對方程(16)進行求解,為了說明所提方法的普適性,考慮三種類型的Hilbert矩陣,分別為:行數(shù)大于列數(shù)、行數(shù)等于列數(shù)和行數(shù)小于列數(shù),不失一般性,分別取系數(shù)矩陣為H30×20、H20×20和H15×20,利用本文方法的計算結果xnew的元素值列于表2中.另外,為了便于比較,表2中也給出了部分最小二乘估計和奇異值截斷解.

由表2的第2和3列可見,系數(shù)矩陣取為H30×20時(行數(shù)大于列數(shù)),即使方程中不含有數(shù)據(jù)噪聲,用最小二乘法計算的結果xlse的元素值也是完全失真的,因為系數(shù)矩陣的病態(tài)性相當嚴重,而采用本文方法所得結果與真值高度一致.因此,對于病態(tài)非常嚴重的方程組,最小二乘法已完全不能使用,故對于系數(shù)矩陣取H20×20和H15×20的情況不再給出最小二乘解.對于后兩種情況,表2的第4和6列分別給出了奇異值截斷法的計算結果xt的元素值,和本文方法計算結果對比發(fā)現(xiàn),奇異值截斷法也能取得較好的精度,但顯然操作步驟上比本文方法更為繁瑣,也增加了額外的計算量.總體而言,本文方法綜合效率最好.

4 結 論

提出了一種統(tǒng)一的奇異值修正公式,以此為基礎提出的奇異值修正法,可用于解決測量平差中的病態(tài)最小二乘問題.和現(xiàn)有的奇異值截斷法或修正法相比,所提新方法更簡便快捷,且能獲得比奇異值截斷法精度更高的計算結果.另外,所提方法普適性強,對方程系數(shù)矩陣的維數(shù)和是否滿秩沒有特殊的要求.以兩個病態(tài)方程為例對所提方法進行了驗證,并將計算結果與最小二乘解和奇異值截斷解進行了比較,結果說明了所提方法的有效性和優(yōu)越性.