基于橢球不確定性的平差模型與算法

宋迎春^1,2, 夏玉國^1,2, 謝雪梅^1,2,3     

1. 有色金屬成礦預測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室(中南大學), 湖南 長沙 410083; 
2. 中南大學地球科學與信息物理學院, 湖南 長沙 410083; 
3. 中南林業(yè)科技大學土木工程學院, 湖南 長沙 410004

收稿日期：2017-10-31；修回日期：2018-07-15

基金項目：國家自然科學基金（41574006；41674009；41674012）

第一作者簡介：宋迎春(1963-), 男, 博士, 教授, 研究方向為測量數(shù)據(jù)處理理論與方法。E-mail:csusyc@csu.edu.cn csusyc@qq.com

通信作者：謝雪梅 E-mail：xiexuemei_2003@126.com

摘要：測量平差模型中的參數(shù)通常存在一些不確定的附加信息或先驗信息，充分利用它們可以對部分參數(shù)進行約束，從而保證參數(shù)解的唯一性和穩(wěn)定性。本文利用橢球集合描述不確定性，建立了一個新的帶有橢球不確定性的平差模型。以兩個橢球交集的外接橢球的特征矩陣的跡最小平差準則，分析了不確定度的傳播規(guī)律，給出了帶有橢球不確定性的平差方法。最后，通過算例驗證了算法的有效性，說明了平差解與帶權(quán)混合估計的關系。

關鍵詞：不確定性    橢球約束    平差模型    病態(tài)問題    集員估計    

Adjustment model and algorithm based on ellipsoid uncertainty

SONG Yingchun^1,2, XIA Yuguo^1,2, XIE Xuemei^1,2,3     

1. Key Laboratory of Metallogenic Prediction of Nonferrous Metals and Geological Environment Monitoring(Central South University), Ministry of Education, Changsha 410083, China;
2. School of Geosciences and Info-Physics, Central South University, Changsha 410083, China;
3. School of Civil Engineering, Central South University of Forestry and Technology, Changsha 410004, China

Foundation support: The National Natural Science Foundation of China (Nos. 41574006; 41674009; 41674012)

First author: SONG Yingchun(1963—), male, PhD, professor, majors in theory and method of measuring data processing.E-mail:csusyc@csu.edu.cn csusyc@qq.com.

Corresponding author: XIE Xuemei E-mail: xiexuemei_2003@126.com.

Abstract: In surveying adjustment models, there usually is some uncertain additional information or prior information on parameters, which can constraint on the parameters, and guarantee uniqueness and stability of parameters solution.In this paper, ellipsoidal sets are used to describe uncertainty, so an adjustment model with ellipsoidal uncertainty is established. The minimization in matrix trace of circumscribed ellipsoid with two ellipsoid intersections is regarded as a proposed adjustment criterion, the propagation law of uncertainty is analyzed, and the adjustment method with ellipsoid uncertainty is given. Finally, a numerical example is given to test and verify the effectiveness of the proposed algorithm, and the relation between the adjustment result and the weighted mixed estimation is illustrated.

Key words: uncertainty    ellipsoid constraint    adjustment model    ill-posed problem    set membership estimation    

不確定性是一種廣義的誤差，它包含可度量的數(shù)值誤差和無法用數(shù)值度量的誤差。不確定性不再是一個具體數(shù)值，它在一定的實數(shù)區(qū)間內(nèi)變動，或者僅是一個模糊數(shù)。抑制不確定性的影響，現(xiàn)有的平差理論還存在局限性。最近有許多學者研究了一種新的不確定性假設“未知但有界(unknown-but-bound，UBB)噪聲”在測量數(shù)據(jù)處理中的應用^[1-3]。由于UBB噪聲不需要太多的先驗條件，只要求噪聲滿足有界假設，這一點在實際測量中容易得到保證。如果能夠找到由所有與觀測數(shù)據(jù)、模型結(jié)構(gòu)和噪聲的有界假設相容的參數(shù)組成的集合^[4-7]，那么，此集合中的任何元素都可以成為參數(shù)解，此集合一般被稱為參數(shù)的可行集(feasible solution set)。在一定的條件下，隨著樣本容量增大，成員集所包含的范圍逐漸縮小，最后成員集最終收斂于系統(tǒng)的真實參數(shù)^[8]。這種基于有界不確定性噪聲(UBB噪聲)的參數(shù)估計方法稱為集員估計(set membership estimation)方法^[4, 9-13]。2005年Mathematical and Computer Modeling of Dynamical Systems雜志出了一期?？榻B集員估計理論與方法的研究成果^[14]。Schweppe(1968)是早期研究橢球集員估計算法的學者。他采用橢球近似描述狀態(tài)可行集^[15]。文獻[16]首先提出基于UBB噪聲的參數(shù)的集員估計方法，其集合的Chebyshev中心可作為參數(shù)真實值的一個點估計。文獻[17]又進一步改善了其算法，將橢球引入?yún)?shù)可行集的近似描述中來，提出了橢球集員估計算法。最近幾年，橢球集員估計算法得到了迅速的發(fā)展^[18-20]。用橢球集合來描述不確定性，實際上就是測量平差中用誤差橢圓來描述點位誤差的擴展，目前測量數(shù)據(jù)處理中，已有一些針對于橢圓和區(qū)間集合的簡單算法^[21-23]。用一個集合來描述不確定性，然后再用集合的特征(如體積)，來度量不確定性是對誤差概念的一種較好的擴展。本文將在橢球集合描述不確定性的基礎上建立一個新的不確定性平差模型。通過定界集合的運算，以兩個集合的交集來研究不確定度的傳遞過程?；跈E球集合特征矩陣的跡最小化建立最小不確定度平差準則，并尋找在此準則下的最優(yōu)解。

1 有界橢球不確定性平差模型

平差模型為

 (1)

式中，A為m×n(m≥n)維設計矩陣，且列滿秩；L為n維觀測向量；X=[x₁ x₂ … x_n]^T為n維參數(shù)向量；e為m維有界觀測噪聲，其有界性用范數(shù)形式來描述，如||e||≤γ，也可用區(qū)間來描述，如e∈[e_l, e_u]，這些描述方法，都可以統(tǒng)一用一個橢球集合來表示

 (2)

式中，P為n階正定矩陣。它是橢球的特征矩陣，用來刻畫橢球的形狀特征，類似于e的協(xié)方差陣描述e的特征。有許多學者研究了矩陣P的構(gòu)造^[24-25]。在幾何上，橢球的扁平程度以及橢球的體積是由橢球的特征矩陣來確定的。由于本文研究的不確定性是一種有界約束(橢球約束)，這個有界性是通過矩陣P來刻畫的，它相當于e的協(xié)方差陣來刻畫e的特征一樣。

若L=AX是相容方程組，取X₀使得L=AX₀。當L=AX不相容時，取X₀=X_LS=(A^TA)^-1A^TL，這時，L≈AX₀，利用式(1)有

e的有界不確定性也可以近似地表示為

 (3)

若X帶有橢球約束先驗信息，X的可行空間可以用下面的橢球集合來表示

 (4)

式中，c是橢球的中心；Q是橢球的特征矩陣，用來刻畫橢球的形狀特征。對于下面的平差模型

 (5)

稱式(5)為帶有橢球不確定性的平差模型。利用式(3)，式(5)的約束條件可以寫成

故式(5)也可以表示為

 (6)

E=E(e)∩E(c, Q)是參數(shù)向量的可行解集。

2 帶有橢球不確定性約束的集員估計

首先，建立一個不確定性橢球最小化的準則來確定式(6)的集員估計解。設E(z₁, P₁)和E(z₂, P₂)為兩個橢球，它們分別定義為

它們的交定義為

 (7)

顯然，兩個橢球的交集不一定是一個橢球，它可以用一個外包橢球來近似^[19]，如圖 1所示。設其外包橢球族為E(z, P)，有

 (8)

 (9)

 (10)

 (11)

 (12)

圖 1 兩個橢球交的最小外包橢球Fig. 1 The minimum circumscribed ellipsoid with two ellipsoid intersections

式中，0 < a < 1。上面的交的運算式(8)—式(12)的推導可以參考文獻[26]。通過優(yōu)化系數(shù)a，可以得到最小跡或最小體積外包橢球^[19, 26]。如圖 1，現(xiàn)在假設式(6)的可行解集E=E(e)∩E(c, Q)的外包橢球為

由式(9)與式(10)，有

 (13)

 (14)

由文獻[16]可知，是式(5)的一個解。從式(14)可以看出，是嶺估計算法的一種推廣，形式上就是文獻[27]提出的加權(quán)混合估計。如果能確定a，不僅可以得到帶有橢球不確定性的參數(shù)估計方法，還可以從不確定性的角度，得到一個加權(quán)混合估計的權(quán)確定方法。

3 ρ和a的計算方法

利用(F-CG^-1D)^-1=F^-1+F^-1C(G-DF^-1C)^-1DF^-1，式(13)和式(14)可以化為

利用上面P_U的計算，可以得到

令

 (15)

有

 (16)

 (17)

式中，。由式(12)，顧及式(15)，有

因此

 (18)

不同的a可以得到不同的測量更新橢球。為了保證橢球交的外包橢球的最小性，可以通過優(yōu)化系數(shù)a，得到最小跡外包橢球。

 (19)

許多文獻給出了式(19)中a的計算方法，但通常較為復雜，本文方法是直接搜索。因為，0 < a < 1，0≤ρ < 1，可知a必須使得

 (20)

由式(13)可知，P_U是一個正定矩陣，因此a還必須滿足

 (21)

a在滿足式(20)和式(21)的條件下，直接利用搜索算法(a從0開始到1止，通過增量Δa，逐步搜索得到使tr(P_U)達到最小的a)求出a的值，從而求出參數(shù)估計值和P_U。

4 算例分析

算例1  為了便于畫出橢圓進行分析說明，特設計如下的平差模型

 (22)

式中，X的真值為[4 7]^T

誤差向量e₁的不確定性和參數(shù)向量X的橢圓約束信息分別定義如下

 (23)

 (24)

式中

計算采用Matlab的隨機函數(shù)生成的滿足式(23)的隨機數(shù)，e₁=[0.033 90.012 9]^T，利用式(19)算得：a=0.059 2，利用式(16)和式(17)計算得到

圖 2給出了3個相應的不確定性橢圓，其中，E(e)通過(3)變換以后的橢圓{X:(X-X₀)^TA^TP^-1A(X-X₀)^T≤1}，其中X₀=A₁^-1L₁=[3.943 87.041 7]^T，E(c, Q)為約束橢圓。這兩個橢圓的交集中的每一個點都可以作為參數(shù)估計的解，交集的最小跡外包橢球為可以作為解的不確定度，通?？梢詫⒖闯墒菐в袡E球不確定性平差模型的解，它的不確定度由確定。

圖 2 算例1中的誤差橢圓，X的約束橢圓及解的不確定性橢圓Fig. 2 Error ellipse, constrained ellipse of X and the uncertainty ellipse of solution in example 1

算例2 在算例1中，為了驗算病態(tài)模型下算法的效率，取

此算例中，A₂的病態(tài)性相較于算例1有了適當?shù)脑黾?為了圖形的效果，沒有進行更大的增加，對于較嚴重的病態(tài)情形，參見算例3)。計算中采用Matlab的隨機函數(shù)生成的滿足式(23)的隨機數(shù)，e₂=[0.021 60.039 3]^T，利用式(19)算得：a=0.057 0，利用式(16)和式(17)計算得到

圖 3給出了3個相應的不確定性橢圓，其中，E(e)通過式(3)變換以后的橢圓{X:(X-X₀)^TA^TP^-1A(X-X₀)^T≤1}，這里X₀=A₂^-1L₂=[3.823 07.131 1]^T，E(c, Q)為約束橢圓。由于橢圓變換中牽涉到了逆矩陣運算，使得變換后的矩陣形狀有了較大的變化。交集的最小跡外包橢球為可以作為解的不確定度，可以看成是帶有橢球不確定性平差模型的解，它的不確定度由確定。

圖 3 算例2中的誤差橢圓，X的約束橢圓及解的不確定性橢圓Fig. 3 Error ellipse, constrained ellipse of X and the uncertainty ellipse of solution in example 2

算例3 設有一測邊網(wǎng)，P₁、P₂為已知點，其坐標分別為(48 580.285 m, 600 500.496 m)和(48 570.013 m, 60 555.845 m)。為了便于分析比較，算例中的點P₃、P₄、P₅、P₆的真實坐標假設為已知(表 1)，邊長的觀測值是利用真實坐標計算，再加上誤差得到的，觀測邊長視為同精度(表 2)。

表 1 真實坐標與近似坐標Tab. 1 The true coordinates and approximate coordinates

表 2 邊長觀測值Tab. 2 The observations on edge length

為了便于分析，假設由前期的觀測得到了P₃、P₄、P₅、P₆的近似坐標(表 1)，以及它們相應的點位精度。相對于近似坐標的改正數(shù)構(gòu)成的未知向量為X=[x₃ y₃ x₄ y₄ x₅ y₅ x₆ y₆]^T，可以得到相應的平差模型的系數(shù)矩陣A和觀測向量L。由于已知點P₂的坐標非?？拷?em>P₁點，導致算法中的系數(shù)矩陣病態(tài)。

相應的平差模型為

 (25)

誤差向量e的不確定性和參數(shù)向量X的橢圓約束信息分別定義如下

 (26)

 (27)

式中

由于在平差算法中沒有直接解算帶有橢球約束的平差方法，在數(shù)學上通常是拉格朗日函數(shù)求極值的方法，轉(zhuǎn)化成嶺估計方法。如文獻[25]，將帶橢球約束的線性模型估計寫成

其拉格朗日函數(shù)為

求得橢球約束下的廣義嶺估計

文獻[25]給出了其嶺參數(shù)的計算方法，但同時說明“當設計陣病態(tài)時，使用這種類似于兩步估計的做法應該格外小心，理論上已經(jīng)表明此時不宜采用廣義最小二乘估計”。因此，在本算例中，采用通常的嶺參數(shù)計算方法求出嶺估計再與本文的方法進行比較。

利用式(19)求得a=0.052 8，利用式(16)和式(17)可以得到

綜上，對實例解算有如下說明與分析：

(1) 可看作是帶有橢球約束不確定信息平差模型式(25)的解，此解包含的不確定度可以用橢球

來進行計算。也可以看作是對解的精度評估。

(2) 從加權(quán)混合估計的角度來看，因為計算得到a=0.052 8，說明參數(shù)約束先驗信息式(27)在參數(shù)估計中的作用更大。這也正好說明當模型出現(xiàn)病態(tài)時，利用參數(shù)先驗信息可以改善其病態(tài)性。

(3) 令

式中，X_true是參數(shù)的真值。從表 3可以看出，最小二乘算法的目標函數(shù)值是最小的。但因為系數(shù)矩陣病態(tài)，使得得到的最小二乘解嚴重偏離了真值，因此M(X)值是最大的。本文給出的算法中，是利用兩種不確定信息進行了加權(quán)融合，因此，結(jié)果接近真值，M(X)=0.129 7。

表 3 法矩陣病態(tài)時算法的結(jié)果比較Tab. 3 Algorithm comparison when normal equation matrix is ill-posed

	真值	最小二乘方法	截斷奇異值法	嶺估計法(L曲線法)	本文算法
	-0.523 0	-1.347 3	-0.535 0	-0.552 7	-0.510 7
	-2.758 0	6.162 5	-2.372 9	-2.263 9	-2.677 4
	0.902 0	10.443 9	1.358 5	1.382 6	1.034 0
	-0.536 0	-0.364 5	-0.530 7	-0.497 5	-0.540 6
	-1.560 0	2.869 2	-1.331 3	-1.328 4	-1.480 0
	2.503 0	-5.908 4	2.067 1	2.080 8	2.368 2
	2.334 0	-5.331 9	1.907 1	1.935 5	2.184 4
	-2.665 0	-16.214 1	-3.387 3	-3.328 9	-2.907 2
	0.004 3	1.686 1×10-5	0.001 2	0.004 4	0.306 8
m	0	504.044 1	1.303 2	1.308 9	0.129 7

(4) 本文算法中不確定度的最小化是通過求橢球最小特征矩陣的跡來實現(xiàn)的，也可以通過最小化的橢球體積(對應的是特征矩陣的行列式最小)，相關的算法可參看文獻[8]。

(5) 算法中，a=0.052 8是一個近擬值。a從0開始，通過增量Δa=0.000 1，逐步搜索得到使tr(P_U)達到最小的a。

(6) 對于病態(tài)模型的其他算法，如表 3中的截斷奇異值算法和嶺估計算法，它們是利用數(shù)學原理來處理病態(tài)系數(shù)矩陣，不能有效地利用先驗信息，計算的結(jié)果不如本文的算法。更重要的是，本文算法不僅能給出參數(shù)估計的值，而且還能對參數(shù)估計的不確定度進行估計。

5 結(jié)束語

不論是觀測過程，還是未知參數(shù)本身，都存在不確定性噪聲的干擾。然而，不確定性噪聲非常復雜，難以確切了解諸如噪聲的分布或均值和方差等統(tǒng)計特性。本文在橢球集合描述不確定性的基礎上建立一個新的不確定性平差模型，以兩個集合的交集來研究不確定度的傳遞?；跈E球集合特征矩陣的跡最小化建立最小不確定度平差準則，得到了在最優(yōu)準則下的最優(yōu)解。它與文獻[27]提出的加權(quán)混合估計在形式上是一致的，都是對先驗信息的利用。a的確定方法可以看作加權(quán)混合估計的權(quán)的一種新的確定方法。不確定性因素會以不確定度的形式反映在測繪數(shù)據(jù)中，其統(tǒng)計特性難以準確獲得，需要在平差解算的同時，盡量使不確定度達到最小。不同于誤差用數(shù)值來描述和度量不確定性，本文嘗試用一個橢球來描述不確定性，然后用橢球特征矩陣的跡來度量不確定性的大小，這是對于不確定性描述與度量的一種嘗試。

【引文格式】宋迎春, 夏玉國, 謝雪梅. 基于橢球不確定性的平差模型與算法. 測繪學報，2019，48(5)：555-562. DOI: 10.11947/j.AGCS.2019.20170611