翁燁，邵德盛

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

0 引 言

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

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

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

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

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

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

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

令

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

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

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

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

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

D和G的表達(dá)式為

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ì)算法，根據(jù)Tikhonov正則化原理，依據(jù)嶺估計(jì)的估計(jì)準(zhǔn)則可以表示為

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

令

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

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

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

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

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

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

3)最后按照式

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

3 算例及分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié) 語

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