病態(tài)加權(quán)總體最小二乘的廣義嶺估計解法

翁燁，邵德盛

（ 1. 云南省地震局, 昆明 650200；2. 昆明理工大學(xué) 國土資源與工程學(xué)院, 昆明 650093 ）

0 引 言

最小二乘估計[1]理論和經(jīng)典平差模型在大地測量數(shù)據(jù)解算以及工程領(lǐng)域中應(yīng)用非常廣泛，其中高斯-赫爾默特模型(GHM)為經(jīng)典平差模型的常用形式[2]. 在實際應(yīng)用中，坐標轉(zhuǎn)換、大地測量反演、控制網(wǎng)平差、精密單點定位(PPP)等平差模型中的系數(shù)矩陣也含有隨機誤差，經(jīng)典的GHM模型擴展為帶有隨機系數(shù)矩陣的變量誤差模型[3]. 基于這種模型的解算，ADCOCK[4]首次提出了觀測向量和系數(shù)矩陣具有隨機誤差的整體最小二乘估計解法，也叫總體最小二乘法[4].

處理病態(tài)總體最小二乘法問題的常用方法有截斷奇異值法[5]、嶺估計法[6-7]、正則化法[8]以及c-K法[9]等. 對于總體最小二乘法嶺估計的解算方法主要在于嶺參數(shù)的選取，其中嶺參數(shù)的選取方法有：嶺跡法[10]、L-曲線法[11]、U-曲線法[12]以及廣義交叉檢核法[13]等.在等權(quán)條件下總體最小二乘法的解算通常采用正則化方法研究的迭代式，迭代方法較為復(fù)雜，計算量偏大. 針對這種問題，王樂洋等[14]推導(dǎo)了病態(tài)加權(quán)總體最小二乘法平差的嶺估計解式以及精度評定方法，對多種選取嶺參數(shù)方法進行了對比分析. 基于此，本文考慮將病態(tài)加權(quán)總體最小二乘法中的嶺估計解換成廣義嶺估計解，通過傳統(tǒng)平差模型推導(dǎo)了未知參數(shù)的求解公式，利用協(xié)方差傳播率導(dǎo)出參數(shù)解的方差-協(xié)方差矩陣、廣義嶺參數(shù)矩陣K中對角線元素ki(i=1,2,···,m)的迭代解式和均方誤差.

1 加權(quán)總體最小二乘的廣義嶺估計模型

經(jīng)典線性模型為

式中：L為觀測向量；B為設(shè)計矩陣且R(B)=m；X為未 知參 數(shù)向 量； Δ 是 隨 機 向 量 誤 差；P是 權(quán) 矩 陣；Σ=σ2P-1是 協(xié)方差矩陣； σ2為單位權(quán)方差. 設(shè)參數(shù)的估值為觀測值誤差Δ的估計值為觀測值的個數(shù)為n個，未知參數(shù)估計值的個數(shù)為m個，總體最小二乘在考慮觀測向量的誤差之外還要考慮系數(shù)矩陣的誤差，設(shè)EB為系數(shù)矩陣B的誤差向量矩陣，即B=B0+EB，因此經(jīng)典線性模型擴展為變量誤差(EIV)模型[1]

式中：eL∈Rn×1、EB∈Rn×m是觀測向量L0和設(shè)計矩陣B0的誤差矩陣，隨機模型為

式中：Q為 (nm+n)×(nm+n) 階 分塊對角矩陣；vec(·)為矩陣的拉直運算，這里表示矩陣按行拉直得到的列向量；Q為觀測值和系數(shù)矩陣元素之間的協(xié)因數(shù)矩陣為n×n階 方陣，為nm×nm階方陣. 將式(2)展開并且舍去二次項，得

令

加權(quán)總體最小二乘準則為

由于eB是EB通過行拉直所得到的列向量，排列順序由左至右，根據(jù)式(5)的模型即有準則

在系數(shù)矩陣B病態(tài)時，矩陣存在的復(fù)共線性導(dǎo)致法矩陣的求逆變得極不穩(wěn)定，王樂洋等[14]引入嶺參數(shù)k，利用嶺估計的優(yōu)點來處理加權(quán)總體最小二乘的病態(tài)性問題. 但是嶺估計的改進方法對于所有的參數(shù)無區(qū)別統(tǒng)一添加嶺參數(shù)影響，對一部分參數(shù)估計的改進效果不明顯. 在考慮廣義嶺估計的基礎(chǔ)上，根據(jù)Tikhonov正則化原理[15]，得到式(5)的總體最小二乘廣義嶺估計準則為

tr(*) 表示求矩陣的跡，K=diag{k1,k2,···,km} ，結(jié)合式(5)和式(7)的模型準則，通過建立拉格朗日目標函數(shù)，分別對其求一階偏導(dǎo)數(shù)，可以得到加權(quán)總體最小二乘的廣義嶺估計解為

下面借鑒文獻[14]中求協(xié)因數(shù)矩陣的方法求Qˉ ，方法如下，由V=eL-EBX0，e=L0+B0X0可以得

D和G的表達式為

G矩陣是由eL矩陣和EB矩陣組成的增廣矩陣，列向量 v ec(G) 中的所有元素eL1,···,eLn,EB1,1,···,EBn,m之間的互協(xié)因數(shù)及協(xié)因數(shù)組成的矩陣為：

從(3)式中協(xié)因數(shù)Q得 到有Qvec(G)=Q，因此根據(jù)協(xié)因數(shù)傳播率可以得到

2 廣義嶺參數(shù)ki的確定

2.1 嶺參數(shù)的確定——L-曲線法

在1970年Mille首先采用L-曲線法求解病態(tài)方程，其后Hansen[16]又對該方法作了深入研究，L-曲線法是一種在均方誤差有意義下的有偏估計算法，根據(jù)Tikhonov正則化原理，依據(jù)嶺估計的估計準則可以表示為

式中： ‖ *‖ 為歐式2-范數(shù)； Ω (X) 是穩(wěn)定泛函；式(11)在最小二乘估計準則基礎(chǔ)上添加了kXTX項，由于這一項的參與，法方程的病態(tài)性得到了抑制，(BTB+kXTX)求逆變得正常起來，所以相比較最小二乘估計可以得到更加可靠的估計值. 嶺估計結(jié)果優(yōu)于最小二乘估計的原因就在于嶺參數(shù)k，若選擇不同的嶺參數(shù)，得到的估計結(jié)果就會有所不同，所以嶺參數(shù)的選取也至關(guān)重要，接下來介紹關(guān)于嶺參數(shù)的選取方法之一L-曲線法.都是嶺參數(shù)k的函數(shù)，選擇不同的嶺參數(shù)，以為 橫坐標，以為縱坐標畫圖，得到許多點再利用多點曲線擬合的方式擬合出一條曲線，因為該曲線通常形狀都像“L”，所以稱為L-曲線.L-曲線的關(guān)鍵就在于選擇曲線上曲率最大的那個點，其對應(yīng)的嶺參數(shù)即為所求. 接下來以對數(shù)形式推導(dǎo)式如下[11]：

令

兩邊取其對數(shù)形式，得出

則L-曲線是由許多點 ( ξ/2,?/2) 擬合成的曲線，用 ξ′,ξ′′,?′,?′′分別代表 ξ 和 ? 的一階和二階導(dǎo)數(shù)，所以 ξ ,?,ξ′,ξ′′,?′,?′′均是嶺參數(shù)k的函數(shù)，那么，在L-曲線上的點的曲率 η 的計算式為

關(guān)于 ξ′、ξ′′、?′、?′′的計算式可以參考文獻[4-5]，對式(14)求取最大值，得到最大曲率 ηmax， ηmax所對應(yīng)的點的嶺參數(shù)即為所求，這樣就定位出L-曲線上曲率最大的點.

從式(11)可以看出，在應(yīng)用L-曲線法選擇嶺參數(shù)的合理性在于平衡數(shù)據(jù)擬合度部分和解部分這種平衡是通過嶺參數(shù)來實現(xiàn)的，需要指出的是，用L-曲線法求得的嶺參數(shù)不是最優(yōu)選擇，只是近似最優(yōu). 以上的推導(dǎo)過程是在權(quán)矩陣為單位陣的時候推導(dǎo)的，如果權(quán)矩陣不是單位矩陣，由于其都是正定型矩陣，可以先將權(quán)矩陣進行單位化，同理可以計算出L-曲線方法.

2.2 廣義嶺參數(shù) ki 的確定

式中：J為正交方陣，是特征向量矩陣； Λ 是特征值對角方陣， Λ =diag{λ1,λ2,···,λm} ， 按照λ1≥λ2≥···≥λm降序排列；Z=B0J， α =JTx為典則參數(shù). 式(5)模型的典則形式為

再根據(jù)協(xié)方差傳播率，式(8)可以得到其協(xié)因數(shù)矩陣為

均方誤差的定義為

將式(8)帶入到式(18)中，得到參數(shù)改正數(shù)估計值的均方誤差為

可由駐點求得滿足均方誤差最小化情況下，廣義嶺估計的嶺參數(shù)為

根據(jù)式(22)確定的ki可 使 M SE(x?)=min ，但由于未知，因此無法由式(22)直接得到K中的每一個嶺參數(shù)值，可用迭代法求得：

3)最后按照式

由式(22)可知，當(dāng)m個廣義嶺參數(shù)相等且為一常數(shù)，即k1=k2=···=km=k時，就退化成嶺估計，廣義嶺估計相比較嶺估計的不同之處在于通過對法方程系數(shù)陣的主對角元素增加不同的k值來改善其病態(tài)性，而嶺估計是在于法方程系數(shù)矩陣的主對角線元素上加上相同的常數(shù)k來達到改善病態(tài)性的目的. 廣義嶺估計是有偏的最小二乘的線性組合估計值，其主要作用在于使得引入的偏差值小于最小二乘估計的方差值，從而達到提高參數(shù)的估計精度. 嶺估計是對最小二乘估計的改進，而廣義嶺估計是對嶺估計的改進，不過廣義嶺估計與嶺估計有著相同的性質(zhì).

3 算例及分析

采用文獻[17-18]中的數(shù)值算例，在文獻[14]算例的基礎(chǔ)上進行改化，含有5個未知數(shù)，病態(tài)設(shè)計矩陣和觀測向量的真值如下所示[16-17]：

0=[0.8，0.8，0.8，0.8，0.8]T. 分別對觀測值真值和設(shè)計矩陣真值添加模擬誤差eL和eB，其是 設(shè)計矩陣B中的誤差向量矩陣按行拉直后得到的列向量元素矩陣. 法矩陣N=B0TB0的條件數(shù)為cond(N)=2.084 × 104，嚴重病態(tài)性. 這里假設(shè)Pvec(EB)矩陣是一個對角線元素為0.5，其他元素為0.2的方陣，且Pvec(EB)∈R50×50. 為了顯示本文解法的可行性和優(yōu)越性，設(shè)計了多個方案.

方案一：不考慮設(shè)計矩陣B誤差的最小二乘解法；

方案二：考慮設(shè)計矩陣B誤差的總體最小二乘解法；

方案三：不考慮設(shè)計矩陣B誤差的嶺估計解法；嶺參數(shù)采用L-曲線法求得；

方案四：考慮設(shè)計矩陣B誤差的總體最小二乘嶺估計解法；嶺參數(shù)采用L-曲線法求得；

方案五：不考慮設(shè)計矩陣B誤差的總體最小二乘廣義嶺估計解法；

方案六：考慮設(shè)計矩陣B誤差的總體最小二乘廣義嶺估計解法.

圖1 方案三嶺參數(shù)（L-曲線）

圖2 方案四嶺參數(shù)（L-曲線）

表1 不同方法的解算結(jié)果

通過計算對比分析，可以得出以下結(jié)論：

1)模型參數(shù)計算中，最小二乘的嶺估計解優(yōu)于普通的最小二乘解，更加貼近于真值，計算得到方案一和方案三的范數(shù)值分 別為1 .3063 和 0 .8547 .可見嶺估計可以改善最小二乘估計.

2)總體最小二乘嶺估計估計相比較總體最小二乘估計，差值范數(shù)更小，由表1可知方案二和方案四的值 分別為 6 .7322 和 0 .5534 .

總體最小二乘廣義嶺估計優(yōu)于總體最小二乘嶺估計，其差值范數(shù)更小，原因在于廣義嶺估計的嶺參數(shù)選取更加全面，更能削弱系數(shù)矩陣的病態(tài)性. 但是在不考慮設(shè)計矩陣誤差的情況下，考慮設(shè)計矩陣誤差的嶺估計解優(yōu)于普通最小二乘的廣義嶺估計解. 可以看出，設(shè)計矩陣自身的誤差對參數(shù)估計影響較大，因此，考慮設(shè)計矩陣的誤差項不可避免.

3)根據(jù)式(17)得出本文方法的協(xié)因數(shù)矩陣如表2所示. 方案六在考慮設(shè)計矩陣誤差的情況比方案五得出的參數(shù)估計值更加貼近于真值，方案五和方案六的差值范數(shù)分 別為 0 .6248 和 0 .4476 .

表2 方案六的協(xié)因數(shù)矩陣

4 結(jié) 語

EIV模型同時顧及了設(shè)計矩陣的誤差項和觀測向量的誤差項，適用解法為總體最小二乘估計，總體最小二乘估計優(yōu)于最小二乘估計. 但是系數(shù)矩陣存在復(fù)共線性時，一般的總體最小二乘估計也難以得到穩(wěn)定的估計值，且估計值的均方誤差也不理想，因此必須要針對設(shè)計矩陣的病態(tài)性進行處理. 在考慮設(shè)計矩陣自身誤差項的基礎(chǔ)上，對法方程系數(shù)矩陣添加改正項，借鑒廣義嶺估計的優(yōu)勢，推導(dǎo)出在病態(tài)加權(quán)下總體最小二乘的廣義嶺估計直接解法，廣義嶺參數(shù)采取迭代求解法，優(yōu)于總體最小二乘的嶺估計解，在處理病態(tài)加權(quán)的變量含誤差模型有一定的借鑒作用.