潮汐调和分析的一种算法_宋志尧
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1. 856 0. 796 283 0. 237 349 0. 144 307 0. 106 238 1. 236 0. 0411
0. 148 0. 157
常规算法 1980 年
1. 752 0. 706 279 0. 246 7 0. 138 311 0. 080 240 1. 117 0. 0167
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
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) 的系数矩阵是对称的, 可用一般的列主元消去法求解。
第 15 卷第 3 期 1997 年 8 月
海 洋 工 程 T HE OCEA N ENG IN EERIN G
Vo l. 15, N o . 3 A ug ust, 1997
潮汐调和分析的一种算法X
宋志尧 严以新 茅丽华
( 河海大学海工所 江苏省南京市 210098)
摘 要 本文在常 规潮汐调和分析方 法的基础上, 提出 了一种由高 低潮资料 进行潮汐分 析的算
的样本值, 拟合曲线除满足( 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
2 基本原理和方法
由潮汐理论知[ 2] , 海洋中的潮汐主要是由月球和太阳对地球的引潮力所致。所谓调和分
析, 就是把月球等天体引起的复杂的潮汐现象, 分离成由许多假想天体相对于地球作均速圆
周运动而产生的潮汐( 即分潮) 之和。即对某一潮位站, 在任意 t 时刻的潮高 h( t) 可表为
∞
∑ h( t) = H 0 + f kH kcos[ Rk t + ( v0 + u) k - gk ]
方 法 文献[ 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 月 本文算法
44
海 洋 工 程 第 15 卷
算法省时省力又有效。
图 3 烟台站不同情况下的潮位预报曲线
分潮的选取需根据记录长度 L 及样本间隔 $t 决定。一般地, 任意两个分潮角速度之差 应大于 2P/ L , 而且每个分潮的角速度 R≤P/ $ t, 否则就会产生混淆而导致错误的分析结果。 具体计算表明, 对高低潮资料, 由于样本间隔 $t i 不等, 故取 $t = min{ $ ti} 即可。
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 为样本总数) , 按最小二乘法原理,
对总体样本必须使
N
∑ $ =
1 N
[ h( ti )
i= 1
-
hi ] 2
( 3)
最小。
一般地有 N >
4 结果分析与比较
为进一步检验该算法的有效性, 体现其特色, 我们进行了以下工作。
首先, 我们选择烟台站 1980 年 2～7 月( 半年) 实测资料, 分别利用常规算法和本文算法
进行计算, 具体结果见表 1。表中可见, 就拟合精度而言, 长期资料的结果两者十分一致, 短
N
∑ 期资料亦基本一致;
而平均误差 1 N
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 月 1 日 4 1. 752
本文算法 ～15 日
1. 764 0. 750 277 0. 211 357 0. 141 308 0. 083 231 1. 142 0. 0196
0. 104 0. 109
图 2 烟台站不同算法潮位拟合与误差曲线
的比较, 零值线附近各是其误差曲线, 显示了本文算法和常规算法的一致性。从图表中还可 以看出, 当记录越长, 所选分潮个数越多时, 不仅两者的结果更趋一致, 精度更为相当, 而且 潮位拟合和预报精度更高。图 3 为烟台站不同情况下的潮位预报曲线与实测值的比较, 从中 亦同样可知应用本文算法要较准确地预报潮位, 同常规算法一样, 应尽可能地利用长期记 录, 越长越好。但本文算法所处理的样本量却少于常规算法的三分之一, 记录越长越显示本
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 ) ]
43
数较少( 因为分潮的选取与记录长度及间隔有关) ; ( c) 与分潮相关的天文因子在选取时存在 着不可避免的差异( 我们发现, 文献[ 3] 和[ 4] 中所列的相同潮位站的调和常数亦略有不同) 。 尽管如此, 两者的调和常数却十分一致。图 2 为半年记录 47 个分潮两者拟合曲线与实测值
表 1 调和分析结果与比较
既繁又常不尽如人意。数学上对任意样本间隔的调和分析是可行的, 但由于所取样本的任意
性, 增加了相关问题研究的复杂性, 因而国内外学者鲜有涉及。我们曾尝试以高低潮资料取
样, 依上述原理进行计算, 一般都能得出与常规算法相一致的结果。但在个别潮位站出现了 很不一致的情况。如北隍城岛海洋站, 由于高低潮间隔变化甚微( 6 小时左右) , 易使周期为
法。具体计算表明, 该算法与现今通用的常规算法( 即等间隔最小二乘法) 相比, 既可大大减少所需 原始样本 量, 在 相同记录长 度内, 是常规算法 的三分之 一还少, 同时减少 样本处理的 前期工作 量; 又能保证潮位拟合的精度和预报的可信度, 两者精度相当, 结果一致。该算法原理简单、实用有效, 对于局部样本缺损较易处理, 具有实际应用价值。
i
然后, 我们分别对大连、基隆、彭湖和那坝四个海洋站 1995 年全年高低潮样本进行调和 分析, 并将各主要分潮的结果与文献[ 4] 所列资料作了比较( 见表 2) , 两者较为一致, 令人满 意。
站 名 大 连 基 隆 澎 湖 那 坝
方 法
文献[ 4] 本文算法 文献[ 4] 本文算法 文献[ 4] 本文算法 文献[ 4] 本文算法
ûh( t i) -
i= 1
hi û两者相差均不超过
1cm ; 其它如平均海平
面及理论基面亦都小于 5cm 。但两者的调和常数与文献[ 3] 中所列值有一定误差, 我们认为
可能是下列原因引起的: ( a) 资料较短( 精确调和常数的获得需 18. 6 年的资料) ; ( b) 分潮个
第 3 期 潮汐调和分析的一种算法
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 ) , 进而可知各分潮的调
和常数。
等样本间隔时由于可使所生成的线性方程组简化为两个互相独立的线性方程组, 使计 算量大为减少而被广泛使用, 但常常需要处理大量的原始样本, 而且对缺损样本的近似处理
12 小时左右的各分潮之间产生混淆现象, 而该站半日分潮又是绝对占优势的, 因而干扰了
方程组的系数矩阵, 导致调和分析结果的严重失真、拟合曲线失态现象( 图 1 中虚线) 。
图 1 北隍城岛潮位拟合曲线
42
海 洋 工 程 第 15 卷
3 算法的改进
为寻求一个普适的由高低潮样本资料进行调和分析的算法, 我们以为针对高低潮这样
关键词 潮汐 调和分析 最小二乘法 高低潮 中图法分类号 P 731. 23
1 前 言
无论是海岸及海洋工程的研究, 还是潮位站编制永久潮汐表等, 都必须对实测潮汐资料 进行调和分析, 得出当地各分潮的调和常数, 以获取相关的潮汐特征值, 籍此进行工程设计 和潮汐预报。目前, 普遍采用的调和分析仍是以一百年前达尔文( Darw in) 和杜德森( Doodson) 等人提出的方法为基础进行, 算法上也基本沿袭最小二乘法原理而展开, 如逐次回归法 等[ 1] 。尽管观测设施和计算手段更为先进, 使得分析结果更趋精细和便捷, 但对原始样本( 即 实测值) 必要的前期工作( 包括平滑、舍弃、缺损处理和录入) 更繁杂, 要求也更高。即使一次 中期( 以 30 天计) 资料, 按间隔 1 小时的常规算法处理, 也有 721 个样本值; 若长期( 1 年以 上) 资料或样本间隔更短, 则工作量更巨大。鉴于这些问题, 本文提出的算法, 可减少样本数 三分之二多, 而且易处理局部样本的缺损, 并与常规算法具有一致的精度, 可广泛应用于相 关工程的潮汐分析和预报计算。
从理论上说分潮的个数可以是无限的, 但实际分析计算时只取有限个。潮汐分析的目的
就是计算出根据当地实际情况而精选的各分潮的调和常数。
在实际计算时, 式( 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 为分潮个
表 2 主要分潮调和常数比较
M2
S2
K1
H
g
H
g
H
g
0. 988
300
0. 291
348
0. 270
14
1. 039
288
( 1)
k= 1
X 宋志尧 男 33 岁 工程师 硕士; 严以新 男 49 岁 教授 博士 博导; 茅丽华 女 42 岁 工程师
第 3 期 潮汐调和分析的一种算法
41
其中 H 0 为平均海平面, f k 为节点因子, Rk 为分潮角速度, ( v0 + u) k 为分潮初相, H k 、gk 即 为分潮调和常数。
-
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 )
0. 148 0. 157
常规算法 1980 年
1. 752 0. 706 279 0. 246 7 0. 138 311 0. 080 240 1. 117 0. 0167
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
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) 的系数矩阵是对称的, 可用一般的列主元消去法求解。
第 15 卷第 3 期 1997 年 8 月
海 洋 工 程 T HE OCEA N ENG IN EERIN G
Vo l. 15, N o . 3 A ug ust, 1997
潮汐调和分析的一种算法X
宋志尧 严以新 茅丽华
( 河海大学海工所 江苏省南京市 210098)
摘 要 本文在常 规潮汐调和分析方 法的基础上, 提出 了一种由高 低潮资料 进行潮汐分 析的算
的样本值, 拟合曲线除满足( 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
2 基本原理和方法
由潮汐理论知[ 2] , 海洋中的潮汐主要是由月球和太阳对地球的引潮力所致。所谓调和分
析, 就是把月球等天体引起的复杂的潮汐现象, 分离成由许多假想天体相对于地球作均速圆
周运动而产生的潮汐( 即分潮) 之和。即对某一潮位站, 在任意 t 时刻的潮高 h( t) 可表为
∞
∑ h( t) = H 0 + f kH kcos[ Rk t + ( v0 + u) k - gk ]
方 法 文献[ 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 月 本文算法
44
海 洋 工 程 第 15 卷
算法省时省力又有效。
图 3 烟台站不同情况下的潮位预报曲线
分潮的选取需根据记录长度 L 及样本间隔 $t 决定。一般地, 任意两个分潮角速度之差 应大于 2P/ L , 而且每个分潮的角速度 R≤P/ $ t, 否则就会产生混淆而导致错误的分析结果。 具体计算表明, 对高低潮资料, 由于样本间隔 $t i 不等, 故取 $t = min{ $ ti} 即可。
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 为样本总数) , 按最小二乘法原理,
对总体样本必须使
N
∑ $ =
1 N
[ h( ti )
i= 1
-
hi ] 2
( 3)
最小。
一般地有 N >
4 结果分析与比较
为进一步检验该算法的有效性, 体现其特色, 我们进行了以下工作。
首先, 我们选择烟台站 1980 年 2～7 月( 半年) 实测资料, 分别利用常规算法和本文算法
进行计算, 具体结果见表 1。表中可见, 就拟合精度而言, 长期资料的结果两者十分一致, 短
N
∑ 期资料亦基本一致;
而平均误差 1 N
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 月 1 日 4 1. 752
本文算法 ～15 日
1. 764 0. 750 277 0. 211 357 0. 141 308 0. 083 231 1. 142 0. 0196
0. 104 0. 109
图 2 烟台站不同算法潮位拟合与误差曲线
的比较, 零值线附近各是其误差曲线, 显示了本文算法和常规算法的一致性。从图表中还可 以看出, 当记录越长, 所选分潮个数越多时, 不仅两者的结果更趋一致, 精度更为相当, 而且 潮位拟合和预报精度更高。图 3 为烟台站不同情况下的潮位预报曲线与实测值的比较, 从中 亦同样可知应用本文算法要较准确地预报潮位, 同常规算法一样, 应尽可能地利用长期记 录, 越长越好。但本文算法所处理的样本量却少于常规算法的三分之一, 记录越长越显示本
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 ) ]
43
数较少( 因为分潮的选取与记录长度及间隔有关) ; ( c) 与分潮相关的天文因子在选取时存在 着不可避免的差异( 我们发现, 文献[ 3] 和[ 4] 中所列的相同潮位站的调和常数亦略有不同) 。 尽管如此, 两者的调和常数却十分一致。图 2 为半年记录 47 个分潮两者拟合曲线与实测值
表 1 调和分析结果与比较
既繁又常不尽如人意。数学上对任意样本间隔的调和分析是可行的, 但由于所取样本的任意
性, 增加了相关问题研究的复杂性, 因而国内外学者鲜有涉及。我们曾尝试以高低潮资料取
样, 依上述原理进行计算, 一般都能得出与常规算法相一致的结果。但在个别潮位站出现了 很不一致的情况。如北隍城岛海洋站, 由于高低潮间隔变化甚微( 6 小时左右) , 易使周期为
法。具体计算表明, 该算法与现今通用的常规算法( 即等间隔最小二乘法) 相比, 既可大大减少所需 原始样本 量, 在 相同记录长 度内, 是常规算法 的三分之 一还少, 同时减少 样本处理的 前期工作 量; 又能保证潮位拟合的精度和预报的可信度, 两者精度相当, 结果一致。该算法原理简单、实用有效, 对于局部样本缺损较易处理, 具有实际应用价值。
i
然后, 我们分别对大连、基隆、彭湖和那坝四个海洋站 1995 年全年高低潮样本进行调和 分析, 并将各主要分潮的结果与文献[ 4] 所列资料作了比较( 见表 2) , 两者较为一致, 令人满 意。
站 名 大 连 基 隆 澎 湖 那 坝
方 法
文献[ 4] 本文算法 文献[ 4] 本文算法 文献[ 4] 本文算法 文献[ 4] 本文算法
ûh( t i) -
i= 1
hi û两者相差均不超过
1cm ; 其它如平均海平
面及理论基面亦都小于 5cm 。但两者的调和常数与文献[ 3] 中所列值有一定误差, 我们认为
可能是下列原因引起的: ( a) 资料较短( 精确调和常数的获得需 18. 6 年的资料) ; ( b) 分潮个
第 3 期 潮汐调和分析的一种算法
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 ) , 进而可知各分潮的调
和常数。
等样本间隔时由于可使所生成的线性方程组简化为两个互相独立的线性方程组, 使计 算量大为减少而被广泛使用, 但常常需要处理大量的原始样本, 而且对缺损样本的近似处理
12 小时左右的各分潮之间产生混淆现象, 而该站半日分潮又是绝对占优势的, 因而干扰了
方程组的系数矩阵, 导致调和分析结果的严重失真、拟合曲线失态现象( 图 1 中虚线) 。
图 1 北隍城岛潮位拟合曲线
42
海 洋 工 程 第 15 卷
3 算法的改进
为寻求一个普适的由高低潮样本资料进行调和分析的算法, 我们以为针对高低潮这样
关键词 潮汐 调和分析 最小二乘法 高低潮 中图法分类号 P 731. 23
1 前 言
无论是海岸及海洋工程的研究, 还是潮位站编制永久潮汐表等, 都必须对实测潮汐资料 进行调和分析, 得出当地各分潮的调和常数, 以获取相关的潮汐特征值, 籍此进行工程设计 和潮汐预报。目前, 普遍采用的调和分析仍是以一百年前达尔文( Darw in) 和杜德森( Doodson) 等人提出的方法为基础进行, 算法上也基本沿袭最小二乘法原理而展开, 如逐次回归法 等[ 1] 。尽管观测设施和计算手段更为先进, 使得分析结果更趋精细和便捷, 但对原始样本( 即 实测值) 必要的前期工作( 包括平滑、舍弃、缺损处理和录入) 更繁杂, 要求也更高。即使一次 中期( 以 30 天计) 资料, 按间隔 1 小时的常规算法处理, 也有 721 个样本值; 若长期( 1 年以 上) 资料或样本间隔更短, 则工作量更巨大。鉴于这些问题, 本文提出的算法, 可减少样本数 三分之二多, 而且易处理局部样本的缺损, 并与常规算法具有一致的精度, 可广泛应用于相 关工程的潮汐分析和预报计算。
从理论上说分潮的个数可以是无限的, 但实际分析计算时只取有限个。潮汐分析的目的
就是计算出根据当地实际情况而精选的各分潮的调和常数。
在实际计算时, 式( 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 为分潮个
表 2 主要分潮调和常数比较
M2
S2
K1
H
g
H
g
H
g
0. 988
300
0. 291
348
0. 270
14
1. 039
288
( 1)
k= 1
X 宋志尧 男 33 岁 工程师 硕士; 严以新 男 49 岁 教授 博士 博导; 茅丽华 女 42 岁 工程师
第 3 期 潮汐调和分析的一种算法
41
其中 H 0 为平均海平面, f k 为节点因子, Rk 为分潮角速度, ( v0 + u) k 为分潮初相, H k 、gk 即 为分潮调和常数。
-
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 )