海洋要素计算与预报............................................................................................ 1 第一部分 数据预处理与统计分析方法.............................................................. 1
第一章 数据预处理...................................................................................... 1
一、数据质量控制................................................................................ 1
1、异常数据的认定和排除.......................................................... 1 2、数据系统性偏差的检查和修正.............................................. 1 二、不规则空间分布数据网格化........................................................ 1
1、 数学插值法 ............................................................................. 1 2、网格统计法.............................................................................. 2 三、要素统计特征................................................................................ 3
1、要素数据标示.......................................................................... 3 2、均值与距平.............................................................................. 3 3、平均差...................................................................................... 3 4、方差.......................................................................................... 3 5、协方差与相关系数.................................................................. 3 6、自协方差与自相关系数.......................................................... 3 7、落后协方差与相关系数.......................................................... 4 8、经验分布.................................................................................. 4
第二章 谱分析.............................................................................................. 5
一、Fourier变换与谱分析 ................................................................... 5 二、功率谱估计.................................................................................... 6 三、交叉谱分析.................................................................................... 7 第三章 经验模态分解.................................................................................. 8
一、前言................................................................................................ 8 二、EMD计算方法与IMF分量 ........................................................ 9 三、EMD方法中存在的问题 ............................................................ 11
1、EMD方法在处理间歇信号时的不可分问题和产生的模态混合问题.................................................................................................. 11
2、EMD分解方法的边界问题 .................................................. 15 四、应用实例...................................................................................... 17
1、SST资料处理 ........................................................................ 17 2、海平面数据处理.................................................................... 17
第四章 回归分析........................................................................................ 18
一、一元线性回归.............................................................................. 19
1、一元线性回归模型................................................................ 19 2、一元线性回归的方差分析.................................................... 19
~ 1 ~
海洋要素计算与预报
3、回归方程的显著性检验........................................................ 20 4、预报值的置信区间................................................................ 20 二、多元线性回归.............................................................................. 21
1、多元线性回归模型................................................................ 21 2、回归方程显著性检验............................................................ 22 3、预报值的置信区间................................................................ 22 三、非线性回归.................................................................................. 23
1、曲线函数线性化.................................................................... 23 2、多项式回归............................................................................ 23
第五章 经验正交函数分解........................................................................ 23
一、主成分的定义.............................................................................. 24
1、两个变量的主成分定义........................................................ 24 2、多变量的主成分定义............................................................ 25 二、主成分的导出.............................................................................. 26 三、主成分的性质.............................................................................. 27 四、主成分的计算.............................................................................. 28 五、经验正交函数分解(EOF) ........................................................... 28 六、时空转换...................................................................................... 29 第六章 最小二乘法潮汐调和分析与潮汐特征值.................................... 30
一、分潮与潮汐调和常数.................................................................. 30 二、最小二乘法潮汐调和分析方法.................................................. 32
1、任意时间间隔观测序列的方程组导出................................ 32 2、等时间间隔观测序列的方程组系数.................................... 34 3、Fourier系数的计算 ............................................................... 35 4、天文变量与调和常数计算.................................................... 36 三、潮流调和常数与潮流椭圆要素.................................................. 42 四、潮汐性质与潮汐特征值.............................................................. 43
1、潮汐性质................................................................................ 43 2、潮汐特征值............................................................................ 43 3、平均海面、平均海平面与陆地高程,海图深度基准面与海图水深.................................................................................................. 45
(4)海图深度基准面与海图水深............................................ 45
第七章 海浪数据分析................................................................................ 48
一、去倾向和去均值处理.................................................................. 48 二、从波面高度序列中读取海浪的波高和周期.............................. 48
1、跨零点波高、周期定义........................................................ 48
~ 2 ~
海洋要素计算与预报
2、极值点波高、周期定义........................................................ 49 三、波面高度分布、波高和周期的分布,波高和周期的联合分布...................................................................................................................... 49
1、波面高度分布........................................................................ 49 2、波高和周期的分布................................................................ 50 四、各种波高计算.............................................................................. 51 五、海浪谱估计.................................................................................. 52
1、海浪谱估计方法.................................................................... 52 2、谱矩的计算............................................................................ 52 3、谱的零阶矩与各种波高的关系............................................ 52 4、海浪谱的谱宽度计算............................................................ 52 5、谱峰频率与周期的关系........................................................ 53
第二部分 海洋数值预报.................................................... 错误!未定义书签。
第一章 二维浅水方程数值模式(潮汐、风暴潮数值计算模式)错误!未定义书签。
一、基本方程组.................................................. 错误!未定义书签。 二、ADI方法 ..................................................... 错误!未定义书签。 第二章 三维温盐流预报模式.................................... 错误!未定义书签。
一、POM模式 .................................................... 错误!未定义书签。 二、MOM模式................................................... 错误!未定义书签。 第三章 海浪谱预报模式............................................ 错误!未定义书签。
一、文模式.......................................................... 错误!未定义书签。 二、WAM模式 ................................................... 错误!未定义书签。 三、SWAN模式 ................................................. 错误!未定义书签。
~ 3 ~
海洋要素计算与预报
海洋要素计算与预报
第一部分 数据预处理与统计分析方法
第一章 数据预处理 一、数据质量控制 1、异常数据的认定和排除
观测和实验测量数据在大多数情况下,都含有一定数量的异常数据,这些数据往往由于人为错误(比如:输入错误、记录笔误,多见于需人工记录的仪器观测和历史资料的数字化)和仪器工作异常(由短时强电磁干扰、仪器位置非永久性短时变动引起)而产生。此种异常数据具有以下特征:单个或短时连续、随机出现、违反数据序列的整体变化规律。
可采用下列方法判别:
1 阈值判别法:根据水文要素可能的取值范围设定相应的阈值,超过阈值○
的数据即可确认为异常数据或列为可疑异常数据。
2 变率判别法:水文要素的时间变化率应在一定的范围内,并具有某种规○
律性。可利用这一特性对异常观测值进行判别。
异常数据的处理:根据水文要素序列的特征和分析要求,采取简单剔除或进行插值处理。
2、数据系统性偏差的检查和修正
由于观测仪器本身设计原因,在观测过程中仪器的零点发生自动漂移,由此而产生观测数据的系统性偏差,此种偏差一般具有单调性,一般可采用线性修正方法处理。
由于仪器安装位置的不可恢复性缓慢变化或跳动,而产生观测数据的系统偏差。 对于仪器位置缓慢变化引起的偏差,其表现形式类似于仪器的零点漂移,但情况会更加复杂,可能会具有较强的阶段性,其修正过程也更加复杂。对于仪器位置跳动所引起的偏差,一般对数据进行零点修正即可。
数据系统性偏差比较难于判别,尤其当变动值较小时,只有清晰了解要素的变化规律才可能加以判别。 二、不规则空间分布数据网格化 1、 数学插值法
采用数学方法根据插值点周围数据计算出插值点上的数据,从而完成数据的网格化。在已知数据覆盖范围的边缘区域,插值的可信度可能会显著下降。适用于同步和准同步大面积观测数据。
~ 1 ~
海洋要素计算与预报
1 距离权重插值 ○
2 线性曲面插值 ○
3 样条曲面插值 ○
11109876543024681012 2、网格统计法
在数据区域按等间距划分网格,落入网格的数据(0~N个)取平均值,同时计算均方差。均方差和数据个数可作为概点统计数据的质量指标。此方法适用于较长时间段内的时间平均场的计算。
3231.53130.53029.52928.52827.527130130.2130.4130.6130.8131131.2131.4131.6131.8132
~ 2 ~
海洋要素计算与预报
三、要素统计特征 1、要素数据标示
某个水文要素记为x,其观测值的一个时间序列,记为:
2,„,n)。 x1,x2,„,xn或{xt},(t1,2、均值与距平
1n均值:xXt,反映要素观测值的平均状况。
nt1距平:xtxtx, (t1,2,„n),也称为离均差,反映个别数据偏离平均值的大小,用于分析要素的变异规律。 3、平均差
1n平均绝对差:va|xtx|,简称平均差,表示要素的变化情况。
nt14、方差
1n方差:s(xix)2
ni12均方差:s1n(xix)2,亦称标准差。方差和均方差描述要素观测值围ni1绕其平均值的变化情况,即要素观测值的离散程度。 5、协方差与相关系数
2,„,n), 设两个水文要素同步观测样本序列:{xt}、{yt},(t1,1n协方差:sxy(xix)(yiy),描述了两变量之间的相互依赖关系。
ni11nxtxyty相关系数:rxy()(),1r1,即标准化变量的协方差,
nt1sxsy可用来度量两变量之间的相关关系。
xtx},标准化变量具有均值为0,方差为1的s{xt}的标准化变量表示为{特性。因此,经过标准化后的水文要素便于相互比较。 6、自协方差与自相关系数
衡量水文要素不同时刻之间的关系密切程度的量是自协方差和自相关系数。对于{xt},则时间间隔为(t2t1,t2t1)的
1自协方差: s()n
(xi1nix)(xix),
~ 3 ~
海洋要素计算与预报
1自相关系数:r()n(t1nxtxxtx)()。 ss称为落后相关步长。为正整数时,称为落后相关系数;反之则为超前相
关系数。
7、落后协方差与相关系数
衡量两个要素不同时刻之间的相关程度的量,常用落后交叉协方差和落后交叉相关系数表示。设{xt}及{yt}分别为两个水文要素的观测时间序列,则时间间隔为(t2t1,t2t1)的
1落后交叉协方差:sxy()n(xi1nix)(yiy),
1落后相关系数:rxy()n8、经验分布
(t1nxtxyty)()。 sxsy水文要素的变化规律还可以用其概率分布来描述它的各种观测值出现的可能性大小,即可计算给出每一个观测值出现的概率。
(1) 频率分布
若{xt}为水文要素的连续观测时间序列,由连续型随机变量的概率性质知道,不应统计各个观测值出现的频率,而应统计它在某—范围内出现的频率。为此,可把观测值的变化范围划分为许多连续但不重叠的区间(每一个区间称为组),统计观测值落在各区间内的频率,即可得到要素的经验频率分布:
fi(x)100mi nxixi1称为组中值或区间中值,代表该频率的观测2其中,mi为xixtxi1的个数;n为观测样本总数;xi为区间端点亦称上下组限;xi1xi称为组距;
值。可以设想,若观测值无限增多,组距逐渐缩小,那么频率密度曲线将趋向概率密度曲线。
(2) 累积频率分布
设{xt}是水文要素的观测值时间序列,将它们按由小到大顺序排列为(X1,X2,,Xn),若以F(x)P(x)表示上述观测值中气象要素的值小于数x的频率,则函数:
0,mF(x),n1,
xX1XmxXm1 称为观测值{xt}的累积频率分布函数,也称为xXn~ 4 ~
海洋要素计算与预报
样本累积频率分布函数,简称累积频率分布,在实际使用中也常称为不及累积频率分布。
如将{xt}由大到小顺序排列,同样可得到函数:
0,mP(x)P(x),n1,率。
xX1XmxXm1称为超越累积频率分布,常称为保证xXn两种累积频率分布满足以下关系:
P(x)F(x)1。
另外,由经验频率分布按观测值进行由大到小或由小到大的频率积分,即可相应得到超越累积频率分布或不及累积频率分布。需注意的是计算超越累积频率分布时观测值应取组下限,计算不及累积频率分布时观测值应取组上限。
第二章 谱分析
一、Fourier变换与谱分析
对于在周期T内的时间函数f(t)[{xi}(i1,2,N)为其时间间隔为1的N个观测值],在一般情况下可以展成如下形式的Fourier级数:
f(t)a0(akcosktbksinkt) (1)
k1其中,k2k/T(k2k/N)
k为谐波的波数,k为第k各波的角速度,Tk为周期,fk为频率。fkk,2TkNT12(Tk)。 kkfkk
(1)式中的a0、ak和bk称为Fourier系数:
1a0N2akNNx,
tt1Nxtcost12k(t1), N2bkN
xtsint1N2k(t1)。 N其中,k1,2,,N/2(只取整数部分)。
~ 5 ~
海洋要素计算与预报
(1)式可改写为以下形式:
f(t)kCeikkt,
Ck称为f(t)的复谱,它是一个复数。
CkAkeik,
Ak12akbk2, 2ktg1bk ak通常Ak称为振幅谱,k称为位相谱,它们随波数k变化,每个波数的谐波又对应于相应的频率,故可用k为横坐标,以Ak、k为纵坐标作图,这种图称为谱图。
Ck与f(t)具有如下转换关系:
1CkTT/2T/2f(t)eikdt
f(t)a0
kCkeikt
Fourier变换与反变换的一般数学描述:
(j1)(k1)Fourier变换:X(k)x(j)N
j1N1Fourier反变换x(j)NX(k)k1N(j1)(k1)N
其中:Ne(2i)/N,x(j)为时间序列,N为时间序列总长度。
一般情况下采用快速Fourier变换计算方法进行Fourier变换与反变换,快速Fourier变换简写为FFT(Fast Fourier Transform),反变换写为IFFT(Inverse Fast Fourier Transform)。
二、功率谱估计
按照一般的物跟模型,若其电阻为1个单位,瞬时电压用f(t)表示,则它的瞬时功率为f2(t),它的总能量为:
f2(t)dt
从统计观点来看,上式就表示以数学期望为零的f(t)的总方差。对于这种总
~ 6 ~
海洋要素计算与预报
能量或总方差,可以分解为各波动的部分能量或方差之和,功率谱分析就是在这种概念下产生的。
时间函数f(t)表示为:
f(t)iktCek k平均功率为:
1T/22f(t)|Ck|2 TT/2k12bk2)。 其中,|Ck|2(ak42由于|Ck2|与|Ck|相同,因此,通常离散功率谱图仅绘出正半轴的谱线
(k0),故令
Sk22|Ck|212(akbk2), 2则平均功率可写为:
1T/222f(t)2|Ck|Sk2。 TT/2k1k1实际计算时,首先对时间序列{xi}(i1,2,N)进行Fourier变换,从而得到Fourier系数ak、bk,k1,2,N/2(只取整数),然后计算Sk2。利用波数与周期的关系Tk功率谱图。
三、交叉谱分析
交叉谱分析是揭露两个时间序列在不同频率上互关系的一种分析方法。 设有时间序列{x1(t)}和{x2(t)},对应的复谱为F1()、F2(),则交叉谱为:
N计算出每一个组成波的周期,即可以Tk为横轴,Sk2为纵轴画出kS12()F2()F1*()
其中,F1*()为F1()的共扼复谱。它是一个复谱,可写成实部与虚部之和的形式:
S12()P12()iQ12()
其中:
P12(k)Q12(k)1(a1ka2kb1kb2k),称为协谱或共谱; 41(a1kb2ka2kb1k),称为正交谱或余谱。 4 ~ 7 ~
海洋要素计算与预报
a1k、b1k、a2k、b2k分别为{x1(t)}和{x2(t)}的Fourier系数。
协谱反映两个时间函数在某一频率上的同位相相关程度。
正交谱反映两个时间函数在某一频率上的位相相差90度时的相关关系。
凝聚谱:反映两个时间函数在某一频率上的总体相关程度。 类似于相关系数,定义凝聚谱如下:
|S12()|2, R()S11()S22()212也可写成为
22P12()Q12()。 R()P11()P22()212
位相谱:反映两个时间函数各种频率波动的位相差关系,因此又称为位相差谱。
Q12() P()1212()tg1若12()0,表示{x2(t)}在频率上位相大于{x1(t)}的位相,即{x1(t)}落后于{x2(t)}变化;反之,若12()0,表示在频率为的振动上{x1(t)}超前于
{x2(t)}变化。
与位相谱相对应,定义落后时间长度谱如下:
12()T 2其中,T为对应于的周期,L的单位与周期相同为时间。
本章参考书:《气象统计原理与方法》、《气象统计分析与预报方法》、《时间
L(T)序列及其谱分析》
第三章 经验模态分解 一、前言
1998年美籍华人N.E. Huang [1998]等人提出了一种新的信号处理方法——
~ 8 ~
海洋要素计算与预报
经验模态分解方法(Empirical Mode Decomposition),简称EMD方法。次年,Huang [1999]又将该方法进行了一些改进。该方法从本质上讲是对一个信号进行平稳化处理,其结果是将信号中真实存在的不同尺度波动或趋势逐级分解开来,产生一系列具有不同特征尺度的数据序列。每一个序列称为一个本征模函数(Intrinsic Mode Function),简称IMF分量。Wang [2000]等人对该方法作了进一步的改进,使EMD方法的尺度分解精度大大提高。
由EMD方法得到的最低频率IMF分量通常情况下代表原始信号的趋势或均值。从这个意义上讲,EMD分解方法可以有效地提取一个信号的趋势或去掉该信号的均值。进行EMD分解的主要目的是为了进一步对信号进行处理分析。
EMD方法自推出以来已经成功地应用在湍流研究、地震研究以及非线性研究等许多领域。可以预期,在不远的将来,该方法必将在更多的研究领域中发挥巨大的作用。
二、EMD计算方法与IMF分量
EMD分解方法是一种平稳化处理过程。其基本思想是:假如一个原始数据序列X(t)的极大值或极小值数目比上跨零点(或下跨零点)的数目多2个(或2个以上),则该数据序列就被认为是非平稳的,需要进行平稳化处理。这种非平稳特性是由不同尺度波动叠加造成的,EMD方法的目的就是将这些不同的尺度分离。图2.1所示的是由两个不同频率波动叠加而形成的非平稳信号,该信号(sig1)的跨零点数目比极值点数目少两个以上。其中:
22X(t)cost0.6sint2004021 (2.1)
X(t)0-1-20100200300t400500600图2.1 非平稳的信号 EMD方法的具体处理过程是:找出X(t)所有的极大值点并将其用三次样条函数拟合成原数据序列的上包络线;找出所有的极小值点并将其用三次样条函数拟合成原数据序列的下包络线;上下包络线的均值为原数据序列的平均包络线
m1;将原数据序列减去该平均包络后即可得到一个去掉低频的新数据序列h1
~ 9 ~
海洋要素计算与预报
X(t)m1h1 (2.2)
图2.2给出了上述处理过程的结果。图2.2a中细实线为原信号,点线为信号的上包络线,点划线为信号的下包络线,粗实线为平均包络线m1。图2.2b给出的是新数据序列h1的图形。
一般来讲,为此需要对它重复上述处理过程。h1不一定是一个平稳数据序列,假如h1序列的平均包络线为m11,去除该包络所代表的低频成分之后的数据序列为h11,则有
h1m11h11 (2.3)
重复进行上述处理过程k次,直到所得到的平均包络趋于零为止,这样就得到了第一个IMF分量C1
210-1-210.5a b h10-0.5-10100200300t400500600 图2.2 进行一次EMD分解之后得到的结果
h1(k1)m1kh1k C1h1k
(2.4)
第一个IMF分量代表原始数据序列中最高频的组分。将原始数据序列X(t)减去第一个IMF分量C1,可以得到一个去掉高频的新数据序列r1,对r1进行上述平稳化处理过程可以得到第二个IMF分量C2,如此重复下去直到最后一个差值序列rn不再可被分解为止,此时rn代表原始数据序列的均值或趋势
r1C2r2,,rn1Cnrn (2.5)
Huang将这样的处理过程形象地比喻为“筛”过程。最后,原始的数据序列即可由这些IMF分量以及一个均值或趋势项表示。
~ 10 ~
海洋要素计算与预报
nX(t)Cjrnj1 (2.6)
由于每一个IMF分量代表一个特征尺度的数据序列,因此“筛”过程实际上将原始数据序列分解为各种不同特征波动序列的叠加。
图2.3给出了对图2.1中信号进行EMD分解得到的最终结果,其中只有两个IMF分量C1和C2,分别对应于公式(2.1)所给定的两个频率。图中虚线代表由公式(2.1)给出的解析信号,因为EMD方法的分解精度很高,因此对于这种简单的信号分解得到的两个模态与原有模态几乎没有差别,而且只需一次分解就可以达到。
10.5C1C20-0.5-110.50-0.5-10100200300t400500600 图2.3 EMD分解得到的最终结果
每一个IMF分量既可以是线性的也可以是非线性的。
以上给出的是EMD分解方法的简单描述,在Huang的文章中同时还给出了该方法的完备性和正交性讨论,在此我们就不作介绍。 三、EMD方法中存在的问题
1、EMD方法在处理间歇信号时的不可分问题和产生的模态混合问题
虽然EMD方法能够将两个不同尺度的波动进行有效的分解,但是它所需要的条件是其中一个尺度的信号不能是强间歇的。
下面的表达式所给出的是一个含有强间歇高频信号的数据序列(sig2)。
2cost,ast1264ort1296640X(t)22cost0.02sint,as1264t129664032
~ 11 ~
海洋要素计算与预报
由该表达式生成的波面如图3.11所示。其中只在1264t1296处存在一个频率为大波20倍,振幅为大波
1的高频波动(在图中圆圈内)。该高频波动的50波面如图3.12所示。这种形式的波面结构在实际信号中经常存在,如骑行在大波之上的毛细重力波等。
1.510.5X(t)0-0.5-1-1.5050010001500t200025003000 图3.11 含有高频强间歇信号 的波面
0.030.020.010-0.01-0.02-0.03120012201240126012801300t13201340136013801400 图3.12 叠加在大波之上的小波波形
叠加在大波之上的高频小振幅波动通常情况下只会使得大波的波面发生微小的形变,而不会改变大波原有的平稳特性。在这种情况下EMD方法将会无视小波的存在,使得尺度分离过程无法进行。只有当小波正好处于大波的波峰或波谷处,才会对大波的平稳特性产生影响,而EMD方法才有可能将其进行分解。
~ 12 ~
海洋要素计算与预报
21C1C2C3r0-1-210.50-0.5-10.20.10-0.1-0.20.10.050-0.05-0.1050010001500t200025003000 图3.13 经EMD分解得到的各IMF分量以及趋势项
图3.13是经EMD分解得到的各IMF分量。可以看出在第一个IMF分量C1中存在尺度相差很大的两种波动,同时,由于在两种不同尺度波动相接的地方是光滑过渡的,因此该范围内的波动属于虚假波动。同时应当注意的是,这种模态混合的后果还引入了虚假的低频组分C3以及虚假的趋势项r。Huang称这种情况为模态混合(Huang 1999)。他给出的解决方案是在每一个IMF分量中选定一种尺度作为标准尺度,大于该尺度的波动被强行置为零。这种做法至少因为以下两个原因使之不可取:在实际信号处理过程中,给出一个合适的尺度标准并不是一件容易并合理的事情;由于在不同尺度波动之间存在有平滑过渡带来的虚假信号,而这种虚假信号一旦产生就是成双的(在这个IMF分量中引入虚假信号,在另一个IMF分量中就会有一个与之相反的虚假信号存在)。
通过在原有信号上叠加一个已知高频强信号(我们称之为中介信号)的方法
~ 13 ~
海洋要素计算与预报
强行将平均包络线靠近大波的波面。
所加的高频强信号为:
2Xtmp(t)atmpcosTtmpt (4.3)
其中atmp的选取原则是,它应足够大,以至于在原有的大波波面上形成一系列极大值和极小值。Ttmp的选取原则是,它应该远离主波的频率。实际应用中,atmp和
Ttmp的选取范围是相当大的,即该方法的可操作性或易用性很好。一般来讲,比中介信号频率高或接近的波动可以利用EMD方法有效地提取出来。值得一提的是,当高于中介信号频率的组分不存在时,低于中介信号频率的最高频率波动也同样会被有效地提取出来。
对于图3.11中的信号,我们选取atmp0.2,Ttmp40作为高频中介信号。该信号于原有信号叠加的结果如图4.9所示。
1.510.5X(t)+Xtmp(t)0-0.5-1-1.5050010001500t200025003000 图4.9 原始信号于中介信号叠加后的波面
应用EMD分解方法将图4.9中的信号进行分解后得到一个高频的IMF分量,该IMF分量与中介信号的差即是在中介信号的作用下所提取出的原信号中的高频组分。所得结果如图4.10所示,其中实线表示从原信号中分离出的高频强间歇信号,虚线表示该信号的真实波面。
~ 14 ~
海洋要素计算与预报
0.030.020.010-0.01-0.02-0.0312001220124012601280130013201340136013801400 图4.10 利用中介信号分离出的高频强间歇信号与真实信号的对比
这种分离的效果可以说是非常令人满意的,而且据我们所知,现在还没有其他方法能够达到如此高精度的分离效果。
2、EMD分解方法的边界问题
作为进行EMD分解的一个重要步骤,准确地得到一个信号的上下包络线是十分关键的。目前,得到包络线的方法是利用三次样条函数对极大值序列和极小值序列进行插值。以信号(sig3)
222X(t)cost0.6cost0.5sint,t[5,95]
5025200为例,利用三次样条函数对已有的极值序列差值得到的上下包络线如图3.16所示,其中短虚线表示信号的上包络线,长虚线表示信号的下包络线。
420X(t)-2-4-6-8SignalUpper Lower 10 20 30 40 50 t60 70 80 90 图3.16 利用三次样条插值函数得到上下包络线
~ 15 ~
海洋要素计算与预报
可以看出,由于该信号只存在3个极大值点和4个极小值点,因此所得到的上下包络线都出现了失真,尤其是上包络线在数据的两端出现了巨大的失真。同时,由于信号很短,由端点处造成的包络误差已经“污染”到整个数据序列,这种偏差如果得不到有效地抑制,所得结果的真实性就无从谈起。另一方面,利用三次样条进行曲线拟合出现这样的问题也是很自然的,由于三次样条插值时需要用到前后各两个邻近点,解决问题的唯一途径是在数据序列的两端各增加两个极大值点和两个极小值点。
Huang [1998]给出一种象征性的解决方案(当然不是他们申请专利时所用的那一种),即在数据序列的端点处利用一个合适的特征波直接向外延拓。但是,利用特征波方法延拓图3.16中的信号存在相当的难度,因为该信号中含有两个特征尺度,不同的起始位置和终止位置所需依据的特征波是不相同的。在实际信号中往往同时存在多个特征尺度,简单地利用特征尺度波延拓就更加不可靠。
图3.17给出的是端点处利用特征波延拓进行EMD分解得到各IMF分量与真实信号的比较。虽然前两个IMF分量C1、C2的数据段中部与真实信号还比较符合,但是在两个端点区域信号的差别是相当大的,其结果是使得第三个IMF分量C3与真实信号形成巨大的偏差。对这种比较规则信号处理的结果尚如此,对于实际过程中复杂信号的处理结果就更加无法预测。如此结果从信号处理的角度来讲是不能接受的。这也就是为什么Huang能够在其发表的文章中给出这种端点延拓方法的同时又将另一种更好的方法申请专利的根本原因。
1C1C2C30-120-210.50515253545t5565758595 ~ 16 ~
海洋要素计算与预报
图3.17 由EMD分解得到的各IMF分量与表达式(2.12)中各模态的
对比结果,其中实线为IMF分量,虚线为真实模态。
解决方法:模糊预测、神经网络预测 四、应用实例 1、SST资料处理
nino 1区SST数据分解。1980年10月29日——2001年8月18日7天平均。
420-2-410-17.247.257.267.277.287.297.37.31x 1057.24420-2-47.2420-27.240.50-0.5-103/30/827.257.267.277.287.297.37.31x 1057.257.267.277.287.297.37.31x 1057.257.267.277.287.297.37.31x 10512/24/8409/20/8706/16/9003/12/9312/07/9509/02/9805/29/01 2、海平面数据处理
图5.7给出的是青岛验潮站1949年至1997年的月平均潮位。
50250-25-50194919541959196419691974Year19791984198919941999 图5.7 青岛验潮站1949年至1997年的月平均潮位
利用EMD改进方法,对其分解得到的3个长周期IMF分量,如图5.8所示。这种分解结果的意义重大,因为它不仅能够准确地给出一个波动的波高与周期,同时还能够给出准确的相位,也就是说这种分解给出的是波动的完整信息。另一
~ 17 ~
海洋要素计算与预报
方面,由于长周期运动往往代表着全球的变化规律,因此利用这些长周期运动,可以进行全球环境的中长期预报。另外利用长江流域月平均流量资料或许可以准确地预测出特大洪水的发生年份;在人工地震方法探测油田的领域利用该方法也会得到比其他方法准确的结果等等。
31.50-1.5-331.50-1.5-331.50-1.5-3194910 YearC116 YearC230 YearC319591969Year197919891999 图5.8 由月平均潮位资料中得到的3个最长周期的波动
本章主要参考文献:
1、邓拥军,王伟,钱成春,王忠,戴德君, 2001,EMD方法及Hilbert变换中边值问题的处理。科学通报,Vol.46,No.3,pp.157-263。
2、Huang, NE. et. al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. A. Vol. 454, pp. 903-995。
3、邓拥军,2000,《EMD方法的改进》,青岛海洋大学硕士学位论文。
第四章 回归分析
回归分析是数理统计中和多元统计分析中的一个重要内容,它是研究变量与变量之间相关关系的数学方法。回归分析在海洋和大气的统计预报中有着广泛的应用,用它来分析某一个变量(称为预报量)与一个或多个自变量(称为预报因子)之间的统计关系,并建立预报量与预报因子制件的回归方程,然后根据这个方程做出对预报量的预测。回归分析有着广泛的用途,如:地震预测、农业产量预测,甚至股票涨落的预测。
回归分析是一种统计模型,它包括线性回归和非线性回归。常用的是线性回归,因为它的模型比较简单,在理论上也比较严谨,在海洋和气象中不少要素之间可以近似地存在这种关系,因而可以利用这种方法作出比较符合实际的分析与预报。非线性回归则考虑预报量与因子之间的非线性关系,在海洋和气象中也有
~ 18 ~
海洋要素计算与预报
不少应用。 一、一元线性回归 1、一元线性回归模型
一元回归处理的是两个变量之间的关系,即—个预报量与一个因子子之间的关系。
一般来说,对抽取容量为N的预报量y与预报因子x的一组样本,如认为y
ˆ与x有如下关系: 与x是一种线性关系,那么预报量y的估计量yˆiabxi,i1,2,3,,N (1) y利用最小二乘法原理,可计算出a0和a。
实测值与预报值之差的平方和(预报残差平方和)表示为:
ˆi)2 (2) Q(a0,a)(yiyi1N预报残差平方和的极小值条件为:
QQ0 0,aa0将式(1)、(2)代入上式并整理可得方程组如下:
nabxiyi
i1i12iNNaxibxxiyi
i1i1i1NNN由方程组可解出:
aybx
b1其中,xN1,xyiNi1NNlxylxxi
(xix),lxy(xix)(yiy)。
2i1i1NNy,li1xx2、一元线性回归的方差分析
从回归方程中可以看到,用预报因子x预报y仅是对预报量y的一个估计,
ˆ值还是有差别的,如何评价这种预报关系的好坏,也即毕竞实际的y与估计的y如何衡量所建立的回归方程的优劣,这就需要进行方差分析。
ˆi,预报量的方差: 预报误差可以表示为:eiyiy ~ 19 ~
海洋要素计算与预报
1ns(yiy)2ni121nˆieiy)2(yni11n1n22ˆiy)ei(yni1ni1上式可写为:
222SySyˆSe
22上式表明预报量的方差Sy可以表示成为回归估计值的方差Syˆ和误差方差
Se2之和,回归估计值的方差又称为回归方差,误差方差又称为残差方差。
ˆi)2称为残差平ˆiy)称为回归平方和;Qe(yiy此外,U(y22ii1i1i1nnn方和。
上面方差分析表明,预报量y的变化可以看成由前期因子x的变化所引起,同时加上随机因素e变化的影响。这种前期因子x的变化影响可以归为一种简单的线性关系,这部分关系的变化可以用回归方差的大小来衡量。显然,如果回归方差大,则表明用这种线性关系来解释y与x的关系比较符合实际情况,这种回归模型就比较好。否则,说明这时的回归模型,或者回归方程比较差。 3、回归方程的显著性检验
从回归中的方差分析可知,当预报量方差不变时,回归方差愈大,残差方差愈小,那么我们就可以用它们之比值来衡量回归方程的效果好坏。
可以证明统计量:
2Syˆ/1FS/(n2)2e
遵从分子自出度为1,分母自由度为(n2)的F分布。在给定显著水平((1)为置信度)下,当FF时,认为回归方程是显著的。 4、预报值的置信区间
假设y服从正态分布N(,),则yi的置信度为(1)置信区间为:
yiK/2
ˆi作为估计值,可用估计量: 其中,K2可查正态分布表得到,yi可用yˆei1n2in2 来估计。因此,预报值的置信区间可近似为:
~ 20 ~
海洋要素计算与预报
ˆ。 yiK/2例如:对于95%的置信度,预报值的置信区间为: ˆ yi1.96即预报量得实测值落在此区间内的概率为95%。
ˆ,其置信度为95.4%;yi3ˆ,置信度为99.7% 对于置信区间yi2二、多元线性回归 1、多元线性回归模型
某一海洋或大气要素的变化可能与多个要素的前期变化有关,此时可建立多元线性回归模型,得到多元线性回归方程,从而进行预报。多元线性回归分析的原理与一元线性回归分析完全相同。
假定预报量y与p个因子的关系是线性的,为研究它们之间的联系作n次抽样,每一次抽样可能发生的预报量之值为:
y1,y2,y3,,yn,
相应的p个因子的观测值为:
。 xi1,xi2,xi3,,xin(i1,2,,p)
多元回归模型表示为:
ˆib0b1xi1b2xi2bpxip (i1,2,,n), y可写成矩阵形式如下:
YXB,
ˆ11y1yˆ2其中,Y,X1ˆynx11x21xn1x1pb0bx2p,B1。 bxnpp与一元回归相同,仍然可采用最小二乘法确定系数矩阵B,从而可以得到多元回归方程:
ˆb0b1x1b2x2bpxp。 y最小二乘法要求,对于n次观测样本,满足下式:
ˆi)20 Q(yiyi1n由极值原理有:
~ 21 ~
海洋要素计算与预报
Qb00Qb0, 1Q0bp由此可得到线性方程组,写成矩阵形式为:
SBSY,
其中,
S11S21SSp1S12S22Sp2S1pS1yb1bSS2p,S2y,B2;
YbSSppppynSi,j(xitxi)(xjtxj)t1(i,j1,2,,n) nS(xx)(yyˆ)iyititt1在一般情况下,系数矩阵S为非奇异矩阵,存在其唯一逆矩阵S1,因此,方程组的解为:
BS1SY,
b0ybixi。
i1p2、回归方程显著性检验
与一元回归相似,仍利用回归平方和U及剩余平方和Q来构成统计量进行显著性检验。统计量:
FU/p
Q/(np1)是服从自由度为p和(np1)的F分布。根据统计检验方法,在给定显著水平下,若FF,则回归方程效果显著。 3、预报值的置信区间
多元回归预报的置信区间估计方法与一元回归完全相同,只是用下式估计:
~ 22 ~
海洋要素计算与预报
ˆ三、非线性回归
ei1n2inp1。
海洋和气象要素之间的关系有一些是非线性的,对于这种情况,就需要进行非线性回归。非线性回归是线性回归方法基础上进行的,主要有曲线函数线性化和多项式展开两种方法。 1、曲线函数线性化
此方法是将曲线函数经过适当的变换化为线性函数,然后进行回归计算。 例如,函数:
1ba yx令,y11、x,则上式即变为:
xyyabx,
从而可利用线性回归方法。 2、多项式回归
当预报量与预报因子之间的关系不是很明确的函数关系时,可利用p阶多项式展开。
yb0b1xb2x2bpxp
px,x2x3,„,xx2,x3令,x1px,则上式改写为:
b2x2bpxyb0b1x1p
此时,即转化为多元线性回归问题。
本章参考书:《气象统计原理与方法》、《气象统计分析与预报方法》
第五章 经验正交函数分解
经验正交函数分解(Empirical Othorgnal Function,EOF)又称为自然正交函数分解、主分量分析、主成分分析(Principal Component Analysis,PCA)或特征向量分析,是海洋、气象上多变量分析中常用方法之一,是一种从多个变量化为少数变量的统计方法。我们知道,多个变量之间是相互影响的,研究它们的关系是非常复杂的。为简化分析而又不损失或少损失信息,并能提取它们之间相互关系的主要特征,我们可以利用多个变量之间的相互关系构造一些新的变量,这些新变
~ 23 ~
海洋要素计算与预报
量不仅能综合反映原来多个变量的信息,而且彼此之间是相互独立的,同时按方差贡献大小排列的。这样,就可以取少数变量来代替原来多变量的研究,从而达到降维分析的目的。
经验正交函数分解也能够把随时间变化的海洋要素场分解为空间函数部分和时间函数(主分量)部分。空间函数部分概括要素场的空间分布特点,这部分是不随时间变化的;而时间函数部分则由空间点(变量)的线性组合所构成,称为主分量,这些主分量的头几个占有原空间点(变量)的总方差的很大部分。研究主分量随时间变化的规律就可以代替对场的随时间变化的研究。 一、主成分的定义 1、两个变量的主成分定义
设x1(青岛港高潮数据)和x2(乳山港高潮数据)两个变量,其n次观测的样本数据排列为:
x11,x12,x1n, x21,x22,x2n。
两变量的点聚图如下:
150100500-50-100-150-100-50050100150 (pca2)
由于x1和x2具有良好的相关关系,使n个数据点近似呈直线分布,显然沿直线方向反映了数据的主要变化趋势,而在直线的垂直方向反映了数据变化的次要趋
~ 24 ~
海洋要素计算与预报
势。因此,如把直线方向作为新变量y1,与直线垂直方向作为新变量y2,则新变量y1和y2具有如下特征:
(1)y1和y2相关关系很小,即相互独立;
(2)n个数据点在y1方向上离散度最大,而在y2方向上离散度最小。 因此,新变量y1和y2综合反映了原变量x1和x2的信息,是相互独立的,并且按方差贡献大小排列的。所以,y1称为x1和x2的第一主分量,y2称为x1和x2的第二主分量。这种变换的结果相当于在原坐标系中旋转一个角度,使新坐标y1在延直线方向(数据离散度最大方向),新坐标y2在垂直直线方向(数据离散度最小方向)。即使原变量x1和x2的离散度在新坐标系下重新分配,y1占绝大部分,而y2占小部分,但是其总和是相等的:
nnnn(xi12 x)(xx)(yy)(yy)1i12i21i12i2222i1i1i1(53.4%) (46%) (99.4%) (0.6%)
由此可知,若只取第一主成分y1来代替变量x1和x2进行分析,即可达到99%的精度,从而达到降维分析的目的,这就是主成分分析的意义。
新变量y1和y2与x1和x2之间的关系如下:
y1x1cos()x2sin()y2x1sin()x2cos()或写成
y1l11x1l21x2L1X y2l12x1l22x2L2X其中:
l11cosl12sinL1,L2 lsinlcos2122X(x1,x2)
L1、L2称为主成分系数向量,并且具有如下性质:
1ij (i,j1,2) LiLj0ij上式说明主成分的系数向量是单位向量,主成分彼此间是无关的。 2、多变量的主成分定义
参照上述两个变量的主成分定义,给出多变量的主成分定义。
设X(x1,x2,,xm)是由m个随机变量组成的随机向量,为讨论方便,设X的数学期望E(X)0,协方差矩阵为SE(XX),则X的第i个的主成分定义为
~ 25 ~
海洋要素计算与预报
, yiLXi (LiLi1,i1,2,,m)
其中,L(L1,L2,,Lm)为主成分系数矩阵;同时满足下列条件:
(1) 第1主成分y1是在yiLiX中,使y的方差达到最大者;
(2) 第2主成分y2是在yiLiX中与y1无关,且使y的方差达到最大者;
(3) 第k主成分yk(km)是在yiLiX中与y1,y2,,yk1无关,且使y的
方差达到最大者。
类似构成二维变量的主成分,m维变量的主成分可写成m个变量的线性函数:
y1l11x1l21x2lm1xmL1XXy2l12x1l22x2lm2xmL2Xyml1mx1l2mx2lmmxmLm写成矩阵形式:
YLX
式中:
y1l11l12yll22212Y,L(L1,L2,,Lm)ylm1lm2ml1ml2m, lmmL称为主成分的系数矩阵。 二、主成分的导出
确定主成分的关键在于求解主成分系数。根据主成分的定义,使原变量转换成主成分,在相互独立的条件下,必须具有最大的方差,并且要求系数向量的长度为1,即LiLj1。因此,问题归结于在满足LiLj1(ij)和LiLj0(ij)的条件下使yiLiX的方差达到最大来确定主成分的系数。
1nS(yiy)2极大,
ni12y其中:yLx,x(x1,x2,xm)为m个变量的平均值向量。新变量
yiLxi
式中xi(x1i,x2i,,xmi)为m个变量第i个样品观测向量。
~ 26 ~
海洋要素计算与预报
1nS(LxiLx)(LxiLx)ni12y1L(xix)(xix)Lni1上式中括号内的部分为m个变量的协方差矩阵,即:
1nS(xix)(xix)
ni1n
于是新变量的方差可表示为
2SyLSL
在条件LL1下的极值问题,利用拉格朗日乘数法转化为 QLSL(LL1) 的求极值问题,即
Q2SL2L0 L整理后得
(SI)L0
要使L有非零解,必须满足
SI0
上式就是矩阵S的特征多项式。因此问题就转化为求矩阵S的特征值及其对应的特征向量。由于S为mm协方差矩阵,设其为非奇异矩阵,则它的m个非零特征值(记为12m0)及其对应的m个特征向量L1,L1,,Lm,均为满
2足Sy最大条件的解。因此,可以得到m个主分量,表示为:
x (k1,2,,m) ykLk(l1k,l2k,,lmk)为第k个特征值k所对应的转置特征向量。 其中,Lk三、主成分的性质
1、主成分Y(y1,y2,,ym)的协方差矩阵为对角阵,即
100002, E(YY)00m说明各主成分之间是相互独立的,对角阵是由m个主成分的方差组成的。
2、m个主成分的方差总和与原m个变量的方差总和相等,即
mmSkk1i1ii
3、各主成分的方差贡献大小是按矩阵S的特征值由大到小顺序排列的,即
~ 27 ~
海洋要素计算与预报
12m0
由此性质可得到两个常用的指标: (1)第k个主成分的方差贡献率为:
Ukk/i (k1,2,,m)
i1m(2)前p个主成分的累积方差贡献率为:
G(p)i/i (pm)
i1i1pm实际应用中,经常用前几个主要成分的累积方差贡献率的标准(如: G(p)80~90)来确定前几个主成分的个数%四、主成分的计算
1、将m个变量的n次观测值排列成mn的资料矩阵; 2、计算资料矩阵的协方差矩阵S;
3、求解协方差矩阵S的特征值1,2,,m和对应的特征向量L1,L2,,Lm; 4、计算每个主成分的累积方差贡献率G(p);
5、根据给定的累积贡献率G0,确定前p个特征向量L1,L2,,Lp (pm),
并根据yiLiX (i1,2,,p)计算所有主成分的值。
五、经验正交函数分解(EOF)
经验正交函数分解是针对分析海洋或气象要素场的应用提出来的。它的基本原理和计算方法与主成分分析完全相同。
设在
m个观测点,进行了n次观测,把m个变量的n次
x11xX21xm1x12x22xm2x1nx2n xmn观测排列成mn的资料矩阵:
我们希望能将xij分解成空间函数与时间函数的线性组合。
xijlikykjli1y1jli2y2jlimymj (i1,2,,m;j1,2,,n)
k1m其中lik称为空间函数,它不随时间j变化,仅依赖于空间点i变化;而ykj称为时间系数,仅依赖于时间j,而与空间点无关。
写成矩阵形式:
XLY (1)
其中:
~ 28 ~
海洋要素计算与预报
l11l12llL2122lm1lm2y11yY21ym1y12y22ym2l1ml2m(l)
ijmmlmmy1ny2n(y) mijnymnL表示空间函数矩阵,Y表示时间函数矩阵。在分解过程中,要求不同的空间函数之间具有正交性,即满足:
1ij LiLj0ij因此空间函数具有性质:LLI(单位矩阵)。
用L左乘(1)式,则有:
YLX
可见上式与主成分的表达式完全相同。因此,EOF的时间系数就是m个空间点的主成分,空间函数就是主成分的系数。因此,EOF的计算也归结为求解m个变量的协方差矩阵的特征值和特征向量,其特征向量就是空间函数。
根据主成分的性质,可以用头几个时间函数(即主成分)与空间函数的线性组合来对原始场做出估计与解释。这就是EOF的主要目的。 六、时空转换
当变量数大于样本数,即mn时,变量的协方差矩阵SXX的阶数较大,对求解S的特征值和特征向量就会增加难度。在这种情况下,可以先求矩阵XX的特征值和特征向量,再求XX的特征向量。因为XX与XX的秩是相同的,即它们的非零特征值相同,因此,利用它们的关系即可求得XX的特征向量。从而完成对原m个空间点的EOF分解。这种方法称为时空转换。
设XX(mm)的特征值所对应的特征向量为LR,XX(nn)的特征值所对应的特征向量为LQ,则:
LRXLQ 其中,为特征值。
本章参考书:《气象统计原理与方法》、《气象统计分析与预报方法》
~ 29 ~
海洋要素计算与预报
第六章 最小二乘法潮汐调和分析与潮汐特征值 一、分潮与潮汐调和常数
引潮力可以展开成为许多余弦振动之和,其中的每一个振动称为一个分潮。显然不能期望实际海洋的潮汐与由平衡潮理论所得出的结果一样。但是可以期望,在某一个频率为f2的周期性变化的引潮力分潮(包括垂直和水平引潮力)的作用下,海洋也要产生这一频率的振荡。在某一定地点,海面高度变化应包含有这个频率的成分,可以写作Hcos,它代表了实际潮汐的一个分潮,其中振幅H对一定地点为常量,位相以速率均匀增加。
既然实际分潮位相的增加速率与引溯力相角的增加速率相同,则将与引潮力相角保持着不变的位相差。当然,可以将与垂直引潮力相应分潮的相角进行比较,也可以与水平引溯力或平衡潮的相应分潮的相角进行比较。在潮汐分析和预报工作中,习惯于与垂直引潮力或平衡潮的相角(这两者位相相同)进行比较,并规定:
(1)
为地方迟角,其中为当地垂直引潮力的位相。之所以被称为迟角是因为如它为正,则当垂直引潮力分潮达最大,即当0时,为,需要再经过这
样一段时间,才能达到0°,亦即实际潮汐分潮才能达到最大。所以反映了实际分期相对于天文分期的位相落后。
它是由把当地实际潮汐分潮位相和当地天文分潮的相角相比较而得出的,的物理意义较明显,早先的潮汐文献中多采用这种迟角。但在实际工作中采用它常常会感到不方便,因为这样做必须对每个不同经度的地方计算该处的天文相角。现在在实际潮汐分析和预报中,更多地采用另一种迟角,其规定是这样的:实际潮汐观测所用的时间是区时,譬如说用东N区时,在这个时间系统的tN时刻,设实际分潮的位相为N,而在格林威治的tGtN时刻,若格林威治天文相角为G,则取迟角:
gGN (2)
时刻的实际分潮位相N,则这样,以后作预报时,如欲知道观测站在N区时tN减去g即得。tN的天文分潮相角G可用格林威治tG这样就不必推算当地的天文
tN相角,格林威治天文相角可以直接用于任何地点。在这里,我们注意,等式tG只是数字上相等,它们在实质上不是同一时刻,在格林威治到达tG时刻要比在东
N时区到达tN时刻迟N小时。
为了建立g和之间的关系,我们要弄清式(1)和(2)之间意义上的差别。在式(1)中,实际际分潮和天文分潮的位相的比较是在同一个时间系统中进行的,也就是说,是天文分潮和实际分潮在同一瞬间的位相之差。所谓同一瞬间意思
~ 30 ~
海洋要素计算与预报
就是在同一时间系统中的同一个时间值。由于位相差不随时间变化,所以时间系统倒是可以随便选用,不一定要用地方时,所要强调的是对天文分期和实际分潮要用同一个时间系统。而式(2)中,天文和实际分期的位相则分别采用不同的时间系统:实际分潮采用区时,天文分潮采用世界时。还不只此,天文分潮用的不是相对于观测地点的相角,而是用相对于格林威治的相角。这样一来,由于格林威治tG时刻要比在数值上与之相等的东N区时的tN时刻落后N小时,故在格林威治tG时刻,实际分潮的位相已不是N,而是NN。而于格林威治tG时刻,位于东经L处的观测地点的天文相角则是G1L。如前所述,迟角与所用的时间系统无关,现在NN和G1L都是在同一时间系统(世界时)中同一时刻(tG)的观测地点的实际分潮和天文分潮的位相,故两者之差为:
G1L(NN)GN1LN
将(2)式代入,则得与g之间的换算关系:
g1LN
注意式中L和N系指东经和东N时区,若是西经或西N时区,则其值为负;1为分潮的第一个Doodson数。
g被称为格林威治迟角,此名称不太确切。
由于迟角与所用的时间系统无关,故在资料中不需要注明采用的时间系统。而对迟角g情况就不同。凡是资料中用了迟角g,就必须注明所用的时区。如果某处的一组调和常数,观测时采用了东付N时区,算得的迟角为g,现在想把它化到东N时区下的迟角g,则可用
gg(NN)
计算。
H和g()叫做实际潮汐分潮的调和常数,它们反映了海洋对这一频率外力的响应。这种响应决定于海洋本身的动力学性质。由于海洋环境的变化十分缓慢,对一般海区,调和常数具有极大的稳定性,在不特别长的时期内,可充分近似地认为是常数。只有在那些短期内地形有重大变化的地点,调和常数才显示出随时间有明显可感觉到的变化。
事实上引潮力是由许多不同周期的振动迭加起来的,因而实际海面升降也是由许多不同周期的振动迭加而成的。故应当用下式表示潮位高度:
hS0Hicos(ii)S0Hicos(it0ii)
ii但更常用的是用
hS0Hicos(igi)S0Hicos(it0igi)
ii来表示潮高。这个式中要理解为格林威治天文相角,为了省略,这里和以后都
~ 31 ~
海洋要素计算与预报
不加下标。0为t0时刻的值。上两式中多S0为长期平均水位高度。 迟角的gi代表了实际分潮相对于引潮力分潮的位相落后,而实际分潮对引潮力分潮的振幅比HiCi(Ci为与分潮相对应的引潮力的相对大小)代表了实际分潮相对于引潮力分潮的放大率,它们一起代表了观测地点对频率fii2的周期性外力的响应特性。把它们看作频率的函数,就叫做响应函数。实际潮汐的响应函数通常是比较平滑的。 二、最小二乘法潮汐调和分析方法
实际水位可以看作是许多调和分潮迭加的结果。不过在实际分析中只能选取其中有限个较主要的分潮。假设我们选取了J个分潮,可将潮位的表达式写作:
JSHhjcos(jtj), 0j1为第j个分潮的初位相,g0j,0j为t0时刻的值。这里,用j其中jj、H、gj等表示该地点平均水位和调和常数等的真值,以区别于由观测资料S0j分析得到的包含误差的分析结果S0、Hj和gj等。
实际观测到的潮位可写作
JrS)r, Hjthhjcos(0jj1其中r为余差或噪声,它包括由气象等因素引起的不规则扰动、观测中存在的误差、数据处理中的差错和截断误差、被忽略的分潮等,有时也简单地称为观测误差。上式也可写作
JJxjsin(jt)r, hSjcos(jt)y0j1j1。 ,ysincosjHjH式中:xjjjj1、任意时间间隔观测序列的方程组导出
如果在N个时刻tt1,t2,,tN,有N个潮高观测值hh1,h2,,hN,那么,就可以建立如下由N个方程构成的方程组:
JJjcos(jt1)yjsin(jt1)h1r1S0xj1j1JJjcos(jt2)yjsin(jt2)h2r2S0x j1j1JJjcos(jtN)yjsin(jtN)hNrNS0xj1j1 ~ 32 ~
海洋要素计算与预报
在实际分析中,不可能事先知道余差值,因而只能对下列方程求解:
JJS0xjcos(jt1)yjsin(jt1)h1j1j1JJS0xjcos(jt2)yjsin(jt2)h2 (1) j1j1JJS0xjcos(jtN)yjsin(jtN)hNj1j1上式中,由于忽略了噪声,能够得到是是具有一定误差的S0、xj和yj值。式中所有观测时间t1,t2,,tN和角速率j(j1,2,,J)都为已知,且hN也是已知的观测值,因此上式为含有2J1个未知量的线性方程组。潮汐分析的任务就在于从上述N个方程中得出这些未知量的尽可能可靠的数值。
因为未知量总共是2J1个,故一般来说如进行2J1次潮位观测便可列出
2J1个方程,从而可解出这2J1个未知量。但是由于观测值h并不正好等于
,而是还包含噪声r,故计算结果就会受到噪声的影响。因此,在实际工作中,h总是希望尽可能多地进行一些观测,以使得计算得到的S0、xj和yj值尽可能接
、xj和yj。这样,在实际潮汐分析中,总有付N2J1。在这种情况近真值S0下,一般不存在一组由2J1个数组成的数组,将它代入方程组之后能够满足所有N个方程。这样的方程组叫矛盾方程组。对于这类方程组可以来用最小二乘法处理。
将方程组(1)写成如下形式:
a00x0a01x1a02x2a0MxMh1axaxaxaxh1001111221MM2 (2) aN0x0aN1x1aN2x2aNMxMhN其中M2J,NM。在NM的情况下,没有一组解(x0,x1,x2,,xM)能满足方程组(2)中的全部方程,但可以找到一组解,使下式
(an0x0an1x1aaMxMhn)2 (3)
n1N达到最小。上式达到最小的条件是:
0 x0x1xM将(3)式代入上式,即可得到具有M1个未知数的M1个方程组成的方程组:
~ 33 ~
海洋要素计算与预报
c00x0c01x1c02x2c0MxMcxcxcxcx1001111221MMcM0x0cM1x1cM2x2cMMxMf1f2fM (4)
其中:
cijanianj(i,j0,1,2,,M) (5)
n1Nfmanmhn(m0,1,2,,M)
n1N(5)式表明方程组的系数矩阵是对称的。
方程组(4)叫做矛盾方程组(2)的法方程或正规方程。法方程中方程的个数已 和未知数的个数相等,可按普通的线性方程组解出,从而计算出S0、xj和yj的值,进一步计算得到Hj和j。 2、等时间间隔观测序列的方程组系数
如果观测时间没有一定的规律,则只能按上一小节的方法计算法方程的各元素。但是经常遇到的是观测时间间隔为等间距的,即tntn1t与n无关。这时计算可大为简化。选取观测中间时刻(它等于第一个记录的观测时刻加上
1(N1)t为时间原点,则可以使计算更简单一些。观测记录个数N可以是奇数2也可以是偶数;不过既然选取中间时刻为t0,则令N为奇数,在下面书写中要更加方便一些。这时,令N1(N1),则观测时刻为Nt,(N1)t,„,2t,0,t,„,(N1)t,Nt。相应的实测水位值亦记作hN,hN1,„,
h1,h0,h1,„,hN1,hN。这时法方程化为:
A00S0A01x1A02x2A0JxJF0ASAxAxAxF1001111221JJ1 (6) AJ0S0AJ1x1AJ2x2AJJxJFJ和
B11y1B12y2B1JyJF1ByByByF2112222JJ2 (7) BJ1y1BJ2y2BJJyJFJ这两个方程组的系数矩阵仍然是对称的。其中:
~ 34 ~
海洋要素计算与预报
A00NNsinjt2A0jAj0(j1,2,,J)1sinjt2 A1NsinNjt(j1,2,,J)jj2sintjNNsin()tsin()tijij122AA(i,j1,2,,J,ij)ijji112sin()tsin()tijij22sinNjt1BjjN(j1,2,,J)2sinjt NNsin()tsin()tijij122BB(i,j1,2,,J,ij)jiij2sin1()tsin1()tijij22NF0hnnN NFhncosnit(i1,2,,J)inNFinNhNnsinnit(i1,2,,J)
由方程组(6)、(7)可解出S0、xj和yj,利用以下两式:
xjHjcosj,yjHjsinj
计算出Hj和j。 3、Fourier系数的计算
Fi和Fi是Fourier系数。如直接用这两个公式来计算,必须求出NJ个余弦和正弦值。而在电子计算机中,每个余弦或正弦值要由一个子程序计算,这就耗费相当长的机器时间。为了节省计算时间,戈尔策(Goertzel,1958)提出了一个迭代算法,这个方法把绝大部分三角函数的计算都取消了,从而大大提高了计算速度。这里给出的戈尔策算法,其证明过程参见《潮汐和潮流的分析与预报》一书。
Kfkcos(k)U2K1cos(K)U2Kcos[(K1)]kK Kfksin(k)U2K1sin(K)U2Ksin[(K1)]kK
~ 35 ~
海洋要素计算与预报
U00UfK1U2(2cos)U1U0fK1 U3(2cos)U2U1fK2U2K1(2cos)U2KU2K1fK4、天文变量与调和常数计算
实际潮汐的分潮从其来源看可分为以下四种:天文分潮、气象分潮、天文-气象分潮和浅水分潮。
从分潮的频率分布来看,分潮在频率上的分布是极不均匀的,而是分成族、群和亚群。在Doodson展开中,按Doodson数1区分潮族,按2区分群,按3区分亚群。在潮族中一般分为长周期分潮族(10)、全日分潮族(11)、半日分潮族(12)、三分日分潮族(13)直到十二分日分潮族(112),共13个潮族。在每一个潮族中,具有不同数量的群和亚群。
在亚群中的各个分潮的角速度是非常接近的,彼此之间只有微小的差异。因此,在资料长度有限的情况下(如一年的逐时观测数据),亚群中的各个分潮是无法区分的。因此,在实际的潮汐分析中,往往将一个亚群合成为一个分潮,此时这一分潮的振幅和迟角不再是常数,而是随着升交点的黄经十分缓慢地变化,一般在相当长(例如一年)的时间内可近似看作不变。这样的分潮实质上是准调和的,但习惯上仍叫做调和分潮。
由于以上原因,潮位表达式实际写成:
hS0fjhjcos(jujgj)S0fjhjcos(tt0jujgj),
jj其中,fj称作交点因子,uj称作交点订正角。此式可用于代替表达式:
hS0Hicos(igi)S0Hicos(it0igi)
ii由以上两式可知,在从方程组解出Hj和j后,应按以下两式计算分潮的调和常数:
hjHjfj
gj0jujj
在上式中,有0j、fj和uj三个量需要计算,这三个量都与天文变量有关。 分潮角速度:
12s3h4p5N6p
其中:
~ 36 ~
海洋要素计算与预报
14.49205211s0.549016530.04106864h(单位:度/平太阳时) p0.004641830.00220641N0.00000196p分潮初位相:
Y年M月D日t时刻(实际计算中是观测数据的起始时间)的天文初相角:
012s3h4p5N6p090
ts277.02129.3848(Y1900)13.1764(ni)24h280.190.2387(Y1900)0.9857(nit)24tp334.3940.6625(Y1900)0.1114(ni)24 tN100.8419.3282(Y1900)0.0530(ni)24tp281.220.0172(Y1900)0.00005(ni)2415tshY1901);n为从Y年1月1日开始式中i为1900年至Y年的闰年数,iINT(4计算的累积日期序数,注意1月1日的日期序数为0;t为时间。以上各式中的单位是度。
fj和uj的计算:
由于fj和uj随时间变化非常缓慢,一般情况下取资料序列的中间时刻计算。各分潮的fj、uj的具体计算公式如下:
Mmmfcos(u)cos(pm45N)m1 Mfsin(u)sin(mpmN)m45m1mm、5和Doodson数见以下表格。 m、4对于所有分潮,其中Mm、Mf、O1、P1、K1、J1、OO1、M2、L2、k2分潮的f和u依照上式计算,其他分潮由这些分潮组合计算。
例如P1分潮:
fcos(u)0.0008cos(2N)0.0112cos(N)10.0015cos(2p)0.0003cos(2pN) fsin(u)0.0008sin(2N)0.0112cos(N)00.0015sin(2p)0.0003sin(2pN)
~ 37 ~
海洋要素计算与预报
~ 38 ~
海洋要素计算与预报
~ 39 ~
海洋要素计算与预报
~ 40 ~
海洋要素计算与预报
~ 41 ~
海洋要素计算与预报
三、潮流调和常数与潮流椭圆要素
由潮汐运动的所引起的海水水平运动称为潮流。一般情况下,潮流观测数据以流速(w)、流向()记录;流向向北为0°,向东为90°,向南为180°,向西为270°。
为进行潮流的调和分析,将潮流矢量分解成东分量和北分量:
uwsin vwcos对于分解以后的潮流东分量和北分量,可以使用与潮位调和分析相同的方法进行调和分析,获得潮流的调和常数。
对于每一个分潮流,其矢量端点随时间变化,其轨迹是一个椭圆。不过这个椭圆的轴不正好在u、v轴上,而是偏转了一个角度。分潮流椭圆长半轴和短半轴即是这个分潮流速可能达到的最大值和最小值,分别叫做该分潮的最大和最小潮流,记作W和w。最小潮流与最大潮流之比叫旋转率,记为k,并且规定,若这个分潮流矢随着时间的增加按逆时针方向旋转,k算作正,否则,算作负。最大分潮流流速W、方向、发生的时间以及旋转率k决定了分潮流椭圆的基本特征,叫做分潮流的椭圆要素。注意,在一个分潮周期内,分潮流两次达到其最大值和最小值。两个最大潮流大小相同,方向相反,且发生时间相差半个分潮周期。最小潮流方向与最大潮流垂直,发生时间相差1/4分潮周期。因此一个潮流椭圆有两个和两个,知道了其中一对和也就知道另一对,但是一个和一个相对应,另一个和另一个相对应,不可混淆。
分潮流的椭圆要素,可用潮流调和常数计算。具体计算方法参见《潮汐和潮流的分析和预报》
v u
~ 42 ~
海洋要素计算与预报
四、潮汐性质与潮汐特征值 1、潮汐性质
潮汐性质的数值指标由下式计算和区分:
HK1HO10.5半日潮HM2HHO10.5K12.0不规则半日潮混合潮HM2 HHO12.0K14.0不规则日潮混合潮HM2HK1HO14.0全日潮HM22、潮汐特征值
平均半潮面:所有高潮位和低潮位的平均值。由于潮汐是具有不同频率的许多振动的叠加,平均半潮面与平均海平面通常是不重合的
平均高潮位:所有高潮位的平均值。平均高潮位事实上只在半日潮港有实际意义。在日潮较大或日潮为主的港口,回归潮时期和分点潮时期的潮汐特征差别很大,同时由于存在较大的日潮不等,故一天中两个高潮的高度也可以差别很大,或者只有一个高潮。在这种情况下,平均高潮位不能用来描述高潮位的基本状况。与平均高潮位相对应的是平均低潮位。
平均潮差:平均高潮位与平均低潮位之差。
大潮平均高潮位:半日潮大潮期间高潮位的平均值。为了减小偶然误差的影响,通常先在朔望日附近取潮差最大的连续三天(在我国它们大都发生在朔望之后)高潮位计算其平均值,并将其作为一次大潮的高潮位,然后计算所有大潮高潮位的平均值。显然,只在半日潮为主的港口需要计算大潮平均高潮位。与大潮平均高潮位相对应的是大潮平均低潮位。
小潮平均高潮位:半日潮小潮期间高潮位的平均值。类似于大潮,在这里应计算上、下弦日附近潮差最小的连续三天(在我国它们常在上、下弦之后)的高潮位平均值。与小潮平均高潮位相对应的是小潮平均低潮位。
平均大潮差:大潮平均高潮位与大潮平均低潮位之差。
平均小潮差:小潮平均高潮位与小潮平均低潮位之差。平均大潮差与平均小潮差之间的差别反映了潮汐振幅随月相变化的平均幅度。
平均高高潮位:所有高高潮位的平均值。在日潮较大的港口,通常存在日潮不等现象。日潮不等可分为三种类型:高潮不等、低潮不等和高潮、低潮都不等。所谓高潮(或低潮)不等,就是在一个太阴日中两个高潮(或低潮)的高度不同。日潮不等随着日潮振幅的增大而加强,所以在回归潮期间日潮不等达到极大值。习
~ 43 ~
海洋要素计算与预报
惯上称一个太阴日内两个高潮中较高的一个为高高潮,较低的一个为低高潮;称两个低潮中较低的一个为低低潮,较高的一个为高低潮。与平均高高潮位相对应的是平均低低潮位。对低高潮和高低潮的类似统计可得平均低高潮位和平均高低潮位。在不规则日潮港和规则日潮港,不统计平均低高潮位和平均高低潮位。平均高高潮位与平均低高潮位之差反映了高潮不等的量值,称作平均高潮不等;平均高低潮位与平均低低潮位之差称平均低潮不等。
平均大的潮差:平均高高潮位与平均低低潮位之差。、 平均小的潮差:平均低高潮位与平均高低潮位之差。
回归潮平均高高潮位:回归潮期间高高潮位的平均值。类似于大潮平均高潮位,通常是先在月赤纬最大日期附近取潮差最大的连续三天高高潮位计算其平均值,然后计算所有这些平均高高潮位的平均值。与回归潮平均高高潮位相对应的是回归潮平均低低潮位。对低高潮和高低潮进行类似统计,可得回归潮平均低高潮位和回归潮平均高低潮位。在日潮为主的港口不计算后两种潮位。回归潮的平均高高潮位与平均低高潮位之差称回归潮平均高潮不等;回归潮的平均高低潮位与平均低低潮位之差称回归潮哥平均低潮不等。
回归潮平均大的潮差:回归潮的平均高高潮位与平均低低潮位之差。 回归潮平均小的潮差:回归潮的平均低高潮位与平均高低潮位之差。 分点潮平均高潮位:分点潮期间所有高潮位的平均值。与回归潮不同,在分点潮期间日潮通常很小,一般可不考虑日潮不等,即在统计时不再区分高高潮和低高潮。与平均高潮位相对应的是分点潮平均低潮位。
平均高潮间隙:高潮间隙的平均值。为了统一,我们约定月球中天时刻是月球在格林威治上、下中天的时刻(格林威治时间),而高潮时刻则是地方用的标准时间(例如我国为东八时)。以半日潮为主的港口,相邻两个高潮的时间差约为半个太阳日,故统计时可以不区分月球的上、下中天。即如果第一个高潮的间隙是从月球上中天时刻算起,相继的第二个高潮的间隙就可以从月球下中天时刻算起。与平均高潮间隙相对应的是平均低潮间隙。低潮间隙是从月球中天时刻到发生第一个低潮的时间间隔。
平均高(低)期间隙除适用于半日潮为主的港口外,还适用于日潮较大或日潮为主港口的分点潮时期。在这些港口,分点潮时期潮汐变化多呈半日潮性质,统计时可在每个分点潮时期取连续三天的观测值计算高(低)潮间隙的平均值。
回归潮平均高高潮、低低潮、低高潮和和高低潮间隙:对以半日潮为主但有显著日潮不等的港口,除了平均高(低)潮间隙外,还应分别统计以上四种不同的平均间隙。统计这些间隙时,除了月中天时刻外,还必须考虑月球赤纬是北还是南。为了统一,我们约定这些间隙都是月球北赤纬时,高、低潮发生时刻与月球上中天时刻的时间差。当月球为南赤纬时,所有这些间隙都应加或减12.42小时。
在日潮为主的港口,只统计回归潮的高高潮和低低潮间隙。
~ 44 ~
海洋要素计算与预报
上述潮汐特征值原则上都应使用实测潮汐资料统计得到,而且所使用的资料长度应不少于一年。在缺乏足够实测资料的地方,可以使用调和常数来计算其中的一部分特征值,具体计算方法见参考书。 3、平均海面、平均海平面与陆地高程
(1)平均海面
由验潮仪观测记录的水位的平均值即是平均海面。一般情况下,是对观测记录每一整小时的潮高进行统计平均求得。可以有日平均海面、月平均海面、年平均海面和多年平均海面。
影响平均海面的主要因素如下:
长周期分潮(周期超过一个太阴日的分潮,主要有:三分之一月、二分之一月、月、半年、年、18.61年)。长周期分潮对全球海平面总的影响,使得赤道高,越向两极越低,最大差值估计约为90厘米
天气变化。气压、风。 海水温、盐变化和海流变化 (2)平均海平面
水位高度等于平均海面高度的平静的理想海面称为平均海平面。一般情况下,采用多于19年资料计算出的平均海面高度作为平均海平面。平均海海平面可以认为是消除了各种随机振动和短周期波动、长周期波动之后的理想海面;它的高度对固定地点是确定的。
(3)陆地高程
陆地高程(高度)是从平均海平面算起的陆地高度(山高)。通常情况下,对于某一个国家或地区,采用一个处于地壳条件比较稳定地区的验潮站的平均海平面,作为高程的基准;不同的国家和地区往往采用不同的基准,目前在全世界还没有统一。
我国以黄海(青岛验潮站)平均海平面作为计算我国陆地山高的标准,青岛验潮站的平均海平面在验潮站零点之上243cm(称作1985国家高程);令外在1985年之前我国曾经采用1956年之前的青岛验潮站的平均海平面(在验潮站零点之上239cm)作为国家高程基准,称作1956国家高程。
美国以波特兰验潮站,日本以东京灵岸岛验潮站,欧洲地区以阿姆斯特丹验潮站的多年平均海平面,分别为国家或该地区的高程基准面。
(4)海图深度基准面与海图水深
海图深度基准面就是算水深时的起算水面,海图上的水深就是该水面至海底的深度。海图深度基准面一般也是潮汐表的潮高起算面,因此,也叫潮高基准面。
在水深测量和编制海图时,通常采用低于平均海平面的一个面作为海图深度
~ 45 ~
海洋要素计算与预报
基准面,它与平均海平面的距离视当地的潮差大小而定。这个面在绝大部分时间内都应在水内,但它又不是表示最小深度的面,在某些很低的低潮时还能露出,忽视这一点,船只就可能有搁浅的危险。如果深度基准面定得高,船只则较易出事故,如基准面定得低,一般不易干出的在海图上干出了,就会降低航道的使用率。同时确定方法应力求统一,以保持连续性。
世界各国所采用的海图深度基准面很不一致,就是同一个国家出版的海图,深度基准面的算法也不一致。如英国在过去是以测量员选的面为深度基准面。海图深度基准面的确定通常无一定准则,故偏高偏低是可能的。近年来,英国为了统一全国的基准面,采用了所谓“最低天文潮面”,即取潮汐预报中出现的最低水位为基准面。
我国过去采用的海图深度基准面很不一致,较常用的有“略最低低潮面”和“理论深度基准面”,现在统一采用“理论深度基准面”。“理论深度基准面”就是理论上可能的最低潮面,它是按照苏联弗拉基米尔斯基方法计算的,具体计算方法见国家标准《海道测量规范》(GB1237-1998)。
~ 46 ~
海洋要素计算与预报
千年一遇百年一遇五十年一遇最高天文潮位设计高水位477cm352cm318cm285cm226cm187cm139cm129cm27cm4.0cm厂址年平均海面黄海平均海水面-1.4cm-12cm-116cm厂址多年平均海面年最低月平均海面年平均高潮位多年平均高潮位年最高月平均海面年最高潮位多年平均低潮位年平均低潮位-165cm设计低水位最低天文潮位-216cm-118cm-219cm-229cm海图深度基准面年最低潮位五十年一遇百年一遇千年一遇-295cm-307cm-346cm
本章参考书目:陈宗镛《潮汐学》,方国洪、郑文振等《潮汐和潮流的分析和预报》
~ 47 ~
海洋要素计算与预报
第七章 海浪数据分析
现代海浪观测一般采用数字化电子仪器,其主要的观测记录数据是波面高度;对于长时间连续观测,也可以记录海浪谱,以缩减数据量。一般采样频率在10Hz以上,记录时间长度一般为20~30分钟,每小时或每3小时观测一次。 一、去倾向和去均值处理
在潮差较大的近岸海域,当在涨急和落急时刻(一般在半潮面附近,即低(高)潮与高(低)潮之间的中间时刻)进行海浪观测,在观测时段内,海面会发生10~20厘米以上的升高或降低,由此导致波面观测数据带有倾向性。观测数据的倾向性将导致,波浪要素的计算产生误差,因此在进行波浪要素分析之前,应将数据的倾向性去掉。可以采用潮汐分析预报方法、EMD分解、高通滤波等方法,进行去倾向性处理。
去均值处理,即实测波面高度数据减去数据的平均值,以便于后续的数据处理。
二、从波面高度序列中读取海浪的波高和周期 1、跨零点波高、周期定义
跨零点:波面自上而下跨过横轴的交点称为下跨零点(图中的红色圆圈),自下而上跨过横轴的交点称为上跨零点(图中的黑色圆圈)。
400300200100h 0-100-200-300-400050100150200250300350400450500 波高:相邻的两个上跨零点(或下跨零点)间的最大(最小)值点,分别定义
~ 48 ~
海洋要素计算与预报
为波峰和波谷,二者之高度差定义为波高。
周期:相邻的两个上跨零点(或下跨零点)间的时间间隔定义为周期T,周期T的平均值于称为平均周期。也可以取两个相邻波峰(或波谷)间的时间间隔为周期,这种定义的周期与前面定义的周期是有差别的,但如果取足够长的时间(例如10—30分钟)的纪录,上述两种方法得出的平均周期则近似地相等。 2、极值点波高、周期定义
波高:相邻的极大值与极小值分别定义为波峰与波谷,二者之高度差定义为波高。
周期:两个相邻波峰(或波谷)间的时间间隔为周期。
此种波高、周期定义方法,常要求先对波面数据进行一定的平滑,以去掉可能由观测误差等因素引起的极小振幅的高频波动。
使用以上两种波高、周期定义方法,将得到不同的平均波高与周期。跨零点定义方法更加常用。
三、波面高度分布、波高和周期的分布,波高和周期的联合分布 1、波面高度分布
把波面高度观测值的变化范围划分为许多连续但不重叠的区间,统计波面高度落在各区间内的频率,即可得到波面高度的经验分布。
21.81.61.41.210.80.60.40.20-0.8-0.6-0.4-0.200.20.40.60.81
~ 49 ~
海洋要素计算与预报
当海浪被视为正态过程,并且认为是不同频率线性波动的线性叠加时,波面高度符合正态分布,对于任意时刻t的波面高度(t),对应的概率为:
2f()exp(2)
(22)21其中,为波面高度,2为波面高度的方差。
大量的观测数据证明波面高度近似地遵从正态分布,在工程和预报等应用问题中,这种假定造成的误差一般是可以忽略的。但对于风与浪之间的能量传递等的研究,以及其他要求高精度的问题中,波面高度的准确统计分布形式是更重要的。
由于真正的海浪是非线性的,而且这种非线性程度随着水深的变浅而加剧,使波面高度分布偏离正态分布。Longuet-Higgins于1963年根据非线性作用给出了更准确的波面高度的理论分布。 2、波高和周期的分布
把波高或周期的变化范围划分为许多连续但不重叠的区间,统计波高或周期落在各区间内的频率,即可得到波高或周期的经验分布。
波高分布:
1.81.61.41.210.80.60.40.2000.20.40.60.811.21.4 基于窄谱和波面高度服从正态分布的假定,得到的波高分布为瑞利分布:
~ 50 ~
海洋要素计算与预报
H2Hf(H)22exp
H0H01841(H0)(H0H)(H0Hrms)
其中:为波高序列的方差,H为平均波高,Hrms为均方根波高。
波高的瑞利分布得到广泛的应用。但70年代后大量的观测证明,波高分布符合Weibull分布:
1Hf(H)1/2m0周期分布:
f()1expH1/2m0a 其中:2.126,8.42,m0为谱的零阶距。
12(12)3/2
其中:
TT,为谱宽度参量。 T四、各种波高计算
平均波高H:资料中所有波高值的平均值。 均方根波高Hrms:波高序列的均方根。
有效波高(1/3大波平均波高)H1/3:将波高序列(不少于100个波)按由大到小顺序排列,取前1/3序列长度个波高的平均值,称为有效波高。
1/10大波平均波高H1/10:将波高序列按由大到小顺序排列,取前1/10序列长度个波高的平均值,称为1/10大波波高。
平均周期T:与波高对应的所有波周期的平均值。
1/3大波对应的周期T1/3:与1/3大波对应的波周期的平均值。 1/10大波对应的周期T1/10:与1/10大波对应的波周期的平均值。 各种波高和周期之间具有如下关系:
Hrms2H1.129H
H1/31.598H
~ 51 ~
海洋要素计算与预报
H1/102.032H T1/31.15T T1/101.14T1/3
T1/101.31T
五、海浪谱估计 1、海浪谱估计方法
相关函数法 周期图法
快速Fourier变换(FFT)
目前常用的是FFT方法。注意由FFT算法得到的只是粗谱,需经过平滑处理。
2、谱矩的计算
谱的r(r0,1,,n)阶距的定义为:
mrrS()d
03、谱的零阶矩与各种波高的关系
Hrms8m0 H1/34.005m0 H1/105.091m0 H1/1006.672m0 4、海浪谱的谱宽度计算
谱宽度反应了谱内能量的集中程度,宽谱表明能量分散,窄谱表明能量集中。谱宽定义如下:
2m2, 1m0m42通常,实际观测资料中,谱宽介于0.2~0.9。
另一种表征谱宽度的参量由Longuet-Higgins提出:
其中,
2,
020m0
~ 52 ~
海洋要素计算与预报
m0m2m12。 2m02在窄谱条件下:
m0m21 m1212
5、谱峰频率与周期的关系
m平均周期:T20m2m平均频率:2
m01/21/2
有效周期:T1/30.937Tmax,Tmax为谱极值对应的周期。
本章参考书目:《海浪理论与计算原理》文圣常、余宙文,《海浪及其预报》钱志春,《物理海洋学》叶安乐、李凤岐。
~ 53 ~
因篇幅问题不能全部显示,请点此查看更多更全内容