翁 燁，邵德盛，2

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

病態(tài)問題廣泛存在于自然科學(xué)和社會科學(xué)領(lǐng)域中，對于這些領(lǐng)域中的求解問題，常用Gauss-Markov模型或是離散動態(tài)模型進(jìn)行建模，而最小二乘估計(Least Squares Estimate, LS)和卡爾曼(Kalman)濾波分別是兩種經(jīng)典模型解算方法。測量系統(tǒng)的病態(tài)性通常表現(xiàn)為誤差方程系數(shù)矩陣具有復(fù)共線性[1]，LS估計和Kalman濾波估計都存在估計方差膨脹導(dǎo)致參數(shù)估值不穩(wěn)定的現(xiàn)象。自Gauss(高斯)提出最小二乘法開始，最小二乘估計理論在測量數(shù)據(jù)處理中得到廣泛的應(yīng)用，最小二乘估計是測量平差的基礎(chǔ)方法，顧及觀測向量的隨機(jī)誤差，其平差準(zhǔn)則是觀測向量的殘差平方和達(dá)到最小[2-3]。對于系數(shù)矩陣也存在誤差的函數(shù)模型通常稱為EIV(Errors-In-Variables)模型，則需要采用總體最小二乘估計(Total Least Squares, TLS)，其平差準(zhǔn)則是觀測向量和系數(shù)矩陣兩類殘差的平方和最小，平差理論更為合理[4-6]。對于EIV模型的解算研究，國外學(xué)者Gene提出在EIV模型下總體最小二乘問題的奇異值分解(Singular Value Decomposition，SVD)解法，推導(dǎo)出在系數(shù)矩陣的譜分解下參數(shù)的求和解式，缺陷在于并沒有消除或者削弱設(shè)計矩陣的病態(tài)性，只是對參數(shù)估計進(jìn)行分開按照求和公式進(jìn)行求解[4];Burkhard基于拉格朗日乘數(shù)法導(dǎo)出加權(quán)總體最小二乘的迭代求解公式，迭代計算過程復(fù)雜，計算量偏大，不利于編程實現(xiàn)[7]，同時基于Gauss-Helemert模型，通過對非線性觀測方程中隨機(jī)誤差和待估參數(shù)進(jìn)行線性化，由經(jīng)典最小二乘法導(dǎo)出未知參數(shù)的迭代解式[8]。學(xué)者沈云中引入非線性最小二乘的Newton-Gauss法，導(dǎo)出加權(quán)總體最小二乘的迭代解式的同時還提出一種評定估值精度的數(shù)值算法[9]。

在大地測量和地球物理等領(lǐng)域的測量數(shù)據(jù)處理模型中，參數(shù)估計往往伴隨著先驗信息，此時傳統(tǒng)的參數(shù)估計問題即擴(kuò)展為附有約束條件的參數(shù)估計，利用參數(shù)間一些精確已知的先驗等式約束信息也可以顯著提高解的準(zhǔn)確性和精度，因此在病態(tài)總體最小二乘問題中，附加正確的等式約束也能夠提高其解的準(zhǔn)確性和精度[10-11]。對于附有等式約束的EIV模型，國外學(xué)者Burkhard和Yaron研究了附有線性等式約束和二次型約束的等權(quán)總體最小二乘問題并基于拉格朗日乘數(shù)法導(dǎo)出參數(shù)的迭代計算公式[13-14]。但是，約束信息有可能會是相關(guān)的，約束方程就會出現(xiàn)病態(tài)，使得約束矩陣的微小擾動就會導(dǎo)致最終的反演結(jié)果產(chǎn)生極大變化，從而導(dǎo)致反演結(jié)果極不穩(wěn)定，學(xué)者王樂洋等人討論了在具有病態(tài)性質(zhì)約束矩陣的等式約束反演問題，并利用SVD方法來處理約束矩陣的病態(tài)性，有效解決約束矩陣病態(tài)造成的反演結(jié)果極不穩(wěn)定的問題[15]，同時還提出了在附有病態(tài)約束矩陣的等式約束反演問題的嶺估計解法，利用嶺參數(shù)來改變系數(shù)矩陣的病態(tài)性，提高參數(shù)解的可信度，使得反演效果良好，同時也帶有嶺估計統(tǒng)一改正的不足[16]。馬洋、歐吉坤等人驗證了降秩法與奇異值截斷法的基本等效性，針對附有病態(tài)約束的反演問題，提出了主模型與約束條件聯(lián)合解算的新方法，這種方法可以解決主模型病態(tài)甚至秩虧的問題，但丟失了部分原始成分信息[17]。為了保留先驗信息的完整性，文中在總體最小二乘平差準(zhǔn)則的基礎(chǔ)上，考慮精確線性約束條件，利用修正奇異值法對病態(tài)設(shè)計矩陣的奇異值進(jìn)行分區(qū)別修正，導(dǎo)出在精確約束條件下的病態(tài)總體最小二乘迭代解式，并進(jìn)行參數(shù)估值的精度評定。

1 總體最小二乘的SVD模型

總體最小二乘的EIV模型為[5]：

L+el=(B+EB)X.

(1)

式中：L為觀測向量;B為設(shè)計矩陣且R(B)=m

(2)

其中,

eB=vec(EB)，

L=BX+e.

(3)

根據(jù)總體最小二乘平差準(zhǔn)則，模型(1)的求解式為：

在觀測向量和設(shè)計矩陣元素的權(quán)陣均為單位陣Pl=Im×n,PB=In×m時，求解式變?yōu)?

(4)

當(dāng)N=BTB病態(tài)情況不容樂觀，如在條件數(shù)rand(N)≥1 000時，最小二乘估計的參數(shù)值變得不可靠。為了參數(shù)估計值的穩(wěn)定性，設(shè)法削弱法矩陣的病態(tài)性用以獲得穩(wěn)定、可靠的法矩陣逆陣，采用截斷奇異值法來處理總體最小二乘估計。

將B矩陣進(jìn)行奇異值分解，得到[15]:

B=UΣVT.

(5)

式中:U為n×m階矩陣;V為m×m階矩陣;均為正交矩陣。Σ為n×m階對角陣;對角元素為B的奇異值，按照降序排列，Σ=diag(λ1,λ2,…,λm)。將式(5)代入到LS參數(shù)估計中，并做譜展開有:

(6)

式中:vi和ui分別是V和U的第i列向量;λi是B的第i個奇異值。通過式(6)看出在奇異值較小的時候，最小二乘估計解值被嚴(yán)重放大，截斷奇異值主要在于將小的奇異值截斷，避免了小奇異值對參數(shù)估計造成大的影響，若去掉后(n-k)個奇異值，則總體最小二乘的截斷奇異值解為:

(7)

2 等式約束條件總體最小二乘模型

2.1 等式約束總體最小二乘的修正奇異值解

等式約束下的EIV病態(tài)模型為[4-6]：

L+el=(B+EB)X，

h=HX.

(8)

式中:H為p×m階約束矩陣，h為p×1階約束常數(shù)項，且有rank(H=p

(9)

其中，μ為p×1階拉格朗日因子，令

表示為：

(10)

根據(jù)文獻(xiàn)[5]可知：

(11)

組合式(10)、式(11)和式(8)可表示為：

(12)

奇異值的修正要選取最小奇異值門限，設(shè)k為選取參數(shù)，將式(5)中的U,Σ,V進(jìn)行矩陣分塊，得到:

(13)

在進(jìn)行奇異值分解后，相對較大的奇異值及其特征向量在參數(shù)模型解算中更可靠，相對較小的奇異值及其特征向量在解算過程中不太可靠。截斷奇異值法在于去除參數(shù)模型中的不太可靠部分，即刪掉相對較小的奇異值，從而減少解的方差，同時這會損害解的分辨率。若是奇異值呈現(xiàn)出如圖1所示的階梯狀時，截斷奇異值法比較適用參數(shù)估計，確定可靠成分和不可靠部分差別很大，截掉不可靠部分的奇異值對解的分辨率影響不大。Hansen從理論上證明了當(dāng)奇異值分布呈現(xiàn)階梯狀時，截斷奇異值的方法與嶺估計的方法基本上是等價的。奇異值如圖1所示呈現(xiàn)的階梯狀且均勻分布，不可靠部分和確定成分的界限不好確定，難以確定截斷參數(shù)，建議采用奇異值修正法。

圖1 奇異值分類示意

奇異值的修正方法為:設(shè)t為截斷奇異值保留最小奇異值的門限，T為對應(yīng)的小于t的奇異值個數(shù)，進(jìn)行奇異值的修正[18]：

(14)

修改后的奇異值變成

(15)

由于式(12)中系數(shù)矩陣和特征向量均存在未知估計參數(shù)，因此需要利用迭代法求解，求解過程如下：

2)根據(jù)式(12)得到:

其中，“:=”意思是“定義為”，將式(15)奇異值修正結(jié)果重新代入到模型式(6)中，再將LS估計值作為初始值參與迭代計算，迭代的次數(shù)就是修正奇異值參與的次數(shù)，就是病態(tài)精確約束總體最小二乘的修正奇異值解。

2.2 最小奇異值門限的確定

在病態(tài)方程中的法矩陣N的條件數(shù)大于1 000時，病態(tài)性較為嚴(yán)重，根據(jù)病態(tài)性情況結(jié)合條件數(shù)的定義[18]：

(16)

在cond(N)2>1 000時，這里的范圍數(shù)可根據(jù)實際情況調(diào)整，即奇異值λk(k=1,2,…,n-1)滿足：

選取的奇異值最小門限為：

t=λk.

(17)

2.3 精度評定

在精確約束條件下病態(tài)總體最小二乘的修正奇異值解的單位權(quán)方差為[5]：

(18)

需要說明的是，由于在附有精確約束下的病態(tài)總體最小二乘問題中，未知參數(shù)的估計值采用的是迭代求解，未知參數(shù)與觀測向量無法分離，加上設(shè)計矩陣也是隨機(jī)矩陣，本身也存在誤差，因此在這種情況下的單位權(quán)方差無偏估計計算相對復(fù)雜，借助有偏殘差估計的單位權(quán)方差，即式(18)的單位權(quán)方差是有偏的。根據(jù)協(xié)方差的傳播率得出參數(shù)的協(xié)方差矩陣為：

(19)

分塊矩陣形式為:

(20)

3 算例及分析

3.1 整周模糊度算例

采用GPS動態(tài)定位的解算，取歷元間隔時間為2 s，觀測了5顆衛(wèi)星，用4個歷元解算整周模糊度，誤差方程的系數(shù)矩陣為[19]:

(21)

(22)

表1 不同方法參數(shù)估計結(jié)果比較

表2 參數(shù)估計值的方差-協(xié)方差陣

圖2 U曲線法確定正則化參數(shù)α

3.2 測邊網(wǎng)算例

表3 真實坐標(biāo)與近似坐標(biāo) m

表4 邊長觀測值(觀測向量和設(shè)計矩陣，

圖3 三角形網(wǎng)

圖4 U曲線法確定正則化參數(shù)α

表5 不同方法參數(shù)估計結(jié)果比較

表6 參數(shù)估計值的方差-協(xié)方差陣

3.3 病態(tài)網(wǎng)算例

表7 不同方法參數(shù)估計結(jié)果比較

表8 參數(shù)估計值的方差-協(xié)方差陣

圖5 空間測邊網(wǎng)的平面分布圖

圖6 U曲線法確定正則化參數(shù)

參數(shù)估計的精確度容易受到誤差方程的系數(shù)矩陣病態(tài)性影響，在病態(tài)性較為嚴(yán)重時，LS估計值近乎失真，其與真值的差值范數(shù)也較大；根據(jù)3個算例的計算結(jié)果得知，TLS在病態(tài)性較為嚴(yán)重時估計的結(jié)果也不理想，此時可以考慮總體最小二乘的截斷奇異值法或總體最小二乘的正則化法。CTLS方法比較依賴于等式約束條件的系數(shù)矩陣，同時又受到誤差方程系數(shù)矩陣的影響。修正奇異值法的直觀展現(xiàn)在于進(jìn)行一次奇異值修正時，參數(shù)估計值與真值的差值范數(shù)為:

算例1：

算例2：

算例3：

進(jìn)而使得參與迭代計算的初始值經(jīng)過奇異值修正后有所改善。

4 結(jié)束語

參數(shù)估計的先驗信息能夠提高參數(shù)估計解的準(zhǔn)確性和精度。在等式約束下，提出了病態(tài)總體最小二乘估計的修正奇異值解法，利用總體最小二乘準(zhǔn)則和協(xié)方差傳播率推導(dǎo)出在等式約束下的總體最小二乘迭代式；并借助有偏殘差估計的單位權(quán)方差，通過協(xié)方差的傳播率得出參數(shù)的近似協(xié)方差矩陣。用數(shù)值分析算例和測邊網(wǎng)算例驗證文中方法，并與TSVD、TLS和CTLS等參數(shù)估計方法作比較，證明文中方法在奇異值按照降序排列且成階梯均勻分布時，用于對參數(shù)估計精度提升最有利，文中方法在于保持參數(shù)估計解的分辨率同時顧及解的精度提升。文中只做了在等式約束下總體最小二乘奇異值修正，后續(xù)將進(jìn)一步研究在隨機(jī)約束、不等式約束下的總體最小二乘奇異值修正方法。