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

宋迎春，夏玉國，謝雪梅,3

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

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

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

平差模型為

L=AX+e

(1)

E(e)={e:eTP-1e≤1}

(2)

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

若L=AX是相容方程組，取X0使得L=AX0。當(dāng)L=AX不相容時(shí)，取X0=XLS=(ATA)-1ATL，這時(shí)，L≈AX0，利用式(1)有

eTP-1e=(L-AX)TP-1(L-AX)=

(L-AX)TP-1(L-AX)≈

(X-X0)TATP-1A(X-X0)

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

E(e)={X:(X-X0)TATP-1A(X-X0)≤1}

(3)

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

E(c,Q)={X:(X-c)TQ-1(X-c)≤1}

(4)

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

L=AX+es.t.e∈E(e)，X∈E(c,Q)

(5)

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

X∈E(e)∩E(c,Q)

故式(5)也可以表示為

L=AX+e， s.t.X∈E(e)∩E(c,Q)

(6)

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

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

首先，建立一個(gè)不確定性橢球最小化的準(zhǔn)則來確定式(6)的集員估計(jì)解。設(shè)E(z1,P1)和E(z2,P2)為兩個(gè)橢球，它們分別定義為

它們的交定義為

E(z1,P1)∩E(z2,P2)={X:X∈E(z1,P1),

X∈E(z2,P2)}

(7)

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

E(z,P)?E(z1,P1)∩E(z2,P2)

(8)

(9)

(10)

(11)

(12)

式中，0

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

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

(13)

[aATP-1AX0+(1-a)Q-1c]

(14)

3 ρ和a的計(jì)算方法

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

利用上面PU的計(jì)算，可以得到

令

(15)

有

(16)

PU=β(I-KA)Q

(17)

因此

(18)

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

a=argmin tr(PU)

(19)

許多文獻(xiàn)給出了式(19)中a的計(jì)算方法，但通常較為復(fù)雜，本文方法是直接搜索。因?yàn)椋?

(20)

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

tr[(I-KA)Q]>0

(21)

4 算例分析

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

L1=A1X+e1

(22)

式中，X的真值為[47]T

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

e1∈E(e)={e:eTP-1e≤1}

(23)

X∈E(c,Q)={X:(X-c)TQ-1(X-c)≤1}

(24)

式中

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

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

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

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

算例3設(shè)有一測(cè)邊網(wǎng)，P1、P2為已知點(diǎn)，其坐標(biāo)分別為(48 580.285 m,600 500.496 m)和(48 570.013 m,60 555.845 m)。為了便于分析比較，算例中的點(diǎn)P3、P4、P5、P6的真實(shí)坐標(biāo)假設(shè)為已知(表1)，邊長的觀測(cè)值是利用真實(shí)坐標(biāo)計(jì)算，再加上誤差得到的，觀測(cè)邊長視為同精度(表2)。

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

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

表2 邊長觀測(cè)值

為了便于分析，假設(shè)由前期的觀測(cè)得到了P3、P4、P5、P6的近似坐標(biāo)(表1)，以及它們相應(yīng)的點(diǎn)位精度。相對(duì)于近似坐標(biāo)的改正數(shù)構(gòu)成的未知向量為X=[x3y3x4y4x5y5x6y6]T，可以得到相應(yīng)的平差模型的系數(shù)矩陣A和觀測(cè)向量L。由于已知點(diǎn)P2的坐標(biāo)非常靠近P1點(diǎn)，導(dǎo)致算法中的系數(shù)矩陣病態(tài)。

相應(yīng)的平差模型為

L=AX+e

(25)

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

e∈E(e)={e:eTP-1e≤1}

(26)

X∈E(c,Q)={X:(X-c)TQ-1(X-c)≤1}

(27)

式中

P=0.005I9

Q=0.01I8c=[-1.523 0-2.758 01.902 0-1.530 0

-1.560 03.503 03.334 0-3.665 0]T

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

min (L-AX)TP-1(L-AX)

s.t. (X-c)TQ-1(X-c)≤1

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

f(X,λ)=(L-AX)TP-1(L-AX)+

λ((X-c)TQ-1(X-c)-1)

求得橢球約束下的廣義嶺估計(jì)

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

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

-1.480 02.368 22.184 4-2.907 2]T

綜上，對(duì)實(shí)例解算有如下說明與分析：

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

(3) 令

f(X)=(L-AX)TP-1(L-AX)

M(X)=(X-Xtrue)TP-1(X-Xtrue)

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

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

(6) 對(duì)于病態(tài)模型的其他算法，如表3中的截?cái)嗥娈愔邓惴ê蛶X估計(jì)算法，它們是利用數(shù)學(xué)原理來處理病態(tài)系數(shù)矩陣，不能有效地利用先驗(yàn)信息，計(jì)算的結(jié)果不如本文的算法。更重要的是，本文算法不僅能給出參數(shù)估計(jì)的值，而且還能對(duì)參數(shù)估計(jì)的不確定度進(jìn)行估計(jì)。

5 結(jié)束語

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

表3 法矩陣病態(tài)時(shí)算法的結(jié)果比較