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

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

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

0 引 言

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

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

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

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

1 通用EIV模型的線性化估計

GH模型為

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

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

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

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

依據LS最優化準則，將式(2)定義為

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

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

2 通用EIV平差模型的嶺估計解

2.1 嶺估計法

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

依據LS原理有

通過建立相應的目標函數

式中， μ 為f×1 階Lagrange向量，將上式分別對其函數元素求一階偏導數，并令其為0有：

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

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

將式(16)作為改正數方程，于是基礎方程變為

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

得到參數估計的協因數矩陣和方差-協方差矩陣為:

μ?作為關聯系數，可以不必計算，直接將式(20)帶入到式(16)，得到

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

2.2 L 曲線法

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

與LS相比較，嶺估計中多添加了kI一項，由于這一項的參與，法方程的病態性得到了抑制，(N+kI)求逆變得正常起來，所以相比較LS能得到更加可靠的估計值. 嶺估計的結果優于LS的部分就在于嶺參數k，若選擇不同的嶺參數，得到的估計結果就會有所不同，所以嶺參數的選取也至關重要. 在式(26)中都是嶺參數k的函數，選擇不同的嶺參數，以為 橫坐標，為縱坐標畫圖，得到許多點再利用多點曲線擬合的方式擬合出一條曲線，因為該曲線通常形狀都像“L”，所以稱為L曲線.L曲線的關鍵就在于選擇曲線上曲率最大的那個點，其對應的嶺參數即為所求.

接下來以對數形式推導公式[3]，令：

兩邊取其對數形式，得出：

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

2.3U 曲線法

關于嶺參數的選取，U曲線法與L曲線法有所異同，相同的部分在于都是尋求曲率最大的點對應的嶺參數即為所求. 不同部分在于U曲線是根據U(k) 函數得到一條U(k)-k曲線，曲線左側曲率的最大點k為U曲線法所確定的嶺參數.

在通用EIV平差模型的LTLS解法基礎上，轉變為GH模型后，根據U曲線的函數定義為

式中：

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

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

式(32)～(33)中，x(k) 和y(k) 的 一 階、二 階 導 數x(k)′、y(k)′、x(k)′′、y(k)′′的 計 算 過 程 可 以 參 考 文獻[17].

采用U曲線法確定嶺參數時，對式(31)在區間(?,*) 內 求出最大值 ηmax，定位出U曲線左側曲率的最大值點，該點所對應的k值就是我們所要求的嶺參數.U曲線不需要曲線擬合，計算量要小于L曲線的計算量，也不會出現曲線發散而沒有取值的現象，相比較L曲線法，U曲線法更加適合選取嶺參數.

3 算例與分析

3.1 算例1

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

表1 含有隨機誤差的模擬數據矩陣

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

圖1 方差分量變化圖

表2 參數解估計值及其方差估計值

式中：r為多余觀測數，r=n-t；n為 觀測值個數；t為必要觀測值個數；r與經典GH模型的多余觀測數相同.

3.2 算例2

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

表3 坐標觀測值及相應權值

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

表4 參數及其方差估計值、迭代次數

圖2 方差分量變化圖

3.3 算例分析

從兩個算例結果以及文中的公式推導可以得出以下結論：

1)由現有的嶺估計方法研究現狀可以發現，嶺估計方法可以用來解決不適定問題即病態問題，使得參數估計提升較為明顯，嶺估計方法具有很好地改善病態矩陣的效果；當觀測數據構成的設計矩陣并非病態矩陣時，嶺估計方法同樣具有一定的改善作用. 在眾多學者對嶺估計驗證的基礎上，文中采取的算例發現，系數矩陣的病態性不明顯，嶺估計方法同樣適用.區別在于嶺參數選取的數值可能較小，迭代次數偏多一點.

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

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

4 結 語

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