段海東　張獻(xiàn)州　馬　龍　李祖來張　拯　張正國　羅　奕

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

GPS動態(tài)變形分析中的傅里葉級數(shù)恒星日濾波

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

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

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

關(guān)鍵詞：GPS動態(tài)變形監(jiān)測;恒星日濾波;廣播星歷法;傅里葉級數(shù);曲線擬合

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

1廣播星歷法計算衛(wèi)星軌道周期的統(tǒng)計規(guī)律

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

1.1廣播星歷法計算衛(wèi)星重復(fù)出現(xiàn)提前時間的正確率

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

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

表1　廣播星歷法計算衛(wèi)星重復(fù)出現(xiàn)提前時間的正確率統(tǒng)計

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

1.2衛(wèi)星提前時間的分布

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

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

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

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

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

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

2衛(wèi)星軌道周期提前時間曲線的函數(shù)擬合

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

(1)

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

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

參數(shù)參數(shù)值參數(shù)參數(shù)值參數(shù)參數(shù)值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 。將擬合參數(shù)代入式(1)，得到衛(wèi)星運行周期提前時間的擬合函數(shù)式，提前時間具有一個周期為3 057.5 d的重復(fù)。

3實例分析

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

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

表3　傅里葉級數(shù)恒星日濾波前后坐標(biāo)殘差序列的均方根值統(tǒng)計

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

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

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

4結(jié)語

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

參考文獻(xiàn)

[1]伊廷華,張永恒,李宏男,等.多路徑效應(yīng)在GPS結(jié)構(gòu)健康監(jiān)測中的研究進(jìn)展[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]黃丁發(fā),陳永奇,丁曉利,等.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多路徑重復(fù)性試驗研究[J].測繪科學(xué),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]殷海濤，甘衛(wèi)軍，肖根如.恒星日濾波的修正以及對高頻GPS定位的影響研究[J].武漢大學(xué)學(xué)報:信息科學(xué)版,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測站坐標(biāo)殘差序列的傅里葉模型建立與分析[J].測繪技術(shù)裝備,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]劉維仲,牟景文,戴積成,等.傅里葉級數(shù)的曲線仿真應(yīng)用及算法[J].天津大學(xué)學(xué)報，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].大地測量與地球動力學(xué)，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多路徑效應(yīng)規(guī)律研究[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

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

DOI：10.14075/j.jgg.2016.07.011

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

中圖分類號：P228

文獻(xiàn)標(biāo)識碼：A

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