基于椭球不确定性的平差模型与算法

宋迎春，夏玉国，谢雪梅,3

1. 有色金属成矿预测与地质环境监测教育部重点实验室(中南大学)，湖南 长沙 410083； 2. 中南大学地球科学与信息物理学院，湖南 长沙 410083； 3. 中南林业科技大学土木工程学院，湖南 长沙 410004

不确定性是一种广义的误差，它包含可度量的数值误差和无法用数值度量的误差。不确定性不再是一个具体数值，它在一定的实数区间内变动，或者仅是一个模糊数。抑制不确定性的影响，现有的平差理论还存在局限性。最近有许多学者研究了一种新的不确定性假设“未知但有界(unknown-but-bound，UBB)噪声”在测量数据处理中的应用[1-3]。由于UBB噪声不需要太多的先验条件，只要求噪声满足有界假设，这一点在实际测量中容易得到保证。如果能够找到由所有与观测数据、模型结构和噪声的有界假设相容的参数组成的集合[4-7]，那么，此集合中的任何元素都可以成为参数解，此集合一般被称为参数的可行集(feasible solution set)。在一定的条件下，随着样本容量增大，成员集所包含的范围逐渐缩小，最后成员集最终收敛于系统的真实参数[8]。这种基于有界不确定性噪声(UBB噪声)的参数估计方法称为集员估计(set membership estimation)方法[4,9-13]。2005年Mathematical and Computer Modeling of Dynamical Systems杂志出了一期专刊介绍集员估计理论与方法的研究成果[14]。Schweppe(1968)是早期研究椭球集员估计算法的学者。他采用椭球近似描述状态可行集[15]。文献[16]首先提出基于UBB噪声的参数的集员估计方法，其集合的Chebyshev中心可作为参数真实值的一个点估计。文献[17]又进一步改善了其算法，将椭球引入参数可行集的近似描述中来，提出了椭球集员估计算法。最近几年，椭球集员估计算法得到了迅速的发展[18-20]。用椭球集合来描述不确定性，实际上就是测量平差中用误差椭圆来描述点位误差的扩展，目前测量数据处理中，已有一些针对于椭圆和区间集合的简单算法[21-23]。用一个集合来描述不确定性，然后再用集合的特征(如体积)，来度量不确定性是对误差概念的一种较好的扩展。本文将在椭球集合描述不确定性的基础上建立一个新的不确定性平差模型。通过定界集合的运算，以两个集合的交集来研究不确定度的传递过程。基于椭球集合特征矩阵的迹最小化建立最小不确定度平差准则，并寻找在此准则下的最优解。

1 有界椭球不确定性平差模型

平差模型为

L=AX+e

(1)

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

(2)

式中，P为n阶正定矩阵。它是椭球的特征矩阵，用来刻画椭球的形状特征，类似于e的协方差阵描述e的特征。有许多学者研究了矩阵P的构造[24-25]。在几何上，椭球的扁平程度以及椭球的体积是由椭球的特征矩阵来确定的。由于本文研究的不确定性是一种有界约束(椭球约束)，这个有界性是通过矩阵P来刻画的，它相当于e的协方差阵来刻画e的特征一样。

若L=AX是相容方程组，取X0使得L=AX0。当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)是参数向量的可行解集。

2 带有椭球不确定性约束的集员估计

首先，建立一个不确定性椭球最小化的准则来确定式(6)的集员估计解。设E(z1,P1)和E(z2,P2)为两个椭球，它们分别定义为

它们的交定义为

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

X∈E(z2,P2)}

(7)

显然，两个椭球的交集不一定是一个椭球，它可以用一个外包椭球来近似[19]，如图1所示。设其外包椭球族为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可以得到不同的测量更新椭球。为了保证椭球交的外包椭球的最小性，可以通过优化系数a，得到最小迹外包椭球。

a=argmin tr(PU)

(19)

许多文献给出了式(19)中a的计算方法，但通常较为复杂，本文方法是直接搜索。因为，0

(20)

由式(13)可知，PU是一个正定矩阵，因此a还必须满足

tr[(I-KA)Q]>0

(21)

4 算例分析

算例1为了便于画出椭圆进行分析说明，特设计如下的平差模型

L1=A1X+e1

(22)

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

误差向量e1的不确定性和参数向量X的椭圆约束信息分别定义如下

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

(23)

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

(24)

式中

计算采用Matlab的随机函数生成的满足式(23)的随机数，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中，为了验算病态模型下算法的效率，取

此算例中，A2的病态性相较于算例1有了适当的增加(为了图形的效果，没有进行更大的增加，对于较严重的病态情形，参见算例3)。计算中采用Matlab的随机函数生成的满足式(23)的随机数，e2=[0.021 60.039 3]T，利用式(19)算得：a=0.057 0，利用式(16)和式(17)计算得到

算例3设有一测边网，P1、P2为已知点，其坐标分别为(48 580.285 m,600 500.496 m)和(48 570.013 m,60 555.845 m)。为了便于分析比较，算例中的点P3、P4、P5、P6的真实坐标假设为已知(表1)，边长的观测值是利用真实坐标计算，再加上误差得到的，观测边长视为同精度(表2)。

图3 算例2中的误差椭圆，X的约束椭圆及解的不确定性椭圆Fig.3 Error ellipse,constrained ellipse of X and the uncertainty ellipse of solution in example 2

表1 真实坐标与近似坐标

表2 边长观测值

为了便于分析，假设由前期的观测得到了P3、P4、P5、P6的近似坐标(表1)，以及它们相应的点位精度。相对于近似坐标的改正数构成的未知向量为X=[x3y3x4y4x5y5x6y6]T，可以得到相应的平差模型的系数矩阵A和观测向量L。由于已知点P2的坐标非常靠近P1点，导致算法中的系数矩阵病态。

相应的平差模型为

L=AX+e

(25)

误差向量e的不确定性和参数向量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

由于在平差算法中没有直接解算带有椭球约束的平差方法，在数学上通常是拉格朗日函数求极值的方法，转化成岭估计方法。如文献[25]，将带椭球约束的线性模型估计写成

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

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

其拉格朗日函数为

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

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

求得椭球约束下的广义岭估计

文献[25]给出了其岭参数的计算方法，但同时说明“当设计阵病态时，使用这种类似于两步估计的做法应该格外小心，理论上已经表明此时不宜采用广义最小二乘估计”。因此，在本算例中，采用通常的岭参数计算方法求出岭估计再与本文的方法进行比较。

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

-1.480 02.368 22.184 4-2.907 2]T

综上，对实例解算有如下说明与分析：

(2) 从加权混合估计的角度来看，因为计算得到a=0.052 8，说明参数约束先验信息式(27)在参数估计中的作用更大。这也正好说明当模型出现病态时，利用参数先验信息可以改善其病态性。

(3) 令

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

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

(4) 本文算法中不确定度的最小化是通过求椭球最小特征矩阵的迹来实现的，也可以通过最小化的椭球体积(对应的是特征矩阵的行列式最小)，相关的算法可参看文献[8]。

(5) 算法中，a=0.052 8是一个近拟值。a从0开始，通过增量Δa=0.000 1，逐步搜索得到使tr(PU)达到最小的a。

(6) 对于病态模型的其他算法，如表3中的截断奇异值算法和岭估计算法，它们是利用数学原理来处理病态系数矩阵，不能有效地利用先验信息，计算的结果不如本文的算法。更重要的是，本文算法不仅能给出参数估计的值，而且还能对参数估计的不确定度进行估计。

5 结束语

不论是观测过程，还是未知参数本身，都存在不确定性噪声的干扰。然而，不确定性噪声非常复杂，难以确切了解诸如噪声的分布或均值和方差等统计特性。本文在椭球集合描述不确定性的基础上建立一个新的不确定性平差模型，以两个集合的交集来研究不确定度的传递。基于椭球集合特征矩阵的迹最小化建立最小不确定度平差准则，得到了在最优准则下的最优解。它与文献[27]提出的加权混合估计在形式上是一致的，都是对先验信息的利用。a的确定方法可以看作加权混合估计的权的一种新的确定方法。不确定性因素会以不确定度的形式反映在测绘数据中，其统计特性难以准确获得，需要在平差解算的同时，尽量使不确定度达到最小。不同于误差用数值来描述和度量不确定性，本文尝试用一个椭球来描述不确定性，然后用椭球特征矩阵的迹来度量不确定性的大小，这是对于不确定性描述与度量的一种尝试。

表3 法矩阵病态时算法的结果比较