缺失數(shù)據(jù)下基于EM算法的GPS高程擬合*

林東方 宋迎春 杜 琨 肖琴琴

(中南大學(xué)地球科學(xué)與信息物理學(xué)院測(cè)繪與國(guó)土信息工程系，長(zhǎng)沙 410083)

缺失數(shù)據(jù)下基于EM算法的GPS高程擬合*

林東方 宋迎春 杜 琨 肖琴琴

(中南大學(xué)地球科學(xué)與信息物理學(xué)院測(cè)繪與國(guó)土信息工程系，長(zhǎng)沙 410083)

GPS高程擬合常因觀測(cè)條件限制，部分重要擬合點(diǎn)的數(shù)據(jù)無(wú)法獲取，致擬合數(shù)據(jù)缺失或數(shù)據(jù)量不充足。采用常規(guī)方法擬合，會(huì)降低似大地水準(zhǔn)面的擬合精度。應(yīng)用EM算法，添加了有益于高程擬合的“潛在數(shù)據(jù)”，即缺失數(shù)據(jù)的條件期望，有效提高了缺失數(shù)據(jù)下的GPS高程擬合精度。

缺失數(shù)據(jù);EM算法;高程擬合;平面相關(guān)擬合;二次曲面擬合

1 引言

高程系統(tǒng)包括大地高、正高和正常高系統(tǒng)。大地高系統(tǒng)的基準(zhǔn)面為參考橢球面，正高系統(tǒng)的基準(zhǔn)面為大地水準(zhǔn)面，而正常高系統(tǒng)以似大地水準(zhǔn)面為基準(zhǔn)面。我國(guó)規(guī)定采用正常高高程系統(tǒng)作為我國(guó)高程的統(tǒng)一系統(tǒng)［1］。GPS高程觀測(cè)數(shù)據(jù)為大地高系統(tǒng)下的高程數(shù)據(jù)，由于基準(zhǔn)面的不同，大地高系統(tǒng)與正常高系統(tǒng)存在一定的差距稱為高程異常。GPS高程數(shù)據(jù)需要進(jìn)行高程擬合，將大地高高程轉(zhuǎn)換為正常高高程。

平面相關(guān)擬合法和二次曲面擬合法是GPS高程擬合的主要方法，平面相關(guān)擬合法認(rèn)為水準(zhǔn)面為一平面適用于小區(qū)域的擬合，二次曲面擬合法認(rèn)為水準(zhǔn)面為二次曲面，觀測(cè)區(qū)域較大時(shí)，擬合效果明顯優(yōu)于平面相關(guān)擬合法［2］。兩種擬合方法同時(shí)也受到擬合點(diǎn)數(shù)據(jù)量和點(diǎn)位分布的影響［3］，若擬合點(diǎn)足夠多且分布均勻，則有益于水準(zhǔn)面的擬合。觀測(cè)條件的限制(如GPS無(wú)信號(hào)、環(huán)境惡劣水準(zhǔn)測(cè)量無(wú)法進(jìn)行、控制點(diǎn)被破壞或找不到)會(huì)造成部分重要的擬合點(diǎn)無(wú)法觀測(cè)。缺失數(shù)據(jù)則直接影響水準(zhǔn)面的擬合效果，降低GPS高程測(cè)量精度。對(duì)于缺失數(shù)據(jù)的處理，國(guó)內(nèi)外學(xué)者通常選擇的方法是只應(yīng)用剩余的數(shù)據(jù)進(jìn)行擬合，這樣降低了水準(zhǔn)面擬合的準(zhǔn)確度，或者采用內(nèi)插方法補(bǔ)充數(shù)據(jù)，結(jié)果并無(wú)太大改進(jìn)［4］。EM(Expectation-Maximization)算法［5-8］是用于非完全數(shù)據(jù)參數(shù)估計(jì)的一種有效方法，它是一種數(shù)據(jù)添加算法，即在觀測(cè)數(shù)據(jù)的基礎(chǔ)上添加一些“潛在數(shù)據(jù)”，利用對(duì)數(shù)據(jù)處理有益的潛在信息，得到數(shù)據(jù)缺失下參數(shù)的最優(yōu)估計(jì)。當(dāng)觀測(cè)數(shù)據(jù)不完全時(shí)，可以利用觀測(cè)數(shù)據(jù)所服從的一些規(guī)律，對(duì)缺失數(shù)據(jù)的取值范圍加以限制，即采用EM算法，利用缺失數(shù)據(jù)在現(xiàn)有條件下的期望值，對(duì)不完全數(shù)據(jù)進(jìn)行處理計(jì)算，最終可以得到一個(gè)很好的參數(shù)估計(jì)結(jié)果。

2 EM算法及其在缺失數(shù)據(jù)下的GPS高程擬合方法

2.1 EM算法

EM算法是一種迭代算法，它的每一步迭代由兩步組成:E步(求期望)和M步(極大化)。E步在給定己觀測(cè)到的數(shù)據(jù)和現(xiàn)有參數(shù)下，求“缺失數(shù)據(jù)”的條件期望;M步計(jì)算參數(shù)的MLE估計(jì)，這與己知似然求參數(shù)的MLE估計(jì)的計(jì)算方法一致。

具體地講，我們以P(θ/Y)表示θ的基于觀測(cè)數(shù)據(jù)的后驗(yàn)分布密度，稱為觀測(cè)數(shù)據(jù)后驗(yàn)分布;以P (θ/Y，Z)表示添加數(shù)據(jù)Z(缺失數(shù)據(jù))后得到的關(guān)于θ的后驗(yàn)分布密度函數(shù)，稱為完全數(shù)據(jù)后驗(yàn)分布;P (Z/θ，Y)表示在給定θ和觀測(cè)數(shù)據(jù)Y下潛在數(shù)據(jù)Z (即缺失數(shù)據(jù))的條件分布密度函數(shù)。我們的目的是計(jì)算觀測(cè)后驗(yàn)分布P(θ/Y)的參數(shù)，于是記θi為第i+1次迭代開(kāi)始時(shí)后驗(yàn)參數(shù)的估計(jì)值，則第i+1次迭代為:

E步:將P(θ/Y，Z)或logP(θ/Y，Z)關(guān)于Z的條件分布求期望，從而把Z積掉，即

M步:將Q(θ/θi，Y)極大化，即找一個(gè)點(diǎn)θi+1，使

如此形成了一次迭代θi→θi+1，θi+1∈M(θi)，這里M(θi)是在整個(gè)參數(shù)空間Ω內(nèi)使得Q(θi+1/θi，Y)最大的θ的值所組成的集合。將上述E步和M步進(jìn)行迭代直至 θi+1-θi或者 Q(θi+1/θi，Y)-Q(θi/θi，Y)充分小時(shí)為止。

EM算法在每一次迭代后均提高觀測(cè)極大似然密度函數(shù)值，具有良好的全局收斂性［8，9］。

2.2 GPS高程擬合數(shù)學(xué)模型

設(shè)GPS高程擬合數(shù)學(xué)模型為

式中，ζ為高程異常，a0，a1，…，a5為擬合參數(shù)，x、y為擬合點(diǎn)坐標(biāo)。等式右邊取前三項(xiàng)為平面擬合法(三參數(shù)法)，取前四項(xiàng)為平面相關(guān)擬合法(四參數(shù)法)，取六項(xiàng)則為二次曲面擬合法(六參數(shù)法)。

若測(cè)區(qū)內(nèi)有n(n≥6)個(gè)控制點(diǎn)，且經(jīng)過(guò)GPS觀測(cè)已知它們的高程異常，根據(jù)式(3)列出誤差方程為［10］:

式中，

利用最小二乘原理［10，11］計(jì)算參數(shù)X的解為:

二次曲面擬合參數(shù)確定后可確定局部區(qū)域的高程異常曲面(似大地水準(zhǔn)面)，繼而推算出其他待定點(diǎn)的高程異常，根據(jù)GPS高程觀測(cè)值計(jì)算出待定點(diǎn)的正常高。

2.3 缺失數(shù)據(jù)下基于EM算法的GPS高程擬合方法

由于某一區(qū)域內(nèi)GPS無(wú)信號(hào);環(huán)境惡劣水準(zhǔn)測(cè)量無(wú)法進(jìn)行;控制點(diǎn)被破壞或找不到等因素，造成部分重要的擬合點(diǎn)數(shù)據(jù)缺失，采用二次曲面擬合時(shí)將會(huì)影響高程異常曲面的準(zhǔn)確度。當(dāng)缺失數(shù)據(jù)過(guò)多時(shí)(n≤6)，導(dǎo)致二次曲面擬合法無(wú)法使用，僅能采用平面擬合法或平面相關(guān)擬合法;測(cè)區(qū)過(guò)大時(shí)，采用這兩種方法無(wú)法保證GPS高程的測(cè)量精度。EM算法可有效解決上述問(wèn)題，二次曲面擬合法的誤差方程中誤差向量V服從高斯正態(tài)分布，應(yīng)用EM算法可建立似然函數(shù)方程P(θ/Y，Z)，

假設(shè)擬合數(shù)據(jù)ζt缺失，式中θ為待估參數(shù)(a0，a1，…，a5，σ)，Y為不完全觀測(cè)數(shù)據(jù)，Z為缺失數(shù)據(jù)ζt。

誤差向量V服從高斯正態(tài)分布，因此可以得到缺失數(shù)據(jù)ζt的條件分布概率密度函數(shù)為:

式中，θi為經(jīng)過(guò)i次迭代后的參數(shù)值。由式(6)、(7)得到EM算法的期望步:

EM算法的極大化步，將Q(θ/θi，Y)極大化，積分后對(duì)各參數(shù)求偏導(dǎo)數(shù)，計(jì)算得到參數(shù)θi+1，將θi+1代入 E步進(jìn)行迭代，循環(huán)直至 θi+1-θi或者Q(θi+1/θi，Y)-Q(θi/θi，Y)充分小時(shí)為止。

3 算例分析

以杭州灣跨海大橋的GPS高程擬合數(shù)據(jù)為例，GPS控制網(wǎng)點(diǎn)位分布見(jiàn)圖1。

圖1 GPS控制網(wǎng)點(diǎn)位略圖Fig.1 GPS control points network

圖1中，GP01至GP07為分布均勻的GPS控制網(wǎng)點(diǎn)，其點(diǎn)位坐標(biāo)和高程異常均已知，可作為GPS高程擬合點(diǎn)使用。

現(xiàn)假設(shè)GP03、GP06點(diǎn)因觀測(cè)條件限制數(shù)據(jù)缺失，擬合點(diǎn)僅有GP01、GP02、GP04、GP05和GP07 5個(gè)數(shù)據(jù)。由于達(dá)不到二次曲線擬合法所需擬合點(diǎn)n≥6的條件，因此采用平面相關(guān)擬合法。應(yīng)用Matlab編程計(jì)算得平面相關(guān)擬合結(jié)果見(jiàn)表1。

采用EM算法進(jìn)行數(shù)據(jù)處理，數(shù)學(xué)模型為二次曲面擬合模型，θ為參數(shù)(a0，a1，…，a5，σ)期望步計(jì)算GP06點(diǎn)的條件期望，以GP05、GP07點(diǎn)的高程異常值為積分區(qū)間，應(yīng)用 Matlab編程迭代計(jì)算，當(dāng)θi+1-θi的值小于0.000 01時(shí)停止迭代，擬合結(jié)果見(jiàn)表2。

完全數(shù)據(jù)下采用二次曲面擬合法，以GP01、GP02、GP04、GP05、GP06和 GP07為擬合點(diǎn)，計(jì)算結(jié)果與平面相關(guān)擬合法和EM算法比較結(jié)果見(jiàn)表3。

表1 平面相關(guān)擬合結(jié)果(單位：cm)Tab.1 Results of related plane fitting(unit：cm)

表2 EM算法擬合結(jié)果(單位：cm)Tab.2 Results from fitting with EM algorithm(unit：cm)

表3 擬合結(jié)果統(tǒng)計(jì)(單位：cm)Tab.3 Statistics of fitting results(unit：cm)

表3表明，在缺失數(shù)據(jù)的情況下，使用平面相關(guān)擬合法很難保證高程擬合的精度，而采用EM算法，引入了缺失數(shù)據(jù)的條件期望，有效地提高了高程擬合精度。與完全數(shù)據(jù)下的二次曲面擬合法相比較，內(nèi)符合精度相差很小，外符合精度提高，表明了EM算法的有效性。

4 結(jié)論

擬合似大地水準(zhǔn)面的精度直接影響著GPS高程測(cè)量的精度，高程擬合時(shí)擬合點(diǎn)分布均勻且數(shù)據(jù)充足，才能較好地?cái)M合出似大地水準(zhǔn)面。而在擬合點(diǎn)觀測(cè)中難免會(huì)遇到部分重要的擬合點(diǎn)因觀測(cè)條件限制而無(wú)法觀測(cè)，造成數(shù)據(jù)缺失。采用不完全的擬合點(diǎn)數(shù)據(jù)擬合似大地水準(zhǔn)面將降低似大地水準(zhǔn)面的擬合準(zhǔn)確度。應(yīng)用EM算法，將觀測(cè)數(shù)據(jù)的誤差分布信息引入數(shù)據(jù)計(jì)算中，引入缺失數(shù)據(jù)的條件期望，添加了“潛在數(shù)據(jù)”，較好地改善了缺失數(shù)據(jù)下的擬合似大地水準(zhǔn)面精度，與完全數(shù)據(jù)下的擬合相比，提高了外符合精度。因此，在缺失數(shù)據(jù)的情況下，應(yīng)用EM算法進(jìn)行GPS高程擬合是一種行之有效的方法。

1 孔祥元，郭際明，劉宗泉.大地測(cè)量學(xué)基礎(chǔ)［M］.武漢:武漢大學(xué)出版社，2002.(Kong Xiangyuan，Guo Jiming and Liu Zhongquan.The base of geodesy［M］.Wuhan;Wuhan University Press，2002)

2 魏立峰，何建國(guó).GPS高程擬合似大地水準(zhǔn)面的方法［J］.地理空間信息，2010，8(4):72-76.(Wei Lifeng and He Jianguo.Methods for GPS height fitting quasi-geoid［J］.Geospatial Information，2010，8(4):72-76)

3 姚吉利，曲國(guó)慶，劉科利.擬合點(diǎn)分布與GPS水準(zhǔn)面擬合精度關(guān)系研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，2008，(5): 50-54.(Yao Jili，Qu Guoqing and Liu Keli.Research on accuracy of quasi-geoid with fitting methods under different distribution of GPS points［J］.Journal of Geodesy and Geodynamics，2008，(5):50-54)

4 姚喜，欒學(xué)科，王志博.GPS水準(zhǔn)擬合方法的研究［J］.測(cè)繪科學(xué)，2010，35:42-43.(Yao Xi，Luan Xueke and Wang Zhibo.Research on GPS level fitting methods［J］.Science of Surveying and Mapping，2010，35:42-43)

5 曾傳璜，等.EM算法的研究［J］.軟件導(dǎo)刊，2008，7(9): 97-98.(Zeng Chuanhuang，et al.Research of EM algorithm［J］.Software Guide，2008，7(9):97-98)

6 Dempster A P，Laird N M and Rubin D B.Maximum likelihood from incomplete data via the EM algorithm［J］.Journal of the Royal Statistifcal Society B，1977，39:1-38.

7 錢俊，舒寧.基于EM算法和單幅雷達(dá)圖像陰影的控制點(diǎn)坡度校正［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2004，29 (12):1 089-1 092.(Qian Jun and Shu Ning.Correction of control point slope based on EM algorithm and shading of single SAR image［J］.Geomatics and Information Science of Wuhan University，2004，29(12):1 089-1 092)

8 Graham C G and Juan C A.Approximate EM algorithms for parameter and state estimation in nonlinear stochastic models［A］.Proceedings of the 44th IEEE conference on decision and control，and the European Control Conference［C］.2005，12-15.

9 王兆軍.EM算法收斂的必要條件［J］.南開(kāi)大學(xué)學(xué)報(bào)(自然科學(xué)版)，1994，(2):85-88.(Wang Zhaojun.The necessary condition on the convergence of the EM algorithm［J］.Acta Scientiarum Naturalium Universitatis Nankaiensis，1994，(2):85-88)

10 武漢大學(xué)測(cè)繪學(xué)院測(cè)量平差學(xué)科組.誤差理論與測(cè)量平差基礎(chǔ)［M］.武漢:武漢大學(xué)出版社，2003.(Research group of Surveying Adjustment，School of Surveying and Mapping，Wuhan University.The base of errors theory and surveying adjustment［M］.Wuhan;Wuhan University Press，2003)

11 楊元喜，張麗萍.中國(guó)大地測(cè)量數(shù)據(jù)處理60年重要進(jìn)展［J］.地理空間信息，2010，8(1):1-6.(Yang Yuanxi and Zhang Liping.Progress of geodetic data processing for 60 years in China［J］.Geospatial Information，2010，8 (1):1-6)

GPS ELEVATION FITTING BASED ON EM ALGORITHM FOR LACK OF DATA

Lin Dongfang，Song Yingchun，Du Kun and Xiao Qinqin
(School of Geosciences and Info-Physics，Central South University，Changsha 410083)

Restricted by the observing conditions，some important data of GPS elevation fitting are often missed or not enough.Using conventional methods will reduce the fitting accuracy of the quasi-geoid.By use of the EM algorithm with the addition of“potential data”，conditional expectation of the missing data，it can effectively improve the accuracy of GPS elevation fitting for lack of data.

missing data;EM algorithm;elevation fitting;related plane fitting;quadratic surface fitting

1671-5942(2012)02-0139-04

2011-11-12

國(guó)家自然科學(xué)基金(40874005);教育部博士點(diǎn)基金(200805331086)

林東方，男，1986年生，碩士研究生，主要研究方向?yàn)?現(xiàn)代測(cè)量數(shù)據(jù)處理.E-mail:lindongfang223@163.com

P207

A