胡小雷,劉小喜
(中交上海航道勘察設(shè)計(jì)研究院有限公司,上海 200120)
在港口與航道工程中,潮位資料是重要的參考資料[1]。為獲取臨時(shí)潮位資料,目前常用的方法有潮位推算、拋設(shè)壓力驗(yàn)潮儀、GNSS驗(yàn)潮等。潮位推算法受限于精度問題,不能直接用于生產(chǎn)[2]。拋設(shè)壓力驗(yàn)潮儀面臨儀器丟失導(dǎo)致數(shù)據(jù)空白、零點(diǎn)漂移、驗(yàn)潮儀在水下滑動(dòng)致使數(shù)據(jù)失真等問題[3-4]。GNSS驗(yàn)潮可分為RTK驗(yàn)潮、PPK驗(yàn)潮和PPP驗(yàn)潮等,其原理均是從GNSS觀測(cè)的瞬間水面高程中獲取潮位數(shù)據(jù),其精度主要受GNSS高程精度控制,而GNSS高程隨時(shí)間呈非平穩(wěn)性、明顯的周期性變化特征。為去除GNSS高程中的長(zhǎng)期趨勢(shì)項(xiàng)及各類噪聲,需要對(duì)GNSS高程進(jìn)行濾波處理。目前,常用的濾波方法有傅里葉變換、小波變換,但傅里葉變換僅適用于平穩(wěn)信號(hào),小波變換則需要預(yù)先設(shè)定小波基函數(shù),缺乏自適應(yīng)性,二者均不適合處理具有非平穩(wěn)、非線性的GNSS高程數(shù)據(jù)[5]。經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)完全拋開了基函數(shù)的束縛,僅僅依據(jù)數(shù)據(jù)自身的時(shí)間尺度特征來進(jìn)行信號(hào)分解,具備自適應(yīng)性[6]。鑒于此,本文利用經(jīng)驗(yàn)?zāi)B(tài)分解法對(duì)PPK潮位數(shù)據(jù)進(jìn)行分解,將去除噪聲后的潮位資料進(jìn)行數(shù)據(jù)重構(gòu),并與同期長(zhǎng)期潮位站數(shù)據(jù)進(jìn)行比對(duì)分析,以探討該方法在PPK潮位測(cè)量中的適用性。
2018年4月16—19日分別在連云港車牛山站、油碼頭站、開山島站等3個(gè)長(zhǎng)期驗(yàn)潮站附近水域開展PPK潮位測(cè)量,基準(zhǔn)站架設(shè)于東西連島(L090),各驗(yàn)潮站PPK潮位測(cè)量時(shí)間段見表1,各驗(yàn)潮站位置及基準(zhǔn)站位置見圖1。PPK潮位觀測(cè)系統(tǒng)由Trimble SPS 855雙頻接收機(jī)、TSS姿態(tài)測(cè)量系統(tǒng)組成,在測(cè)量船重心上方位置分別固定GNSS接收機(jī)和TSS傳感器,測(cè)定GNSS天線在船體坐標(biāo)系下到水面的垂直距離。GNSS的采樣頻率為1 Hz,衛(wèi)星截止高度角為10°,TSS傳感器的采樣頻率為5 Hz。PPK潮位測(cè)量期間,同步收集車牛山站、油碼頭站、開山島站長(zhǎng)期潮位站潮位資料。PPK潮位測(cè)量期間天氣狀況良好。
表1 各驗(yàn)潮站PPK潮位測(cè)量時(shí)段Table 1 Tidal level measurement period of PPK at each tide gauge station
圖1 驗(yàn)潮站位置及基準(zhǔn)站位置示意圖Fig.1 Sketch map of the location of tide gauge station and reference station
經(jīng)驗(yàn)?zāi)B(tài)分解法(EMD)對(duì)原始時(shí)間序列進(jìn)行分解,將時(shí)間序列信號(hào)分解為一系列具有不同時(shí)間尺度的IMF分量(Intrinsic Mode Function)和1個(gè)最終余項(xiàng),每個(gè)IMF分量需要滿足2個(gè)條件:1)在全體數(shù)據(jù)序列中,極值點(diǎn)和過零點(diǎn)的個(gè)數(shù)必須相同或者最多相差1個(gè);2)在任意一點(diǎn),由局部極小值點(diǎn)組成的下包絡(luò)線和局部極大值點(diǎn)確定的上包絡(luò)線的算術(shù)平均值為0[7]。分解過程通過一個(gè)稱為“篩選”的步驟來完成,對(duì)原始信號(hào)x(t)進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解的基本步驟為:
1)求得原始信號(hào)x(t)的局部極大值和局部極小值,分別擬合上、下包絡(luò)線;
2)分別計(jì)算上、下包絡(luò)線的均值;
3)分別計(jì)算原始信號(hào)x(t)減去上、下包絡(luò)線的均值;
4)判斷是否滿足IMF分量的2個(gè)篩選條件,若滿足,則進(jìn)入下一步;若不滿足,則將步驟3)中結(jié)果返回步驟1)重新計(jì)算;
5)分解結(jié)果作為第i個(gè)IMF分量;
6)判斷是否滿足EMD分解終止條件,若滿足轉(zhuǎn)入步驟8),若不滿足則轉(zhuǎn)入步驟7);
7)從原始信號(hào)中減去該層IMF分量作為新的原始信號(hào)返回步驟2),計(jì)算第i+1層IMF分量;
8)得到n個(gè)IMF分量和剩余分量。經(jīng)分解得到一系列IMF分量后可表示為:
式中:rn(t)代表最終余項(xiàng),表示信號(hào)的趨勢(shì)成分;imfi(t)表示頻率從高到低排列的n個(gè)IMF分量。
根據(jù)經(jīng)驗(yàn)?zāi)B(tài)分解過程可知,分解獲得的n個(gè)IMF分量中,噪聲信號(hào)對(duì)每個(gè)IMF分量的影響逐漸減小,信號(hào)對(duì)IMF分量的影響則不斷增大。因此,為了獲取準(zhǔn)確的潮位數(shù)據(jù),需要在m個(gè)IMF分量中確定潮位數(shù)據(jù)和噪聲信號(hào)的分界IMF分量,以達(dá)到去除噪聲信號(hào)的目的。基本步驟如下:首先,求取m個(gè)IMF分量與原序列的相關(guān)系數(shù);其次,搜索求得的一系列相關(guān)系數(shù)中第1個(gè)取局部極小值的相關(guān)系數(shù)對(duì)應(yīng)的IMF分量,其為噪聲與信號(hào)的分界imfk分量;最后,將k及k之前的IMF分量重構(gòu)為噪聲序列,將k之后的IMF分量重構(gòu)為周期項(xiàng)序列,通過重構(gòu)被分解序列可以有效提取原序列的周期項(xiàng)[8]。
PPK數(shù)據(jù)后處理采用TBC解算軟件,計(jì)算過程中均使用精密星歷。由于PPK處理后獲得水面高程的基面為GNSS橢球面,而潮位數(shù)據(jù)的基面為理論深度基準(zhǔn)面,二者之間需要轉(zhuǎn)換。利用車牛山站、油碼頭站、開山島站轉(zhuǎn)換關(guān)系得到各站基于理論深度基準(zhǔn)面的PPK潮位,并和各長(zhǎng)期站同步潮位數(shù)據(jù)共同繪制潮位過程線圖。以車牛山站為例,從圖2可以看出,車牛山站PPK驗(yàn)潮數(shù)據(jù)和車牛山站長(zhǎng)期驗(yàn)潮站同期潮位變化趨勢(shì)基本一致,但PPK驗(yàn)潮數(shù)據(jù)呈明顯波動(dòng)特征,這說明車牛山站PPK潮位數(shù)據(jù)中存在噪聲序列。因此,為獲取準(zhǔn)確的潮位數(shù)據(jù),必須對(duì)PPK潮位數(shù)據(jù)進(jìn)行濾波處理,以去除潮位數(shù)據(jù)中的噪聲序列。
圖2 車牛山站長(zhǎng)期潮位與PPK潮位過程線圖Fig.2 Chart of long-term tide level and PPK tide level at Cheniushan station
通過經(jīng)驗(yàn)?zāi)B(tài)分解車牛山站、油碼頭站、開山島站PPK潮位,并以車牛山站為例說明經(jīng)驗(yàn)?zāi)B(tài)分解效果,車牛山站2018年4月16—17日PPK潮位數(shù)據(jù)分解結(jié)果如圖3所示。
圖3 車牛山站PPK潮位經(jīng)驗(yàn)?zāi)B(tài)分解結(jié)果Fig.3 Empirical mode decomposition results of PPK tidal level at Cheniushan station
從圖3可以看出,基于潮位數(shù)據(jù)自身特點(diǎn)將信號(hào)自適應(yīng)地分解為12個(gè)高頻分量imf1~imf12、1個(gè)低頻分量imf13和1個(gè)剩余分量。
分別求取各站分解后IMF分量和原始潮位信號(hào)的相關(guān)系數(shù),計(jì)算結(jié)果見圖4。從圖4可以看出,車牛山站、開山島站和油碼頭站imf1~imf12分量和原始潮位信號(hào)的相關(guān)性均小于0.1,相關(guān)性很差;車牛山站、開山島站和油碼頭站imf13分量的相關(guān)系數(shù)分別為0.987、0.996和0.998,相關(guān)性很好。因此,將imf1~imf12分量作為噪聲序列予以去除,將imf13分量作為周期項(xiàng)進(jìn)行重構(gòu)。潮位數(shù)據(jù)重構(gòu)后,分別繪制各站PPK潮位和長(zhǎng)期站同步潮位過程線,并以車牛山站為例說明,車牛山站潮位過程線見圖5。由圖5可知,去噪后的PPK潮位和長(zhǎng)期站潮位趨勢(shì)一致,且潮位信號(hào)波形清晰,這表明經(jīng)驗(yàn)?zāi)B(tài)分解法能夠有效濾除PPK潮位數(shù)據(jù)中噪聲序列。
圖4 各站IMF分量和原始數(shù)據(jù)相關(guān)系數(shù)圖Fig.4 Graph of correlation coefficients between IMF components and raw data at each station
圖5 車牛山站長(zhǎng)期潮位和經(jīng)驗(yàn)?zāi)B(tài)分解前后PPK潮位過程線圖Fig.5 PPK tide level process diagram before and after long-term tidal level and empirical mode decomposition at Cheniushan Station
為定量說明經(jīng)驗(yàn)?zāi)B(tài)分解法的去噪效果,分別將去噪前后的PPK潮位和同期長(zhǎng)期站潮位的最大誤差、平均絕對(duì)誤差和均方根誤差進(jìn)行對(duì)比統(tǒng)計(jì)。由于車牛山站長(zhǎng)期潮位時(shí)間間隔為5 min,油碼頭站和開山島站均為10 min,因此按各長(zhǎng)期潮位站的時(shí)間間隔提取重構(gòu)后的PPK潮位數(shù)據(jù)。各潮位站的潮位精度統(tǒng)計(jì)結(jié)果見表2。
表2 經(jīng)驗(yàn)?zāi)B(tài)分解前后PPK潮位和長(zhǎng)期站潮位精度統(tǒng)計(jì)表Table 2 Statistical table of tidal level accuracy of PPK tidal level and long-term station before and after empirical mode decomposition cm
從統(tǒng)計(jì)結(jié)果可以看出,各驗(yàn)潮站最大誤差、平均絕對(duì)誤差和均方根誤差均由大變小,如油碼頭站經(jīng)驗(yàn)?zāi)B(tài)分解前后的最大誤差分別為15.5 cm和9.7 cm、平均絕對(duì)誤差分別為6.3 cm和2.9 cm、均方根誤差分別為7.1 cm和4.6 cm。由此可見,經(jīng)驗(yàn)?zāi)B(tài)分解能夠有效濾除PPK潮位中的噪聲信號(hào),可將潮位精度從10 cm量級(jí)提高至5 cm左右。
此外,距基準(zhǔn)站40 km的車牛山站、開山島站經(jīng)驗(yàn)?zāi)B(tài)分解后的平均絕對(duì)誤差為5.5 cm、4.1 cm,而距基準(zhǔn)站3.5 km的油碼頭站經(jīng)驗(yàn)?zāi)B(tài)分解后的平均絕對(duì)誤差僅為2.9 cm。由此可見,PPK潮位測(cè)量的精度和距基準(zhǔn)站的距離呈負(fù)相關(guān)關(guān)系,即隨基準(zhǔn)站和驗(yàn)潮站距離的增大而減小。
為提高PPK潮位測(cè)量精度,本文利用經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)對(duì)車牛山站、油碼頭站、開山島站獲得的PPK潮位時(shí)間序列進(jìn)行分析、重構(gòu),剔除PPK潮位測(cè)量中的噪聲序列,并將去噪后的PPK潮位和同期長(zhǎng)期驗(yàn)潮站潮位進(jìn)行對(duì)比。結(jié)果表明:
1)經(jīng)驗(yàn)?zāi)B(tài)分解法能有效剔除PPK潮位數(shù)據(jù)中的高頻噪聲序列,有助于獲取真實(shí)潮位數(shù)據(jù)。
2)PPK潮位和長(zhǎng)期站潮位的平均絕對(duì)誤差在10 cm以內(nèi),經(jīng)驗(yàn)?zāi)B(tài)分解、重構(gòu)后平均絕對(duì)誤差在5 cm左右,PPK潮位數(shù)據(jù)可滿足港口與航道工程的測(cè)量精度要求。
3)PPK潮位精度和基準(zhǔn)站與驗(yàn)潮站之間距離呈負(fù)相關(guān)關(guān)系,即驗(yàn)潮站距基準(zhǔn)站越遠(yuǎn),PPK解算精度越低。