馬 俊，周曉慧，朱兆涵

(武漢大學測繪學院，湖北 武漢 430079)

聯(lián)合區(qū)域疊加濾波法與小波變換去除GPS站坐標時間序列噪聲

馬 俊，周曉慧，朱兆涵

(武漢大學測繪學院，湖北 武漢 430079)

GPS站坐標時間序列中存在的周期性與非周期性誤差嚴重影響了對測站運動特征的分析及其非線性變化的物理機制解釋。因此，為削弱噪聲的影響，本文首先利用區(qū)域疊加濾波法去除了南加利福尼亞地區(qū)16個測站時間序列的共模誤差，以此削弱時間序列中存在的包括周年和半周年誤差在內(nèi)的周期性誤差。為去除濾波后殘留的噪聲，對濾波后的信號進行靜態(tài)離散小波變換，提取了周期為半周年以上的信號。結(jié)果表明，聯(lián)合區(qū)域疊加濾波法與小波變換對GPS站坐標時間序列進行處理，既能夠削弱周期性誤差對信號的影響，又能較好地提取測站的非線性運動信號。

GPS站坐標時間序列；共模誤差；小波變換；信噪分離

積累多年的GPS站坐標時間序列可以反映出測站的線性與非線性的運動特征，為地殼形變監(jiān)測及一些地球物理現(xiàn)象的解釋提供了寶貴的基礎數(shù)據(jù)。然而有研究表明，GPS坐標時間序列中除含有非周期性誤差外，還包含周期性誤差[1]，稱之為虛假的周期性信號。這些周期性誤差是由未模型化的系統(tǒng)性誤差(包括由地球物理因素引起的誤差、系統(tǒng)誤差、GPS數(shù)據(jù)處理中采用的計算方法的誤差)傳播到GPS坐標時間序列中生成，其周期范圍為兩個星期至一年[2-3]。這些周期性及非周期性的誤差嚴重影響著對測站運動狀態(tài)的分析，且會導致對引起測站非線性變化的地球物理因素的錯誤解釋[3-4]，因此有必要采取措施削弱其對GPS站坐標時間序列分析的影響。

由于共模誤差(common mode error,CME)被認為是GPS站坐標時間序列誤差的主要來源之一[5]，因此許多研究通過去除CME來削弱噪聲的強度，提高坐標時間序列的信噪比[6-9]。此外，由于小波變換可以通過對信號的高頻分解系數(shù)進行閾值量化去除噪聲，因此也被用于GPS站坐標時間序列分析，具有較好的去噪效果[12-13]。以上方法有效地削弱了噪聲對GPS坐標時間序列的影響，但仍存在一些不足。首先，小波變換無法分辨和剔除GPS坐標時間序列中包括周年及半周年在內(nèi)的長周期性誤差。其次，去除CME僅可以剔除區(qū)域GPS站坐標時間序列所受的共同的誤差影響，但在每個測站剩余的時間序列中依然存在較強的噪聲。由于CME很可能是由未模型化或未完善的衛(wèi)星軌道誤差、參考框架誤差、地球定向參數(shù)誤差造成的，且未完善的接收機和衛(wèi)星天線相位中心誤差以及大尺度的大氣影響也是共模誤差的潛在來源[5]，因此CME與周期性誤差具有共同的誤差源。根據(jù)以上分析，本文結(jié)合這兩種方法的優(yōu)點，首先通過區(qū)域疊加濾波法去除GPS站坐標時間序列中的CME，以此削弱周期性誤差；然后對濾波后的時間序列進行小波分解，并提取周年及半周年信號，最后獲得具有較高信噪聲比的時間序列。

1 時間序列分析方法

1.1 區(qū)域疊加濾波法

提取GPS站坐標時間序列中的共模誤差可以去除測站相同分量時間序列之間的相關(guān)性，提高時間序列的信噪比[7]。由于本文所選擇的GPS測站之間的距離不超過500 km，因此利用區(qū)域疊加濾波法去除共模誤差[14]。該算法假設共模誤差在某一區(qū)域是均勻分布的，將單日解的誤差作為權(quán)重因子，利用簡單的算數(shù)平均計算共模誤差[7]

(1)

式中，εi為第i個歷元的共模誤差值；v為殘差坐標時間序列；σi,S為單日坐標解的誤差；S為參與計算共模誤差的站臺數(shù)。其中，v通過去除坐標時間序列的線性趨勢以及周年和半周年諧波獲得。

1.2 靜態(tài)離散小波變換

小波變換是一種時頻局部化方法，具有多分辨率分析的特點，可以探測信號中的瞬態(tài)成分及頻率成分，因此被稱為數(shù)學顯微鏡。小波變換主要包括連續(xù)小波變換和離散小波變換。在實際問題的數(shù)值計算中常采用離散形式，即離散小波變換(discrete wavelet transform,DWT)。小波分析相關(guān)公式如下

(2)

式中,φ(t)為小波母函數(shù)，對其進行平移和伸縮得到小波基函數(shù)集φa,b(t);a為尺度因子;b為平移因子;R為實數(shù)集。通常取

(3)

代入式(2)可得

(4)

此時離散小波變換公式為

(5)

式中，fDWT(a,b)為連續(xù)小波變換系數(shù)；*表示共軛。一維離散小波變換實現(xiàn)的算法一般是Mallat算法，即先對較大的信號進行小波變換，獲得近似系數(shù)和細節(jié)系數(shù)，再選取近似系數(shù)在原尺度的1/2尺度上進行小波變換。其中每次小波分解都進行下采樣，只保留偶數(shù)項，因此每次離散小波變換所獲得的小波系數(shù)的長度為上一次小波變換系數(shù)長度的一半，且丟棄了奇數(shù)項所含的時移信息。根據(jù)采樣定理，為獲得包含周年及半周年信號的小波系數(shù)，本文需對GPS站坐標時間序列至少分解為8層。因此，若利用離散小波變換，坐標時間序列的長度至少為28。為解除離散小波變換層數(shù)對坐標時間序列長度的限制，獲得包含全部時移信息的小波系數(shù)，本文采用靜態(tài)離散小波變換(static discrete wavelet transform,SWT)。在對信號的分解過程中，SWT不對分解的系數(shù)進行下采樣，因此每次分解得到的系數(shù)長度都與原信號相同。在信號重建過程中，近似系數(shù)和細節(jié)系數(shù)分別作用重建低通和高通濾波后，直接相加就可以得到重建的上一層次的近似系數(shù)(或原始信號)。

2 實例分析

2.1 GPS數(shù)據(jù)

本文采用南加利福尼亞地區(qū)16個GPS測站2010—2016年的坐標時間序列(ftp:∥sideshow.jpl.nasa.gov/pub/usrs/mbh/point/)。這些測站是PBO(Plate Boundary Observation)核心網(wǎng)絡的組成部分，具體位置如圖1所示。其坐標時間序列由JPL(Jet Propulsion Laboratory)利用PPP數(shù)據(jù)使用GIPSY軟件生成的日解組成。JPL給出了剔除粗差后干凈的時間序列，確定了階躍并估計了測站的運動速度。在分析GPS站坐標時間序列之前，筆者根據(jù)JPL給出的階躍估值對剔除粗差的時間序列進行改正。在日解處理過程中采用的改正模型如下：地球、太陽、月亮及其他行星的引力模型；DE421行星星歷模型；IAU06歲差和章動模型；IERS21010潮汐模型；FES2004海洋負載模型；IGS衛(wèi)星以及接收機天線相位中心模型。

圖1 GPS測站分布

2.2 試驗結(jié)果

本文首先采用區(qū)域疊加濾波法去除GPS站坐標時間序列的共模誤差。由于垂直方向時間序列的誤差強度最大，受篇幅所限，本文僅列出BBRY、BSRY，以及TABL站垂向時間序列的分析結(jié)果。圖2為這3個測站垂直方向原始坐標時間序列以及去除共模誤差后時間序列的功率譜。

分析圖2可知，去除共模誤差之前，BBRY與BSRY站垂向時間序列中皆含有明顯的周年及半周年信號，而TABL站僅含有能量較強的半周年信號。此外，從圖2中還可以看出，原始時間序列中還含有眾多能量弱于周年半周年的周期性誤差。這些誤差很可能是由地球物理效應及日解中的GPS系統(tǒng)誤差造成的虛假周期信號。研究表明，未模型化的半周日和周日海潮負荷位移信號在GPS站坐標時間序列中產(chǎn)生了兩星期、半年和周年周期信號[14]。理論計算表明，這種虛假的“潮汐分量”會對GPS時間序列產(chǎn)生垂向0.5～1 mm、水平向0.1～0.2 mm的偏差。同時GPS衛(wèi)星的重復軌道周期也會引起混疊效應，在GPS站坐標時間序列中產(chǎn)生虛假的周期信號[15]。以上周期性誤差會導致對地球物理信息的錯誤解釋。

對比濾波前后信號的功率譜可知，濾波后BBRY與BSRY站垂向時間序列的周年和半周年信號的功率均有較大程度的下降，且BBRY站頻率為3 cpy的周期信號基本上被消除。TABL站的半周年信號功率幾乎減少一半，但其周年信號功率卻大幅度增加。此外，頻率大于2 cpy的周期性誤差的功率也被大幅度削弱。由于共模誤差主要包含與GPS相關(guān)的技術(shù)誤差(如高階電離層延遲、對流層延遲)，以及包含少量的由其他地球物理因素(如海洋潮汐、大氣潮汐及大氣負載)造成的誤差，且以上因素都能夠?qū)PS站坐標時間序列中的周年及半周年信號造成不同程度的影響[16-17]。因此，共模誤差序列中存在年周期信號、半年周期信號和其他周期性信號。通過去除共模誤差削弱了GPS站坐標時間序列中包括周年及半周年誤差在內(nèi)的周期性誤差的強度，致使GPS站坐標時間序列的周年和半周年信號的功率降低，還原了被掩蓋的真實的周年及半周年信號的功率，提高了時間序列的信噪比。另外，從圖2還可以看出，濾波前后這3個測站半周年的功率皆大于周年，尤其BSRY站周年半周年的功率濾波后相差最大。根據(jù)JPL提供的這3個測站在垂直方向周年與半周年的正余弦系數(shù)可知，BSRY與TABL站在垂直方向的半周年振幅大于周年振幅，而BBRY站的周年振幅大于半周年，因此BSRY和TABL站垂向時間序列中半周年信號的功率也大于周年信號。對于BBRY站而言，在1998—2016年間其垂向時間序列中周年信號的振幅大于半周年信號，而在2010—2016年內(nèi)半周年信號的振幅很可能大于周年信號，因此圖2中BBRY站半周年信號的功率大于周年信號。從圖2還可以看出，以上測站濾波后垂向時間序列的功率譜中仍存在未剔除干凈的高頻的周期性誤差，如在TABL站的功率譜中，頻率為3、10 cpy處仍存在功率為100 db的周期性誤差。

圖2 垂向原始時間序列去除共模誤差前后的功率譜

周年、半周年信號是當今IGS全球基準站最主要的信號特征，而濾波后的時間序列中仍存在未剔除干凈的較強的噪聲。因此本文通過小波變換，將濾波后的時間序列中信號的不同頻率成分分解到互不重疊的頻帶上，以此將周期為半周年以上的信號和噪聲分開。首先剔除濾波后時間序列的性趨勢，然后選用小波函數(shù)“coif5”，利用靜態(tài)離散小波變換將其分解為8層。圖3—圖5為以上測站垂向時間序列各層的小波分解信號及其功率譜圖。其中，d1—d8為對信號進行小波分解后得到的高頻信號，包含了信號的細節(jié)信息；a8為第8層分解得到的低頻信號，反映了信號的整體非線性運動特征。根據(jù)采樣定理，d1—d8以及a8信號的主要頻帶范圍分別是91.25～182.5、45.625～91.25、22.81～45.625、11.39～22.81、5.69～11.39、2.85～5.69、1.42～2.85、0.73～1.42及0～0.73。因此d7與d8層信號分別包含了周期為半周年與周年的周期信號。

圖4 BSRY站小波分解及各層功率譜

圖5 TABL站小波分解及各層功率譜

從各層的重建信號及其相應的功率譜可以看出，在d1—d2層的28～30 cpy頻帶上存在功率不足1 db的周期性誤差；d3、d4層分別在20～30 cpy以及10～24 cpy頻帶存在功率不超過20 db的周期性誤差；d5、d6層分別在6～12 cpy以及2～6 cpy頻帶存在最大功率不超過40 db的周期性誤差；在d7～d8層信號中存在較為明顯的振幅均小于5 mm的周年及半周年信號，其功率分別為280和410 db，其中TABL與BSRY的周年振幅均小于2 mm。此外，在a8的功率譜中還存在功率較小但周期約為兩年的信號。由于a8包含了原始信號中除線性趨勢之外的總體趨勢信息，考慮到信號的擬合度，本文對其予以保留。根據(jù)以上分析，本文利用小波逆變換提取含有周期大于半周年的d7、d8以及a8信號，舍棄包含高頻誤差的d1—d6信號，最后將提取出的信號與最初從原始時間序列中剔除的線性趨勢信號相加，獲得具有較高信噪比的時間序列。

為進一步比較濾波與小波變換后曲線擬合的效果，本文假設只存在白噪聲的情況下，利用最小二乘分別對經(jīng)過濾波和同時經(jīng)過濾波與小波變換后的時間序列去除線性趨勢、周年半周年信號，獲得殘差時間序列的其標準差，具體結(jié)果見表1。其中，Ⅰ列數(shù)據(jù)為原始時間序列經(jīng)過區(qū)域疊加濾波后擬合殘差的標準差；Ⅱ列數(shù)據(jù)為對時間序列進行區(qū)域疊加濾波以及小波變換后擬合殘差的標準差。比較Ⅰ列與Ⅱ列數(shù)據(jù)可知，經(jīng)過濾波以及小波處理后的時間序列所含噪聲比濾波后的時間序列大幅度減小。其中水平方向減小幅度約為40%，垂直方向減小幅度超過50%。

3 結(jié) 語

本文首先利用區(qū)域疊加濾波法去除了南加利福尼亞地區(qū)16個測站坐標時間序列的共模誤差，在此基礎上通過靜態(tài)離散小波變換去除頻率大于2.85 cpy的誤差。結(jié)果表明，通過區(qū)域疊加濾波法去除共模誤差，可以削弱周期性誤差尤其是周年與半周年誤差對原始坐標時間序列周期信號的影響。對去除共模誤差后的時間序列再進行靜態(tài)離散小波變換，可以從濾波后的坐標時間序列中提取周期為半周年以上的長周期信號，進一步剔除剩余誤差的影響，最終實現(xiàn)信噪的有效分離。這有助于提高GPS測站的位置精度及其非線性運動的地球物理因素的解釋，對于高精度的地殼形變監(jiān)測具有積極的意義。

表1 濾波及小波處理后測站各分量時間序列擬合殘差的標準差 mm

[1] DONG D, FANG P, BOCK Y, et al. Anatomy of Apparent Seasonal Variations from GPS-derived Site Position Time Series[J]. Journal of Geophysical Research Atmospheres, 2002, 107(B4): 9-16.

[2] STEWART M P, PENNA N T, LICHTI D D, et al. Investigating the Propagation Mechanism of Unmodelled Systematic Errors on Coordinate Time Series Estimated Using Least Squares[J]. Journal of Geodesy, 2005, 79(8):479-489.

[3] PENNA N T, KING M A, STEWART M P. GPS Height Time Series: Short-period Origins of Spurious Long-period Signals[J]. Journal of Geophysical Research, 2007, 112(B2):1074-1086.

[4] PENNA N T, STEWART M P. Aliased Tidal Signatures in Continuous GPS Height Time Series[J]. Geophysical Research Letters, 2003, 30(23):69-73.

[5] DONG D, FANG P, BOCK Y, et al. Spatiotemporal Filtering Using Principal Component Analysis and Karhunen-Loeve Expansion Approaches for Regional GPS Network Analysis[J]. Journal of Geophysical Research Atmospheres, 2006, 111(B3):1581-1600.

[6] WDOWINSKI S, BOCK Y, ZHANG J, et al. Southern California Permanent GPS Geodetic Array: Spatial Filtering of Daily Positions for Estimating Coseismic and Postseismic Displacements Induced by the 1992 Landers Earthquake[J]. Journal of Geophysical Research Atmospheres, 1997, 1021(1021):18057-18070.

[7] NIKOLAIDIS R. Observation of Geodetic and Seismic Deformation with the Global Positioning System[M]. [S.l.]: American Association for Cancer Research, 2002.

[8] JACKSON D A, CHEN Y. Robust Principal Component Analysis and Outlier Detection with Ecological Data[J]. Environmetrics, 2010, 15(2):129-139.

[9] 楊博, 周偉, 陳阜超,等. GPS連續(xù)站水平位置時間序列共模白噪聲識別與估計的歐拉-濾波法[J]. 山東科技大學學報(自然科學版), 2010, 29(3):26-31.

[10] 范玉磊, 王賀, 黃聲享，等. 利用小波分析IGS跟蹤站的非線性運動特征[J]. 測繪工程, 2014, 23(9): 5-8.

[11] 田亮, 孫付平.基于GPS測站坐標殘差序列的小波工具應用與分析[J]. 測繪工程, 2013, 22(1): 44-47.

[12] 李尊建，臧斌. 非平穩(wěn)大地測量信號特征信息小波識別[J]. 山東理工大學學報(自然科學版), 2009, 23(4): 58-61.

[14] 盛傳貞.中國大陸非構(gòu)造負荷地殼形變的區(qū)域性特征與改正模型[D]. 北京:中國地震局地質(zhì)研究所, 2013.

[15] 袁林果，丁曉利，陳武,等．香港GPS基準站坐標序列特征分析[J]．地球物理學報，2008，51(5)：1372-1384.

[16] 姜衛(wèi)平, 李昭, 鄧連生, 等. 高階電離層延遲對GPS 坐標時間序列的影響分析[J]. 科學通報, 2014, 59(10): 913-923.

[17] 唐江森, 曲國慶, 袁興明. 區(qū)域GPS網(wǎng)共模誤差的提取與分析[J]. 山東理工大學學報(自然科學版), 2016, 30(6): 48-52.

[18] 賀小星.GPS坐標序列噪聲模型估計方法研究[D]. 武漢:武漢大學, 2016.

[19] WU H, LI K, SHI W, et al. A Wavelet-based Hybrid Approach to Remove the Flicker Noise and the White Noise from GPS Coordinate Time Series[J]. GPS Solutions, 2015, 19(4):511-523.

CombineRegionStackFilteringandWaveletTransformationtoReduceNoiseinGPSCoordinateTimeSeries

MA Jun，ZHOU Xiaohui，ZHU Zhaohan

(School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China)

The periodic and nonperiodic errors in coordinate time series of GPS stations have seriously affected the analysis of the characteristics of the stations and the explanation of the physical mechanism of the nonlinear motion. Therefore, in order to weaken the influence of noise, the region stack filtering method is used to remove the common mode error of the time series of 16 stations in Southern California, which weakens the periodic errors including annual and semi-annual errors in time series. In order to remove the residual noise after filtering, the filtered signals are subjected to static discrete wavelet transform processing, and the signals with the periods of more than half a year are extracted. The results show that the combined region superposition filtering method and wavelet transform can deal with the time series of GPS station coordinate time series, which can weaken the influence of periodic errors of the signal and extract the nonlinear motion signals of the station as well.

GPS station coordinate time series；common mode error； wavelet transformation； signal noise separation

2017-04-10

國家杰出青年科學基金(41525014)；國家自然科學基金(41374033)

馬 俊(1985—)，男，碩士，主要研究方向為GPS站坐標時間序列噪聲模型。E-mail:794598508@qq.com

馬俊，周曉慧，朱兆涵.聯(lián)合區(qū)域疊加濾波法與小波變換去除GPS站坐標時間序列噪聲[J].測繪通報，2017(12)：6-11.

10.13474/j.cnki.11-2246.2017.0369.

P228

A

0494-0911(2017)12-0006-06