張鳳燁, 魏澤勛 , 王新怡, 王永剛, 方國洪

(1. 海洋環(huán)境科學(xué)與數(shù)值模擬國家海洋局重點實驗室, 山東 青島 266061; 2. 國家海洋局 第一海洋研究所,山東 青島266061)

潮汐調(diào)和分析方法的探討

張鳳燁1,2, 魏澤勛1,2, 王新怡1,2, 王永剛1,2, 方國洪1,2

(1. 海洋環(huán)境科學(xué)與數(shù)值模擬國家海洋局重點實驗室, 山東 青島 266061; 2. 國家海洋局 第一海洋研究所,山東 青島266061)

為實現(xiàn)對不等時距潮汐資料的分析, 基于 Matlab內(nèi)部函數(shù)功能, 提出了一種調(diào)和分析方法?；谶@種方法, 分別對大連、北海兩個站位1985年的全年等時間間距取樣的資料和非等時間間距取樣的資料進(jìn)行了調(diào)和分析, 結(jié)果顯示, 由等時間距資料和非等時間距資料計算的調(diào)和常數(shù)基本吻合。對大連、北海兩個站位的全年資料進(jìn)行多個不同時間間距取樣分析, 發(fā)現(xiàn)當(dāng)分潮頻率大于取樣頻率的二分之一時, 分潮發(fā)生頻率混淆。若分潮周期明顯大于樣品長度, 該分潮的分析結(jié)果產(chǎn)生很大誤差。最后得出結(jié)論為： 此調(diào)和分析方法, 適合對非等時間間距、非連續(xù)潮汐潮流資料進(jìn)行調(diào)和分析, 并且能夠獲得與傳統(tǒng)方法精度相當(dāng)?shù)慕Y(jié)果。

Matlab; 潮汐; 調(diào)和分析; 取樣

在實際的海洋調(diào)查中, 尤其是長時間的觀測,由于受到儀器故障、惡劣天氣、地理位置制約, 觀測方式等因素的影響, 很難得到從觀測初始時刻到結(jié)束時刻這段時間內(nèi)完整的高質(zhì)量數(shù)據(jù)資料, 也可以說獲得的海洋資料是不完整的、不連續(xù)的。對于潮汐資料, 調(diào)和分析方法作為一種主要的分析方法,發(fā)展至今, 理論基礎(chǔ)已經(jīng)很成熟, 但是傳統(tǒng)的調(diào)和分析大多是對觀測時間間隔為等距的資料進(jìn)行調(diào)和分析, 這是由于非等時序列數(shù)據(jù)的法方程系數(shù)矩陣很大, 所需的計算量也非常大, 通常要求觀測數(shù)據(jù)等時間間隔以減少計算量, 而且常用的觀測資料也是以觀測時間間隔為等距離為主。若資料中出現(xiàn)數(shù)據(jù)缺測漏測, 或者因為異常天氣和其他原因使得部分實測數(shù)據(jù)相對正常潮汐存在大的誤差, 就首先要對數(shù)據(jù)進(jìn)行前期處理[1], 然后應(yīng)用最小二乘法對潮汐進(jìn)行調(diào)和分析, 這種方法在計算時, 首先將原始矛盾方程進(jìn)行求導(dǎo)處理, 從而得到矛盾方程的法方程, 這樣大大增加了計算的難度[2]。如果我們從線性代數(shù)這門學(xué)科的角度來考慮, 矛盾方程組只是一個線性代數(shù)中的超定方程組, 只需要直接解這個超定方程組即可, 基于 Matlab內(nèi)部函數(shù)功能求解線性方程組能很好地解決這一問題, 并且方程組所得解即為最小二乘解[3-4], 從而簡化了調(diào)和分析的算法。有多少樣本, 就有多少個方程, 更能體現(xiàn)最小二乘法的靈活性, 完全不受等時或非等時取樣的限制。

本文采用基于Matlab內(nèi)部函數(shù)功能的最小二乘法實現(xiàn)潮汐資料的調(diào)和分析, 并應(yīng)用于大連、北海兩站位的水位觀測資料分析, 通過對兩站位全年資料的調(diào)和分析, 從非理論推導(dǎo)的角度, 分析基于Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法是否可應(yīng)用于等時間間隔和不等時間間隔的觀測資料的調(diào)和分析,以及觀測資料的取樣間隔對于調(diào)和分析精度的影響,然后通過對兩站位 7月份資料的調(diào)和分析, 來分析采樣長度對潮汐調(diào)和分析精度的影響。

1 改進(jìn)的潮汐調(diào)和分析方法介紹

調(diào)和分析是建立在分潮的概念上[5], 將潮汐看成是以不同頻率傳播的各種潮波疊加產(chǎn)生的現(xiàn)象,調(diào)和分析的目的就是求出各個分潮的振幅和遲角,遲角是指： 某時刻、某一地點實際的分潮的相角與理論上該時刻的分潮相角的差值, 遲角分為3種： 有區(qū)時遲角、地方遲角、格林威治遲角?，F(xiàn)在國際上通用的是格林威治遲角, 用字母g表示, 本文所用的遲角即為格林威治遲角[6]。

潮位表示式如下：

式中A0為平均水位高度,i為分潮序列號,ζ表示分潮潮高,ω為分潮的角速度,f,u分別表示由于月球軌道18.61 a的周期變化引進(jìn)的對平均振幅H和相角的訂正值, 即交點因子和交點訂正角,V0為格林威治初相角,H和g為分潮的調(diào)和常數(shù)。

本文利用潮位公式, 建立一種基于 Matlab內(nèi)部函數(shù)功能的潮汐調(diào)和分析方法, 其區(qū)別于傳統(tǒng)的調(diào)和分析方法, 其特殊之處在于, 這種基于 Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法可以對不等時間間隔的潮汐、潮流資料進(jìn)行調(diào)和分析, 而且計算易于實現(xiàn)。下面將介紹這種基于Matlab內(nèi)部函數(shù)功能的潮汐調(diào)和分析的方法和原理。

設(shè)有觀測資料：ζ(j),j=1,2,3…N,N代表數(shù)據(jù)個數(shù), 該資料每?t小時為一次記錄, 此處的?t既可以是整數(shù), 也可以是小數(shù)或者是一個不定值。取數(shù)據(jù)的開始時刻為時間原點。

設(shè)潮汐為M個分潮疊加而成：

由于每個分潮的ω是一個定值, 所以每個時刻的cosωt, sinωt是一個常數(shù), 這里A0是平均水位,我們可以把它看做是ω為0的一個分潮。這樣就形成了一個有2M+1個未知量的方程個數(shù)為N的線性方程組, 然后應(yīng)用基于 Matlab的內(nèi)部函數(shù)功能的最小二乘法, 直接對線性方程組進(jìn)行求解, 從而得最小二乘解, 即每個對應(yīng)分潮的A,B。

每個分潮的ω根據(jù)下式求得：

通過 Matlab可以求出一個平均水位、N/2個A值和N/2個B值。根據(jù)得出各分潮的振幅和相位。最后只要依據(jù)天文要素計算出初始時刻的f,u和V0, 此處f,u和V0為格林威治時刻值。根據(jù)即可求出調(diào)和常數(shù), 其中f,u的計算方法參照引潮力的第四展開式[1,8]。

2 站點水位的潮汐調(diào)和分析

下面將采用基于Matlab內(nèi)部函數(shù)功能的潮汐調(diào)和分析方法對大連、北海兩個站點等間隔資料、非等間隔資料以及不同時距的等時間隔進(jìn)行調(diào)和分析,討論本方法的適用性以及適用條件。

2.1 1 a資料的調(diào)和分析

利用基于Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法對兩個站點的1 a的水位資料進(jìn)行調(diào)和分析, 從潮族中選取具有代表意義的分潮, 選取了 30個分潮, 分別為 K1, O1, Q1, P1, MP1, J1、SO1, OO1, M1, M2, S2, N2,K2, OQ2, L2, O3, MO3, K3, MK3, M4, MS4, Mk4, SN4,2MK5, MSK5, M6, MSk6, 2Mk6, Sa, Ssa。

分別對大連、北海1 a(1985年)的觀測資料進(jìn)行調(diào)和分析, 這里原始觀測數(shù)據(jù)取樣的時間間隔是1 h, 而對于非等時的觀測資料的分析, 首先對原始數(shù)據(jù)進(jìn)行處理, 去掉 3月份的所有的觀測資料, 然后隨機刪除一些水位觀測資料, 使其達(dá)到一種非等時距采樣的效果, 然后對其進(jìn)行調(diào)和分析, 從圖 1可以看到, 等時調(diào)和分析預(yù)報的結(jié)果和不等時調(diào)和分析預(yù)報的結(jié)果與原始的觀測資料基本上是重合的, 從表1、表2可以看出, 等時間序列調(diào)和分析得到的調(diào)和常數(shù)與非等時間序列的調(diào)和常數(shù)的值非常接近, 表中提到的原始值代表的是利用傳統(tǒng)的調(diào)和分析方法對潮汐觀測數(shù)據(jù)進(jìn)行調(diào)和分析所得的調(diào)和常數(shù), 由于本文中提到的原始調(diào)和常數(shù)是對26 a數(shù)據(jù)進(jìn)行調(diào)和分析由135個分潮所得, 所以一些分潮的調(diào)和常數(shù)與本文計算的調(diào)和常數(shù)有一定的差別, 如因為M1分潮的杜德森數(shù)自身的問題, 所以其調(diào)和常數(shù)不予考慮, 但是主要分潮的調(diào)和常數(shù)與本文計算所得分潮調(diào)和常數(shù)相當(dāng)接近。這說明基于 Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法適用于潮汐資料的調(diào)和分析。

下面的實驗主要是探討數(shù)據(jù)的取樣間隔對該方法所得調(diào)和常數(shù)的影響, 為了從結(jié)果上討論取樣頻率是否對基于Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法有影響, 本文分別對大連、北海兩站位時長為1 a、取樣時間間隔分別為1, 2, 3, 4, 5, 6 h的觀測數(shù)據(jù), 應(yīng)用基于Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法進(jìn)行調(diào)和分析, 然后對調(diào)和常數(shù)進(jìn)行對比, 結(jié)果見表 3、表4、表 5、表 6。

圖1 1985年1月大連、北海潮位觀測數(shù)據(jù)、等時和非等時后報值圖Fig. 1 The chart of tide values, showing the observational and the hindcast values of harmonic analysis of data sampled at equal time intervals and non-equal time intervals in Dalian and Beihai in January 1985

表1 大連原始、等時與不等時資料調(diào)和常數(shù)比較Tab. 1 The comparison of harmonic analysis of data sampled at equal time intervals and non-equal time intervals in Dalian

表2 北海原始、等時與不等時資料調(diào)和常數(shù)比較Tab. 2 The comparison of harmonic analysis of data sampled at equal time intervals and non-equal time intervals in Beihai

在表 3、表 4、表 5、表 6中, 只列出了個別分潮, 本文引入了無單位量綱的負(fù)相關(guān)系數(shù)R, 表示潮汐預(yù)報值與原始觀測數(shù)據(jù)的吻合程度[7], 0≤R≤1,其值越大, 代表回歸效果越顯著, 計算公式如下：

上式中R為負(fù)相關(guān)系數(shù),U為后報值的離差平方和, 稱為回歸方差,S為觀測資料的總方差,為后報值,為后報值的平均數(shù),yi為觀測數(shù)據(jù),為觀測數(shù)據(jù)平均數(shù)。

從表7中能看出, 當(dāng)取樣間隔為6 h的時候, 變量與原始觀測數(shù)據(jù)的負(fù)相關(guān)系數(shù)為： 大連為 0.9468;北海為 0.9708, 同其余各取樣頻率結(jié)果的負(fù)相關(guān)系數(shù)相比, 相差不是很大。從表4、表6中, 當(dāng)取樣間隔為6 h時, S2分潮遲角一直是60°, 從中分析主要原因為 S2的圓頻率為 30°, 為了進(jìn)一步驗證這個假設(shè),當(dāng)取樣間隔為12, 18 h時, 其遲角都是60°/ h, 尤其是其取樣間隔為12 h的整數(shù)倍時, 此程序無法回報潮位, 也就是說當(dāng)圓頻率與取樣間隔的乘積為 180°的整數(shù)倍時, 這種基于 Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法將是無效的。從表3、表4、表5、表6中不同時間間隔的各分潮的調(diào)和常數(shù)可以總結(jié)出, 隨著取樣時間的增加, 高頻分潮的調(diào)和常數(shù)誤差變大,這就涉及取樣間隔與分潮分辨問題, 如果采樣間隔?t取得大, 則采樣頻率fs(fs=1/?t)低, 當(dāng)所分析信號的最高頻率fmax大于采樣頻率fs的二分之一時, 就會引起“頻率混淆”現(xiàn)象, 使得原信號中的頻率成分出現(xiàn)在數(shù)字信號中完全不同的頻率處, 造成信號的失真。

表3 大連原始振幅與各取樣間隔不同的觀測資料調(diào)和分析振幅Tab. 3 The original amplitudes and amplitudes of the harmonic analysis of data sampled at various time intervals in Dalian

表4 大連原始遲角與各取樣間隔不同的觀測資料調(diào)和分析遲角Tab. 4 The original phases and phases of the harmonic analysis of data sampled at various time intervals in Dalian

表5 北海原始振幅與各取樣間隔不同的觀測資料調(diào)和分析振幅Tab. 5 The original amplitudes and amplitudes of the harmonic analysis of data sampled at various time intervals in Beihai

表6 北海原始遲角與各取樣間隔不同的觀測資料調(diào)和分析遲角Tab. 6 The original phases and phases of the harmonic analysis of data sampled at various time intervals in Beihai

表7 大連、北海兩站點負(fù)相關(guān)系數(shù)Tab. 7 The negative correlation coefficients in Dalian and Beihai

取樣間隔為1 h, 能分辨出圓頻率小于180°/s的分潮; 取樣間隔為2 h, 能分辨出圓頻率小于90°/s的分潮; 取樣間隔為3 h, 能分辨出圓頻率小于60°/s的分潮。依照σ≤π/?t公式, 以此類推。

因為取樣間隔為1 h和2 h時, 都能分辨出圓頻率小于90°/s的分潮, 由于本文用到的30個分潮的圓頻率都小于90°/s, 從表 3、表4、表5、表6中可以看到, 當(dāng)間隔是1 h和2 h時, 分潮的調(diào)和常數(shù)基本上是一樣的, 但隨著取樣間隔的增大, 像M6, MS4這樣的高頻分潮就會出現(xiàn)誤差, 即發(fā)生頻率混淆, 高頻分潮在調(diào)和分析的過程中不能被分辨。

2.2 1個月資料的調(diào)和分析

對1個月資料的調(diào)和分析, 樣品取樣間隔為1 h,同樣用30個分潮： K1, O1, Q1, P1, MP1, J1, SO1, OO1,M1, M2, S2, N2, K2, OQ2, L2, O3, MO3, K3, MK3, M4,MS4, MK4, SN4, 2MK5, MSK5, M6, MSK6, 2Mk6, Sa,Ssa。選取的月份為 7月份, 選擇的站點同樣是大連和北海。

應(yīng)用該方法對潮位進(jìn)行調(diào)和分析時, 結(jié)果顯示,在大連, Sa的振幅和遲角的值, 以及Ssa的振幅和遲角的值同樣不正確, 變量與原始觀測數(shù)據(jù)的負(fù)相關(guān)系數(shù)為0.9911。在北海, Sa的振幅和遲角的值, 以及Ssa的振幅和遲角的值明顯不正確, 但是變量與原始觀測數(shù)據(jù)的負(fù)相關(guān)系數(shù)為0.9723。該結(jié)果說明, 雖然此種方法應(yīng)用于潮汐的調(diào)和分析時, Ssa, Sa兩個分潮的調(diào)和常數(shù)明顯不正確, 但是對于原始潮位資料與后報潮位的擬合, 該方法還是可以實現(xiàn)的。為了說明這個問題, 我們在調(diào)和分析的過程中, 去掉 Sa,Ssa兩個分潮, 利用其余的 28個分潮對其進(jìn)行調(diào)和分析, 結(jié)果發(fā)現(xiàn), 在大連, 變量與原始觀測數(shù)據(jù)的負(fù)相關(guān)系數(shù)為0.9910。在北海, 變量與原始觀測數(shù)據(jù)的負(fù)相關(guān)系數(shù)為0.9926; 從相關(guān)系數(shù)上判斷, 28個分潮與30個分潮調(diào)和分析后報值基本上差不多, 甚至30個分潮擬合的結(jié)果更加精確(圖2)。從Sa, Ssa振幅可以看出, 該方法沒有能分辨出這兩個分潮, 其主要原因是因為Sa, Ssa是周期約為1 a的長周期分潮, 1個月的資料根本就不能得到他們的調(diào)和常數(shù), 也就是說,在調(diào)和分析的過程中, 分潮的選取與取樣的長度有很大的關(guān)系, 最根本的條件就是取樣的時間長度必須大于所選分潮的周期。

圖2 1985年7月大連, 北海1個月潮位觀測值、28個分潮和30個分潮調(diào)和后報值圖Fig. 2 The chart of monthly observation values and hindcast values in Dalian and Beihai in July 1985

3 結(jié)論

應(yīng)用基于Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法對已知觀測點水位進(jìn)行調(diào)和分析時, 不需要尋找觀測資料時間序列的中間時刻來實現(xiàn)對計算的簡化,從初始時刻計算即可。其在回報時, 可以采用各個相應(yīng)時刻的f和u, 使后報結(jié)果稍加精確?；贛atlab內(nèi)部函數(shù)功能的調(diào)和分析方法與傳統(tǒng)調(diào)和分析方法相比較, 前者更加簡便, 其更顯著的優(yōu)點是此方法可以對任意的觀測時間間隔的潮汐水位進(jìn)行調(diào)和分析以及可處理數(shù)據(jù)漏測或間斷情況下的潮汐資料。

通過對大連、北海1985年全年資料的調(diào)和分析結(jié)果顯示, 基于 Matlab內(nèi)部函數(shù)功能的調(diào)和分析方法可以應(yīng)用于等時間序列和非等時間序列的潮汐資料的調(diào)和分析, 而且得到的調(diào)和常數(shù)的精確性比較高, 其預(yù)報所得的潮位與實際觀測潮位也比較吻合,這說明此種方法完全是可行的。對于不同的取樣間隔的觀測資料, 隨著取樣間隔的增加, 對高頻的分潮的分辨能力會越來越低, 出現(xiàn)頻率混淆的現(xiàn)象。對于 1個月的觀測資料的調(diào)和分析, 調(diào)和分析得到的Ssa, Sa等長周期分潮的調(diào)和常數(shù)異常, 說明1個月的資料不能分析出周期大于 1個月的分潮, 即觀測數(shù)據(jù)的取樣長度必須要大于分潮的1個周期。

本文只是對基于Matlab內(nèi)部函數(shù)的調(diào)和分析方法是否可以應(yīng)用于潮汐調(diào)和分析上的一個初步討論,接下來的工作將是針對不同長度、不同取樣間隔的潮汐資料對分潮的選取上做進(jìn)一步的深入了解, 從而使其調(diào)和分析的精度更加精確, 使潮位的回報值更加接近原始資料。

[1]方國洪, 鄭文振, 陳宗鏞, 等. 潮汐和潮流的分析和預(yù)報[M]. 北京：海洋出版,1986： 90-93.

[2]王如龍, 童章龍, 陳耀登, 等. 基于連續(xù)函數(shù)最小二乘法的潮汐迭代調(diào)和分析方法[J]. 中國水利, 2007,7(3)： 116-118.

[3]同濟(jì)大學(xué)應(yīng)用數(shù)學(xué)系. 線性代數(shù)[M]. 北京： 高等教育出版社, 2004.

[4]David C L. Linear Algebra and Its Applications[M]. 北京： 機械工業(yè)出版社,2005.

[5]馮士笮,李鳳岐,李少菁, 等. 海洋科學(xué)導(dǎo)論[M].北京： 高等教育出版社, 1999： 215-222.

[6]陳宗鏞. 潮汐學(xué)[[M]. 北京： 科學(xué)出版社, 1980：127-144.

[7]施能. 氣象科研與預(yù)報中的多元分析方法[M]. 北京：氣象出版社, 2002. 2： 34-37.

[8]黃祖珂, 黃磊. 潮汐原理與計算[M]. 青島： 中國海洋大學(xué)出版社, 2005： 58-63.

Tidal harmonic analysis

ZHANG Feng-ye1,2, WEI Ze-xun1,2, WANG Xin-yi1,2, WANG Yong-gang1,2, FANG Guo-hong1,2
(1. Key Laboratory of Marine Science and Numerical Modeling, State Oceanic Administration, Qingdao 266061, China, 2. First Institute of Oceanography, State Oceanic Administration, Qingdao 266061, China)

Jan., 24, 2011

Matlab; tide; harmonic analysis; sample

To analyze observational data sampled at non-equal time intervals, we devised a least squares method of tidal harmonic analysis, based on Matlab’s function. By this method, the harmonic constants were almost the same with observational data sampled at either equal time intervals or and non-equal time intervals at Dalian and Beihai in 1985. The harmonic analysis with data sampled at different time intervals showed that frequency aliasing occured when the tidal frequency was more than 1/2 of the sampling frequency, and there were great error when the period of the constituent was significantly larger than the sample length. It was shown that our method could obtain results of the same precision of the traditional method in analyzing discontinuous data or data sampled at non-equal time intervals.

P731.23

A

1000-3096(2011)06-0068-08

2011-01-24;

2011-04-21

國家自然科學(xué)基金項目(40976016); 國家海洋局第一海洋研究所基本科研業(yè)務(wù)費專項資金項目(GY02-2010G07); 國家高技術(shù)研究發(fā)展計劃(863計劃)重點項目(2009AA121402)

張鳳燁(1985-), 男, 河北唐山人, 碩士研究生, 主要從事物理海洋學(xué)研究, E-mail： zhangfy@fio.org.cn; 魏澤勛, 通信作者, 研究員, E-mail： weizx@fio.org.cn

劉珊珊)