GIS 叠置图层方差分量的极大似然估计

第30卷第10期2005年10月武汉大学学报・信息科学版

G eomatics and Information Science of Wuhan University Vol. 30No. 10

Oct. 2005

文章编号:167128860(2005) 1020892205文献标志码:A

GIS 叠置图层方差分量的极大似然估计

王新洲1,2　游扬声1,3　刘　星3　史文中4

(1　武汉大学测绘学院, 武汉市珞喻路129号,430079)

(2　武汉大学灾害监测与防治研究中心, 武汉市珞喻路129号,430079)

(3　重庆大学土木工程学院, 重庆市沙坪坝区沙正街,400024) (4　香港理工大学土地测量与地理资讯学系, 香港九龙红磡)

摘　要:针对GIS 叠置中的同名点, 以维希特分布密度为似然函数, 方法。该方法不依赖残差, 关键词:叠置; 同名点; 方差分量估计; 极大似然估计; 中图法分类号:P207. 2; P208

　　GIS , , 其。点元的位置不确定性是研究线元、面元不确定性的基础。在研究线元、面元等的位置不确定性和GIS 操作的误差传播时, 通常假设点元的位置不确定度已知。实际上, GIS 的数字化数据很少带有位置不

确定度信息(如协方差阵) , 这是目前不确定性的理论成果难以实际应用的重要原因。有关具体研究见文献[1～5]。为了更严密地估计GIS 叠置操作中各图层权同名点元的方差, 本文首先在顾及统计量的相关性的基础上推导出统计量的分布, 然后应用极大似然估计原理导出了同名点元各图层的方差分量的严密估计公式。

cov (x i , j , x k , l ) =0, cov (y i , j , y k , l ) =0

　　(i , j , k , l =1, 2, …, n , i ≠k 或者j ≠l ) cov (x i , j , y k , l ) =0(i , j , k , l =1, 2, …, n )

(1)

即假设位置数据误差服从正态分布, 同一图层上的坐标数据具有相同的方差, 且各图层之间以及纵横坐标之间误差独立。

要对叠置后输出的坐标值作出最优估计, 需要知道各叠置图层坐标的方差, 并根据方差定权。由于各图层的单位权因子不同, 所以需要进行方差分量估计。1. 2　方差分量估计

文献[6～12]首先对未知参数作出估计, 利用估计得到的残差构造统计量, 然后根据统计量与方差分量的函数关系计算出方差分量估值。估计出方差分量后, 重新对观测值定权, 并重新估计未知参数。如此反复迭代, 最终给出未知参数的最优估计值和方差分量的估计结果。

由于上述传统的方差分量估计法的计算工作量非常大, 因此, 本文设法避开传统方法, 先估计方差分量, 再据此方差分量定权, 不用迭代直接解出未知参数的最优估值的新方法。即针对GIS 中同名点的误差特性, 直接以真误差的统计量的密度函数作为似然函数,

给出方差分量的极大似

1　同名点元方差分量的极大似然估

计

1. 1　基本假设

设在对m (m ≥2) 幅地图进行叠置分析中获得了n 个同名点的坐标数据, 令x i , j 、y i , j 分别为第

i 层第j 个同名点的纵横坐标, x j 、y j 分别为第j

个同名点的坐标真值。设

x i , j ～N ( x j , σy j , σi ) , y i , j ～N ( i )

收稿日期:2005207226。

项目来源:国家自然科学基金资助项目(40474003) ; 香港政府研究资助局资助项目(3_ZB40) 。

　第30卷第10期王新洲等:GIS 叠置图层方差分量的极大似然估计

893

然估计, 然后根据估计得到的方差定权, 直接给出未知参数的最优估计。1. 3　似然函数

令L =(x 1y 1x 2y 2… x n y n ) , L i =(x i 1y i 1

2n ×1

中心化维希特分布, 记作S ～W m -1(2n , V ) ,

其密度函数为(证明见文献[13]) :

(det S ) n-1-2exp -f (S ) =

(-1)

tr V S 2

m-1

x i 2y i 2…x in y in ) , i =1, 2, …, m , Er i =L -L i ,

2n ×1

Er i 是真误差序列, 服从期望为0、方差阵为σI

2n ×1

2i 2n

(m-1)

(det V ) n

n +∏2j =1

的正态分布(I 2n 为2n 阶单位阵) 。

令d i -1=L i -L 1, 显然, d i -1=Er 1-2n ×1

2n ×1

Er i 。

　　　S 为正定阵

0, 　　　其他

(5)

根据方差2协方差传播律知, d k (k =1, …, m -1)

服从期望为0、方差阵为(σ1+σk +1) I 2n 的正态分布。令

(m-1) 2n ×

　　当m =2时, 易得D =(l 2, 1-l 1, 1, l 2, 2-l 1, 2,

…, l 2, 2n -l 1, 2n ) , S =

i =1

∑(l

2, i

-l 1, i ) 2, 退化为一

D =(d 1d 2, …, d m-1) =D D

1T 2…D

T 2T

个非负实数, V =(θ1+θ2) , 退化为一个正实数。再以m =2代入式(5) , 得:

S ) =

n-1

则D i 是(m -1) 维随机向量, 且D i

是真误差的线

性函数, 仍服从正态分布。易知, E (D i ) =0, 其分布密度函数为:

T -1-D i V i f (D i ) =exp 222(2π) (det V )

22σ1σ211

(m-1) ×(m-1)

21)

, S >0

, (6)

可见, 它正是(θ1+θ2) χ2n 的分布密度。因此, 维希

特分布是χ分布的推广。有人将导出维希特分

V =

1σ+σσ

2123σ…22

…σ1+σm

…2σ1

(3)

布的时间作为多元分析诞生的时间, 可见维希特

分布的重要性[13]。在测量数据处理中, χ分布在假设检验、单位权方差估计等方面应用非常广, 但目前维希特分布应用得还很少。

S 是真误差的统计量, 与未知参数无关。由式(4) 知, S 的密度值由θ确定, 于是得到θ的似然函数为:

θθθL (1, 2, …, m , S ) =f (S )

1. 4　方差分量的极大似然估计

(7)

rank (V ) =m -1

为记述方便, 令θj =σj , 将V 作如下分解:

11…111

(m-1) ×(m-1)

V 1=

(m-1) ×(m-1)

V 2=

…111000…0

0000

……

, 00

θ的极大似然估计θθ^满足:L (^, S ) =max L

(θ, S ) 。考虑到函数ln x 严格单增, max L (θ, S )

θ, S ) , 易得:Ζmaxln L (

θθln L (, S ) =

n -1-……

　　…

(m-1) ×(m-1)

lndet S -4

(-1)

tr V S -2

n (m -1) ln2--

π-n lndet V ln

m-1

V m =

…0

…

则

V =

j =1

(8) n +∑2j =1

略去与θ无关的量, 化简得:

θ, S ) Ζmin [t r (V -1S ) +2n lndet V ]max ln L (

θθ

(4)

) , 利用　　为求θ^(为了简便, 有时也将θ^记为θ

θV ∑

　　式(3) 表明, d 1, d 2, …, d m -1是不独立的。由

式(1) 知, D 1, D 2, …, D 2n 是独立的, 均服从期望为0、协方差阵为V 的m -1维正态分布。根据统计理论, 当2n >(m -1) 时, () S () =D T D 服从

m -1×m -1

矩阵及其行列式微商运算规则, 可令

-1=

-t r (V

-1

V j V

-1

S ) +2n tr (V V j ) =0(9)

式(9) 是由m 个方程组成的方程组, 有m 个未知

894

武汉大学学报・信息科学版2005年10月

数θj , 则方程组(9) 的解θ^就是方差分量θ的极大似然估计。

V =V

-1

θθ1+θ21

θθ11+θ3

1+θ3θ

-θ1

2　方差分量极大似然估计的解

2. 1　低阶符号矩阵求逆解法

θ1+θ2θθθ　　令d e =det V =θ12+θ13+θ23, 则式(9) 可写为:22

θθ2nd e (θ2+θ3) -θ3s 1, 1-223s 1, 2-θ2s 2, 2=0

d e

θ2nd e (θ1+θ3) -(1+θ3) s 1, 1+2

θθ　　21(1+θ3) s 1, 2-θ1s 2, 2=02

θθ2nd e (θ1+θ2) -θ1s 1, 1+21(1+θ2) s 1, 2-2

　　(θ1+θ2) s 2, 2=0

-θ1

当m 很小时, 容易给出符号矩阵求逆公式,

若将式(9) 化为θj 的高次多项式方程组, 这样就给求解带来很大的方便。

当m =2, 即仅两个图层叠置分析时, 方差阵

V 只有一个数θ1+θ2, S 只有一个数d 1d 1(简记为

s ) , 维希特分布也退化成χ分布。方程组(9) 容易写成:

(12)

θ(10) 1+θ2=s/2n

线性方程组(10) 秩亏, 必须增加约束才有惟一解。如按文献[2], 依图层比例尺1∶M 定权, 等价于

θ令θ2=M 21/M 1, 则可解得:

θθ=11+θ31

θ() 1=M 1s/2n M 1+M θθθ111+θ4(θ(2=M 2s/n 12

θθθθθθθθ故d e =det V =θ123+θ124+θ134+θ234

, 2=s/n 当m :

θθθθθθ-θ-θ13+θ14+θ341414-1

θθθθθθ-θ-θV =1412+θ14+θ2412

d e

式(12) 的解不惟一。注意到θ是方差分量, θ0, j ≥

则式(12) 解中的全正实数解向量即为θ的极大似然估计θ^。对于式(12) 的解法, [14]。

当m 4时,

θθ121

θ-θ13θθ12

θE (1) =

θθθ12+θ13+θ23

n (m -1) (m -2) θ1=θ1

n (m -1) (m -2)

则式(9) 仍可写成与式(12) 类似的多项式方程组。可见, 当叠置层数m =4时, 符号矩阵求逆的办法已经不很方便了; 当m ≥5时, 符号矩阵求逆相当困难。

2. 2　式(9) 的牛顿解法

因为符号矩阵求逆比较困难, 故在实际计算中采用的是牛顿解法[14]。由式(9) 易得:

tr (V -1V j V -1S ) =2n t r (V -1V j )

上式左右同乘θj , 再对j =1, 2, …, m -1求和, 顾及式(4) , 得:

t r (V -1S ) =2n (m -1)

m-1

即θ1是θ1的无偏估计。同理可知, 其他的θj 也是

000

θj 的无偏估计。如果θθj

得到初值θ后, 将式(9) 在θ处按泰勒级数展开, 取至一次项, 得:

-1-1-1

-tr (V V j V S ) +2n t r (V V j ) -

(13)

-1-10δ+|θ=θk k k =1

-10δ=02n |θ=θk

5k k =1

式中, V 是以θ代替θ构成的方差阵。

若令V =

S /2n , 式(13) 成立。注意到式(3) , 可令θ01=

[s i , j -t r (S ) ],

2n (m -1) (m -2) i , j =1

(15)

∑

θ0j =-θ1, j =2, …, m

(14)

即用S 的非主对角元素来估计θ1, 用S 的主对角

) , Q 为V 令d v =det (V 的代数余子式构成的

-1

矩阵, 则V =Q , 则有:

d v

元素来估计其余的θj 。由统计理论知, s i , j /v i , j 服

χ从自由度为2n 的χ2分布。分布具有可加性, 0

其期望为自由度。考虑到S 是对称阵, 则估计θ1

-1-1(0=|θ=θQV k QV j QS ) 3t r k d v -1(0=-|θ=θ2t r QV k QV j ) k d v

令

的自由度为n (m -1) (m -2) , 有:

　第30卷第10期王新洲等:GIS 叠置图层方差分量的极大似然估计

895

f j =2nd v t r (QV j ) -t r (QV j QS ) a j , k =2n t r (QV k QV j ) -m

d v

tr (QV k QV j QS )

3　算　例

本算例用模拟的真误差来考察本文提出的估

计方法, 为节省篇幅, 直接列出5个叠置层、10个点的真误差序列于表1(单位:0. 1mm ) 。事实上, 它们是由计算机生成的服从期望为0的正态

分布的随机数序列, 其理论方差列于σ行, 用真误差计算各层同名点元的方差估值列于θ行, 如

θ1=Er 1Er 1/20。

按式(14) 计算得初值, 列于θ行, 按式(15)

111

～式(17) 计算得θ, 列于θ行; 再以θ为初值, 反

复计算得到θ, …, 最后得到方差分量的极大似然估计θ^, 列于θ^行。为便于比较则式(13) 可写为:

k =1

∑a

j , k

δk =f j , 其矩阵形式为:

m ×1

m ×m m ×1

A δ=F

(16)

　　注意到V 是数值矩阵, 则A 、F 都是数值矩

ε>0, 若m >2, 则A 满阵, 可以计算出。设F ≥

秩, 易解得δ=A -1F , 则有:

10θ(17) =θ+δ式(14) ～式(17) 构成了用牛顿法解算式(9) 的公

式。再令θ为初值, 反复计算, 直至

^x i =

j =1

y θ/θ, ^

j =1

θ/θ

j =1

, (18)

表　Tab. 1of of Common Points

编号

2345678910

2σ

11y 22y Er 3x Er 3y Er 4x Er 4y Er 5x Er 5y

384-196-0. 2051. 194-0. 131-1. 9920. 567-0. 2270. 881. 262

0. 766-1. 104-0. 514-0. 8180. 570. 1220. 741-0. 444-0. 0830. 50. 7000. 6030. 33080. 37940. 39410. 39400. 39380. 3938470. 39395

519-0. 7730. 1990. 939-0. 585-0. 0781. 578-1. 078-0. 111-0. 272

1. 0

0. 465-0. 924-1. 3260. 2080. 7970. 7371. 0240. 399-1. 1870. 720. 7550. 3840. 55320. 63980. 67730. 68260. 68280. 6827650. 68276

-0. 542-0. 062-0. 106-1. 4010. 515-0. 61-1. 571. 461-0. 7761. 552

1. 5

-0. 3280. 4931. 319-2. 09

1. 6710. 0050. 3510. 551

1. 777-0. 513-0. 785-0. 231-3. 7463. 363-0. 678-0. 877-0. 063. 7732. 02. 9391. 9652. 19482. 30092. 30642. 30532. 30522. 3052482. 3047

1. 905-2. 159-0. 902-2. 8951. 468-0. 0690. 0781. 614-2. 066-2. 212

2. 5

-2. 074-0. 4681. 2060. 7790. 6870. 2570. 2610. 091-1. 625-0. 2612. 0663. 3162. 02362. 44462. 68762. 74542. 74792. 7479442. 74792

0. 822-1. 58

-0. 908-1. 353-0. 279-1. 099-1. 1140. 275

0. 1392. 008

-0. 811-1. 041

0θ1θ2θ3θ4θ5θ

1. 0271. 2631. 4341. 58441. 62351. 62621. 62621. 6262411.

62641

θ^

表2　叠置输出图同名点坐标误差及方差

Tab. 2　Error and Variance of the Coordinate of Common Points for the Overlapped Coverage

编号

12345

V ^x

2) 误差(σ

(H )

(θ^)

编号

678910

2) 误差(σ(H )

(θ^)

0. 877-0. 784-0. 022

0. 4440. 0380. 194

0. 791-0. 58

0. 7170. 13

0. 791-0. 580

0. 7170. 130

-0. 9610. 081

0. 7310. 143

-1. 1370. 431-0. 1120. 4620. 296

0. 6500. 2660. 2650. 3480. 272

-1. 1370. 431-0. 1120. 4620. 296

0. 6500. 2660. 2650. 3480. 272

-0. 163-0. 4980. 212-0. 937

-0. 115-0. 8020. 488-0. 001-0. 018-0. 659

0. 200-0. 0060. 134-0. 5820. 1410. 255

0. 5350. 246

0. 316-0. 6120. 316-0. 612

896

武汉大学学报・信息科学版2005年10月

　　依式(18) , 采用σ估计^x 的方差为0. 219, 采用θ^(或者H ) 估计^x 的方差为0. 185, 而采用^x 的真误差估计得^x 的方差为0. 29左右。这几个数之间的差值比预期的大。更多的模拟计算表明, 按式(18) 估计^x 的方差, 由于没有考虑到θ^的误差, 一般会较真实的方差小。

从表1可以看出, θ^的值与σ或者θ的值相

6　崔希璋, 於宗俦, 陶本藻, 等. 广义测量平差(新版) .

武汉:武汉测绘科技大学出版社, 2001

7　於宗俦. 方差2协方差分量极大似然估计的通用公式.

测绘学报, 1994,23(2) :6～13

8　王新洲, 刘经南, 陶本藻. GPS 基线向量网方差和协方

差分量的MINQU E 估计. 武汉测绘科技大学学报,

1993,18(3) :57～66

9　王新洲, 于正林. 方差分量估计的快速算法. 武测科

差不大, 与H 的值基本相等; 从表2可以看出, 用赫尔默特法和本文提出的方法估计得到的叠置输出图上的同名点坐标是完全一致的。这验证了本文所提方法的正确性和可行性。

参　考　文　献

1　刘文宝. GIS 空间数据的不确定理论:[博士论文].武

技, 1994

10　王新洲. 变二次估计为线性估计的拉直变换法. 武

汉测绘科技大学学报, 1994,19(3) :239～243

11　王新洲, 刘经南, 陶本藻. 稳健最优不变二次无偏估

计. 武汉测绘科技大学学报, 1995,20(3) :240～245

12　彭军还, 陶本藻, 黄　

, 等. 基于M 残差的方差分

量估计. 见:朱耀仲, 孙和平. 大地测量与地球动力学进展. 武汉:, 2004. 824～

836

13. . 北京:科学出

汉:武汉测绘科技大学,1995

2　刘文宝, 史文中. GIS 叠置前后同名点元的方差估计.

武汉测绘科技大学学报, 1996,21(4) :355～360

3　戴洪磊, 徐泮林, 卢秀山. GIS , . . 武汉:武

方差估计. ,23(1) :14～17

4　陶本藻, . 汉大学出版社,2002

第一作者简介:王新洲, 教授, 博士生导师。主要从事空间信息处理理论与应用方面的研究。

E 2mail :[email protected]. edu. cn

方法. , 2001, 26(2) :101～

104

5　游扬声, 王新洲, 汪海洪. 利用GIS 叠置中的同名点

元估计方差. 测绘工程,2002,11(4) :18～21

Maximum Likelihood Estimation for V ariance of GIS

Overlapped Coverage

A N G X i nz hou

1, 2

　YOU Yangs heng 1, 3　L I U X i ng 3　S H I W enz hong 4

(1　School of Geodesy and Geomatics , Wuhan University , 129Luoyu Road , Wuhan 430079, China )

(2　Hazard Monitoring and Prevention Research Center ,Wuhan University , 129L uoyu Road ,Wuhan 430079, China ) (3　School of Civil Engineering , Chongqing University , Shazheng Street , Shapingba District , Chongqing 400024, China ) (4　Depart ment of Land Surveying and Geo 2informatics , Hong K ong Polytechnic University , K owloon , Hong K ong )

Abstract :Taking on Wishart dist ribution f unction as likelihood f unction , t his paper repre 2sent s a maximum likelihood estimatio n of variance component s of every coverage aimed at common point s in t he overlapped coverage. The met hod can estimate t he unknown parame 2ters and variance component s independent of remnant errors and iterative comp utation. K ey w ords :overlapped 2coverage ; common 2point s ; variance component s estimation ; maxi 2mum likelihood estimation ; Wishart dist ribution

About the f irst author :WANG X inzhou , professor , Ph. D supervisor. He is concentrated on the research and education in the theory and ap 2plication of spatial data processing. E 2mail :[email protected]. edu. cn

第30卷第10期2005年10月武汉大学学报・信息科学版

G eomatics and Information Science of Wuhan University Vol. 30No. 10

Oct. 2005

文章编号:167128860(2005) 1020892205文献标志码:A

GIS 叠置图层方差分量的极大似然估计

王新洲1,2　游扬声1,3　刘　星3　史文中4

(1　武汉大学测绘学院, 武汉市珞喻路129号,430079)

(2　武汉大学灾害监测与防治研究中心, 武汉市珞喻路129号,430079)

(3　重庆大学土木工程学院, 重庆市沙坪坝区沙正街,400024) (4　香港理工大学土地测量与地理资讯学系, 香港九龙红磡)

cov (x i , j , x k , l ) =0, cov (y i , j , y k , l ) =0

　　(i , j , k , l =1, 2, …, n , i ≠k 或者j ≠l ) cov (x i , j , y k , l ) =0(i , j , k , l =1, 2, …, n )

(1)

即假设位置数据误差服从正态分布, 同一图层上的坐标数据具有相同的方差, 且各图层之间以及纵横坐标之间误差独立。

给出方差分量的极大似

1　同名点元方差分量的极大似然估

计

1. 1　基本假设

设在对m (m ≥2) 幅地图进行叠置分析中获得了n 个同名点的坐标数据, 令x i , j 、y i , j 分别为第

i 层第j 个同名点的纵横坐标, x j 、y j 分别为第j

个同名点的坐标真值。设

x i , j ～N ( x j , σy j , σi ) , y i , j ～N ( i )

收稿日期:2005207226。

项目来源:国家自然科学基金资助项目(40474003) ; 香港政府研究资助局资助项目(3_ZB40) 。

　第30卷第10期王新洲等:GIS 叠置图层方差分量的极大似然估计

893

然估计, 然后根据估计得到的方差定权, 直接给出未知参数的最优估计。1. 3　似然函数

令L =(x 1y 1x 2y 2… x n y n ) , L i =(x i 1y i 1

2n ×1

中心化维希特分布, 记作S ～W m -1(2n , V ) ,

其密度函数为(证明见文献[13]) :

(det S ) n-1-2exp -f (S ) =

(-1)

tr V S 2

m-1

x i 2y i 2…x in y in ) , i =1, 2, …, m , Er i =L -L i ,

2n ×1

Er i 是真误差序列, 服从期望为0、方差阵为σI

2n ×1

2i 2n

(m-1)

(det V ) n

n +∏2j =1

的正态分布(I 2n 为2n 阶单位阵) 。

令d i -1=L i -L 1, 显然, d i -1=Er 1-2n ×1

2n ×1

Er i 。

　　　S 为正定阵

0, 　　　其他

(5)

根据方差2协方差传播律知, d k (k =1, …, m -1)

服从期望为0、方差阵为(σ1+σk +1) I 2n 的正态分布。令

(m-1) 2n ×

　　当m =2时, 易得D =(l 2, 1-l 1, 1, l 2, 2-l 1, 2,

…, l 2, 2n -l 1, 2n ) , S =

i =1

∑(l

2, i

-l 1, i ) 2, 退化为一

D =(d 1d 2, …, d m-1) =D D

1T 2…D

T 2T

个非负实数, V =(θ1+θ2) , 退化为一个正实数。再以m =2代入式(5) , 得:

S ) =

n-1

则D i 是(m -1) 维随机向量, 且D i

是真误差的线

性函数, 仍服从正态分布。易知, E (D i ) =0, 其分布密度函数为:

T -1-D i V i f (D i ) =exp 222(2π) (det V )

22σ1σ211

(m-1) ×(m-1)

21)

, S >0

, (6)

可见, 它正是(θ1+θ2) χ2n 的分布密度。因此, 维希

特分布是χ分布的推广。有人将导出维希特分

V =

1σ+σσ

2123σ…22

…σ1+σm

…2σ1

(3)

布的时间作为多元分析诞生的时间, 可见维希特

分布的重要性[13]。在测量数据处理中, χ分布在假设检验、单位权方差估计等方面应用非常广, 但目前维希特分布应用得还很少。

S 是真误差的统计量, 与未知参数无关。由式(4) 知, S 的密度值由θ确定, 于是得到θ的似然函数为:

θθθL (1, 2, …, m , S ) =f (S )

1. 4　方差分量的极大似然估计

(7)

rank (V ) =m -1

为记述方便, 令θj =σj , 将V 作如下分解:

11…111

(m-1) ×(m-1)

V 1=

(m-1) ×(m-1)

V 2=

…111000…0

0000

……

, 00

θ的极大似然估计θθ^满足:L (^, S ) =max L

(θ, S ) 。考虑到函数ln x 严格单增, max L (θ, S )

θ, S ) , 易得:Ζmaxln L (

θθln L (, S ) =

n -1-……

　　…

(m-1) ×(m-1)

lndet S -4

(-1)

tr V S -2

n (m -1) ln2--

π-n lndet V ln

m-1

V m =

…0

…

则

V =

j =1

(8) n +∑2j =1

略去与θ无关的量, 化简得:

θ, S ) Ζmin [t r (V -1S ) +2n lndet V ]max ln L (

θθ

(4)

) , 利用　　为求θ^(为了简便, 有时也将θ^记为θ

θV ∑

　　式(3) 表明, d 1, d 2, …, d m -1是不独立的。由

式(1) 知, D 1, D 2, …, D 2n 是独立的, 均服从期望为0、协方差阵为V 的m -1维正态分布。根据统计理论, 当2n >(m -1) 时, () S () =D T D 服从

m -1×m -1

矩阵及其行列式微商运算规则, 可令

-1=

-t r (V

-1

V j V

-1

S ) +2n tr (V V j ) =0(9)

式(9) 是由m 个方程组成的方程组, 有m 个未知

894

武汉大学学报・信息科学版2005年10月

数θj , 则方程组(9) 的解θ^就是方差分量θ的极大似然估计。

V =V

-1

θθ1+θ21

θθ11+θ3

1+θ3θ

-θ1

2　方差分量极大似然估计的解

2. 1　低阶符号矩阵求逆解法

θ1+θ2θθθ　　令d e =det V =θ12+θ13+θ23, 则式(9) 可写为:22

θθ2nd e (θ2+θ3) -θ3s 1, 1-223s 1, 2-θ2s 2, 2=0

d e

θ2nd e (θ1+θ3) -(1+θ3) s 1, 1+2

θθ　　21(1+θ3) s 1, 2-θ1s 2, 2=02

θθ2nd e (θ1+θ2) -θ1s 1, 1+21(1+θ2) s 1, 2-2

　　(θ1+θ2) s 2, 2=0

-θ1

当m 很小时, 容易给出符号矩阵求逆公式,

若将式(9) 化为θj 的高次多项式方程组, 这样就给求解带来很大的方便。

当m =2, 即仅两个图层叠置分析时, 方差阵

V 只有一个数θ1+θ2, S 只有一个数d 1d 1(简记为

s ) , 维希特分布也退化成χ分布。方程组(9) 容易写成:

(12)

θ(10) 1+θ2=s/2n

线性方程组(10) 秩亏, 必须增加约束才有惟一解。如按文献[2], 依图层比例尺1∶M 定权, 等价于

θ令θ2=M 21/M 1, 则可解得:

θθ=11+θ31

θ() 1=M 1s/2n M 1+M θθθ111+θ4(θ(2=M 2s/n 12

θθθθθθθθ故d e =det V =θ123+θ124+θ134+θ234

, 2=s/n 当m :

θθθθθθ-θ-θ13+θ14+θ341414-1

θθθθθθ-θ-θV =1412+θ14+θ2412

d e

式(12) 的解不惟一。注意到θ是方差分量, θ0, j ≥

则式(12) 解中的全正实数解向量即为θ的极大似然估计θ^。对于式(12) 的解法, [14]。

当m 4时,

θθ121

θ-θ13θθ12

θE (1) =

θθθ12+θ13+θ23

n (m -1) (m -2) θ1=θ1

n (m -1) (m -2)

则式(9) 仍可写成与式(12) 类似的多项式方程组。可见, 当叠置层数m =4时, 符号矩阵求逆的办法已经不很方便了; 当m ≥5时, 符号矩阵求逆相当困难。

2. 2　式(9) 的牛顿解法

因为符号矩阵求逆比较困难, 故在实际计算中采用的是牛顿解法[14]。由式(9) 易得:

tr (V -1V j V -1S ) =2n t r (V -1V j )

上式左右同乘θj , 再对j =1, 2, …, m -1求和, 顾及式(4) , 得:

t r (V -1S ) =2n (m -1)

m-1

即θ1是θ1的无偏估计。同理可知, 其他的θj 也是

000

θj 的无偏估计。如果θθj

得到初值θ后, 将式(9) 在θ处按泰勒级数展开, 取至一次项, 得:

-1-1-1

-tr (V V j V S ) +2n t r (V V j ) -

(13)

-1-10δ+|θ=θk k k =1

-10δ=02n |θ=θk

5k k =1

式中, V 是以θ代替θ构成的方差阵。

若令V =

S /2n , 式(13) 成立。注意到式(3) , 可令θ01=

[s i , j -t r (S ) ],

2n (m -1) (m -2) i , j =1

(15)

∑

θ0j =-θ1, j =2, …, m

(14)

即用S 的非主对角元素来估计θ1, 用S 的主对角

) , Q 为V 令d v =det (V 的代数余子式构成的

-1

矩阵, 则V =Q , 则有:

d v

元素来估计其余的θj 。由统计理论知, s i , j /v i , j 服

χ从自由度为2n 的χ2分布。分布具有可加性, 0

其期望为自由度。考虑到S 是对称阵, 则估计θ1

-1-1(0=|θ=θQV k QV j QS ) 3t r k d v -1(0=-|θ=θ2t r QV k QV j ) k d v

令

的自由度为n (m -1) (m -2) , 有:

　第30卷第10期王新洲等:GIS 叠置图层方差分量的极大似然估计

895

f j =2nd v t r (QV j ) -t r (QV j QS ) a j , k =2n t r (QV k QV j ) -m

d v

tr (QV k QV j QS )

3　算　例

本算例用模拟的真误差来考察本文提出的估

计方法, 为节省篇幅, 直接列出5个叠置层、10个点的真误差序列于表1(单位:0. 1mm ) 。事实上, 它们是由计算机生成的服从期望为0的正态

分布的随机数序列, 其理论方差列于σ行, 用真误差计算各层同名点元的方差估值列于θ行, 如

θ1=Er 1Er 1/20。

按式(14) 计算得初值, 列于θ行, 按式(15)

111

～式(17) 计算得θ, 列于θ行; 再以θ为初值, 反

复计算得到θ, …, 最后得到方差分量的极大似然估计θ^, 列于θ^行。为便于比较则式(13) 可写为:

k =1

∑a

j , k

δk =f j , 其矩阵形式为:

m ×1

m ×m m ×1

A δ=F

(16)

　　注意到V 是数值矩阵, 则A 、F 都是数值矩

ε>0, 若m >2, 则A 满阵, 可以计算出。设F ≥

秩, 易解得δ=A -1F , 则有:

10θ(17) =θ+δ式(14) ～式(17) 构成了用牛顿法解算式(9) 的公

式。再令θ为初值, 反复计算, 直至

^x i =

j =1

y θ/θ, ^

j =1

θ/θ

j =1

, (18)

表　Tab. 1of of Common Points

编号

2345678910

2σ

11y 22y Er 3x Er 3y Er 4x Er 4y Er 5x Er 5y

384-196-0. 2051. 194-0. 131-1. 9920. 567-0. 2270. 881. 262

0. 766-1. 104-0. 514-0. 8180. 570. 1220. 741-0. 444-0. 0830. 50. 7000. 6030. 33080. 37940. 39410. 39400. 39380. 3938470. 39395

519-0. 7730. 1990. 939-0. 585-0. 0781. 578-1. 078-0. 111-0. 272

1. 0

0. 465-0. 924-1. 3260. 2080. 7970. 7371. 0240. 399-1. 1870. 720. 7550. 3840. 55320. 63980. 67730. 68260. 68280. 6827650. 68276

-0. 542-0. 062-0. 106-1. 4010. 515-0. 61-1. 571. 461-0. 7761. 552

1. 5

-0. 3280. 4931. 319-2. 09

1. 6710. 0050. 3510. 551

1. 777-0. 513-0. 785-0. 231-3. 7463. 363-0. 678-0. 877-0. 063. 7732. 02. 9391. 9652. 19482. 30092. 30642. 30532. 30522. 3052482. 3047

1. 905-2. 159-0. 902-2. 8951. 468-0. 0690. 0781. 614-2. 066-2. 212

2. 5

-2. 074-0. 4681. 2060. 7790. 6870. 2570. 2610. 091-1. 625-0. 2612. 0663. 3162. 02362. 44462. 68762. 74542. 74792. 7479442. 74792

0. 822-1. 58

-0. 908-1. 353-0. 279-1. 099-1. 1140. 275

0. 1392. 008

-0. 811-1. 041

0θ1θ2θ3θ4θ5θ

1. 0271. 2631. 4341. 58441. 62351. 62621. 62621. 6262411.

62641

θ^

表2　叠置输出图同名点坐标误差及方差

Tab. 2　Error and Variance of the Coordinate of Common Points for the Overlapped Coverage

编号

12345

V ^x

2) 误差(σ

(H )

(θ^)

编号

678910

2) 误差(σ(H )

(θ^)

0. 877-0. 784-0. 022

0. 4440. 0380. 194

0. 791-0. 58

0. 7170. 13

0. 791-0. 580

0. 7170. 130

-0. 9610. 081

0. 7310. 143

-1. 1370. 431-0. 1120. 4620. 296

0. 6500. 2660. 2650. 3480. 272

-1. 1370. 431-0. 1120. 4620. 296

0. 6500. 2660. 2650. 3480. 272

-0. 163-0. 4980. 212-0. 937

-0. 115-0. 8020. 488-0. 001-0. 018-0. 659

0. 200-0. 0060. 134-0. 5820. 1410. 255

0. 5350. 246

0. 316-0. 6120. 316-0. 612

896

武汉大学学报・信息科学版2005年10月

从表1可以看出, θ^的值与σ或者θ的值相

6　崔希璋, 於宗俦, 陶本藻, 等. 广义测量平差(新版) .

武汉:武汉测绘科技大学出版社, 2001

7　於宗俦. 方差2协方差分量极大似然估计的通用公式.

测绘学报, 1994,23(2) :6～13

8　王新洲, 刘经南, 陶本藻. GPS 基线向量网方差和协方

差分量的MINQU E 估计. 武汉测绘科技大学学报,

1993,18(3) :57～66

9　王新洲, 于正林. 方差分量估计的快速算法. 武测科

参　考　文　献

1　刘文宝. GIS 空间数据的不确定理论:[博士论文].武

技, 1994

10　王新洲. 变二次估计为线性估计的拉直变换法. 武

汉测绘科技大学学报, 1994,19(3) :239～243

11　王新洲, 刘经南, 陶本藻. 稳健最优不变二次无偏估

计. 武汉测绘科技大学学报, 1995,20(3) :240～245

12　彭军还, 陶本藻, 黄　

, 等. 基于M 残差的方差分

量估计. 见:朱耀仲, 孙和平. 大地测量与地球动力学进展. 武汉:, 2004. 824～

836

13. . 北京:科学出

汉:武汉测绘科技大学,1995

2　刘文宝, 史文中. GIS 叠置前后同名点元的方差估计.

武汉测绘科技大学学报, 1996,21(4) :355～360

3　戴洪磊, 徐泮林, 卢秀山. GIS , . . 武汉:武

方差估计. ,23(1) :14～17

4　陶本藻, . 汉大学出版社,2002

第一作者简介:王新洲, 教授, 博士生导师。主要从事空间信息处理理论与应用方面的研究。

E 2mail :[email protected]. edu. cn

方法. , 2001, 26(2) :101～

104

5　游扬声, 王新洲, 汪海洪. 利用GIS 叠置中的同名点

元估计方差. 测绘工程,2002,11(4) :18～21

Maximum Likelihood Estimation for V ariance of GIS

Overlapped Coverage

A N G X i nz hou

1, 2

　YOU Yangs heng 1, 3　L I U X i ng 3　S H I W enz hong 4

(1　School of Geodesy and Geomatics , Wuhan University , 129Luoyu Road , Wuhan 430079, China )

GIS 叠置图层方差分量的极大似然估计

相关内容

热门内容

标签