缺失數據下基于EM算法的GPS高程擬合*

林東方 宋迎春 杜 琨 肖琴琴

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

缺失數據下基于EM算法的GPS高程擬合*

林東方 宋迎春 杜 琨 肖琴琴

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

GPS高程擬合常因觀測條件限制，部分重要擬合點的數據無法獲取，致擬合數據缺失或數據量不充足。采用常規方法擬合，會降低似大地水準面的擬合精度。應用EM算法，添加了有益于高程擬合的“潛在數據”，即缺失數據的條件期望，有效提高了缺失數據下的GPS高程擬合精度。

缺失數據;EM算法;高程擬合;平面相關擬合;二次曲面擬合

1 引言

高程系統包括大地高、正高和正常高系統。大地高系統的基準面為參考橢球面，正高系統的基準面為大地水準面，而正常高系統以似大地水準面為基準面。我國規定采用正常高高程系統作為我國高程的統一系統［1］。GPS高程觀測數據為大地高系統下的高程數據，由于基準面的不同，大地高系統與正常高系統存在一定的差距稱為高程異常。GPS高程數據需要進行高程擬合，將大地高高程轉換為正常高高程。

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

2 EM算法及其在缺失數據下的GPS高程擬合方法

2.1 EM算法

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

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

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

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

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

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

2.2 GPS高程擬合數學模型

設GPS高程擬合數學模型為

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

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

式中，

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

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

2.3 缺失數據下基于EM算法的GPS高程擬合方法

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

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

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

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

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

3 算例分析

以杭州灣跨海大橋的GPS高程擬合數據為例，GPS控制網點位分布見圖1。

圖1 GPS控制網點位略圖Fig.1 GPS control points network

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

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

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

完全數據下采用二次曲面擬合法，以GP01、GP02、GP04、GP05、GP06和 GP07為擬合點，計算結果與平面相關擬合法和EM算法比較結果見表3。

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

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

表3 擬合結果統計(單位：cm)Tab.3 Statistics of fitting results(unit：cm)

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

4 結論

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

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

2 魏立峰，何建國.GPS高程擬合似大地水準面的方法［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水準面擬合精度關系研究［J］.大地測量與地球動力學，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 姚喜，欒學科，王志博.GPS水準擬合方法的研究［J］.測繪科學，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］.軟件導刊，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］.武漢大學學報(信息科學版)，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］.南開大學學報(自然科學版)，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 武漢大學測繪學院測量平差學科組.誤差理論與測量平差基礎［M］.武漢:武漢大學出版社，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 楊元喜，張麗萍.中國大地測量數據處理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

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

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

P207

A