馮壯壯,徐 工
(山東理工大學 建筑工程學院,山東 淄博 255049)
病態問題廣泛存在于大地測量控制網平差、航空重力向下延拓和GPS快速定位等領域中,線性方程組設計矩陣具有嚴重的復共線性,數據微小的變化對導數參數估計結果產生較大的影響[1]。雖然線性最小二乘具有線性無偏性和方差最小性,但由于小奇異值影響方差被嚴重放大,得不到所期望的結果。目前處理線性最小二乘病態問題的算法多以正則化理論為基礎,通過修正設計矩陣病態性或者避免矩陣求逆來提高參數估計的解[2-3]。較常用的線性正則化算法有截斷奇異值方法(TSVD)、Tikhonov算法和兩參數LIU估計等[4-5]。正則化算法是通過引入偏差量來降低方差,從而提高參數解的精度和穩定性[6]。
TSVD實質為分割法,該方法將小奇異值及其對應的奇異分量刪除,保留大奇異值及其對應的主成分信息,雖然能夠提高解的穩定性和精度,但降低了解的分辨率[7]。標準Tikhonov方法對線性方程組設計矩陣無差異修正,引入了偏差量,降低了參數估計的可靠性[8]。此外,截斷參數和正則化參數的確定問題是制約線性正則化算法求解精度和穩定性的主要因素,常用的正則化算法有L曲線、方差最小法、和GCV廣義交叉法[9-11]。但由于問題復雜性,各種算法都有其實用性和局限性,即使最為經典的L曲線法,針對某些問題也會出現無法收斂的問題,為此開展新的正則化參數確定方法顯得尤為重要,比如最近基于L曲線法發展起來的U曲線法,克服了L曲線法的不足,能夠更為準確的定位正則化參數[12-13]。
本文針對TSVD方法和標準Tikhonov方法存在的不足,基于兩種算法處理病態問題的思想,根據奇異值分布空間差異性,提出了一種混合正則化方法。該方法主要是將奇異值進行分割,保留大奇異值及其主成分信息,采用正則化參數來修正小奇異值及其奇異分量,以提高解的精度和參數估計的可靠性,最后采用實驗驗證了該方法的有效性。
設觀測向量為L=[L1,L2,…,Ln]T,改正向量為V=[V1,V2,…,Vn]T,則誤差方程為
V=Ax-L
(1)
式中,A為系數矩陣,x為待求參數向量。線性最小二乘目標函數為
φ(x)=VTV
(2)
根據最小二乘定理,按照拉格朗日求極值的思想,可解得線性最小二乘估計解為
xLS=(ATA)-1ATL
(3)
參數估值的協方差陣為
(4)

(5)
式中,λ1>λ2>…>λm為設計矩陣的奇異值。最小二乘估計解在線性估計類具有方差最小性。然而方程組設計矩陣呈病態時,較小的奇異值能夠過度放大參數估值的方差,導致最小二乘求解性質變壞。為準確獲取參數估計的最優解,改善系數矩陣秩虧或病態對參數估值的影響,可引入穩定泛函極小準則,將最小二乘優化問題轉化為
φ(xtik)=VTV+αxTRx
(6)
式中,α為正則化參數,用來平衡解的不穩定性和光滑性;R為正則化矩陣,主要是平滑參數分量和殘差分量。根據拉格朗日求極值的思想,基于穩定泛函極小準則可得Tikhonov正則化解為
xtik=(ATA+αR)-1ATL
(7)
Tikhonov正則化是一類通過修正病態矩陣而降低均方誤差的有偏估計理論。嶺估計是正則化矩陣為單位陣的Tikhonov正則化的特例,稱為標準Tikhonov正則化,其估計解為
xtik=(ATA+αI)-1ATL
(8)
嶺估計在修正病態矩陣的同時,不可避免的引入了偏差量。嶺估計參數估值的協方差陣為
(9)
則參數估值協方差陣的跡為
(10)
嶺估計的偏差計算公式為
bias(xα)=E(xα)-x=
(ATA+αI)-1ATEL-x=
((ATA+αI)-1ATA-I)x=
-α(ATA+αI)-1x
(11)
綜上分析,嶺估計無差別修正雖能夠良好的修正小奇異值的奇異分量,但不可避免的給大奇異值的主成分信息摻入了虛假信息,以致于標準Tikhonov正則化算法有效降低LS估計的均方誤差的同時不可避免的引入了偏差。為此,正則化矩陣的合理構造是提高Tikhonov正則化估計解的可信度的重要手段之一。

(12)
式中,k為截斷參數且1≤k≤m,該值主要是根據奇異值和正則化參數值的相對關系進行確定,該方法能夠將小奇異和大奇異值進行分割處理,這樣既能夠保留大奇異值及其對應的主成分信息,也能夠采用Tikhonov正則化方法對小奇異值及其對應的奇異信息進行修正。
(13)
式(13)稱為混合正則化方法,該方法能夠保證解的分辨率,又能降低偏差量的引入,相比于嶺估計和截斷奇異值法可靠性更高。
正則化估計的作用效果主要取決于正則化參數的選擇。獲取正則化參數較常用的方法,主要有GCV交叉檢驗法、L曲線法和U曲線法,本文將采用U曲線法確定正則化參數。U曲線法基于SVD分解理論,由正則化參數α>0的殘差范數的倒數平方和信號范數的倒數平方之和定義U(α)。以U(α)為縱軸,α為橫軸所表示的曲線形狀U曲線,故稱為U曲線法,定義為
(14)

(15)
式中,U(α)′、U(α)″分別是U(α)的一階導數和二階導數,并且
(16)
采用文獻[14]模擬空間測邊網的案例。設P1,P2,...,P9為9個已知點,且已知點到待測點P10,P11(模擬值分別為(0,0,0) 和(7,10,-5)的觀測距離即表1所示,兩個待測點的觀測距離為13.107 85,觀測圖形如圖1所示,要求根據19個觀測距離確定待測點的坐標。距離為等精度觀測,中誤差為±0.001 m。假定坐標近似值為(0.01,-0.01,0.02)和(7.01,9.99,-5.01)。經計算設計矩陣條件數為cond(N(x))=4.601 5×103,表明方程組具有嚴重的病態性。
表1 測邊網控制點坐標和已知觀測距離
Tab.1 Control point coordinates of trilateration network and known observation distance

點號坐標觀測距離XYZdi,10di,11P123.00010.000 0.01025.078 616.765 17P2- 10.000 9.9900.00014.134 517.719 65P335.00010.010- 0.01036.415 928.442 94P4100.00019.9900.005101.479 493.168 39P5- 36.00010.005- 0.00037.364 243.299 05P60.00010.010- 0.00510.010 08.600 60P756.000 9.9950.01056.996 149.256 18P8- 15.00010.015- 0.01018.035 922.559 66P9-1.700 010.0080.01510.150 610.043 82

圖1 測邊網平面觀測圖形Fig.1 Plane observation graph of trilateration network
圖2給出了U曲線確定正則化參數的變化曲線圖,由圖可知,U曲線確定正則化參數相對于L曲線法不同之點,是U曲線存在兩個極值點,通過兩個極值點去逼近最佳正則化參數,能夠克服L曲線法無法收斂的情況。經確定,該案例最佳正則化參數為0.4130。分別采用LS估計、TSVD、嶺估計和混合正則化進行解算,對比其偏差量bias和均方誤差RMSE,結果詳見表2。

圖2 U曲線確定正則化參數的變化曲線圖Fig.2 The curve diagram of U curve determining regularization parameter
表2給出不同方法的對比結構。由表能夠看出,LS估計由于線性方程組病態的影響,其估計解以偏離真實解。嶺估計和截斷奇異值方法均提高了線性估計解的精度,但嶺估計作用于小奇異值及其奇異值分量時,也會給大奇異值及其主成分信息摻雜虛假信息,降低了解的可靠性。截斷奇異值方法降低了解的分辨率。本文提出的方法,所計算的結果優于嶺估計和截斷奇異值法,相比于兩種方法可靠性更好。這是因為提出的方法能夠降低了偏差的攝入量,減少對解分辨率的影響。
表2 不同方法的解空間對比表
Tab.2 Solution space comparison table for different methods

參數真值LS估計/m嶺估計/mTSVD/m新方法/mX1 0-0.036 8-0.032 8-0.064 8-0.044 6Y100.051 80.008 0-0.017 10.029 5Z109.363 10.295 40.143 20.088 2X277.048 86.980 36.977 76.973 2Y2105.596 39.882 19.687 810.038 1Z2-5-4.648 2-4.914 9-4.960 4-4.871 7Bias010.353 30.331 60.352 90.171 1RMSE04.226 70.135 40.144 10.069 8
線性方程組求解過程中,若設計矩陣具有病態性,即使最小二乘具有線性無偏性和線性無偏估計類的方差最小性,但也會由于奇異值及其奇異分量的影響,方差被過度放大,最小二乘性質變壞。截斷奇異值方法雖能夠提高線性估計解的精度,但會降低解的分辨率,嶺估計無差異修正方式,會給大奇異值及其主成分信息摻雜虛假信息,參數估計產生偏差,降低解的可靠性。本文提出的方法,能夠分割處理大奇異值和小奇異值對應的成分信息,降低偏差量的攝入,提高了參數估計的可靠性和精度,最后實驗驗證了該方法的解算精度明顯優于嶺估計和TSVD方法。