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

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

不確定性是一種廣義的誤差，它包含可度量的數(shù)值誤差和無法用數(shù)值度量的誤差。不確定性不再是一個具體數(shù)值，它在一定的實數(shù)區(qū)間內(nèi)變動，或者僅是一個模糊數(shù)。抑制不確定性的影響，現(xiàn)有的平差理論還存在局限性。最近有許多學(xué)者研究了一種新的不確定性假設(shè)“未知但有界(unknown-but-bound，UBB)噪聲”在測量數(shù)據(jù)處理中的應(yīng)用[1-3]。由于UBB噪聲不需要太多的先驗條件，只要求噪聲滿足有界假設(shè)，這一點(diǎn)在實際測量中容易得到保證。如果能夠找到由所有與觀測數(shù)據(jù)、模型結(jié)構(gòu)和噪聲的有界假設(shè)相容的參數(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雜志出了一期專刊介紹集員估計理論與方法的研究成果[14]。Schweppe(1968)是早期研究橢球集員估計算法的學(xué)者。他采用橢球近似描述狀態(tài)可行集[15]。文獻(xiàn)[16]首先提出基于UBB噪聲的參數(shù)的集員估計方法，其集合的Chebyshev中心可作為參數(shù)真實值的一個點(diǎn)估計。文獻(xiàn)[17]又進(jìn)一步改善了其算法，將橢球引入?yún)?shù)可行集的近似描述中來，提出了橢球集員估計算法。最近幾年，橢球集員估計算法得到了迅速的發(fā)展[18-20]。用橢球集合來描述不確定性，實際上就是測量平差中用誤差橢圓來描述點(diǎn)位誤差的擴(kuò)展，目前測量數(shù)據(jù)處理中，已有一些針對于橢圓和區(qū)間集合的簡單算法[21-23]。用一個集合來描述不確定性，然后再用集合的特征(如體積)，來度量不確定性是對誤差概念的一種較好的擴(kuò)展。本文將在橢球集合描述不確定性的基礎(chǔ)上建立一個新的不確定性平差模型。通過定界集合的運(yùn)算，以兩個集合的交集來研究不確定度的傳遞過程?；跈E球集合特征矩陣的跡最小化建立最小不確定度平差準(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]。在幾何上，橢球的扁平程度以及橢球的體積是由橢球的特征矩陣來確定的。由于本文研究的不確定性是一種有界約束(橢球約束)，這個有界性是通過矩陣P來刻畫的，它相當(dāng)于e的協(xié)方差陣來刻畫e的特征一樣。

若L=AX是相容方程組，取X0使得L=AX0。當(dāng)L=AX不相容時，取X0=XLS=(ATA)-1ATL，這時，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帶有橢球約束先驗信息，X的可行空間可以用下面的橢球集合來表示

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

(4)

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

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 帶有橢球不確定性約束的集員估計

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

它們的交定義為

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

X∈E(z2,P2)}

(7)

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

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

(8)

(9)

(10)

(11)

(12)

式中，0

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

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

(13)

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

(14)

3 ρ和a的計算方法

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

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

令

(15)

有

(16)

PU=β(I-KA)Q

(17)

因此

(18)

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

a=argmin tr(PU)

(19)

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

(20)

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

tr[(I-KA)Q]>0

(21)

4 算例分析

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

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)

式中

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

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

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

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

算例3設(shè)有一測邊網(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的真實坐標(biāo)假設(shè)為已知(表1)，邊長的觀測值是利用真實坐標(biāo)計算，再加上誤差得到的，觀測邊長視為同精度(表2)。

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

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

表2 邊長觀測值

為了便于分析，假設(shè)由前期的觀測得到了P3、P4、P5、P6的近似坐標(biāo)(表1)，以及它們相應(yīng)的點(diǎn)位精度。相對于近似坐標(biāo)的改正數(shù)構(gòu)成的未知向量為X=[x3y3x4y4x5y5x6y6]T，可以得到相應(yīng)的平差模型的系數(shù)矩陣A和觀測向量L。由于已知點(diǎn)P2的坐標(biāo)非?？拷黀1點(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)化成嶺估計方法。如文獻(xiàn)[25]，將帶橢球約束的線性模型估計寫成

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)

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

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

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

-1.480 02.368 22.184 4-2.907 2]T

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

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

(3) 令

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

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

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

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

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

5 結(jié)束語

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

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