線性化通用EIV平差模型的嶺估計(jì)解法

翁燁，邵德盛,2，甘淑,3

( 1. 昆明理工大學(xué) 國土資源工程學(xué)院, 昆明650093；2. 云南省地震局, 昆明 650041；3. 云南省高校高原山地空間信息測(cè)繪技術(shù)應(yīng)用工程研究中心, 昆明 650093 )

0 引 言

EIV(errors-in-variables)的多元非線性模型[1]不僅考慮了觀測(cè)值的誤差，還顧及了系數(shù)矩陣的誤差，大部分模型的系數(shù)矩陣并非常數(shù)陣，而是具有一定誤差的變量矩陣. 總體最小二乘法[2](TLS)自Golub等[3]在1980年引入測(cè)量領(lǐng)域之后，開始應(yīng)用于控制網(wǎng)平差、大地測(cè)量反演、多系統(tǒng)組合定位及航空攝影測(cè)量等方面，是一種常用于EIV模型中參數(shù)的估計(jì)方法.張俏等[4]針對(duì)北斗衛(wèi)星導(dǎo)航系統(tǒng)(BDS)整周模糊度在經(jīng)典最小二乘估計(jì)(LS)平差模型中的解算問題，研究出一種利用TLS快速解算雙差整周模糊度的改進(jìn)LAMBDA算法. 武曙光等[5]關(guān)于高精度GPS超短基線場(chǎng)，運(yùn)用GAMIT/GLOBK軟件，采用不同解算方案進(jìn)行了高精度數(shù)據(jù)處理. 何安良等[6]在高程異常擬合問題中，針對(duì)加權(quán)總體最小二乘(WTLS)法估計(jì)GPS高程擬合中沒有考慮到觀測(cè)數(shù)據(jù)存在粗差的情況，采用穩(wěn)健WTLS求解GPS高程平面擬合參數(shù).部分國產(chǎn)衛(wèi)星影像定位可能缺少實(shí)測(cè)控制數(shù)據(jù)，在區(qū)域網(wǎng)平差時(shí)可以采用TLS 法進(jìn)行解算[7].

TLS主要應(yīng)用于EIV模型的間接平差或附有限制條件的間接平差，即觀測(cè)向量的系數(shù)矩陣為單位陣；間接平差模型作為四種經(jīng)典平差模型之一，并非通用平差模型，經(jīng)過一些學(xué)者的思量，考慮從高斯-赫爾默特模型(GH)變換成通用的EIV模型[8]. 文獻(xiàn)[8-9]中提出了通用EIV模型，將GH模型觀測(cè)向量和參數(shù)向量的系數(shù)矩陣由固定矩陣推廣到隨機(jī)矩陣，包含隨機(jī)系數(shù)矩陣的各類情況. 曾文憲等[8]提出了基于TLS估計(jì)準(zhǔn)則下通用EIV模型的WTLS，并推導(dǎo)出了參數(shù)估計(jì)值精度的近似計(jì)算公式，缺點(diǎn)在于參數(shù)估計(jì)量較多時(shí)迭代計(jì)算量偏大、耗時(shí)長(zhǎng). 文獻(xiàn)[10]中針對(duì)非線性通用EIV模型的計(jì)算量偏大問題，提出了一種利用非線性平差原理，將模型展開的二次項(xiàng)納入平差方程的常數(shù)項(xiàng)，從而將非線性通用EIV模型變成經(jīng)典GH模型，推導(dǎo)出了通用EIV模型的線性化總體最小二乘(LTLS)求解式. 通用EIV模型在任意權(quán)下TLS雖然具有一般性，其優(yōu)點(diǎn)在于結(jié)構(gòu)更為嚴(yán)謹(jǐn)，涵蓋了傳統(tǒng)的四種平差模型形式[8]，缺點(diǎn)在于觀測(cè)向量和參數(shù)向量的系數(shù)矩陣都病態(tài)時(shí)，TLS解的穩(wěn)定度會(huì)產(chǎn)生偏移，估計(jì)值的均方誤差也會(huì)隨之增加.

嶺估計(jì)可以有效地改善病態(tài)矩陣的病態(tài)性，在模型最小化準(zhǔn)則中添加穩(wěn)定泛函，有效減小參數(shù)估計(jì)解的穩(wěn)定度偏移量，得到更為穩(wěn)定的參數(shù)估計(jì)解. 文獻(xiàn)[11]中針對(duì)嶺估計(jì)無法單次計(jì)算使得均方誤差達(dá)到最小問題，提出了關(guān)于嶺估計(jì)的迭代解法；文獻(xiàn)[12]中采用嶺估計(jì)法處理LTLS平差的病態(tài)性問題并推導(dǎo)了相應(yīng)的求解公式及均方誤差評(píng)定精度的方法；文獻(xiàn)[13]中基于抗差嶺估計(jì)模型，介紹了應(yīng)用L曲線法對(duì)系數(shù)陣病態(tài)和觀測(cè)粗差的影響，推導(dǎo)出抗差嶺估計(jì)模型的計(jì)算公式. 嶺參數(shù)的確定至關(guān)重要，嶺參數(shù)起到平衡正則化矩陣與病態(tài)設(shè)計(jì)矩陣的作用，進(jìn)一步反應(yīng)病態(tài)設(shè)計(jì)矩陣的權(quán)重[14]. 關(guān)于嶺參數(shù)的確定方法，國內(nèi)外研究學(xué)者提出了許多有效地確定方法，如嶺跡法、L曲線法[15]、U曲線法[16-17]、廣義交叉核實(shí)法[18]、方差分量估計(jì)法[19]等.

本文在線性化通用EIV平差模型解的基礎(chǔ)上添加嶺估計(jì)準(zhǔn)則，利用L曲線法和U曲線法確定嶺參數(shù)，利用算例驗(yàn)證了嶺估計(jì)的促進(jìn)作用；發(fā)現(xiàn)嶺估計(jì)可以有效減少參數(shù)解的迭代次數(shù)，使得參數(shù)估計(jì)方差分量收斂更快.

1 通用EIV模型的線性化估計(jì)

GH模型為

式中：L0和vL為n×1 階觀測(cè)向量和觀測(cè)向量的改正數(shù)向量；X∈Rm×1是 未知參數(shù)向量；A為f×n階觀測(cè)值向量對(duì)應(yīng)的系數(shù)矩陣；B為f×m階參數(shù)向量對(duì)應(yīng)的設(shè)計(jì)矩陣，且 R (B)=m＜f；e為f×1 階常數(shù)向量. 在經(jīng)典平差函數(shù)模型中，常將A化 作單位矩陣，且認(rèn)為在A和B矩陣中不含有隨機(jī)誤差，是固定矩陣.

由文獻(xiàn)[1]可知，若在參數(shù)向量的系數(shù)矩陣B中含有隨機(jī)誤差，式(1)就變成經(jīng)典的EIV模型；在觀測(cè)向量的系數(shù)矩陣A中也含有隨機(jī)誤差時(shí)，式(1)就變成了一般通用的EIV式[8-9]

式中：vA和vB分 別是f×n階 、f×m階A、B矩陣的改正數(shù)矩陣；其余符號(hào)意思同式(1)；A、L0、B為含有隨機(jī)誤差的數(shù)據(jù)矩陣，因此式(2)可以轉(zhuǎn)變成：

式中： v ec(.) 是矩陣的拉直符號(hào)，這里做列向量拉直；為未知單位權(quán)方差；L、v表 示t×1 階的觀測(cè)向量及其改正數(shù)向量，其中t=fn+fm+m；D(L) 、P和Q表示L的方差-協(xié)方差矩陣、權(quán)矩陣和協(xié)因數(shù)矩陣.

依據(jù)LS最優(yōu)化準(zhǔn)則，將式(2)定義為

WTLS在參數(shù)估計(jì)量較多時(shí)計(jì)算量偏大，關(guān)于通用EIV平差模型的LTLS解式可以參考文獻(xiàn)[10]，利用二次項(xiàng)作為模型誤差帶入到常數(shù)項(xiàng)，線性化后的通用EIV模型解法上等價(jià)于GH模型，可以用經(jīng)典LS得到參數(shù)解，減少了計(jì)算量. 設(shè)參數(shù)的估計(jì)值為

將通用的EIV式(2)展開得出

2 通用EIV平差模型的嶺估計(jì)解

2.1 嶺估計(jì)法

由式(7)可知，線性化后的模型與GH模型形式一致，亦可采用LS估計(jì)方法得到觀測(cè)向量和參數(shù)解值，推導(dǎo)過程如下：

依據(jù)LS原理有

通過建立相應(yīng)的目標(biāo)函數(shù)

式中， μ 為f×1 階Lagrange向量，將上式分別對(duì)其函數(shù)元素求一階偏導(dǎo)數(shù)，并令其為0有：

由式(12)～(13)轉(zhuǎn)置可得：

式(7)、(14)～(15)中，共有(n+m+f)個(gè)方程，待求的未知數(shù)為m個(gè)改正數(shù)、n個(gè) 參數(shù)和f個(gè)聯(lián)系數(shù)，即方程個(gè)數(shù)等于未知參數(shù)的個(gè)數(shù)，所以可以由vTPv=min求得滿足條件的一組唯一解值. 用P-1左乘式(15)得到

將式(16)作為改正數(shù)方程，于是基礎(chǔ)方程變?yōu)?/p>

解算基礎(chǔ)方程，通常在于將其中的改正數(shù)方程帶入原始方程，得到一組包含有拉格朗日乘子向量的對(duì)稱線性方程組，即有

得到參數(shù)估計(jì)的協(xié)因數(shù)矩陣和方差-協(xié)方差矩陣為:

μ?作為關(guān)聯(lián)系數(shù)，可以不必計(jì)算，直接將式(20)帶入到式(16)，得到

通用EIV模型下的LTLS算法步驟如下：

2.2 L 曲線法

1970年Mille首先采用L曲線法求解病態(tài)方程，其后Hansen和Peter R等又對(duì)該方法作了深入研究，并于上世紀(jì)90年代發(fā)表了多篇有關(guān)該方法的論文著作. 至今，L曲線法逐漸被業(yè)界認(rèn)同采用，它是一種在均方誤差有意義下的有偏估計(jì)算法，根據(jù)Tikhonov正則化原理，依據(jù)嶺估計(jì)的估計(jì)準(zhǔn)則可以表示為

與LS相比較，嶺估計(jì)中多添加了kI一項(xiàng)，由于這一項(xiàng)的參與，法方程的病態(tài)性得到了抑制，(N+kI)求逆變得正常起來，所以相比較LS能得到更加可靠的估計(jì)值. 嶺估計(jì)的結(jié)果優(yōu)于LS的部分就在于嶺參數(shù)k，若選擇不同的嶺參數(shù)，得到的估計(jì)結(jié)果就會(huì)有所不同，所以嶺參數(shù)的選取也至關(guān)重要. 在式(26)中都是嶺參數(shù)k的函數(shù)，選擇不同的嶺參數(shù)，以為 橫坐標(biāo)，為縱坐標(biāo)畫圖，得到許多點(diǎn)再利用多點(diǎn)曲線擬合的方式擬合出一條曲線，因?yàn)樵撉€通常形狀都像“L”，所以稱為L(zhǎng)曲線.L曲線的關(guān)鍵就在于選擇曲線上曲率最大的那個(gè)點(diǎn)，其對(duì)應(yīng)的嶺參數(shù)即為所求.

接下來以對(duì)數(shù)形式推導(dǎo)公式[3]，令：

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

關(guān)于 ξ′、ξ′′、?′、?′′的計(jì)算公式可以參考文獻(xiàn)[20]，對(duì)式(29)求取最大值，就可以得到最大曲率 ηmax， ηmax所對(duì)應(yīng)的點(diǎn)的嶺參數(shù)即為所求，這樣就定位出了L曲線上曲率最大的點(diǎn). 由式(26)可以看出，在應(yīng)用L曲線法選擇嶺參數(shù)的合理性在于平衡數(shù)據(jù)擬合度部分和解部分這種平衡是通過嶺參數(shù)來實(shí)現(xiàn)的，需要指出的是，用L曲線法求得的嶺參數(shù)不是最優(yōu)選擇的，只是近似最優(yōu).

2.3U 曲線法

關(guān)于嶺參數(shù)的選取，U曲線法與L曲線法有所異同，相同的部分在于都是尋求曲率最大的點(diǎn)對(duì)應(yīng)的嶺參數(shù)即為所求. 不同部分在于U曲線是根據(jù)U(k) 函數(shù)得到一條U(k)-k曲線，曲線左側(cè)曲率的最大點(diǎn)k為U曲線法所確定的嶺參數(shù).

在通用EIV平差模型的LTLS解法基礎(chǔ)上，轉(zhuǎn)變?yōu)镚H模型后，根據(jù)U曲線的函數(shù)定義為

式中：

f=為矩陣B經(jīng)過奇異值分解后的左側(cè)矩陣； α1≥α2≥···≥αn＞0 為 矩陣B經(jīng)過奇異值分解過后的按照降序排列的n個(gè)奇異值.U曲線和L曲線關(guān)鍵都是定位曲率最大的點(diǎn)，就U曲線而言，其曲率可以表示為

式中，U(k)′和U(k)′′分別為U(k) 的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)，并且有：

式(32)～(33)中，x(k) 和y(k) 的 一 階、二 階 導(dǎo) 數(shù)x(k)′、y(k)′、x(k)′′、y(k)′′的 計(jì) 算 過 程 可 以 參 考 文獻(xiàn)[17].

采用U曲線法確定嶺參數(shù)時(shí)，對(duì)式(31)在區(qū)間(?,*) 內(nèi) 求出最大值 ηmax，定位出U曲線左側(cè)曲率的最大值點(diǎn)，該點(diǎn)所對(duì)應(yīng)的k值就是我們所要求的嶺參數(shù).U曲線不需要曲線擬合，計(jì)算量要小于L曲線的計(jì)算量，也不會(huì)出現(xiàn)曲線發(fā)散而沒有取值的現(xiàn)象，相比較L曲線法，U曲線法更加適合選取嶺參數(shù).

3 算例與分析

3.1 算例1

采用文獻(xiàn)[10]中的模擬數(shù)據(jù)驗(yàn)證本文算法的正確性. 根據(jù)通用EIV模型式(2)設(shè)計(jì)真值矩陣，定義參數(shù)向量X=[5,10]T，設(shè)置A、L0、B矩陣的真值，給A、L0、B分別添加服從正態(tài)分布的誤差得到含有誤差的矩陣如表1所示，假設(shè)A、L0、B中都含有隨機(jī)誤差，隨機(jī)誤差分別服從標(biāo)準(zhǔn)差為0.01、0.02、0.03的正態(tài)分布，3個(gè)矩陣之間相互獨(dú)立，設(shè)置迭代閾值為1 ×10-8.

表1 含有隨機(jī)誤差的模擬數(shù)據(jù)矩陣

分別利用通用EIV平差模型的TLS、線性化、以及兩種嶺估計(jì)算法進(jìn)行解算、估計(jì)值與真值的差值2-范數(shù)參數(shù)估值的方差-協(xié)方差矩陣跡 t r[D(X?)] 計(jì)算結(jié)果、迭代次數(shù)如表2所示；U曲線法方差分量隨迭代次數(shù)變化如圖1所示. 由于LTLS與WTLS估計(jì)值相一致，對(duì)比情況僅選擇LTLS估計(jì)值. 由于文獻(xiàn)[10]中根據(jù)隨機(jī)誤差進(jìn)行過1 000次實(shí)驗(yàn)，得出了LTLS和WTLS方法的求參數(shù)解平均值和協(xié)因數(shù)矩陣都一致，可以直接添加嶺參數(shù)，進(jìn)行計(jì)算即可；通過U曲線法[16-17]確定嶺參數(shù)為α=9.5×10-3，L曲線法[15]確定嶺參數(shù)為α=1.0×10-2. 通用EIV平差模型的LTLS的單位權(quán)方差計(jì)算式為

圖1 方差分量變化圖

表2 參數(shù)解估計(jì)值及其方差估計(jì)值

式中：r為多余觀測(cè)數(shù)，r=n-t；n為 觀測(cè)值個(gè)數(shù)；t為必要觀測(cè)值個(gè)數(shù)；r與經(jīng)典GH模型的多余觀測(cè)數(shù)相同.

3.2 算例2

采用文獻(xiàn)[21]中的直線擬合算例，具體數(shù)據(jù)如表3所示，權(quán)矩陣如參考文獻(xiàn)所示，收斂條件設(shè)置為1×10-8，通用EIV模型中涉及觀測(cè)向量系數(shù)矩陣的誤差，待估參數(shù)為，直線擬合模型為

表3 坐標(biāo)觀測(cè)值及相應(yīng)權(quán)值

TLS、線性化、以及兩種嶺估計(jì)算法解算結(jié)果及其迭代次數(shù)如表4所示；U曲線法方差分量變化如圖2所示. 通過U曲線法[16-17]確定嶺參數(shù)為α=0.075，L曲線法[15]確定嶺參數(shù)為 α =0.080 .

表4 參數(shù)及其方差估計(jì)值、迭代次數(shù)

圖2 方差分量變化圖

3.3 算例分析

從兩個(gè)算例結(jié)果以及文中的公式推導(dǎo)可以得出以下結(jié)論：

1)由現(xiàn)有的嶺估計(jì)方法研究現(xiàn)狀可以發(fā)現(xiàn)，嶺估計(jì)方法可以用來解決不適定問題即病態(tài)問題，使得參數(shù)估計(jì)提升較為明顯，嶺估計(jì)方法具有很好地改善病態(tài)矩陣的效果；當(dāng)觀測(cè)數(shù)據(jù)構(gòu)成的設(shè)計(jì)矩陣并非病態(tài)矩陣時(shí)，嶺估計(jì)方法同樣具有一定的改善作用. 在眾多學(xué)者對(duì)嶺估計(jì)驗(yàn)證的基礎(chǔ)上，文中采取的算例發(fā)現(xiàn)，系數(shù)矩陣的病態(tài)性不明顯，嶺估計(jì)方法同樣適用.區(qū)別在于嶺參數(shù)選取的數(shù)值可能較小，迭代次數(shù)偏多一點(diǎn).

2)由于LTLS算法在參考文獻(xiàn)[10]中就驗(yàn)證過參數(shù)估值精度等同于WTLS算法，因此算例中主要以LTLS算法對(duì)比分析. 從算例1中可以發(fā)現(xiàn)，添加嶺估計(jì)可以提升LTLS算法得精度，L曲線法與U曲線法的計(jì)算結(jié)果相似，但是U曲線法略優(yōu)于L曲線法. LTLS算法、LTLS嶺估計(jì)算法(U曲線法和L曲線法)三種方法的估計(jì)值與真值的差值范數(shù)分別為1.828 853 ×10-4、 1.332 065 × 10-4和 7.687 158 × 10-4；從表2可以看出，三種估計(jì)方法的方差-協(xié)方差矩陣的跡依次減小，分別為1.151 777 ×10-4、8.641 684×10-4和8.077 596 × 10-4. 對(duì)應(yīng)的迭代次數(shù)分別為11、8和8.

3)算例2中，LTLS算法、LTLS嶺估計(jì)算法(U曲線法和L曲線法)三種方法的迭代次數(shù)分別為18、14和13. 由于參數(shù)真值未知，通過比較參數(shù)估值的方差分量變化情況進(jìn)行精度分析. 考慮具有一般性的通用EIV模型的平差計(jì)算結(jié)果，對(duì)其估值進(jìn)行方差分量估計(jì)表達(dá)；由表2、表4可以看出，三種方法計(jì)算得到單個(gè)未知參數(shù)方差數(shù)值接近，說明嶺估計(jì)方法計(jì)算得到參數(shù)估值具有可信度，而在相同迭代閾值約束下，嶺估計(jì)方法迭代計(jì)算量更少，收斂更快. 從圖1～2可以看出，未知參數(shù)在嶺估計(jì)(U曲線法)下估值的方差分量收斂較快.

4 結(jié) 語

關(guān)于TLS模型的解算研究是測(cè)繪數(shù)據(jù)處理研究熱點(diǎn)之一. 通用EIV平差模型作為EIV模型的通用形式，有方差分量估計(jì)解、WTLS以及線性化估計(jì)解法. 為了降低參數(shù)估計(jì)解的方差，進(jìn)而使得均方誤差降低，引入處理病態(tài)問題的嶺估計(jì)方法，主要是利用嶺參數(shù)來平衡正則化矩陣和病態(tài)系數(shù)矩陣之間的關(guān)系，也可以運(yùn)用在TLS中. 本文在通用EIV平差模型的LTLS基礎(chǔ)上，引入嶺估計(jì)原理，推導(dǎo)出通用EIV平差模型的嶺估計(jì)迭代解式，采用L曲線法和U曲線法確定嶺參數(shù)，來修正各奇異值，可有效減小參數(shù)估計(jì)的方差，減少偏差的引入，得到更為可靠的參數(shù)估計(jì). 通過算例驗(yàn)證了在通用EIV平差模型的嶺估計(jì)解可有效減少計(jì)算量，方差分量收斂更快，迭代次數(shù)低于LTLS迭代次數(shù).