邵明政

（廣饒縣水利工程公司，山東 廣饒 257300）

0.引言

20世紀(jì)末我國(guó)開(kāi)始進(jìn)行CORS建設(shè)，CORS的主要作用包括提供位置服務(wù)、推動(dòng)氣象學(xué)與地球動(dòng)力學(xué)等多領(lǐng)域?qū)W科發(fā)展[1，2]。近些年，我國(guó)的城市建設(shè)不斷加快，各種城市化建設(shè)所需的信息需求也在不斷加大，大多數(shù)城市已經(jīng)建設(shè)滿足當(dāng)?shù)匚恢梅?wù)需求的CORS網(wǎng)并投入使用。CORS系統(tǒng)作為一種集成互聯(lián)網(wǎng)技術(shù)、數(shù)據(jù)通信技術(shù)、衛(wèi)星導(dǎo)航技術(shù)的實(shí)時(shí)導(dǎo)航定位服務(wù)系統(tǒng)，具有自動(dòng)化、網(wǎng)絡(luò)化等優(yōu)勢(shì)[3]。目前，CORS系統(tǒng)在基礎(chǔ)測(cè)繪、工程測(cè)量、環(huán)境監(jiān)測(cè)等領(lǐng)域發(fā)揮著重要作用，因此，準(zhǔn)確掌握CORS系統(tǒng)中站點(diǎn)觀測(cè)數(shù)據(jù)特征，評(píng)估站點(diǎn)在各種復(fù)雜環(huán)境下的穩(wěn)定性狀態(tài)，對(duì)于CORS站提供更好、更穩(wěn)定、更持續(xù)的位置服務(wù)具有重要意義。

作為一種對(duì)非線性、非平穩(wěn)性信號(hào)具有良好處理效果的處理方法，小波分析的出現(xiàn)打破了GNSS數(shù)據(jù)處理方法的壁壘，具有的時(shí)頻特征使其在GNSS數(shù)據(jù)領(lǐng)域受到越來(lái)越多的應(yīng)用[4，5]。近些年，不斷有學(xué)者對(duì)小波分析進(jìn)行深入研究，使其不斷優(yōu)化。奇異譜分析是一種對(duì)信號(hào)具有穩(wěn)定識(shí)別與強(qiáng)化功能的信號(hào)處理方法，近些年也被用于GNSS數(shù)據(jù)處理中[6]。本文充分依托小波分析在數(shù)據(jù)處理中的優(yōu)勢(shì)，將小波分析應(yīng)用于GNSS觀測(cè)數(shù)據(jù)不同成分提取中，實(shí)現(xiàn)CORS站點(diǎn)穩(wěn)定性分析。最后將小波分析與奇異譜分析共同應(yīng)用于CORS站點(diǎn)GNSS觀測(cè)數(shù)據(jù)趨勢(shì)項(xiàng)提取中，并對(duì)二者提取效果進(jìn)行了對(duì)比。

1.研究方法

1.1 小波變換分析

小波變換方法自提出起就快速成為信號(hào)領(lǐng)域應(yīng)用最為廣泛的方法之一，相比于傳統(tǒng)的Fourier變換信號(hào)處理方法，小波變換能夠取得更加顯著的效果。

1.1.1連續(xù)小波變換

假設(shè)存在Ψ（t）∈L2（R），其Fourier變換為φ^(ω)，當(dāng)φ^(ω)滿足完全重構(gòu)或者恒等分分辨率條件時(shí)[7]：

這說(shuō)明式（1）是有界的，此時(shí)Ψ（t）為基本小波，伸縮平移基本小波后得到新的小波序列：

式（2）中，a、b分別為小波序列的伸縮因子與平移因子。

連續(xù)小波變換主要有以下幾個(gè)方面的性質(zhì)：

（1）線性性質(zhì)

（2）時(shí)標(biāo)定理

（3）時(shí)移共變性

（4）能量守恒

小波變換過(guò)程中不增加信號(hào)能量，信號(hào)能量始終保持穩(wěn)定。

1.1.2離散小波變換

連續(xù)小波變換中，a與b不斷變化，為減少冗余，可進(jìn)行離散化處理。離散方式是對(duì)a與b作冪級(jí)數(shù)，分別取a=ajo，b=kajob0，其中j∈Z，擴(kuò)展步長(zhǎng)a0>1，離散小波函數(shù)可表示為：

重構(gòu)公式為：

式（9）中的常數(shù)c與信號(hào)無(wú)關(guān)。

1.1.3小波分析分解層次的確定

通常來(lái)說(shuō)，在進(jìn)行信號(hào)小波分析時(shí)，信號(hào)中的噪聲會(huì)隨著分解層次的增加被更多地濾除。然而，分解層次在一定范圍內(nèi)是可行的，超過(guò)一定范圍會(huì)將信號(hào)中真實(shí)信號(hào)剔除，造成信號(hào)缺失。因此，在進(jìn)行小波分析時(shí)分解層次尤為關(guān)鍵，直接影響小波分析去噪效果。信號(hào)的信噪比對(duì)于小波分析分解層次的影響較大，若信噪比較小，信號(hào)中包含的噪聲過(guò)多，此時(shí)若將分解層次設(shè)置過(guò)小會(huì)剔除過(guò)多噪聲。

在對(duì)CORS站監(jiān)測(cè)數(shù)據(jù)進(jìn)行小波去噪時(shí)，信噪比無(wú)法得知，因此若要得到最優(yōu)分解層次須通過(guò)估計(jì)的方式。首先將分解層次定為1，在此基礎(chǔ)上依次累加，計(jì)算在不同分解層次下去噪前后信號(hào)均方根誤差得到最優(yōu)分解層次：

式（10）中，f（n）與fj（n）分別表示去噪前后信號(hào)，計(jì)算分解層次為j+1與分解層次為j時(shí)的均方根誤差比：

當(dāng)存在最接近1的r值時(shí)，表明噪聲已被全部濾除，此時(shí)分解層次可定為j或者j+1。

1.2 奇異譜分析

奇異譜分析是由多元統(tǒng)計(jì)方法、時(shí)間序列分析方法等多種技術(shù)方法組成的，能夠?qū)Ψ欠€(wěn)定性信號(hào)進(jìn)行有效處理，提取得到信號(hào)中包含的振動(dòng)信息與變形信息，實(shí)現(xiàn)變形體變形趨勢(shì)的有效判斷[8]。

實(shí)現(xiàn)奇異譜分析功能的主要途徑為重建成分RC（Recons truction Component），假設(shè)存在某序列xi，該序列的重建成分RC可由T-EOF與T-PC組成，可表示為：

原始信號(hào)可以看成是若干個(gè)重建成分RC的集合，即將所有重建成分RC相加就是原始信號(hào)序列，表示為：

奇異譜分析中，若要判斷重建成分RC是否為趨勢(shì)項(xiàng)，可通過(guò)Kendall非參數(shù)檢驗(yàn)的方式，主要步驟為：

（1）當(dāng)滿足xi，k

（2）若檢驗(yàn)得到xk為非趨勢(shì)項(xiàng)，表明τ服從正態(tài)分布，均值為0、均方差為S：

（3）置信度α取0.05，當(dāng)滿足τ<-1.96S或τ>1.96S時(shí)，表明第k個(gè)重建成分RC為趨勢(shì)項(xiàng)成分。

2.CORS坐標(biāo)序列解算

作為一款GNSS數(shù)據(jù)解算領(lǐng)域較為成熟穩(wěn)定的軟件，GAMIT/GLOBK軟件自問(wèn)世以來(lái)就得到了廣泛應(yīng)用。本文選擇GAMIT/GLOBK軟件對(duì)CORS站點(diǎn)觀測(cè)數(shù)據(jù)進(jìn)行解算，解算模型為雙差模型，通過(guò)雙差模型與衛(wèi)星星歷可達(dá)到毫米級(jí)坐標(biāo)解算精度。

通過(guò)GAMIT和GLOBK這兩個(gè)模塊，GAMIT/GLOBK軟件可快速高精度實(shí)現(xiàn)GNSS觀測(cè)數(shù)據(jù)解算，通過(guò)GAMIT得到單天解h文件并通過(guò)GLOBK對(duì)h文件進(jìn)行后處理，基于GAMIT/GLOBK軟件進(jìn)行GNSS數(shù)據(jù)解算的主要流程，如圖1所示。

圖1基于GAMIT/GLOBK軟件的GNSS數(shù)據(jù)解算

實(shí)驗(yàn)數(shù)據(jù)選擇BRTD站點(diǎn)2015年1月1日至2018年12月31日觀測(cè)GNSS數(shù)據(jù)，通過(guò)GAMIT/GLOBK軟件解算得到實(shí)驗(yàn)數(shù)據(jù)單日松弛解，使用標(biāo)準(zhǔn)化均方根殘差NRMS與基線重復(fù)性對(duì)解算結(jié)果進(jìn)行評(píng)價(jià)。得到NRMS最小值與最大值分別為0.138 31與0.197 22，相關(guān)基線解算規(guī)范中規(guī)定NRMS值小于0.25[9]，結(jié)果表明GAMIT基線解算結(jié)果滿足精度指標(biāo)。

計(jì)算得到基線重復(fù)性常數(shù)部分與比例系數(shù)如表1所示。通過(guò)表1可知：基線解算重復(fù)性滿足精度要求。

表1基線向量重復(fù)性

GAMIT基線解算完成后，通過(guò)GLOBK軟件得到CORS站點(diǎn)在ITRF 2005坐標(biāo)系框架系的坐標(biāo)時(shí)間序列如圖2所示。

圖2 CORS站ITRF 2015坐標(biāo)系下解算坐標(biāo)

通過(guò)圖2可知：N、E方向上坐標(biāo)時(shí)間序列呈現(xiàn)的是線性變化，并且在N方向隨時(shí)間推遲上逐漸遞減，在E方向上隨時(shí)間推遲逐漸遞增，U方向上的變化趨勢(shì)為一種周期性并伴跳躍。

3.結(jié)果與分析

3.1 趨勢(shì)項(xiàng)成分提取

作為表述外界環(huán)境（如，地質(zhì)變化、土地變化）的一種狀態(tài)，趨勢(shì)項(xiàng)變化是GNSS觀測(cè)數(shù)據(jù)中是一個(gè)重要組成部分。通過(guò)提取CORS站監(jiān)測(cè)數(shù)據(jù)中的趨勢(shì)項(xiàng)，可以實(shí)現(xiàn)站點(diǎn)的趨勢(shì)性變形特征。本文通過(guò)小波分析方法提取站點(diǎn)觀測(cè)數(shù)據(jù)中的趨勢(shì)項(xiàng)，首先對(duì)站點(diǎn)3個(gè)方向的坐標(biāo)時(shí)間序列進(jìn)行小波分析，計(jì)算不同分解層次下的均方根誤差比值，當(dāng)分解層次為9與8時(shí)比值為1.06，可將分解層次定為8；以U方向時(shí)間序列為例進(jìn)行小波分解，U方向小波分解逼近信號(hào)與細(xì)節(jié)信號(hào)分別如圖3、圖4所示。

圖3 U方向逼近信號(hào)

圖4 U方向細(xì)節(jié)信號(hào)

通過(guò)圖3、圖4可知：信號(hào)的波動(dòng)隨著分解層次的增加逐漸緩和，本文小波分析方法中選擇db4小波基進(jìn)行8層分解。3個(gè)方向經(jīng)重構(gòu)后得到的重構(gòu)信號(hào)，如圖5所示，其中藍(lán)色細(xì)線表示原始時(shí)間序列，紅色粗線表示經(jīng)小波分析方法提取趨勢(shì)項(xiàng)成分。方向偏移與向北偏移；在U方向上呈現(xiàn)向上偏移。相較于U方向，N、E方向的偏移量更大并且3方向的偏移速率均有所不同。

3.2 周期項(xiàng)成分提取

周期性變化是變形體在一定時(shí)間內(nèi)呈現(xiàn)的規(guī)律性的往返變化，CORS站點(diǎn)GNSS觀測(cè)數(shù)據(jù)在外界環(huán)境影響下也存在穩(wěn)定周期性變化。作為一種能夠?qū)τ^測(cè)時(shí)間序列中的周期項(xiàng)成分進(jìn)行準(zhǔn)確提取的方法，本文選擇使用功率譜估計(jì)方法提取CORS站點(diǎn)3方向周期項(xiàng)成分。該方法包括兩個(gè)部分，先進(jìn)行傅里葉變換，再使用變換后信號(hào)除以樣本量，可表示為[10]：

上式中，N表示時(shí)間序列長(zhǎng)度；f表示時(shí)間序列頻率。

觀測(cè)數(shù)據(jù)經(jīng)過(guò)傅里葉變換后會(huì)表現(xiàn)出周期性變化，將周期性變化制作得到周期圖（PSD），CORS站點(diǎn)3方向上的PSD結(jié)果如圖6所示。

圖5 3方向趨勢(shì)項(xiàng)成分提取結(jié)果

通過(guò)圖5可知：經(jīng)小波分析提取趨勢(shì)項(xiàng)與原始時(shí)間序列的變形趨勢(shì)基本保持一致，能夠?qū)υ甲冃乌厔?shì)進(jìn)行反映。

對(duì)經(jīng)小波分析提取的變形趨勢(shì)進(jìn)行量化，得到變形趨勢(shì)如表2所示。

表2 CORS站點(diǎn)提取變形趨勢(shì)

通過(guò)表2可知：CORS站點(diǎn)在N、E方向上分別呈現(xiàn)向西

圖6 CORS站點(diǎn)3方向周期圖

通過(guò)圖6可知：N、E方向譜峰值均為292 d、512 d、682 d；U方向譜峰值分別為186 d、341 d。作為功率譜估計(jì)中的關(guān)鍵參數(shù)，窗口長(zhǎng)度的選取對(duì)于譜估計(jì)結(jié)果影響較大，一般取在序列長(zhǎng)度三分之一以內(nèi)。

計(jì)算在不同窗口長(zhǎng)度下經(jīng)小波重構(gòu)后信號(hào)與原始信號(hào)均方根誤差的變化如圖7所示。

圖7 CORS站點(diǎn)3方向重構(gòu)信號(hào)與原始信號(hào)均方根誤差隨窗口長(zhǎng)度變化

通過(guò)圖7可知：當(dāng)N、E、U方向的窗口長(zhǎng)度分別為292 d、292 d、341 d時(shí)，均方根變化圖會(huì)產(chǎn)生拐點(diǎn)，此對(duì)應(yīng)均方根誤差分別為0.095 12 mm、0.105 7 mm、0.002 147 mm。統(tǒng)計(jì)3方向上的周期性如表3所示。

表3 CORS站點(diǎn)3個(gè)方向的周期

通過(guò)上述圖表得知：CORS站點(diǎn)觀測(cè)數(shù)據(jù)在U方向的周期性要強(qiáng)于N、E方向。

3.3 與奇異譜分析方法的比較

將奇異譜分析應(yīng)用到CORS站點(diǎn)觀測(cè)時(shí)間序列趨勢(shì)項(xiàng)提取中，并將提取結(jié)果與小波分析提取結(jié)果進(jìn)行對(duì)比，結(jié)果如圖8所示。

圖8 CORS站點(diǎn)3方向小波分析與奇異譜分析提取結(jié)果

通過(guò)圖8可知：小波分析方法與奇異譜分析方法均不同程度將原始序列中包含噪聲剔除，提取得到時(shí)間序列主要變形趨勢(shì)。然而相較于奇異譜分析，小波分析去噪后序列更加平滑，反映出更多的時(shí)間序列特征。

4.結(jié)束語(yǔ)

本文通過(guò)對(duì)CORS站點(diǎn)觀測(cè)GNSS數(shù)據(jù)進(jìn)行解算以及趨勢(shì)項(xiàng)、周期項(xiàng)提取，結(jié)果表明：使用GAMIT/GLOBK軟件可快速高精度實(shí)現(xiàn)GNSS觀測(cè)數(shù)據(jù)解算，解算結(jié)果的標(biāo)準(zhǔn)化均方根殘差NRMS與基線重復(fù)性均滿足規(guī)范精度要求。小波分析提取CORS站點(diǎn)GNSS觀測(cè)數(shù)據(jù)N、E方向上有-0.0326 mm/d的負(fù)向趨勢(shì)運(yùn)動(dòng)、0.0779 mm/d的正向趨勢(shì)運(yùn)動(dòng)；在U方向上有0.0026 mm/d的正向趨勢(shì)運(yùn)動(dòng)。通過(guò)小波分析結(jié)合周期圖提取CORS站點(diǎn)GNSS觀測(cè)數(shù)據(jù)周期項(xiàng)，表明在N、E方向上年周期較為明顯；在U方向上季度、半年、年周期較為明顯。最后將小波分析與奇異譜分析提取CORS站觀測(cè)GNSS數(shù)據(jù)趨勢(shì)項(xiàng)結(jié)果進(jìn)行對(duì)比，表明兩種方法均能將時(shí)間序列主要特征表現(xiàn)出來(lái)，但是小波分析提取序列更加平滑，去噪效果更好。