等式約束病態(tài)總體最小二乘的修正奇異值解及精度評定

翁 燁，邵德盛，2

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

病態(tài)問題廣泛存在于自然科學(xué)和社會科學(xué)領(lǐng)域中，對于這些領(lǐng)域中的求解問題，常用Gauss-Markov模型或是離散動態(tài)模型進(jìn)行建模，而最小二乘估計(jì)(Least Squares Estimate, LS)和卡爾曼(Kalman)濾波分別是兩種經(jīng)典模型解算方法。測量系統(tǒng)的病態(tài)性通常表現(xiàn)為誤差方程系數(shù)矩陣具有復(fù)共線性[1]，LS估計(jì)和Kalman濾波估計(jì)都存在估計(jì)方差膨脹導(dǎo)致參數(shù)估值不穩(wěn)定的現(xiàn)象。自Gauss(高斯)提出最小二乘法開始，最小二乘估計(jì)理論在測量數(shù)據(jù)處理中得到廣泛的應(yīng)用，最小二乘估計(jì)是測量平差的基礎(chǔ)方法，顧及觀測向量的隨機(jī)誤差，其平差準(zhǔn)則是觀測向量的殘差平方和達(dá)到最小[2-3]。對于系數(shù)矩陣也存在誤差的函數(shù)模型通常稱為EIV(Errors-In-Variables)模型，則需要采用總體最小二乘估計(jì)(Total Least Squares, TLS)，其平差準(zhǔn)則是觀測向量和系數(shù)矩陣兩類殘差的平方和最小，平差理論更為合理[4-6]。對于EIV模型的解算研究，國外學(xué)者Gene提出在EIV模型下總體最小二乘問題的奇異值分解(Singular Value Decomposition，SVD)解法，推導(dǎo)出在系數(shù)矩陣的譜分解下參數(shù)的求和解式，缺陷在于并沒有消除或者削弱設(shè)計(jì)矩陣的病態(tài)性，只是對參數(shù)估計(jì)進(jìn)行分開按照求和公式進(jìn)行求解[4];Burkhard基于拉格朗日乘數(shù)法導(dǎo)出加權(quán)總體最小二乘的迭代求解公式，迭代計(jì)算過程復(fù)雜，計(jì)算量偏大，不利于編程實(shí)現(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ù)估計(jì)往往伴隨著先驗(yàn)信息，此時傳統(tǒng)的參數(shù)估計(jì)問題即擴(kuò)展為附有約束條件的參數(shù)估計(jì)，利用參數(shù)間一些精確已知的先驗(yàn)等式約束信息也可以顯著提高解的準(zhǔn)確性和精度，因此在病態(tài)總體最小二乘問題中，附加正確的等式約束也能夠提高其解的準(zhǔn)確性和精度[10-11]。對于附有等式約束的EIV模型，國外學(xué)者Burkhard和Yaron研究了附有線性等式約束和二次型約束的等權(quán)總體最小二乘問題并基于拉格朗日乘數(shù)法導(dǎo)出參數(shù)的迭代計(jì)算公式[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)約束矩陣的等式約束反演問題的嶺估計(jì)解法，利用嶺參數(shù)來改變系數(shù)矩陣的病態(tài)性，提高參數(shù)解的可信度，使得反演效果良好，同時也帶有嶺估計(jì)統(tǒng)一改正的不足[16]。馬洋、歐吉坤等人驗(yàn)證了降秩法與奇異值截?cái)喾ǖ幕镜刃裕槍Ω接胁B(tài)約束的反演問題，提出了主模型與約束條件聯(lián)合解算的新方法，這種方法可以解決主模型病態(tài)甚至秩虧的問題，但丟失了部分原始成分信息[17]。為了保留先驗(yàn)信息的完整性，文中在總體最小二乘平差準(zhǔn)則的基礎(chǔ)上，考慮精確線性約束條件，利用修正奇異值法對病態(tài)設(shè)計(jì)矩陣的奇異值進(jìn)行分區(qū)別修正，導(dǎo)出在精確約束條件下的病態(tài)總體最小二乘迭代解式，并進(jìn)行參數(shù)估值的精度評定。

1 總體最小二乘的SVD模型

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

L+el=(B+EB)X.

(1)

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

(2)

其中,

eB=vec(EB)，

L=BX+e.

(3)

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

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

(4)

當(dāng)N=BTB病態(tài)情況不容樂觀，如在條件數(shù)rand(N)≥1 000時，最小二乘估計(jì)的參數(shù)值變得不可靠。為了參數(shù)估計(jì)值的穩(wěn)定性，設(shè)法削弱法矩陣的病態(tài)性用以獲得穩(wěn)定、可靠的法矩陣逆陣，采用截?cái)嗥娈愔捣▉硖幚砜傮w最小二乘估計(jì)。

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

B=UΣVT.

(5)

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

(6)

式中:vi和ui分別是V和U的第i列向量;λi是B的第i個奇異值。通過式(6)看出在奇異值較小的時候，最小二乘估計(jì)解值被嚴(yán)重放大，截?cái)嗥娈愔抵饕谟趯⑿〉钠娈愔到財(cái)?，避免了小奇異值對參?shù)估計(jì)造成大的影響，若去掉后(n-k)個奇異值，則總體最小二乘的截?cái)嗥娈愔到鉃?

(7)

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

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

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

L+el=(B+EB)X，

h=HX.

(8)

式中:H為p×m階約束矩陣，h為p×1階約束常數(shù)項(xiàng)，且有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ù)模型解算中更可靠，相對較小的奇異值及其特征向量在解算過程中不太可靠。截?cái)嗥娈愔捣ㄔ谟谌コ齾?shù)模型中的不太可靠部分，即刪掉相對較小的奇異值，從而減少解的方差，同時這會損害解的分辨率。若是奇異值呈現(xiàn)出如圖1所示的階梯狀時，截?cái)嗥娈愔捣ū容^適用參數(shù)估計(jì)，確定可靠成分和不可靠部分差別很大，截掉不可靠部分的奇異值對解的分辨率影響不大。Hansen從理論上證明了當(dāng)奇異值分布呈現(xiàn)階梯狀時，截?cái)嗥娈愔档姆椒ㄅc嶺估計(jì)的方法基本上是等價(jià)的。奇異值如圖1所示呈現(xiàn)的階梯狀且均勻分布，不可靠部分和確定成分的界限不好確定，難以確定截?cái)鄥?shù)，建議采用奇異值修正法。

圖1 奇異值分類示意

奇異值的修正方法為:設(shè)t為截?cái)嗥娈愔当Ａ糇钚∑娈愔档拈T限，T為對應(yīng)的小于t的奇異值個數(shù)，進(jìn)行奇異值的修正[18]：

(14)

修改后的奇異值變成

(15)

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

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

其中，“:=”意思是“定義為”，將式(15)奇異值修正結(jié)果重新代入到模型式(6)中，再將LS估計(jì)值作為初始值參與迭代計(jì)算，迭代的次數(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ù)實(shí)際情況調(diào)整，即奇異值λk(k=1,2,…,n-1)滿足：

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

t=λk.

(17)

2.3 精度評定

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

(18)

需要說明的是，由于在附有精確約束下的病態(tài)總體最小二乘問題中，未知參數(shù)的估計(jì)值采用的是迭代求解，未知參數(shù)與觀測向量無法分離，加上設(shè)計(jì)矩陣也是隨機(jī)矩陣，本身也存在誤差，因此在這種情況下的單位權(quán)方差無偏估計(jì)計(jì)算相對復(fù)雜，借助有偏殘差估計(jì)的單位權(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ù)估計(jì)結(jié)果比較

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

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

3.2 測邊網(wǎng)算例

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

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

圖3 三角形網(wǎng)

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

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

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

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

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

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

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

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

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

算例1：

算例2：

算例3：

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

4 結(jié)束語

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