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

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

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

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

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

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

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

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

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

1廣播星歷法計(jì)算衛(wèi)星軌道周期的統(tǒng)計(jì)規(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)，故不予考慮),以廣播星歷法計(jì)算衛(wèi)星重復(fù)出現(xiàn)的提前時(shí)間。

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

以236～254 s作為提前時(shí)間的合理波動(dòng)范圍，檢驗(yàn)廣播星歷法計(jì)算提前時(shí)間的正確率(表1)。

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

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

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

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

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

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

圖1　PRN9號(hào)衛(wèi)星的日均提前時(shí)間散點(diǎn)圖Fig.1　Scatter of daily average preset time of PRN9

圖2　PRN17號(hào)衛(wèi)星的日均提前時(shí)間散點(diǎn)圖Fig.2　Scatter of daily average preset time of PRN17

圖3　PRN25號(hào)衛(wèi)星的日均提前時(shí)間散點(diǎn)圖Fig.3　Scatter of daily average preset time of PRN25

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

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

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

(1)

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

圖4　提前時(shí)間及傅里葉擬合曲線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)星運(yùn)行周期提前時(shí)間的擬合函數(shù)式，提前時(shí)間具有一個(gè)周期為3 057.5 d的重復(fù)。

3實(shí)例分析

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

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

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

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

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

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

4結(jié)語

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

參考文獻(xiàn)

[1]伊廷華,張永恒,李宏男,等.多路徑效應(yīng)在GPS結(jié)構(gòu)健康監(jiān)測(cè)中的研究進(jìn)展[J].振動(dòng)與沖擊，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高層建筑物常荷載振動(dòng)測(cè)試的小波分析[J].振動(dòng)與沖擊,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ù)性試驗(yàn)研究[J].測(cè)繪科學(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)軍，肖根如.恒星日濾波的修正以及對(duì)高頻GPS定位的影響研究[J].武漢大學(xué)學(xué)報(bào):信息科學(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測(cè)站坐標(biāo)殘差序列的傅里葉模型建立與分析[J].測(cè)繪技術(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]劉維仲,牟景文,戴積成,等.傅里葉級(jí)數(shù)的曲線仿真應(yīng)用及算法[J].天津大學(xué)學(xué)報(bào)，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].大地測(cè)量與地球動(dòng)力學(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].測(cè)繪通報(bào),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)介：段海東，碩士生，主要從事精密工程測(cè)量與變形監(jiān)測(cè)研究，E-mail:466481290@qq.com。

DOI：10.14075/j.jgg.2016.07.011

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

中圖分類號(hào)：P228

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

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