馬 俊，周曉慧，朱兆涵

(武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079)

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

馬 俊，周曉慧，朱兆涵

(武漢大學(xué)測(cè)繪學(xué)院，湖北 武漢 430079)

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

GPS站坐標(biāo)時(shí)間序列；共模誤差；小波變換；信噪分離

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

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

1 時(shí)間序列分析方法

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

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

(1)

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

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

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

(2)

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

(3)

代入式(2)可得

(4)

此時(shí)離散小波變換公式為

(5)

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

2 實(shí)例分析

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

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

圖1 GPS測(cè)站分布

2.2 試驗(yàn)結(jié)果

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

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

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

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

周年、半周年信號(hào)是當(dāng)今IGS全球基準(zhǔn)站最主要的信號(hào)特征，而濾波后的時(shí)間序列中仍存在未剔除干凈的較強(qiáng)的噪聲。因此本文通過小波變換，將濾波后的時(shí)間序列中信號(hào)的不同頻率成分分解到互不重疊的頻帶上，以此將周期為半周年以上的信號(hào)和噪聲分開。首先剔除濾波后時(shí)間序列的性趨勢(shì)，然后選用小波函數(shù)“coif5”，利用靜態(tài)離散小波變換將其分解為8層。圖3—圖5為以上測(cè)站垂向時(shí)間序列各層的小波分解信號(hào)及其功率譜圖。其中，d1—d8為對(duì)信號(hào)進(jìn)行小波分解后得到的高頻信號(hào)，包含了信號(hào)的細(xì)節(jié)信息；a8為第8層分解得到的低頻信號(hào)，反映了信號(hào)的整體非線性運(yùn)動(dòng)特征。根據(jù)采樣定理，d1—d8以及a8信號(hào)的主要頻帶范圍分別是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層信號(hào)分別包含了周期為半周年與周年的周期信號(hào)。

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

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

從各層的重建信號(hào)及其相應(yīng)的功率譜可以看出，在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層信號(hào)中存在較為明顯的振幅均小于5 mm的周年及半周年信號(hào)，其功率分別為280和410 db，其中TABL與BSRY的周年振幅均小于2 mm。此外，在a8的功率譜中還存在功率較小但周期約為兩年的信號(hào)。由于a8包含了原始信號(hào)中除線性趨勢(shì)之外的總體趨勢(shì)信息，考慮到信號(hào)的擬合度，本文對(duì)其予以保留。根據(jù)以上分析，本文利用小波逆變換提取含有周期大于半周年的d7、d8以及a8信號(hào)，舍棄包含高頻誤差的d1—d6信號(hào)，最后將提取出的信號(hào)與最初從原始時(shí)間序列中剔除的線性趨勢(shì)信號(hào)相加，獲得具有較高信噪比的時(shí)間序列。

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

3 結(jié) 語(yǔ)

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

表1 濾波及小波處理后測(cè)站各分量時(shí)間序列擬合殘差的標(biāo)準(zhǔn)差 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ù)站水平位置時(shí)間序列共模白噪聲識(shí)別與估計(jì)的歐拉-濾波法[J]. 山東科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2010, 29(3):26-31.

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

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

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

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

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

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

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

[18] 賀小星.GPS坐標(biāo)序列噪聲模型估計(jì)方法研究[D]. 武漢:武漢大學(xué), 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

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

馬 俊(1985—)，男，碩士，主要研究方向?yàn)镚PS站坐標(biāo)時(shí)間序列噪聲模型。E-mail:794598508@qq.com

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

10.13474/j.cnki.11-2246.2017.0369.

P228

A

0494-0911(2017)12-0006-06