劉倩，孫小淇

(青島大學數(shù)學與統(tǒng)計學院，山東 青島 266071)

引言

海洋表面因受天體運動產(chǎn)生的引潮力的作用，海水會出現(xiàn)周期性漲落的現(xiàn)象，稱此現(xiàn)象為潮汐[1]。潮流的變化會影響海上船舶的進出、貨物的進出口作業(yè)和施工計劃的制定等。青島港口涉及到運輸、系泊和沉放等相關業(yè)務，都與潮汐直接相關。

青島港位于黃海西海岸，半島南岸西部的膠州灣口附近，港內水域寬闊，為我國著名的天然良港。對潮汐的調和分析即找出影響某地潮汐的若干主要分潮，從而分析該處潮汐特性，或對該處潮汐變化進行預報。近代海洋潮汐學產(chǎn)生于17世紀后半葉，其任務是研究海洋潮汐的成因，分析潮汐特征及其隨時間的變化規(guī)律，并作出預報[2-5]。平衡潮理論于1687年由牛頓提出。1868年，開爾文設計了用于預報的調和分析法，并發(fā)明潮汐預報機。DARWIN[6]最早采用調和分析法預報潮汐，之后DOODSON[7-8]改進了調和分析法，提高了計算精度。HORN[9]最早用計算機進行潮汐調和分析，并在1960年第一個運用計算機進行潮汐調和分析的計算。近年來，我國學者改進原有理論方法，并應用到本國實際情況中。鄭有任等[10]采用T_TIDE程序對南海北部進行調和分析，總結海域潮汐特征，并利用結果預測潮高，比較全年和分季節(jié)的結果，統(tǒng)計分析預測結果與實測數(shù)據(jù)的殘差。程誠和鄭小萌[11]應用潮汐資料分析T_TIDE程序，對比連云港站點的調和分析結果與實測數(shù)據(jù)的誤差，并根據(jù)計算的調和常數(shù)進行潮位預測分析。尹朝暉等[12]對工程施工海域的潮汐資料進行調和分析，并探究不同分潮個數(shù)對預報精度的影響。

本文對青島港2019年全年逐時潮位資料進行調和分析，之后從2019年逐時潮高的調和分析結果中選擇不同分潮建立調和預報模型，進而對2019年1月的潮高進行預測，通過對相對誤差、判定系數(shù)結果進行分析，確定最優(yōu)調和預報模型。

1 分析原理與數(shù)據(jù)

1.1 調和分析法

潮汐靜力理論(或稱平衡潮理論)[13]是研究海水在引潮力作用下產(chǎn)生潮汐過程的理論，僅分析潮汐現(xiàn)象變化的基本規(guī)律和特點，不能獲得具體的潮汐變化規(guī)律。為研究青島港口潮汐的規(guī)律，據(jù)觀測的實時潮汐水位數(shù)據(jù)進行調和分析，得到其主要分潮的調和常數(shù)，然后進行潮汐預測。

調和分析法[14]是根據(jù)實測潮位數(shù)據(jù)計算各分潮的調和常數(shù)，即分潮振幅和分潮遲角，然后應用兩者預測潮位。為降低調和分析模型預報產(chǎn)生的偏差，研究人員常選取振幅較大的分潮進行預測。根據(jù)潮汐周期性運動的特點，調和分析法將潮汐視為若干個余弦函數(shù)的疊加，故潮汐水位可由下式表示：

(1)

式中：n為分潮數(shù)量，fk為交點因子，Hk為分潮振幅，σk為分潮角速度，vk為分潮初相位角，uk為天文相角的交角訂正角，gk為遲角，H0為海平面，故每個分潮為fHcos(σt+v+u-g)。其中f、σ、(v+u)均已知，故只需求出每個分潮的H和g，則此分潮便可求出。H和g稱為調和常數(shù)。

傳統(tǒng)的潮汐調和分析常采用 FOREMAN and NEUFELD[15]基于FORTRAN程序編寫的軟件進行，但古典潮汐計算理論忽略了對近日點的調制。2002年，PAWLOWICZ et al.[16]在FOREMAN and NEUFELD[15]程序的基礎上進行改進，發(fā)表了T_TIDE潮汐處理軟件包。相比于FOREMAN and NEUFELD[15]的程序，T_TIDE程序可對整年或較短的潮汐資料進行調和分析，還可解釋未確定的潮汐成分，所得分析結果包含置信區(qū)間[17-20]。此外，T_TIDE程序也允許少量的缺測存在。t_predic.m利用t_tide.m得到的調和常數(shù)預測未來某時刻的潮汐水位。本文利用T_TIDE工具箱對青島港口區(qū)域進行潮汐的調和分析。

1.2 資料來源與處理方法

文中收集了2019年1月1日—12月31日青島港驗潮站逐小時的潮位觀測資料。青島港驗潮站位于膠州灣口東北方，座落于青島大港一號碼頭，即36.08°N，120.32°E，港池由防波堤環(huán)抱而成，其入口呈喇叭形向西南敞開，最窄處寬約270 m，港池內水深一般為5～13 m。該港池為大型固定碼頭，此次仿真實驗數(shù)據(jù)來源于中國港口網(wǎng)(http://www.chinaports.com/tidal)。

使用T_TIDE程序包執(zhí)行帶有節(jié)點校正、推理的古典調和分析，選擇傳統(tǒng)的最小平方法求解青島港的調和常數(shù)值，并用Matlab加載和分析時間序列，對青島港每日不同時刻的潮位情況數(shù)據(jù)進行分析。其中t_tide.m用于分析，t_predic.m用于潮汐預報。具體步驟如下：

(1)數(shù)據(jù)預處理。對原始資料的奇異點進行濾波處理，去掉毛刺。在處理極少部分缺測資料時，利用GIL and DE TORO[21]介紹的最小二乘插值法進行插值，并使用格林尼治標準時。

(2)程序實現(xiàn)。首先，讀取原始潮位數(shù)據(jù)。接著，使用t_tide分析逐時的潮位資料，得到輸出參數(shù)分潮的符號(name)、振幅(amp)、遲角(pha)、信噪比(snr)、回報的潮位(xout)等。一般認為信噪比Vsnr>2的分潮是顯著的。

(3)潮汐預報。調和分析結束后，用t_predic進行潮汐的預報，得到不含計算數(shù)據(jù)平均值的預報潮位。

2 結果與分析

采用PAWLOWICZ et al.[16]編寫的T_TIDE軟件包進行驗潮站潮位資料的調和分析和預報。首先，分析青島驗潮站一號碼頭2019年全年潮汐水位數(shù)據(jù)的性質。其次，為研究潮位資料長度與分潮振幅、平均海平面的關系，將全年資料分成不同的時間長度，對比分析仿真結果的變化。之后，選擇不同分潮建立調和預報模型，通過2019年1月潮位的回報檢驗模型可靠性，確定最優(yōu)調和預報模型。最后，用改進的模型預報五號碼頭2019年8月的潮位變化，結合觀測資料檢驗調和預報效果。

2.1 調和分析潮汐的性質

選用青島港驗潮站2019年1月1日—12月31日全年間隔1 h的潮位資料進行調和分析，T_TIDE中設置開始時間為2019年1月1日，緯度為36.08°N，信噪比為2，并添加M10分潮(太陰淺水1/10日分潮)。各觀測站利用所有的逐時觀測節(jié)點，用K1分潮(太陰-太陽赤緯全日分潮)推P1分潮(太陽赤緯全日分潮)，用S2分潮(太陽主要半日分潮)推K2分潮(太陰-太陽赤緯半日分潮)，置信度為95%。而其他參數(shù)按文獻[16]中確定的取值不變，計算出該站點處的調和常數(shù)。使用青島港1 a長時間序列數(shù)據(jù)能夠消除各個主要分潮的相互影響，故能得到代表性良好的調和常數(shù)。最終得到67個分潮的振幅、遲角、信噪比以及年平均海平面(239 cm)，主要分潮的振幅貢獻如圖1所示。

圖1 青島港主要顯著分潮振幅Fig.1 Amplitude of main significant tidal components in Qingdao Port

由圖1可以看出，青島港潮汐以M2分潮(太陰主要半日分潮)為主，其振幅高達136.4 cm。其次是S2分潮和N2分潮(太陰主要橢率半日分潮)，振幅分別為42.0 cm 和 25.9 cm。K1分潮和O1分潮(太陰赤緯全日分潮)的振幅分別為 25.8 cm和 21.6 cm；SA分潮(太陽年周潮)的振幅為 20.1 cm； M4分潮(太陰淺水1/4日分潮)和 MS4分潮(太陰-太陽淺水1/4日分潮)的振幅為 13.6 cm和10.6 cm。

潮汐是各分潮疊加而形成的，實際中常用振幅較大的分潮M2、S2、K1、O1作為潮汐類型的判據(jù)[18]。潮汐類型的判別公式如下[22]：

(2)

(3)

(4)

(5)

其中，HK1、HO1和HM2，分別為分潮K1、O1和M2的平均振幅。在青島港口中，分潮S2的振幅小于M2，故可用此劃分標準。

通過公式(2)—(5)可以得到，該海域潮型數(shù)約為0.35，潮汐特征為規(guī)則半日潮，符合青島港屬于正規(guī)半日潮港的事實。

2.2 不同時段數(shù)據(jù)的調和分析

利用T_TIDE的調和分析程序，研究不同潮位資料長度與分潮振幅、平均海平面的關系。為探究月資料和季節(jié)資料(即時間段長度為3個月)調和分析結果的關系，將2019年的潮汐水位數(shù)據(jù)逐時資料分成12個月和3個季節(jié)(春季：3月1日—5月31日；夏季：6月1日—8月31日；秋季：9月1日—11月30日)。同時，進一步探究季節(jié)資料和整年資料調和分析結果的關系。由Matlab仿真實驗數(shù)據(jù)結果，繪制出12個月、3個季節(jié)資料長度的主要分潮對比圖(圖2、3)以及各月平均海平面高度對比圖(圖4)。此外，還得到了3個季節(jié)資料長度的平均海平面高度，分別為236 cm、257 cm、245 cm。

圖2 青島港各季節(jié)主要顯著分潮振幅比較Fig.2 Comparison of amplitude of main significant tidal components in each season in Qingdao Port

圖3 青島港各月主要顯著分潮振幅比較Fig.3 Comparison of amplitude of main significant tidal components in each month in Qingdao Port

圖4 青島港各月平均海平面比較Fig.4 Comparison of monthly average sea level in each month in Qingdao Port

由圖2和圖3知，3—11月和3個季節(jié)資料長度的仿真結果都主要受半日分潮 M2的影響，其次受分潮S2、N2、K1和O1影響。其中 M2分潮振幅在3月、4月、5月和春季分別為134.33 cm、136.43 cm、138.25 cm和136.04 cm，將3—5月結果與春季對比知，M2分潮振幅的平均差值為1.44 cm，其余分潮振幅的差值均在2.55 cm上下波動，其中存在2個較大差值，即在3月N2分潮處的7.77 cm和在5月S2分潮處的6.83 cm。對比6—8月和夏季的結果，其中 M2分潮振幅在6月、7月、8月和夏季分別為144.97 cm、135.44 cm、136.09 cm和139.53 cm，其余分潮振幅差值的平均值為4.51 cm，同時也出現(xiàn)了2個較大差值，即在6月N2分潮處的9.92 cm和在8月S2分潮處的10.30 cm。 對比9—11月和秋季的結果，在N2分潮處出現(xiàn)了1個較大差值。由圖1、2可知，3個季節(jié)和整年資料長度的結果也主要受分潮M2、S2的影響。對比其分潮振幅結果發(fā)現(xiàn)，其振幅的最大差值均出現(xiàn)在半日分潮S2處。由此可見，不同潮位資料長度對半日分潮S2、N2的振幅影響較大，但差值均小于11 cm，故所選取潮位資料長度對分潮振幅的影響較小。

青島港逐月平均海平面高度如圖4所示，3個季節(jié)資料長度的平均海平面高度分別為236 cm、257 cm、245 cm，年平均海平面高度為239 cm。逐月海平面高度表現(xiàn)為“偏峰型”特征，峰值出現(xiàn)在8月，與季節(jié)海平面高度峰值出現(xiàn)在夏季(6—8月)相契合，而此時正值青島地區(qū)的臺風期。9號臺風“利奇馬”于2019年8月11日17時在青島市南偏西方向約10 km的海面上(35.2°N，120.1°E)登陸，致使外圍最大風力有9級(風速達23 m·s-1)[23]。正是因為受臺風影響，不同潮位資料長度的平均海平面的最大差值出現(xiàn)在8月和夏季。

對比1個月、3個月和整年的分潮振幅及平均海平面高度數(shù)據(jù)可知，仿真結果與所選取潮位資料長度幾乎無關。但為建立最優(yōu)的預報模型，下文使用由整年潮位資料長度得到的仿真結果進行分析預報。

2.3 潮汐調和分析效果檢驗

理論上，在進行潮汐調和預報時選取的分潮個數(shù)越多，預測的精度越高。但實際應用中發(fā)現(xiàn)大部分分潮的振幅很小，因此可以忽略不計。為了定量分析分潮個數(shù)對調和回歸效果的影響，選取由2019年全年調和分析所得的不同分潮個數(shù)，通過分析2019年1月的預測結果，確定最優(yōu)的調和預報模型。

為評價不同預報模型仿真測試結果的性能，引入2個指標：(1)相對誤差(E)，指測量的絕對誤差(即預測潮水位與實測潮水位的差值)與實測的潮水位之比，它能更好地反映測量的可信程度；(2)判定系數(shù)(R2)，指統(tǒng)計學中“相關系數(shù)”的平方，用以表征兩個變量的相關程度。其具體結果由公式(6)和公式(7)計算所得。

(6)

(7)

采用增加、減少分潮數(shù)的方法建立最優(yōu)回歸模型，先將平均振幅從大到小進行排列，然后依次選取分潮數(shù)(5，10，15，……)進行仿真實驗。選取5個分潮時，調和預報相對誤差為 0.205 8，判定系數(shù)為 0.872 5；增加到 25 個分潮之后，調和預報相對誤差降到0.106 6，判定系數(shù)增至0.969 4，預報效果可以得到明顯提高；但是從 25 個分潮增加到50個時，則幾乎沒有改進(表1)。

表1 不同分潮個數(shù)調和預報的性能評價指標

但是，對于量化研究而言，時間序列數(shù)據(jù)中存在異常值，故必須要對相對誤差數(shù)據(jù)作預處理，使仿真數(shù)據(jù)更加客觀，以便于更好地探究潮水位的規(guī)律。本文采取去極值的方法，用以降低極值對相對誤差的影響，提升數(shù)據(jù)結論的準確性。選取5個分潮時，于1月25日15時出現(xiàn)一個極大值，數(shù)值為5.710 5,剔除后的相對誤差為0.155 0；當增加到25個分潮之后，同樣于1月25日15時出現(xiàn)一個極大值，數(shù)值為2.945 5，剔除后的相對誤差為0.089 6；但是從25個分潮增加到50個時，由去極值方法得到的相對誤差波動幅度很小。

為確切地得到最優(yōu)的調和預報模型，接著選取21、22、23和24個分潮進行調和預報，繪制出當分潮數(shù)為21、22、23和24個時的青島港潮汐預測圖(圖5)。由圖5可見，仿真結果和青島港口的實際潮位基本一致。對比其相對誤差，都在0附近上下起伏，但在16 h、40 h和132 h起伏較大。而16 h和40 h分別對應1月25日15時和26日15時，此時青島的風力達到5級，溫度日變量為7 ℃[24]，故潮汐預測精度可能受大風、涌浪等非潮汐因素的影響。

圖5 取21、22、23和24個分潮時的調和預報模型(a、c、e、g)及其相對誤差 (b、d、f、h)Fig.5 Harmonic forecast models (a/c/e/g) with 21, 22, 23, and 24 tidal components and their relative errors (b/d/f/h)

為更加直觀地對比數(shù)據(jù)，將圖5匯總成表2。其中，選取21個分潮時，由圖4可見，于148 h(1月31日03時)出現(xiàn)一個極大值，為14.404 1，剔除后的相對誤差為0.112 5；當增加到 22、23、24個分潮之后，由圖5可知，均于16 h (1月25日15時)出現(xiàn)一個極大值，分別為3.787 6、2.661 6、3.394 3，剔除后的相對誤差分別為0.078 2、0.084 4、0.068 3。

表2 取21、22、23和24個分潮時的調和預報的性能評價指標

由表2可知，當分潮個數(shù)由21個增加至24個時，調和預報模型的判定系數(shù)提高了0.021 9，相對誤差降低了0.109 5，去極值的相對誤差降低了0.044 2。但其由24個進行增加時，判定系數(shù)、相對誤差和去極值的相對誤差值幾乎不變。由此可見，使用24個分潮足以包含絕大部分的分潮信號，可以滿足調和預報的需求。由2019年1月的觀測潮位、調和預報潮位及其相對誤差曲線可知：建立模型時，應選取 24個分潮(即M2、S2、N2、K1、O1、SA、M4、K2、MS4、P1、L2、NU2、MN4、Q1、T2、MK4、MU2、LDA2、SSA、M6、2N2、2MS6、S4、NO1)，其精確的振幅數(shù)值可見圖1。

2.4 潮汐預報效果對比

五號碼頭坐落于青島港一號碼頭的西北方向，由東北向西南伸展共12個順岸泊位，五號碼頭線總長約2 174 m，水深約8 m。同上節(jié)分析，選擇24個分潮建立潮汐調和預報模型，分析預報五號碼頭2019年8月的潮位變化。使用相對誤差和判定系數(shù)來檢驗回報精度，回報結果及相對誤差見圖6。仿真結果(圖6)顯示，潮位的回報精度為0.941 6，相對誤差僅為0.082 7。這表明本文建立的模型可以較好地回報青島港附近海域的潮汐變化情況。

圖6 五號碼頭的調和預報模型(a)及其相對誤差(b) Fig.6 Harmonic forecast model (a) for Pier 5 and its relative error (b)

3 小結

對青島港口潮汐進行調和分析，使用T_TIDE潮汐調和分析預報程序計算其調和常數(shù)，進而對其進行分析和預報。青島港口潮汐性質為：潮汐為正規(guī)半日潮港，主要分潮有M2、S2、N2、K1、O1、M4、K2、MS4和P1，其中M2分潮尤其顯著。對比潮位資料長度為1個月、3個月及整年的調和分析結果，可知各分潮的振幅、平均海表面高度與用于調和分析的潮位資料長度幾乎無關。選擇不同分潮數(shù)會在一定程度上影響調和預報精度，對比實驗發(fā)現(xiàn)，選擇24個分潮的預報潮位和實測潮位的相對誤差為0.088 1，判定系數(shù)為0.969 8，此時為最優(yōu)選擇。

利用上述結果建立調和預測模型，分析五號碼頭2019年8月的潮位。仿真結果顯示，預測的潮位與真實的潮位基本一致，可較好地預報青島港附近海域的潮汐變化。由仿真結果分析知，潮汐預測誤差較大的情況大多是由于風、海浪、風暴潮等高頻自然因素運動造成的，故在潮汐預報時考慮非潮汐和自然等因素是值得研究的方向。