胡小雷，劉小喜

（中交上海航道勘察設(shè)計(jì)研究院有限公司，上海 200120）

0 引言

在港口與航道工程中，潮位資料是重要的參考資料[1]。為獲取臨時(shí)潮位資料，目前常用的方法有潮位推算、拋設(shè)壓力驗(yàn)潮儀、GNSS驗(yàn)潮等。潮位推算法受限于精度問題，不能直接用于生產(chǎn)[2]。拋設(shè)壓力驗(yàn)潮儀面臨儀器丟失導(dǎo)致數(shù)據(jù)空白、零點(diǎn)漂移、驗(yàn)潮儀在水下滑動致使數(shù)據(jù)失真等問題[3-4]。GNSS驗(yàn)潮可分為RTK驗(yàn)潮、PPK驗(yàn)潮和PPP驗(yàn)潮等，其原理均是從GNSS觀測的瞬間水面高程中獲取潮位數(shù)據(jù)，其精度主要受GNSS高程精度控制，而GNSS高程隨時(shí)間呈非平穩(wěn)性、明顯的周期性變化特征。為去除GNSS高程中的長期趨勢項(xiàng)及各類噪聲，需要對GNSS高程進(jìn)行濾波處理。目前，常用的濾波方法有傅里葉變換、小波變換，但傅里葉變換僅適用于平穩(wěn)信號，小波變換則需要預(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)行信號分解，具備自適應(yīng)性[6]。鑒于此，本文利用經(jīng)驗(yàn)?zāi)B(tài)分解法對PPK潮位數(shù)據(jù)進(jìn)行分解，將去除噪聲后的潮位資料進(jìn)行數(shù)據(jù)重構(gòu)，并與同期長期潮位站數(shù)據(jù)進(jìn)行比對分析，以探討該方法在PPK潮位測量中的適用性。

1 數(shù)據(jù)來源

2018年4月16—19日分別在連云港車牛山站、油碼頭站、開山島站等3個(gè)長期驗(yàn)潮站附近水域開展PPK潮位測量，基準(zhǔn)站架設(shè)于東西連島（L090），各驗(yàn)潮站PPK潮位測量時(shí)間段見表1，各驗(yàn)潮站位置及基準(zhǔn)站位置見圖1。PPK潮位觀測系統(tǒng)由Trimble SPS 855雙頻接收機(jī)、TSS姿態(tài)測量系統(tǒng)組成，在測量船重心上方位置分別固定GNSS接收機(jī)和TSS傳感器，測定GNSS天線在船體坐標(biāo)系下到水面的垂直距離。GNSS的采樣頻率為1 Hz，衛(wèi)星截止高度角為10°，TSS傳感器的采樣頻率為5 Hz。PPK潮位測量期間，同步收集車牛山站、油碼頭站、開山島站長期潮位站潮位資料。PPK潮位測量期間天氣狀況良好。

表1 各驗(yàn)潮站PPK潮位測量時(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

2 計(jì)算方法

2.1 經(jīng)驗(yàn)?zāi)B(tài)分解

經(jīng)驗(yàn)?zāi)B(tài)分解法（EMD）對原始時(shí)間序列進(jìn)行分解，將時(shí)間序列信號分解為一系列具有不同時(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è)稱為“篩選”的步驟來完成，對原始信號x（t）進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解的基本步驟為：

1）求得原始信號x（t）的局部極大值和局部極小值，分別擬合上、下包絡(luò)線；

2）分別計(jì)算上、下包絡(luò)線的均值；

3）分別計(jì)算原始信號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）從原始信號中減去該層IMF分量作為新的原始信號返回步驟2），計(jì)算第i+1層IMF分量；

8）得到n個(gè)IMF分量和剩余分量。經(jīng)分解得到一系列IMF分量后可表示為：

式中：rn（t）代表最終余項(xiàng)，表示信號的趨勢成分；imfi（t）表示頻率從高到低排列的n個(gè)IMF分量。

2.2 數(shù)據(jù)重構(gòu)

根據(jù)經(jīng)驗(yàn)?zāi)B(tài)分解過程可知，分解獲得的n個(gè)IMF分量中，噪聲信號對每個(gè)IMF分量的影響逐漸減小，信號對IMF分量的影響則不斷增大。因此，為了獲取準(zhǔn)確的潮位數(shù)據(jù)，需要在m個(gè)IMF分量中確定潮位數(shù)據(jù)和噪聲信號的分界IMF分量，以達(dá)到去除噪聲信號的目的?；静襟E如下：首先，求取m個(gè)IMF分量與原序列的相關(guān)系數(shù)；其次，搜索求得的一系列相關(guān)系數(shù)中第1個(gè)取局部極小值的相關(guān)系數(shù)對應(yīng)的IMF分量，其為噪聲與信號的分界imfk分量；最后，將k及k之前的IMF分量重構(gòu)為噪聲序列，將k之后的IMF分量重構(gòu)為周期項(xiàng)序列，通過重構(gòu)被分解序列可以有效提取原序列的周期項(xiàng)[8]。

3 結(jié)果與分析

3.1 PPK解算結(jié)果

PPK數(shù)據(jù)后處理采用TBC解算軟件，計(jì)算過程中均使用精密星歷。由于PPK處理后獲得水面高程的基面為GNSS橢球面，而潮位數(shù)據(jù)的基面為理論深度基準(zhǔn)面，二者之間需要轉(zhuǎn)換。利用車牛山站、油碼頭站、開山島站轉(zhuǎn)換關(guān)系得到各站基于理論深度基準(zhǔn)面的PPK潮位，并和各長期站同步潮位數(shù)據(jù)共同繪制潮位過程線圖。以車牛山站為例，從圖2可以看出，車牛山站PPK驗(yàn)潮數(shù)據(jù)和車牛山站長期驗(yàn)潮站同期潮位變化趨勢基本一致，但PPK驗(yàn)潮數(shù)據(jù)呈明顯波動特征，這說明車牛山站PPK潮位數(shù)據(jù)中存在噪聲序列。因此，為獲取準(zhǔn)確的潮位數(shù)據(jù)，必須對PPK潮位數(shù)據(jù)進(jìn)行濾波處理，以去除潮位數(shù)據(jù)中的噪聲序列。

圖2 車牛山站長期潮位與PPK潮位過程線圖Fig.2 Chart of long-term tide level and PPK tide level at Cheniushan station

3.2 經(jīng)驗(yàn)?zāi)B(tài)分解與數(shù)據(jù)重構(gòu)

通過經(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)將信號自適應(yīng)地分解為12個(gè)高頻分量imf1～imf12、1個(gè)低頻分量imf13和1個(gè)剩余分量。

分別求取各站分解后IMF分量和原始潮位信號的相關(guān)系數(shù)，計(jì)算結(jié)果見圖4。從圖4可以看出，車牛山站、開山島站和油碼頭站imf1～imf12分量和原始潮位信號的相關(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潮位和長期站同步潮位過程線，并以車牛山站為例說明，車牛山站潮位過程線見圖5。由圖5可知，去噪后的PPK潮位和長期站潮位趨勢一致，且潮位信號波形清晰，這表明經(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 車牛山站長期潮位和經(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

3.3 誤差對比與分析

為定量說明經(jīng)驗(yàn)?zāi)B(tài)分解法的去噪效果，分別將去噪前后的PPK潮位和同期長期站潮位的最大誤差、平均絕對誤差和均方根誤差進(jìn)行對比統(tǒng)計(jì)。由于車牛山站長期潮位時(shí)間間隔為5 min，油碼頭站和開山島站均為10 min，因此按各長期潮位站的時(shí)間間隔提取重構(gòu)后的PPK潮位數(shù)據(jù)。各潮位站的潮位精度統(tǒng)計(jì)結(jié)果見表2。

表2 經(jīng)驗(yàn)?zāi)B(tài)分解前后PPK潮位和長期站潮位精度統(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)潮站最大誤差、平均絕對誤差和均方根誤差均由大變小，如油碼頭站經(jīng)驗(yàn)?zāi)B(tài)分解前后的最大誤差分別為15.5 cm和9.7 cm、平均絕對誤差分別為6.3 cm和2.9 cm、均方根誤差分別為7.1 cm和4.6 cm。由此可見，經(jīng)驗(yàn)?zāi)B(tài)分解能夠有效濾除PPK潮位中的噪聲信號，可將潮位精度從10 cm量級提高至5 cm左右。

此外，距基準(zhǔn)站40 km的車牛山站、開山島站經(jīng)驗(yàn)?zāi)B(tài)分解后的平均絕對誤差為5.5 cm、4.1 cm，而距基準(zhǔn)站3.5 km的油碼頭站經(jīng)驗(yàn)?zāi)B(tài)分解后的平均絕對誤差僅為2.9 cm。由此可見，PPK潮位測量的精度和距基準(zhǔn)站的距離呈負(fù)相關(guān)關(guān)系，即隨基準(zhǔn)站和驗(yàn)潮站距離的增大而減小。

4 結(jié)語

為提高PPK潮位測量精度，本文利用經(jīng)驗(yàn)?zāi)B(tài)分解（EMD）對車牛山站、油碼頭站、開山島站獲得的PPK潮位時(shí)間序列進(jìn)行分析、重構(gòu)，剔除PPK潮位測量中的噪聲序列，并將去噪后的PPK潮位和同期長期驗(yàn)潮站潮位進(jìn)行對比。結(jié)果表明：

1）經(jīng)驗(yàn)?zāi)B(tài)分解法能有效剔除PPK潮位數(shù)據(jù)中的高頻噪聲序列，有助于獲取真實(shí)潮位數(shù)據(jù)。

2）PPK潮位和長期站潮位的平均絕對誤差在10 cm以內(nèi)，經(jīng)驗(yàn)?zāi)B(tài)分解、重構(gòu)后平均絕對誤差在5 cm左右，PPK潮位數(shù)據(jù)可滿足港口與航道工程的測量精度要求。

3）PPK潮位精度和基準(zhǔn)站與驗(yàn)潮站之間距離呈負(fù)相關(guān)關(guān)系，即驗(yàn)潮站距基準(zhǔn)站越遠(yuǎn)，PPK解算精度越低。