林東方 宋迎春 杜 琨 肖琴琴

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

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

林東方 宋迎春 杜 琨 肖琴琴

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

GPS高程擬合常因觀測條件限制，部分重要擬合點的數(shù)據(jù)無法獲取，致擬合數(shù)據(jù)缺失或數(shù)據(jù)量不充足。采用常規(guī)方法擬合，會降低似大地水準(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ī)定采用正常高高程系統(tǒng)作為我國高程的統(tǒng)一系統(tǒng)［1］。GPS高程觀測數(shù)據(jù)為大地高系統(tǒng)下的高程數(shù)據(jù)，由于基準(zhǔn)面的不同，大地高系統(tǒng)與正常高系統(tǒng)存在一定的差距稱為高程異常。GPS高程數(shù)據(jù)需要進行高程擬合，將大地高高程轉(zhuǎn)換為正常高高程。

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

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

2.1 EM算法

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

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

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

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

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

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

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

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

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

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

式中，

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

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

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

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

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

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

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

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

3 算例分析

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

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

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

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

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

完全數(shù)據(jù)下采用二次曲面擬合法，以GP01、GP02、GP04、GP05、GP06和 GP07為擬合點，計算結(jié)果與平面相關(guān)擬合法和EM算法比較結(jié)果見表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)計(單位：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高程測量的精度，高程擬合時擬合點分布均勻且數(shù)據(jù)充足，才能較好地擬合出似大地水準(zhǔn)面。而在擬合點觀測中難免會遇到部分重要的擬合點因觀測條件限制而無法觀測，造成數(shù)據(jù)缺失。采用不完全的擬合點數(shù)據(jù)擬合似大地水準(zhǔn)面將降低似大地水準(zhǔn)面的擬合準(zhǔn)確度。應(yīng)用EM算法，將觀測數(shù)據(jù)的誤差分布信息引入數(shù)據(jù)計算中，引入缺失數(shù)據(jù)的條件期望，添加了“潛在數(shù)據(jù)”，較好地改善了缺失數(shù)據(jù)下的擬合似大地水準(zhǔn)面精度，與完全數(shù)據(jù)下的擬合相比，提高了外符合精度。因此，在缺失數(shù)據(jù)的情況下，應(yīng)用EM算法進行GPS高程擬合是一種行之有效的方法。

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

2 魏立峰，何建國.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 姚吉利，曲國慶，劉科利.擬合點分布與GPS水準(zhǔn)面擬合精度關(guān)系研究［J］.大地測量與地球動力學(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］.測繪科學(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算法和單幅雷達圖像陰影的控制點坡度校正［J］.武漢大學(xué)學(xué)報(信息科學(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］.南開大學(xué)學(xué)報(自然科學(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é)測繪學(xué)院測量平差學(xué)科組.誤差理論與測量平差基礎(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 楊元喜，張麗萍.中國大地測量數(shù)據(jù)處理60年重要進展［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

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

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

P207

A