翁 燁，邵德盛，2

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

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

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

1 總體最小二乘的SVD模型

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

L+el=(B+EB)X.

(1)

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

(2)

其中,

eB=vec(EB)，

L=BX+e.

(3)

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

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

(4)

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

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

B=UΣVT.

(5)

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

(6)

式中:vi和ui分別是V和U的第i列向量;λi是B的第i個(gè)奇異值。通過(guò)式(6)看出在奇異值較小的時(shí)候，最小二乘估計(jì)解值被嚴(yán)重放大，截?cái)嗥娈愔抵饕谟趯⑿〉钠娈愔到財(cái)啵苊饬诵∑娈愔祵?duì)參數(shù)估計(jì)造成大的影響，若去掉后(n-k)個(gè)奇異值，則總體最小二乘的截?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)

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

(13)

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

圖1 奇異值分類示意

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

(14)

修改后的奇異值變成

(15)

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

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

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

2.2 最小奇異值門(mén)限的確定

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

(16)

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

選取的奇異值最小門(mén)限為：

t=λk.

(17)

2.3 精度評(píng)定

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

(18)

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

(19)

分塊矩陣形式為:

(20)

3 算例及分析

3.1 整周模糊度算例

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

(21)

(22)

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

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

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

3.2 測(cè)邊網(wǎng)算例

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

表4 邊長(zhǎng)觀測(cè)值(觀測(cè)向量和設(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 空間測(cè)邊網(wǎng)的平面分布圖

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

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

算例1：

算例2：

算例3：

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

4 結(jié)束語(yǔ)

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