劉倩，孫小淇

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

引言

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

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

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

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

1.1 調(diào)和分析法

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

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

(1)

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

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

1.2 資料來源與處理方法

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

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

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

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

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

2 結(jié)果與分析

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

2.1 調(diào)和分析潮汐的性質(zhì)

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

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

由圖1可以看出，青島港潮汐以M2分潮(太陰主要半日分潮)為主，其振幅高達(dá)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。

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

(2)

(3)

(4)

(5)

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

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

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

利用T_TIDE的調(diào)和分析程序，研究不同潮位資料長度與分潮振幅、平均海平面的關(guān)系。為探究月資料和季節(jié)資料(即時(shí)間段長度為3個(gè)月)調(diào)和分析結(jié)果的關(guān)系，將2019年的潮汐水位數(shù)據(jù)逐時(shí)資料分成12個(gè)月和3個(gè)季節(jié)(春季：3月1日—5月31日；夏季：6月1日—8月31日；秋季：9月1日—11月30日)。同時(shí)，進(jìn)一步探究季節(jié)資料和整年資料調(diào)和分析結(jié)果的關(guān)系。由Matlab仿真實(shí)驗(yàn)數(shù)據(jù)結(jié)果，繪制出12個(gè)月、3個(gè)季節(jié)資料長度的主要分潮對比圖(圖2、3)以及各月平均海平面高度對比圖(圖4)。此外，還得到了3個(gè)季節(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個(gè)季節(jié)資料長度的仿真結(jié)果都主要受半日分潮 M2的影響，其次受分潮S2、N2、K1和O1影響。其中 M2分潮振幅在3月、4月、5月和春季分別為134.33 cm、136.43 cm、138.25 cm和136.04 cm，將3—5月結(jié)果與春季對比知，M2分潮振幅的平均差值為1.44 cm，其余分潮振幅的差值均在2.55 cm上下波動(dòng)，其中存在2個(gè)較大差值，即在3月N2分潮處的7.77 cm和在5月S2分潮處的6.83 cm。對比6—8月和夏季的結(jié)果，其中 M2分潮振幅在6月、7月、8月和夏季分別為144.97 cm、135.44 cm、136.09 cm和139.53 cm，其余分潮振幅差值的平均值為4.51 cm，同時(shí)也出現(xiàn)了2個(gè)較大差值，即在6月N2分潮處的9.92 cm和在8月S2分潮處的10.30 cm。 對比9—11月和秋季的結(jié)果，在N2分潮處出現(xiàn)了1個(gè)較大差值。由圖1、2可知，3個(gè)季節(jié)和整年資料長度的結(jié)果也主要受分潮M2、S2的影響。對比其分潮振幅結(jié)果發(fā)現(xiàn)，其振幅的最大差值均出現(xiàn)在半日分潮S2處。由此可見，不同潮位資料長度對半日分潮S2、N2的振幅影響較大，但差值均小于11 cm，故所選取潮位資料長度對分潮振幅的影響較小。

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

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

2.3 潮汐調(diào)和分析效果檢驗(yàn)

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

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

(6)

(7)

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

表1 不同分潮個(gè)數(shù)調(diào)和預(yù)報(bào)的性能評價(jià)指標(biāo)

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

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

圖5 取21、22、23和24個(gè)分潮時(shí)的調(diào)和預(yù)報(bào)模型(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個(gè)分潮時(shí)，由圖4可見，于148 h(1月31日03時(shí))出現(xiàn)一個(gè)極大值，為14.404 1，剔除后的相對誤差為0.112 5；當(dāng)增加到 22、23、24個(gè)分潮之后，由圖5可知，均于16 h (1月25日15時(shí))出現(xiàn)一個(gè)極大值，分別為3.787 6、2.661 6、3.394 3，剔除后的相對誤差分別為0.078 2、0.084 4、0.068 3。

表2 取21、22、23和24個(gè)分潮時(shí)的調(diào)和預(yù)報(bào)的性能評價(jià)指標(biāo)

由表2可知，當(dāng)分潮個(gè)數(shù)由21個(gè)增加至24個(gè)時(shí)，調(diào)和預(yù)報(bào)模型的判定系數(shù)提高了0.021 9，相對誤差降低了0.109 5，去極值的相對誤差降低了0.044 2。但其由24個(gè)進(jìn)行增加時(shí)，判定系數(shù)、相對誤差和去極值的相對誤差值幾乎不變。由此可見，使用24個(gè)分潮足以包含絕大部分的分潮信號，可以滿足調(diào)和預(yù)報(bào)的需求。由2019年1月的觀測潮位、調(diào)和預(yù)報(bào)潮位及其相對誤差曲線可知：建立模型時(shí)，應(yīng)選取 24個(gè)分潮(即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 潮汐預(yù)報(bào)效果對比

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

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

3 小結(jié)

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

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