胡小雷,劉小喜
(中交上海航道勘察設計研究院有限公司,上海 200120)
在港口與航道工程中,潮位資料是重要的參考資料[1]。為獲取臨時潮位資料,目前常用的方法有潮位推算、拋設壓力驗潮儀、GNSS驗潮等。潮位推算法受限于精度問題,不能直接用于生產[2]。拋設壓力驗潮儀面臨儀器丟失導致數據空白、零點漂移、驗潮儀在水下滑動致使數據失真等問題[3-4]。GNSS驗潮可分為RTK驗潮、PPK驗潮和PPP驗潮等,其原理均是從GNSS觀測的瞬間水面高程中獲取潮位數據,其精度主要受GNSS高程精度控制,而GNSS高程隨時間呈非平穩性、明顯的周期性變化特征。為去除GNSS高程中的長期趨勢項及各類噪聲,需要對GNSS高程進行濾波處理。目前,常用的濾波方法有傅里葉變換、小波變換,但傅里葉變換僅適用于平穩信號,小波變換則需要預先設定小波基函數,缺乏自適應性,二者均不適合處理具有非平穩、非線性的GNSS高程數據[5]。經驗模態分解(EMD)完全拋開了基函數的束縛,僅僅依據數據自身的時間尺度特征來進行信號分解,具備自適應性[6]。鑒于此,本文利用經驗模態分解法對PPK潮位數據進行分解,將去除噪聲后的潮位資料進行數據重構,并與同期長期潮位站數據進行比對分析,以探討該方法在PPK潮位測量中的適用性。
2018年4月16—19日分別在連云港車牛山站、油碼頭站、開山島站等3個長期驗潮站附近水域開展PPK潮位測量,基準站架設于東西連島(L090),各驗潮站PPK潮位測量時間段見表1,各驗潮站位置及基準站位置見圖1。PPK潮位觀測系統由Trimble SPS 855雙頻接收機、TSS姿態測量系統組成,在測量船重心上方位置分別固定GNSS接收機和TSS傳感器,測定GNSS天線在船體坐標系下到水面的垂直距離。GNSS的采樣頻率為1 Hz,衛星截止高度角為10°,TSS傳感器的采樣頻率為5 Hz。PPK潮位測量期間,同步收集車牛山站、油碼頭站、開山島站長期潮位站潮位資料。PPK潮位測量期間天氣狀況良好。

表1 各驗潮站PPK潮位測量時段Table 1 Tidal level measurement period of PPK at each tide gauge station

圖1 驗潮站位置及基準站位置示意圖Fig.1 Sketch map of the location of tide gauge station and reference station
經驗模態分解法(EMD)對原始時間序列進行分解,將時間序列信號分解為一系列具有不同時間尺度的IMF分量(Intrinsic Mode Function)和1個最終余項,每個IMF分量需要滿足2個條件:1)在全體數據序列中,極值點和過零點的個數必須相同或者最多相差1個;2)在任意一點,由局部極小值點組成的下包絡線和局部極大值點確定的上包絡線的算術平均值為0[7]。分解過程通過一個稱為“篩選”的步驟來完成,對原始信號x(t)進行經驗模態分解的基本步驟為:
1)求得原始信號x(t)的局部極大值和局部極小值,分別擬合上、下包絡線;
2)分別計算上、下包絡線的均值;
3)分別計算原始信號x(t)減去上、下包絡線的均值;
4)判斷是否滿足IMF分量的2個篩選條件,若滿足,則進入下一步;若不滿足,則將步驟3)中結果返回步驟1)重新計算;
5)分解結果作為第i個IMF分量;
6)判斷是否滿足EMD分解終止條件,若滿足轉入步驟8),若不滿足則轉入步驟7);
7)從原始信號中減去該層IMF分量作為新的原始信號返回步驟2),計算第i+1層IMF分量;
8)得到n個IMF分量和剩余分量。經分解得到一系列IMF分量后可表示為:

式中:rn(t)代表最終余項,表示信號的趨勢成分;imfi(t)表示頻率從高到低排列的n個IMF分量。
根據經驗模態分解過程可知,分解獲得的n個IMF分量中,噪聲信號對每個IMF分量的影響逐漸減小,信號對IMF分量的影響則不斷增大。因此,為了獲取準確的潮位數據,需要在m個IMF分量中確定潮位數據和噪聲信號的分界IMF分量,以達到去除噪聲信號的目的。基本步驟如下:首先,求取m個IMF分量與原序列的相關系數;其次,搜索求得的一系列相關系數中第1個取局部極小值的相關系數對應的IMF分量,其為噪聲與信號的分界imfk分量;最后,將k及k之前的IMF分量重構為噪聲序列,將k之后的IMF分量重構為周期項序列,通過重構被分解序列可以有效提取原序列的周期項[8]。
PPK數據后處理采用TBC解算軟件,計算過程中均使用精密星歷。由于PPK處理后獲得水面高程的基面為GNSS橢球面,而潮位數據的基面為理論深度基準面,二者之間需要轉換。利用車牛山站、油碼頭站、開山島站轉換關系得到各站基于理論深度基準面的PPK潮位,并和各長期站同步潮位數據共同繪制潮位過程線圖。以車牛山站為例,從圖2可以看出,車牛山站PPK驗潮數據和車牛山站長期驗潮站同期潮位變化趨勢基本一致,但PPK驗潮數據呈明顯波動特征,這說明車牛山站PPK潮位數據中存在噪聲序列。因此,為獲取準確的潮位數據,必須對PPK潮位數據進行濾波處理,以去除潮位數據中的噪聲序列。

圖2 車牛山站長期潮位與PPK潮位過程線圖Fig.2 Chart of long-term tide level and PPK tide level at Cheniushan station
通過經驗模態分解車牛山站、油碼頭站、開山島站PPK潮位,并以車牛山站為例說明經驗模態分解效果,車牛山站2018年4月16—17日PPK潮位數據分解結果如圖3所示。

圖3 車牛山站PPK潮位經驗模態分解結果Fig.3 Empirical mode decomposition results of PPK tidal level at Cheniushan station
從圖3可以看出,基于潮位數據自身特點將信號自適應地分解為12個高頻分量imf1~imf12、1個低頻分量imf13和1個剩余分量。
分別求取各站分解后IMF分量和原始潮位信號的相關系數,計算結果見圖4。從圖4可以看出,車牛山站、開山島站和油碼頭站imf1~imf12分量和原始潮位信號的相關性均小于0.1,相關性很差;車牛山站、開山島站和油碼頭站imf13分量的相關系數分別為0.987、0.996和0.998,相關性很好。因此,將imf1~imf12分量作為噪聲序列予以去除,將imf13分量作為周期項進行重構。潮位數據重構后,分別繪制各站PPK潮位和長期站同步潮位過程線,并以車牛山站為例說明,車牛山站潮位過程線見圖5。由圖5可知,去噪后的PPK潮位和長期站潮位趨勢一致,且潮位信號波形清晰,這表明經驗模態分解法能夠有效濾除PPK潮位數據中噪聲序列。

圖4 各站IMF分量和原始數據相關系數圖Fig.4 Graph of correlation coefficients between IMF components and raw data at each station

圖5 車牛山站長期潮位和經驗模態分解前后PPK潮位過程線圖Fig.5 PPK tide level process diagram before and after long-term tidal level and empirical mode decomposition at Cheniushan Station
為定量說明經驗模態分解法的去噪效果,分別將去噪前后的PPK潮位和同期長期站潮位的最大誤差、平均絕對誤差和均方根誤差進行對比統計。由于車牛山站長期潮位時間間隔為5 min,油碼頭站和開山島站均為10 min,因此按各長期潮位站的時間間隔提取重構后的PPK潮位數據。各潮位站的潮位精度統計結果見表2。

表2 經驗模態分解前后PPK潮位和長期站潮位精度統計表Table 2 Statistical table of tidal level accuracy of PPK tidal level and long-term station before and after empirical mode decomposition cm
從統計結果可以看出,各驗潮站最大誤差、平均絕對誤差和均方根誤差均由大變小,如油碼頭站經驗模態分解前后的最大誤差分別為15.5 cm和9.7 cm、平均絕對誤差分別為6.3 cm和2.9 cm、均方根誤差分別為7.1 cm和4.6 cm。由此可見,經驗模態分解能夠有效濾除PPK潮位中的噪聲信號,可將潮位精度從10 cm量級提高至5 cm左右。
此外,距基準站40 km的車牛山站、開山島站經驗模態分解后的平均絕對誤差為5.5 cm、4.1 cm,而距基準站3.5 km的油碼頭站經驗模態分解后的平均絕對誤差僅為2.9 cm。由此可見,PPK潮位測量的精度和距基準站的距離呈負相關關系,即隨基準站和驗潮站距離的增大而減小。
為提高PPK潮位測量精度,本文利用經驗模態分解(EMD)對車牛山站、油碼頭站、開山島站獲得的PPK潮位時間序列進行分析、重構,剔除PPK潮位測量中的噪聲序列,并將去噪后的PPK潮位和同期長期驗潮站潮位進行對比。結果表明:
1)經驗模態分解法能有效剔除PPK潮位數據中的高頻噪聲序列,有助于獲取真實潮位數據。
2)PPK潮位和長期站潮位的平均絕對誤差在10 cm以內,經驗模態分解、重構后平均絕對誤差在5 cm左右,PPK潮位數據可滿足港口與航道工程的測量精度要求。
3)PPK潮位精度和基準站與驗潮站之間距離呈負相關關系,即驗潮站距基準站越遠,PPK解算精度越低。