GPS動態變形分析中的傅里葉級數恒星日濾波

段海東　張獻州　馬　龍　李祖來張　拯　張正國　羅　奕

1　西南交通大學地球科學與環境工程學院，成都市犀安路999號，611756 2　高速鐵路運營安全空間信息技術國家地方聯合工程實驗室，成都市犀安路999號，611756 3　成都鐵路局成都高鐵工務段，成都市荷花池路3號，610000

GPS動態變形分析中的傅里葉級數恒星日濾波

段海東1,2張獻州1,2馬龍1,2李祖來1,2張拯1,2張正國3羅奕3

1西南交通大學地球科學與環境工程學院，成都市犀安路999號，611756 2高速鐵路運營安全空間信息技術國家地方聯合工程實驗室，成都市犀安路999號，611756 3成都鐵路局成都高鐵工務段，成都市荷花池路3號，610000

摘要：針對恒星日濾波中計算GPS衛星軌道運行周期提前時間的問題，提出一種基于傅里葉級數曲線擬合的提前時間估計方法。根據跨度達8 526 d的廣播星歷計算每天中所有可用衛星運行周期相對于前1 d的平均提前時間，將此平均提前時間作為采樣點，以8階三角函數形式展開的傅里葉級數進行曲線擬合，將衛星運行周期的提前時間擬合成關于時間的單一變量函數，免去了繁瑣的計算步驟。實驗結果表明,傅里葉級數恒星日濾波法縮小了監測點N、E、U三個方向坐標殘差序列的均方根值，有效消除了GPS動態變形監測分析中的多路徑影響，提高了復雜工程環境下動態變形監測中變形信息提取的能力。

關鍵詞：GPS動態變形監測;恒星日濾波;廣播星歷法;傅里葉級數;曲線擬合

近年來，高頻甚至超高頻GPS接收機不斷出現，提升了復雜工程環境下變形監測的實時動態監測能力。多路徑影響是制約GPS精密動態變形監測技術發展的關鍵因素之一，干擾了動態變形信息的提取過程[1-2]。恒星日濾波利用GPS衛星運行的重復性特點，將低通濾波之后的坐標殘差經過時間平移進行差分運算，消除兩天之間具有相關性的誤差，包括多路徑誤差[3-6]。開普勒軌道條件下，衛星重復出現在同一位置的時間應該比前1 d提前236 s，但是美國國防部并未將衛星的軌道周期約束為恒星周期，而是把地面跟蹤站設定為固定站，所以衛星每天的實際軌道并不完全一致[6]。Choi等[1,7]在2004年通過對衛星軌道的分析證明，GPS衛星幾何重復周期略小于恒星日，即衛星較前1 d重復出現在同一位置的時間提前量略大于236 s，約為245 s。本文采用開普勒第三運動定律計算GPS衛星運動周期及其變化情況，將此變化通過三角函數形式的傅里葉級數進行曲線擬合[8-9]，以此對恒星日濾波進行修正，然后將修正后的恒星日濾波應用于GPS動態變形監測中，以便正確提取變形體的變形信息。

1廣播星歷法計算衛星軌道周期的統計規律

本文選取1992-02-12～2015-06-16跨度達8 526 d的GPS廣播星歷，覆蓋了40顆衛星，包括32顆常用衛星(PRN1～PRN32)和8顆非常用衛星(PRN47、PRN65、PRN68、PRN69、PRN73、PRN76、PRN79、PRN80，因這8顆非常用衛星只在1992-08-04、08出現，故不予考慮),以廣播星歷法計算衛星重復出現的提前時間。

1.1廣播星歷法計算衛星重復出現提前時間的正確率

以236～254 s作為提前時間的合理波動范圍，檢驗廣播星歷法計算提前時間的正確率(表1)。

由表1可見，在歷元數目上，PRN3號衛星最多，而PRN32號衛星最少，分別為110 164個和39 765個。相應地，在廣播星歷法計算提前時間的正確數目上，PRN3號衛星也最多，而PRN32號衛星最少，分別為109 883個和37 081。但是，在計算正確率上，PRN19號衛星最高，而PRN29號衛星最低，分別為99.86%和83.63%。廣播星歷法計算32顆衛星重復出現提前時間的正確率平均值達到95.37%。

表1　廣播星歷法計算衛星重復出現提前時間的正確率統計

廣播星歷法計算衛星的重復出現提前時間存在一定的錯誤率，引起這一錯誤的原因包括地球非球形引力場攝動力、潮汐攝動力(包括地球固體潮、海潮、極移潮、負荷潮等)、日月及其他行星引力、太陽光壓、地球反照壓、地球輻射壓、大氣阻力、大氣浮力、衛星本身熱輻射壓、姿態調整及變軌機動力、衛星退役等。當出現這種錯誤時需要使用其他方法進行計算，如最大互相關系數法、最小均方根誤差值法等。

1.2衛星提前時間的分布

以波動在236 s～254 s范圍內的每日平均提前時間作為豎軸、日期(1992-01-01～2016-01-01)為橫軸作散點圖。限于篇幅，只展示了PRN9、PRN17、PRN25共3顆衛星運行的日均提前時間散點分布情況(圖1～3)。

由圖1～3可見，廣播星歷法計算得到的32顆衛星的提前時間量不盡相同，單顆衛星在不同時期的提前時間也不相同。但是也不難看出：1)在大于1 a的時間跨度內，單顆衛星的提前時間在波動范圍內具有一定的周期重復性，但每顆衛星的重復周期不盡相同，單顆衛星在不同時期內的重復周期也不盡相同；2)在1 a以下的時間跨度內，每顆衛星的提前時間具有一定的趨勢性，但每顆衛星的趨勢不盡相同，單顆衛星在不同時期內的趨勢也可能不一樣。

圖1　PRN9號衛星的日均提前時間散點圖Fig.1　Scatter of daily average preset time of PRN9

圖2　PRN17號衛星的日均提前時間散點圖Fig.2　Scatter of daily average preset time of PRN17

圖3　PRN25號衛星的日均提前時間散點圖Fig.3　Scatter of daily average preset time of PRN25

根據在短時間內衛星運行周期提前時間的趨勢性，如遇到上文所述的廣播星歷法計算衛星運行周期提前時間出現錯誤時，可用插值法確定其正確的提前時間(簡單的線性內插基本能夠滿足要求)，也可以在短時間內對衛星運行周期的提前時間作出預測。

2衛星軌道周期提前時間曲線的函數擬合

將同1 d內所有可見衛星的日均提前時間進行均值運算，即為當日用于削弱多路徑影響的坐標殘差序列提前時間，圖4中較細曲線即為日期-提前時間曲線。1992-02-12～2015-06-16每恒星日提前時間基本在244～247 s，提前時間隨時間波動，呈現出一定的周期重復性。將提前時間用傅里葉1～8階三角函數展開式進行最小二乘曲線擬合[8-9]。經計算，擬合殘差在8階時已經趨向收斂。同時，為了不增加過多的計算量，選取擬合殘差最小的8階傅里葉展開式作為擬合函數：

(1)

將1992-02-12作為時間起算原點，每恒星日作為1個采樣間隔，每日的提前時間作為采樣值，用傅里葉8階展開式進行曲線擬合，在95%的置信水平下求得曲線擬合函數的參數，見表2。

圖4　提前時間及傅里葉擬合曲線Fig.4　Preset time and Fourier series curve fitting

參數參數值參數參數值參數參數值a0245.4000ω0.0021a1-0.1045b1-0.0047a2-0.0212b2-0.0238a3-0.0817b3-0.1551a4-0.0461b40.1003a5-0.0521b50.2606a6-0.1182b60.1414a70.0063b7-0.1518a8-0.1097b80.1532

擬合均方根誤差RMSE= 0.496 5 。將擬合參數代入式(1)，得到衛星運行周期提前時間的擬合函數式，提前時間具有一個周期為3 057.5 d的重復。

3實例分析

某隧道地處天山山脈中段南部、塔里木盆地北緣，附近約3 km 范圍內荒無人煙，幾乎無植被，僅有少許耐旱、耐鹽雜草零星散布。隧址區GPS動態變形監測在隧道的兩個洞口外各布設2個監測點，進口端布設GPS1、GPS3點，出口端布設GPS2、GPS4點，GPS3距離洞口較遠，GPS3與GPS1、GPS2、GPS4分別相距約200 m、2.3 km、2.4 km。以GPS3作為固定站進行單歷元差分解算，短距離基線有利于削弱GPS衛星信號傳播路徑中的相關性誤差。

將2014-09-19和2014-09-20兩天的13:50～17:15高頻(1 s采樣率)觀測數據采用GAMIT/GLOBK軟件的Track模塊進行單歷元解算，差分后的坐標殘差序列先用db8小波低通濾波去除高頻觀測噪聲[10]，之后利用上述的傅里葉擬合曲線估計當日衛星運行周期的提前時間，最后用恒星日濾波法削弱多路徑誤差[11]。傅里葉擬合曲線估計得到2014-09-20衛星運行周期的提前時間為245.412 s，將當日的坐標殘差序列平移245個歷元，與19日的坐標殘差序列作差，即可用恒星日濾波法削弱當日的多路徑誤差。由于篇幅原因，本文只對GPS1點坐標殘差中的多路徑誤差削弱情況作出圖形分析。由表3可見，N、E、U3個方向的坐標殘差均方根值由恒星日濾波前的3.79 mm、3.14 mm、8.80 mm分別縮小到2.82 mm、2.27 mm、6.53 mm，分別縮小25.47%、27.86%、25.80%。圖5～7為GPS1點在N、E、U3個方向在恒星日濾波前后的坐標殘差對比。恒星日濾波后的坐標殘差削弱了多路徑影響，坐標殘差序列趨向更加平穩、無規律的小幅波動，主要以隨機噪聲和變形體的變形信息為主。

表3　傅里葉級數恒星日濾波前后坐標殘差序列的均方根值統計

圖5　恒星日濾波后GPS1點N方向坐標殘差Fig.5　N direction residuals before and after sidereal filtering of GPS1 station

圖6　恒星日濾波后GPS1點E方向坐標殘差Fig.6　E direction residuals before and after sidereal filtering of GPS1 station

圖7　恒星日濾波后GPS1點U方向坐標殘差Fig.7　U direction residuals before and after sidereal filtering of GPS1 station

4結語

本文基于三角函數展開形式的傅里葉級數對恒星日濾波方法進行改進，采用超過23 a時長跨度的海量廣播星歷數據，由廣播星歷法計算GPS衛星運行周期的每日平均提前時間，將此提前時間作為擬合采樣點，用8階傅里葉三角函數展開式將GPS衛星運行周期的提前時間擬合成關于時間的函數。在實際應用中，只需基于本文提出的擬合函數，通過日期計算衛星運行周期的提前時間，即可免去繁雜的提前時間計算過程，用改進后的恒星日濾波法消除多路徑影響?；谌呛瘮嫡归_式的傅里葉級數的恒星日改進方法，可以縮小點位N、E、U3個方向的坐標殘差序列均方根值，消除多路徑影響之后的坐標殘差序列主要以隨機噪聲和變形體的變形信息為主。本文提出的基于三角函數展開式的傅里葉級數的恒星日濾波方法，其計算步驟簡單并能有效消除多路徑影響，有利于提高工程動態變形監測中變形信息的提取能力，可廣泛應用于復雜環境下的工程動態變形監測。

參考文獻

[1]伊廷華,張永恒,李宏男,等.多路徑效應在GPS結構健康監測中的研究進展[J].振動與沖擊，2009,28(9):102-108(Yi Tinghua, Zhang Yongheng, Li Hongnan, et al. Study Progress of Multipath Effect in Structural Health Monitoring Based on GPS[J].Journal of Vibration and Shock，2009,28(9):102-108)

[2]黃丁發,陳永奇,丁曉利,等.GPS高層建筑物常荷載振動測試的小波分析[J].振動與沖擊,2001,20(1):12-17( Huang Dingfa, Chen Yongqi, Ding Xiaoli, et al. Wavelet-Based Analysis Technique for the Monitoring of Tall Structure under Normal Loading Using GPS[J].Journal of Vibration and Shock,2001,20(1):12-17)

[3]Choi K, Bilich A, Larson K M, et al. Modified Sidereal Filtering: Implications for High-Rate GPS Positioning[J]. Journal of the Korean Orthopaedic Association, 1999, 15(2):178-198

[4]鐘萍，袁林果，丁曉利.GPS多路徑重復性試驗研究[J].測繪科學,2010,35(2):43-45(Zhong Ping, Yuan Linguo, Ding Xiaoli. Experimental Analysis of GPS Multipath Day-to-Day Repeatability[J]. Science of Surveying and Mapping, 2010,35(2):43-45)

[5]Schwahn W, S?hne W. Modified Sidereal Filtering-Tool for the Analysis of High-Rate GPS Coordinate Time Series[J]. International Association of Geodesy Symposia, 2009, 134:219-224

[6]殷海濤，甘衛軍，肖根如.恒星日濾波的修正以及對高頻GPS定位的影響研究[J].武漢大學學報:信息科學版,2011, 36(5):609-616( Yin Haitao，Gan Weijun，Xiao Genru. Modified Sidereal Filter and Its Effect on High-Rate GPS Positioning[J]. Geomatics and Information Science of Wuhan University,2011, 36(5):609-616)

[7]Larson K M, Bilich A, Axelrad P. Improving the Precision of High-Rate GPS[J]. Journal of Geophysical Research,2007,32(B05 422):1-11

[8]田亮,趙淑芬,趙亞平,等. 基于GPS測站坐標殘差序列的傅里葉模型建立與分析[J].測繪技術裝備,2012,3(14):11-13(Tian Liang, Zhao Shufen, Zhao Yaping, et al. The Establishment and Analysis of the Fourier Model Based on the GPS Measuring Station Coordinate Residual Series[J]. Geomatics Technology and Equipment,2012,3(14):11-13)

[9]劉維仲,牟景文,戴積成,等.傅里葉級數的曲線仿真應用及算法[J].天津大學學報，1996,29(1):108-112(Liu Weizhong, Mou Jingwen, Dai Jicheng,et al. The Application and Algorithm of Curve Simulation Based on Fourier Series[J]. Journal of Tianjin University, 1996,29(1):108-112)

[10]蘇小寧,孟國杰,胡新康,等.高頻GPS單歷元定位精度及噪聲特征[J].大地測量與地球動力學，2014,34(2):150-154(Su Xiaoning, Meng Guojie, Hu Xinkang, et al. Precision and Noise Characteristics of High-Frequency GPS Epoch-by-Epoch Positioning[J]. Journal of Geodesy and Geodynamics，2014,34(2):150-154)

[11]孫常建,楊曉超.GPS多路徑效應規律研究[J].測繪通報,2006(11):12-15(Sun Changjian, Yang Xiaochao. The Research of the Rules of GPS Multipath Effect[J]. Bulletin of Surveying and Mapping,2006(11):12-15)

Foundation support:Program for Changjiang Scholars and Innovative Research Team in University,No. IRT13092;State-Province Joint Engineering Laboratory of Spatial Information Technology for High-Speed Railway Safety;Science and Technology Research and Development Program of China Railway,No.2012G009-C.

About the first author:DUAN Haidong, postgraduate, majors in precise engineering surveying and deformation monitoring,E-mail:466481290@qq.com.

Fourier Series Sidereal Filtering in GPS Dynamic Deformation Analysis

DUANHaidong1,2ZHANGXianzhou1,2MALong1,2LIZulai1,2ZHANGZheng1,2ZHANGZhengguo3LUOYi3

1Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, 999 Xi’an Road, Chengdu 611756, China 2State-Province Joint Engineering Laboratory of Spatial Information Technology for High-Speed Railway Safety,999 Xi’an Road, Chengdu 611756, China 3High-Speed Railway Track Maintenance Division of Chengdu Railway Bureau,3 Hehuachi Road, Chengdu 610000, China

Abstract:High-rate GPS receivers are popular in real-time dynamic monitoring, but the multi-path effect is one of the key obstacles that restricts its development. However, a sidereal filter can eliminate this effect and improve positioning accuracy. The sidereal filter requires the preset time of each satellite compared to the day before. In this research, 8 526 days broadcast ephemerides are used to calculate the preset time. The Fourier series of 8 orders trigonometric functions is used to fit the series of the preset time. It shows a simple and effective way to estimate the preset time, instead of a more complex calculation. Using this method, the root mean square values of the residual coordinate sequences of N, E and U directions are reduced, and the ability of extracting deformation information is improved.This method therefore can be applied in engineering dynamic deformation monitoring.

Key words:GPS dynamic deformation monitoring; sidereal filtering; method based on broadcast ephemeris; Fourier series; curve fitting

收稿日期：2015-07-25

第一作者簡介：段海東，碩士生，主要從事精密工程測量與變形監測研究，E-mail:466481290@qq.com。

DOI：10.14075/j.jgg.2016.07.011

文章編號：1671-5942(2016)07-0609-04

中圖分類號：P228

文獻標識碼：A

項目來源：長江學者和創新團隊發展計劃( IRT13092);高速鐵路運營安全空間信息技術國家地方聯合工程實驗室項目;原鐵道部科技研究開發計劃(2012G009-C)。