从理论上说分潮的个数可以是无限的, 但实际分析计算时只取有限个。潮汐分析的目的
就是计算出根据当地实际情况而精选的各分潮的调和常数。
在实际计算时, 式( 1) 改写为
M
∑ h( t ) = A 0 + ( A kcosRk t + BksinRkt )
( 2)
k= 1
式 中 A 0= H 0, A k = f kH kcos[ gk - ( v0 + u) k ] , Bk = f kH ksin [ gk- ( v 0 + u) k] , M 为分潮个
10 1. 927
本文算法 2～4 月
1. 909 0. 771 284 0. 216 346 0. 151 295 0. 085 235 1. 334 0. 0459
0. 164 0. 170
常规算法 1980 年 本文算法 2 月
1. 824 0. 734 285 0. 255 352 0. 143 311 0. 093 240 1. 191 0. 0355 7 1. 830
2. 044 0. 778 286 0. 205 344 0. 162 291 0. 089 232
1. 313 1. 358 1. 311 1. 321
0. 0329 0. 0345 0. 0526 0. 0555
0. 142 0. 147 0. 185 0. 193
常规算法 1980 年
1. 926 0. 732 286 0. 230 345 0. 147 297 0. 085 238 1. 335 0. 0436
的样本值, 拟合曲线除满足( 3) 式的条件外, 还必须使
∑ $′=
1 N
N
[ h′( ti )
i= 1
-
0] 2
( 4)
最小, 其中 h′( ti ) 是 h( t ) 在 ti 时刻的导数值。而式( 3) 、( 4) 同时为最小的充要条件是
N
∑ $ +
$ ′=
1 N
{ [ h( ti )
i= 1
i
然后, 我们分别对大连、基隆、彭湖和那坝四个海洋站 1995 年全年高低潮样本进行调和 分析, 并将各主要分潮的结果与文献[ 4] 所列资料作了比较( 见表 2) , 两者较为一致, 令人满 意。
站 名 大 连 基 隆 澎 湖 那 坝
方 法
文献[ 4] 本文算法 文献[ 4] 本文算法 文献[ 4] 本文算法 文献[ 4] 本文算法
43
数较少( 因为分潮的选取与记录长度及间隔有关) ; ( c) 与分潮相关的天文因子在选取时存在 着不可避免的差异( 我们发现, 文献[ 3] 和[ 4] 中所列的相同潮位站的调和常数亦略有不同) 。 尽管如此, 两者的调和常数却十分一致。图 2 为半年记录 47 个分潮两者拟合曲线与实测值
表 1 调和分析结果与比较
方 法 文献[ 3]
记录 长度
平均海平面
分
( m)
潮实计
M2
数测算H g
主要分潮调和常数
S2 Hg
K1 Hg
O1 Hg
平均误差
理论基面拟合精度
N
( 平海下,
m)
$
∑ 1
N
ûh( ti)
i= 1
- hiû
0. 74 306 0. 22 6 0. 18 325 0. 13 280
常规算法 本文算法 1980 年 常规算法 2～7 月 本文算法
12 小时左右的各分潮之间产生混淆现象, 而该站半日分潮又是绝对占优势的, 因而干扰了
方程组的系数矩阵, 导致调和分析结果的严重失真、拟合曲线失态现象( 图 1 中虚线) 。
图 1 北隍城岛潮位拟合曲线
42
海 洋 工 程 第 15 卷
3 算法的改进
为寻求一个普适的由高低潮样本资料进行调和分析的算法, 我们以为针对高低潮这样
-
hi] 2 +
h′2( ti ) }
( 5)
为最小。
依最小二乘法原理生成的线性方程组为( 设 R0= 0)
N
∑ ai0A 0 +
( aij A j + bij Bj ) = ei
j= 1
M
∑ ck0A 0 +
( ckj A j + dkj Bj ) = f k
( 6)
j= 1
( i = 0, 1, 2, …, M , k = 1, 2, …, M )
ûh( t i) -
i= 1
hi û两者相差均不超过
1cm ; 其它如平均海平
面及理论基面亦都小于 5cm 。但两者的调和常数与文献[ 3] 中所列值有一定误差, 我们认为
可能是下列原因引起的: ( a) 资料较短( 精确调和常数的获得需 18. 6 年的资料) ; ( b) 分潮个
第 3 期 潮汐调和分析的一种算法
表 2 主要分潮调和常数比较
M2
S2
K1
H
g
H
g
H
g
0. 988
300
0. 291
348
0. 270
14
1. 039
288
M+
1,
为此令55A$0 =
5$ 5A k
=
5$ 5Bk
=
0( k =
1, 2, …, M ) , 从而得出由 2M +
1个
未知量构成的线性方程组。由此计算出 A 0、A k 和 Bk( k= 1, 2, …, M ) , 进而可知各分潮的调
和常数。
等样本间隔时由于可使所生成的线性方程组简化为两个互相独立的线性方程组, 使计 算量大为减少而被广泛使用, 但常常需要处理大量的原始样本, 而且对缺损样本的近似处理
法。具体计算表明, 该算法与现今通用的常规算法( 即等间隔最小二乘法) 相比, 既可大大减少所需 原始样本 量, 在 相同记录长 度内, 是常规算法 的三分之 一还少, 同时减少 样本处理的 前期工作 量; 又能保证潮位拟合的精度和预报的可信度, 两者精度相当, 结果一致。该算法原理简单、实用有效, 对于局部样本缺损较易处理, 具有实际应用价值。
c kj = bj k
N
∑ dkj =
1 N
[ sin( Rkt l) sin( Rj t l) +
l= 1
RkRj cos ( Rkt l) cos ( Rj tl) ]
N
∑ e i=
1 N
[ hl cos( Rit l) ]
l= 1
N
∑ f
k=
1 N
[ hl sin( Rktl ) ]
l= 1
显然方程组( 6) 的系数矩阵是对称的, 可用一般的列主元消去法求解。
( 1)
k= 1
X 宋志尧 男 33 岁 工程师 硕士; 严以新 男 49 岁 教授 博士 博导; 茅丽华 女 42 岁 工程师
第 3 期 潮汐调和分析的一种算法
41
其中 H 0 为平均海平面, f k 为节点因子, Rk 为分潮角速度, ( v0 + u) k 为分潮初相, H k 、gk 即 为分潮调和常数。
既繁又常不尽如人意。数学上对任意样本间隔的调和分析是可行的, 但由于所取样本的任意
性, 增加了相关问题研究的复杂性, 因而国内外学者鲜有涉及。我们曾尝试以高低潮资料取
样, 依上述原理进行计算, 一般都能得出与常规算法相一致的结果。但在个别潮位站出现了 很不一致的情况。如北隍城岛海洋站, 由于高低潮间隔变化甚微( 6 小时左右) , 易使周期为
N
∑ 其中 aij =
1 N
l= 1 [ cos ( Ritl ) co s( Rj tl ) +
Ri Rj sin( Ritl ) sin( Rj tl ) ]
N
∑ bij =
1 N
[ cos( Ritl ) co s( Rj t l) -
l= 1
RiRj sin( Rit l) s in( Rj tl ) ]
关键词 潮汐 调和分析 最小二乘法 高低潮 中图法分类号 P 731. 23
1 前 言
无论是海岸及海洋工程的研究, 还是潮位站编制永久潮汐表等, 都必须对实测潮汐资料 进行调和分析, 得出当地各分潮的调和常数, 以获取相关的潮汐特征值, 籍此进行工程设计 和潮汐预报。目前, 普遍采用的调和分析仍是以一百年前达尔文( Darw in) 和杜德森( Doodson) 等人提出的方法为基础进行, 算法上也基本沿袭最小二乘法原理而展开, 如逐次回归法 等[ 1] 。尽管观测设施和计算手段更为先进, 使得分析结果更趋精细和便捷, 但对原始样本( 即 实测值) 必要的前期工作( 包括平滑、舍弃、缺损处理和录入) 更繁杂, 要求也更高。即使一次 中期( 以 30 天计) 资料, 按间隔 1 小时的常规算法处理, 也有 721 个样本值; 若长期( 1 年以 上) 资料或样本间隔更短, 则工作量更巨大。鉴于这些问题, 本文提出的算法, 可减少样本数 三分之二多, 而且易处理局部样本的缺损, 并与常规算法具有一致的精度, 可广泛应用于相 关工程的潮汐分析和预报计算。
2. 069 0. 739 288 0. 219 345 0. 153 294 0. 085 235 47 2. 067
2. 082 0. 761 286 0. 227 343 0. 163 293 0. 091 238
2. 069 0. 741 288 0. 217 345 0. 152 294 0. 085 235 10 2. 067
数。这样, 求调和常数 H k、gk 就等价于计算 A k、Bk 的值。
若用( 2) 去逼近 ti 时刻的样本值 hi( i = 1, 2, …, N , N 为样本总数) , 按最小二乘法原理,