一个简化的潮汐预报准调和分析方法
王如云1，2,李慧娟1，2，蒋风芝2
（1 河海大学水文水资源与水利工程科学国家重点实验室，南京210098；
2 河海大学海洋学院，南京210098）
摘要：在用现有的浅水港日潮汐准调和分析预报方法进行潮汐分析预报时，发现最小二乘法的法方
程组的系数矩阵条件数很大，数量级在108，因此矩阵是坏条件的(或为病态的)，算法不稳定。根据
潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建
立了一个简化的浅水准调和分析模型。利用连云港的多年实测数据检验，简化的准调和分析模型相
对于原准调和分析模型来讲，最小二乘法的法方程组系数矩阵条件数小很多，因此简化后的模型计
算更为稳定。在实测数据时间较长的情况下，简化前后的模型预报精度相当。但当实测数据较短时，
简化前的原模型却没有传统的调和分析模型的预报结果精度高，而简化后的模型却能保持比传统的
调和分析模型的预报结果有一定的改善。特别是简化后比简化前的模型计算时间减少了68%。
关键词：浅水潮汐；准调和分析；潮汐预报
1引言
在潮汐预报方面，一般采用调和分析方法，在深水区域此方法可以获得很好的预报效果，但在浅水区域尤其是河口区域，由于浅水潮汐的复杂性，采用此方法往往不能获得满意的效果。例如杜德森提出的60个分潮[1]，其结果不能令人满意。为此，杜德森后来又提出了一个直接对高低潮进行浅水改正的方法[2]，该方法虽然使高低潮的预报精度有了提高，但把它应用到逐时潮位预报上则有许多困难和不便之处。在浅水区域由于非线性效应的加大，潮波往往产生畸变。此时，高频振动的作用必须予以充分考虑。为了提高浅水区域潮汐预报的精度，从调和分析方法来讲就必须增加高频的浅水分潮。在水深不太浅的区域，浅水分潮的振幅会随着阶数的增高而迅速减小，所以在一般港口采用较少数目的主要浅水分潮即可满足潮汐预报的要求。但在浅水区，常常需要考虑到六阶甚至更高阶的相互作用，才能满足潮汐预报的要求。上个世纪六十年代，一些潮汐学者试图通过扩充高频分潮的数目以使预报结果获得改进，如Zelter and Cumimngs[3]以及Rossiter and Lennon[4]曾将分潮的数目扩充到110多个，但效果并不理想。方国洪等人认为，不理想的原因在于随着频率的增加，高频分潮的数目极速的增加，不可能从中挑选出少数分潮近似代替所有分潮，难以用有限数目的浅水分潮来体现总的浅水效应。可以认为通过增加浅水分潮以改进潮汐预报，其效果可能是比较有限的。
基于以上分析，方国洪等提出了一个浅水潮汐预报的准调和分析方法[5]，可以用来推算任意时刻的潮高，也可以用来推算高、低潮，效果比传统的调和分析法有了显著的改进。但我们使用此方法对连云港的多年潮位实测数据进行分析预报时，发现最小二乘法的法方程组系数矩阵条件数很大，算法不稳定。为此，我们对浅水准调和分析模型进行了简化，简化后的模型计算更为稳定，计算时间大为减少。
2 准调和分析方法介绍
方国洪等人[5]提出的浅水预报准调和方法思路是把潮高分做两部分，一部分为低频部分，由
基金项目：水文水资源与水利工程科学国家重点实验室开放研究基金(2005407411)；中国教育部科学技术研究重点项目（104104）；中国江苏省普通高等学校高新技术产业发展项目(JH03-010)
作者简介：王如云(1963-)，教授，男，安徽芜湖人，从事计算物理学研究，E-mail:wangry@https://www.doczj.com/doc/495065127.html,
潮族0，1，2组成，主要是天文源潮波，另一部分为高频部分，属于浅水分潮。分析分两步进行，
首先对于低频部分，即潮族0，1，2，用如下表1列出的分潮计算。表1中包含了Doodson [5]
所用的所有属于潮族0，1，2的分潮。对实测潮汐()t ξ进行调和分析，求出各分潮的调和常数。然后用实测水位减去平均水位和0，1，2族的潮位，依据剩余值再作进一步分析。
表1 调和分潮
假如利用一年潮位资料()t ξ，对其进行调和分析，计算各分潮的潮汐调和常数，再依据下式将长周期、全日、半日潮族分别加以组合，得出三个基本准调和分潮，然后把三分日及以上的高频潮表示为这些基本准调和分潮的函数。随着频率的增加，只增加少量的准调和项，每一个准调和项可以近似看作一群频率相近的分潮之和。计算每小时0，1，2族的振幅()R t 和()r t 。各族的振幅和位相的变化已不再是常数，而具有缓慢的变化。称这些量是准调和分量。
5
001
cos cos[()]i i i i i i R r f H t v u g σ==++-∑
5
001sin sin[()]i i i i i i R r f H t v u g σ==++-∑
32
116cos cos[()]i i i i i i R r f H t v u g σ==++-∑ （1）
32
116
sin cos[()]i i i i i i R r f H t v u g σ==++-∑
58
2233cos cos[()]i i i i i i R r f H t v u g σ==++-∑
58
2233
sin cos[()]i i i i i i R r f H t v u g σ==++-∑
从动力学的原因来看，浅水分潮由两种非线性产生。一种是运动方程中的平流项，如u u
x
??，...，和连续方程中的非线性项，如
()
u x
???，...，能够产生高级摄动项。另一种是由运动方程中的摩擦项所产生，如u u ρ（这里u 是流速，?是水位，ρ是摩擦系数）。将上式定义的准调和项视为单一潮波，并用cos()kA a θ-表示由非线性产生的属于高频部分的各阶摄动项，则潮高可表达为：
cos()j j j j h k A a ξθ=+-∑ （2） 式中
58
01
cos[()]i i i i i i h A f H t v u g σ==+++-∑ （3）
式(2)中cos()KA a θ-为浅水分潮部分，用34个准调和项表示浅水效应，式中,K θ为准调和常数，,A a 为已知参数，它们的表达式见表2。预报时，首先由表1中第1到第58个分潮的调和常数按(3)式计算h 。再由(1)式计算各族的R 和r ，根据R 和r 按表2计算,A a 。
利用文[1]中的潮汐预报准调和分析方法对连云港多年的实测潮位资料进行了分析预报，得出的结论是改善效果明显，但浅水系数矩阵的条件数极大（见表4），这意味着方程组病态很严重，该算法不稳定。在进行浅水准调和分析时，选择的浅水分潮并非越多越好，选择的分潮达到一定个数后，浅水系数矩阵奇异严重，预报精度不仅没有改善，反而会由于矩阵奇异给预报可靠度带来负面作用。实际上在观测记录数据资料的过程中，会有恶劣天气或仪器磨损等意外情况导致潮位数据缺测或者具有重大误差等情况发生及观测数据长度等的不同，算法若不稳定，这些原
始误差就会在预报过程中被扩大，导致预报结果的可靠性下降，这对预报是极为不利的。鉴于这些，我们对原准调和分析方法进行了简化。根据潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建立了一个简化的浅水准调和分析模型。
表2 浅水准调和项中
,A a 表达式
j
j A j a
j
j A
j a 1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
21R 21R R 221R R 21R R 2
21R R 321R R 31R 21R R 221R R 23
1R R 22R 32R 2
20R R 221R R 21R R 221R R 321R R
13
r 21r r + 21r r + 122r r - 122r r -
122r r - 14r 212r r + 212r r + 212r r + 22r 22r 22r 213r r + 212r r +
212r r +
212r r + 18
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
2
21R R 321R R 221R R 22R 32R 42R 221R R 321R R
321R R 421R R 32R 42R 421R R 42R 52R 52R 62R
123r r - 123r r - 2122r r + 23r 23r 23r 213r r + 213r r + 124r r - 124r r - 24r 24r 214r r + 25r 25r 26r 26r
3 简化的准调和分析模型
在忽略了摩檫力作用的前提下，研究潮波在一维等深等宽半无限长沟渠中传播的情况，简化
后的方程组00u u u g t x x
u u h u t
x x x ζ
ζζζ????++=???????????+++=??????中，第一个方程的第二项，第二个方程的第三和第四项
与相应方程的主要项相比，其量级即使在浅水地区也小于1，作为零级近似，把它们略去了。因
此倍潮波的振幅与相应源潮波振幅的整数次方成正比例，且幂次与频率的倍数相等；复合潮波的振幅则正比例于源潮波振幅的乘积，这个结果对由无摩檫引起的线性浅水分潮大致上是成立的，但是受摩檫力非线性作用引起的浅水潮波，它们的振幅与源潮波的振幅则不一定遵从以上关系，
实际的振幅变化幅度要小于理论的变幅[6]
。因为摩檫力比例于潮汐振幅的平方，而不是线性关系
[1]
，不能认为引潮力增大若干倍，潮汐振幅也相应增大相同的倍数。由此我们假设浅水分潮的一
般表达式为cos()K a θ-（表达式中a 的含义同原来的模型），同时由于简化后的计算模型避免了对表2中浅水准调和项表达式A 的计算，可节省计算时间。
简化模型求解浅水调和常数时方法同原来的相似。最后进行预报并与原模型进行分析比较。对于表2中a相同的项，我们只须保留一项即可，否则浅水系数矩阵的行列式由于存在完全相同的两行或多行，其行列式为零，从而系数矩阵奇异，方程组无解。简化后的准调和项为17项，a的表达式见表3。
表4为简化前后模型的预报结果比较，仍以1967年实测潮位作为分析资料，来预报其它年份水位。为书写方便，我们称简化前模型为模型1，简化后模型为模型2。
表3 准调和项中a的表达式
表4 模型1和模型2均方差比较表
年份调和分析均方
差（㎝）
模型1均方差
（㎝）
模型2均方差
（㎝）
1962 30.5 28.6 29.5
1963 24.2 21.3 22.6
1964 24.5 21.7 23
1965 24.8 22.1 23.3
1966 24.6 21.7 23
1967 20.8 17.3 18.9
1968 24.2 21.6 22.9
1969 28.3 26.1 27.1
1970 26 23.5 24.7
1971 24.8 22.1 23.5
1972 27.6 25.3 26.4
1973 27.1 24.8 25.9
1974 28.5 26.3 27.3
1975 26.8 24.6 25.5
1982 24.4 21.2 22.6
浅水系数矩阵条件数模型1：35912244 模型2：14
分析表4，预报精度虽没有所提高，但两者均方差相差仅为0.9cm到1.6cm之间，同时模型2的浅水系数矩阵条件数相比原来的算法显然小很多，说明简化后的算法稳定。算法若不稳定，原始误差就会在预报过程中被扩大，导致数据的可靠性下降，这对预报是极为不利的。如果算法稳定，即使用于分析的资料存在误差，在预报时，对这种误差也不会有很明显的放大，这对预报结果的可靠性是有利的。
为了进一步验证模型2稳定性好于模型1，我们改变用于分析的潮位资料的时间长度，比较两种模型的预报结果。
表5显示，当用于准调和分析的数据时间段长度缩减到204日时，模型2的均方差开始小于
传统的调和分析方法及模型1。继续缩减到168日（约半年）之后，模型1的预报均方差开始远远大于以传统的调和分析方法预报的均方差，以致大到难以接受。而模型2较传统方法的均方差有所改进，显然简化后的模型稳定性要好得多。
对简化前后模型进行浅水准调和分析计算（包括将各潮族的分潮进行迭加到求出浅水准调和常数的全部过程）所花的时间做了统计（表6）：
表5 不同分析时间长度两个模型预报精度比较
表6 简化前后模型计算时间比较
年份原模型计算时间
（秒）
改进后模型计算时间（秒）
1962 2.48 0.78 1963 2.41 0.77 1964 2.47 0.77 1965 2.41 0.80 1966 2.44 0.77 1967 2.41 0.78 1968 2.39 0.77 1969 2.41 0.78 1970 2.41 0.75 1971 2.41 0.75 1972 2.39 0.75 1973 2.39 0.77 1974 2.44 0.77
1975 2.45 0.77
1982 2.41 0.75
平均计算时间 2.42 0.77
相对原来模型节省计算时间：68%
由表6可以看出，简化后的模型对浅水准调和分析部分计算时间仅有原来的68%左右，当分析数据量很大时，模型二的优势体现的就更为明显了。
5 结论
利用连云港潮位资料，对简化后的模型进行了分析检验，首先其法方程系数矩阵条件数为14，而简化前的模型高达108，说明简化后的模型其算法更为稳定，而原模型计算不稳定。其次，当分析时段潮位资料时间长度短于204日时，模型1的预报误差较传统调和方法来讲，不但没有减小，反而更大了，而模型2较传统方法来讲，仍然有所改进。由于简化后的模型避免了对浅水准调和项表达式中A的计算，且经简化后准调和项只剩下17项，从而在进行浅水准调和分析求准调和常数这一步骤节省了68%的时间。
参考文献
[1] Doodson,A.T.,1928.The analysis of tidal obserbations. Phil.Trans. Roy. Soc.,Ser. A,227:223-279.
[2] Doodson, A.T.，1957.The analysis and prediction of tides in shallow water. Int, Hydr. Rev. 33(1):
85-126.
[3] Zetler, B. D. and R. A. Cummings, 1967. A harmonic method for predicting shallow-water tides.
J. Mar. Res. 25(1):103-114.
[4] Rossiter J.R.and G. W. Lennon, 1968. An intensive analysis of shallow water tides. Geophys.J.R
Astr.Soc.16:275-293
[5] 方国洪，于克俊．浅水港口潮汐预报的一个方法[J]. 海洋与湖沼，1981,12(5)
[6] 方国洪.潮汐摩檫的非线性效应[J].海洋与湖沼,1980,11(2)．
[7] 方国洪，郑文振，陈宗镛等. 潮汐和潮流的分析和预报[M]. 北京:海洋出版社,1986.
[8] 袁慰平，张令敏，黄新芹等.数值分析[M].东南大学出版社
[9] 陈宗镛. 潮汐学[M]. 北京：科学出版社，1980．
[10] 戈丁，李身铎译，潮汐分析[M]，北京:海洋出版社，1980．
[11] 王如云,童章龙,陈耀登,王锦,申宏伟. 基于连续函数最小二乘法的潮汐迭代调和分析方法[J].中国水运，
2007，7(3):116-118.
潮汐的变化规律 由于太阳与月亮对地球的引力作用，我国大部分沿海地区均有一昼夜各出现海水涨落两次的潮汐现象。每月的农历初一至初五（或农历十六至二十）为大潮汐（当地人称“大活汛”）；农历初六至十二（或农历二十一至农历二十五）为小潮汐（当地人称“死汛”）；而初九或二十四为最小潮（当地人称“死汛底”）。每天的潮汐时间均后延45分钟左右，如此周而复始 有个计算公式共，仅供大家参考。 满潮时间=（农历日—1或16）乘以0.8+10:32 干潮时间=满潮时间加或减6:12 潮汐表编辑 潮汐预报表的简称。它预报沿海某些地点在未来一定时期的每天 潮汐情况。在航运方面，有些水道和港湾须在高潮前后才能航行和进出港；在军事方面，有时为了选择有利的登陆地点和时间，就必须考虑和掌握潮汐的情况；在生产方面，沿海的渔业、水产养殖业、农业、盐业、资源开发、港口工程建设、测量、环境保护和潮汐发电等，都要掌握潮汐变化的规律。潮汐表就是为这些方面服务的。 中文名 潮汐预报表 外文名
Tidal prediction table 作用 预报沿海某些地点潮汐情况 服务行业 航运，军事，生产... 最早文献 《海涛志》 包括 主港逐日预报表，附港差比数等 目录 1简介 2文献来源 3港差比数 4潮汐信息 5简便算法 6潮汐时间 1简介编辑 cháo xī biǎo 潮汐表 tide tables 潮汐表又称潮汐长期预测表，即在正常天气情况下由天文因素影响所
产生的潮汐。 2文献来源编辑 英国开尔文 中国唐代窦叔蒙在《海涛志》一文中提出了根据月相推算高潮时刻的图表法，这是保存下来的介绍潮汐预报方法的最早的文献，大约比英国的《伦敦桥潮候表》早400年。19世纪60年代末，英国开尔文和G.H.达尔文等人提出了潮汐调和分析方法，后来还设计和制造了机械的潮汐推算机，使潮汐表的编算工作得到迅速发展。自20世纪60年代以来，电子计算机已广泛应用在潮汐推算工作中。 潮汐表一般包括主港逐日预报表（通常有高潮和低潮的时间和潮高，有的港还有每小时的潮高）、附港差比数、潮信和任意时刻的潮高计算等内容。 主港逐日预报表 潮汐现象可视为由许多不同周期的分潮叠加而成，故任意时刻的潮高可表示为 图片中A为平均海平面在潮高基准面上的高度，表示分潮的圆频率，为交点因子，d为格林威治开始时的天文相角，H和为分潮的调和常数──振幅和迟角。这样，应用已求出的该港的潮汐调和常数，就能
潮汐原理复习思考题整理 （第四章~第五章） 第四章 1.什么是中期观测资料分析和短期观测资料分析，以及调和常数求解的实际步骤 中期观测资料分析：属于不同群的分潮的会合周期最长为1个月，因此把长度长于一个月但不足一年的观测记录称为中期观测资料 短期观测资料分析：观测的时间长度只有一天或几天 调和常数求解的实际步骤： ?中期观测资料分析(TB P103-107) 1）区分主分潮和随从分潮2）取L 段观测记录，式（4.4）可以写为（4.6） 3）将式（4.6）的余弦函数展开得到（4.7） 4）式（4.7）是包含2（P+Q）+1个未知数的由 () 1 L l l M = ∑ 个方程组成的矛盾方程组 5）通过最小二乘法得到矛盾方程组的法方程(4.10) 6）当L=1时，法方程(4.10)变为TB P106 7）引入Q个随从分潮与相应的主分潮的差比关系后，将给出另外2Q个方程(4.11) 8）进一步求得(4.12) ?短期观测资料分析(TB P116-119) 1）潮汐调和常数的初算2）潮流调和常数的计算 3）噪声方差的估计4）不合理数据的舍弃 5）调和常数和余流的计算6）潮流椭圆要素的计算 2.短期资料观测引入的参数D 和d 代表什么含义，具有什么作用？ 振幅系数D 和迟角订正d 用准调和分潮表达式比用调和分潮表达式要简单的多，不但可以简化许多分析过程，而 且对分析实际潮汐特征也能使得问题变得更容易。 3.什么是准调和分潮，它和调和分潮有什么区别 ?实际准调和分潮的振幅和相角与A 小时前的引潮力准调和分潮相应量有关，与其余时刻，特别是与当时引潮力则没有关系，故A 叫做准调和分潮的潮龄 ?区别 4.了解潮汐和潮流的自报TB P119 第五章 1.潮汐特征值的含义TB P120-121 2.对于不同潮汐类型港口潮汐特征值的计算方法
海洋要素计算作业之二——潮汐（威海2013年五月份） 一．本次潮汐调和分析共选取了十三个分潮： MSf，Q1，O1，K1，P1，K2，N2，M2，S2，MK3，M4，MS4，M6 为使您查看方便，将本次大作业的放在本文件夹各文件内，具体参考如下： 1.原数据为：qd.dat; 2.Fortran编程见该文件夹内：tide.f90文件； 3.求各分潮调和常数H、g的值及其中间过程得到的各值见：qd_tide.dat文件；二．对比回报值和实测值： 1. 回报1968年一月份的水位值见：huibao.dat； 2. 用matlab绘制的潮汐过程曲线见：潮汐过程曲线.bmp 3. 用给定的六个分潮求得的高潮和低潮发生的时刻及潮位值见—：gaodichao.dat; 运行tide.f90后求得威海地区2013年5月份的平均潮差。 由图可知：由于只计算了一个月的潮汐数据，所以回报值和实测值相符的不是很好，如果计算一年的数据，应该会取得比较良好的结果。
三．程序 %% 潮汐过程曲线图 clear,clc %% huibao=load('G:\chaoxi\huibao.dat'); % huibao=fread(fhuibao); shice=load('G:\chaoxi\qd.dat'); % shice=fread(fshice); %huibao_y=zeros(1,12*62); %shice_y=zeros(1,12*62); huibao=double(huibao'); huibao_y=double(huibao(:)); %shice_y=reshape(shice',1,[]) %for i=1:12; % for j=1:62 % huibao_y(i)=huibao(i,j) % shice_y(i)=shice(i,j) %end %end shice=double(shice'); shice_y=double(shice(:)); x=linspace(1,31,length(huibao_y)); plot(x,huibao_y,'r-') hold on plot(x,shice_y,'b-') title('威海(37°31′N ,122°08′E)2013年五月潮汐调和分析图') legend('回报值','实测值') xlabel('时间（2013年五月份）') ylabel('水位（m）')
中期水位观测资料的最小二乘分析报告 摘要： 本次实验中采用了KM站（28.05N,121.17E）1997年8月的逐时潮位数据，运用中期水位观测资料的最小二乘分析方法，通过奇异值剔除、调和分析、逐时潮位回报、高低潮及余水位计算等工作，对此验潮站的数据进行了初步分析，并了解了中期水位资料分析的基本流程。 报告主要内容： （1）11个主要分潮（MSf Q1 O1 K1 N2 M2 S2 MK3 M4 MS4）及两个随从分潮（P1 K2）的调和常数H和g （2）图像和数据文件的基本信息 （3）平均潮差和潮汐类型 （4）余水位特征分析 （5）误差分析 （6）程序的相关说明 （1）潮汐的调和常数： 利用最小二乘原理，通过引入差比关系的方法，我们可以成功得到11个主要分潮和2个随从分潮的调和常数如下： 分潮名称调和常数H 调和常数g MSf 121.2222 -32.38747 Q1 62.95736 233.5120 O1 225.5294 238.7111 K1 266.1612 113.2537 N2 420.5689 167.2492 M2 1922.772 174.8581 S2 679.3940 197.3759 MK3 33.19594 252.1002 M4 32.43390 121.7806 MS4 33.60584 198.8826 M6 3.762754 94.29744 P1 73.46050 109.5160 K2 192.9479 201.4156 程序运行结果如图：
其中H关系到分潮的振幅，g关系到分潮的相位。从表中可以看出，M2分潮的振幅最强，对当地潮位的贡献最大，这与实际情况相符，但K1分潮的调和常数H仅有266.1412，结果偏小。 （2）图像和数据文件的基本信息： 本次报告中包含以下数据文件： 1. KM9708new.dat 数据原始文件。 2. KM9708new_02.dat 经过奇异值订正的数据文件，为方便画图时读取，没有输入数据质量信息。 3. 调和常数.txt 保存了调和常数的相关数据 4.回报值.txt 保存了利用六个主要分潮进行数据回报得到的结果，同样为方便读取，没有输入质量控制信
潮汐知识点汇总： 1.中版潮汐表一共有六册，其中第1-3册是我国海区的，4-6册是国外海区的。潮汐表中的潮时是当地标准时。 2.中版潮汐表的主要内容：主港潮汐预报表、潮流预报表、差比数和潮信表、格林尼治月中天时刻表、梯形图卡。梯形图卡作用是求取任意潮时的潮高和任意潮高的潮时。我国一些重要港口每小时潮高可以在主港潮汐预报表里面查找。 3.差比数是指主附港之间的潮时差、潮差比和改正值；潮差比是附港平均潮差与主港平均潮差之比（对于日潮港来说，潮差比是附港回归潮大的潮差与主港回归潮大的潮差之比）；潮时差是主港和附港高（低）潮时之差，正号表示附港比主港发生的晚，负号表示早。 4.下面因素的急剧变化会导致潮汐变化的反常：台风（增水）、寒潮（减水）、春季气旋入海、降水、气压、结冰。（注意：海啸不会导致潮汐变化的反常） 5.中版潮汐表预报潮高的误差是20-30CM，预报潮时的误差是20-30分钟。 6.利用中版潮汐表查主港潮汐，可从目录查主港所在页码。 7.附港潮时的计算公式：附港潮时=主港潮时+潮时差 附港潮高计算，对于中国海区的1-3册，可以利用下面公式: 当主附港季节改正数较大时（SC>10CM）： 附港潮高=（主港潮高-（主港平均海面+主港平均海面季节改正））*潮差比+（附港平均海面+附港平均海面季节改正） 当主附港季节改正不大时(SC<10CM)：附港潮高=主港潮高*潮差比+改正值 对于外国海区的4-6册，可以利用下面公式: 附港潮高=主港潮高*潮差比+改正数+潮高季节改正数 8.利用潮信资料算潮汐：潮信资料包括平均大（小）潮升，平均高（低）潮间隙，平均海面。当地高（低）潮潮时=格林威治月上（下）中天时+当地高（低）潮间隙 对于半日潮港，可用下列方法近似求月中天时间： 上半月：月上中天时=（农历日期-1）*0.8+1200月下中天时=月上中天时±1225 下半月：月上中天时=（农历日期-16）*0.8月下中天时=月上中天时±1225 平均大潮高潮高(MHWS)=大潮升(SR) 平均大潮低潮高MLWS=2*MSL-SR 平均小潮高潮高（MHWN）=小潮升(SN) 平均小潮低潮高MLWN=2*MSL-SN 大潮日指的是初三、十八；小潮日是十一、二十五。 其他日潮高估算： 高潮潮高=大潮升-（大潮升-小潮升）/7.5*（所求日与大潮日间隔天数） 低潮潮高=2*MSL-高潮高 9.求任意潮时的潮高和任意潮高的潮时 T=涨落潮时间间隔t是任意时与高低潮的时间间隔 10.实际水深=海图水深+潮高+（CD-TD）=测水仪测的水深+船舶吃水 过浅滩最小潮高：海图水深+潮高+（CD-TD）>=船舶吃水+富余水深 过高架最大潮高：船舶水面是最大高度+潮高+安全余量>=大潮升+净空高度 其中船舶水面最大高度=主大桅高+型深—吃水
一个简化的潮汐预报准调和分析方法 王如云1，2,李慧娟1，2，蒋风芝2 （1 河海大学水文水资源与水利工程科学国家重点实验室，南京210098； 2 河海大学海洋学院，南京210098） 摘要：在用现有的浅水港日潮汐准调和分析预报方法进行潮汐分析预报时，发现最小二乘法的法方 程组的系数矩阵条件数很大，数量级在108，因此矩阵是坏条件的(或为病态的)，算法不稳定。根据 潮汐动力学寻找高频潮族与低频潮族之间可能的相互作用关系，在只考虑相角的变化率情况下，建 立了一个简化的浅水准调和分析模型。利用连云港的多年实测数据检验，简化的准调和分析模型相 对于原准调和分析模型来讲，最小二乘法的法方程组系数矩阵条件数小很多，因此简化后的模型计 算更为稳定。在实测数据时间较长的情况下，简化前后的模型预报精度相当。但当实测数据较短时， 简化前的原模型却没有传统的调和分析模型的预报结果精度高，而简化后的模型却能保持比传统的 调和分析模型的预报结果有一定的改善。特别是简化后比简化前的模型计算时间减少了68%。 关键词：浅水潮汐；准调和分析；潮汐预报 1引言 在潮汐预报方面，一般采用调和分析方法，在深水区域此方法可以获得很好的预报效果，但在浅水区域尤其是河口区域，由于浅水潮汐的复杂性，采用此方法往往不能获得满意的效果。例如杜德森提出的60个分潮[1]，其结果不能令人满意。为此，杜德森后来又提出了一个直接对高低潮进行浅水改正的方法[2]，该方法虽然使高低潮的预报精度有了提高，但把它应用到逐时潮位预报上则有许多困难和不便之处。在浅水区域由于非线性效应的加大，潮波往往产生畸变。此时，高频振动的作用必须予以充分考虑。为了提高浅水区域潮汐预报的精度，从调和分析方法来讲就必须增加高频的浅水分潮。在水深不太浅的区域，浅水分潮的振幅会随着阶数的增高而迅速减小，所以在一般港口采用较少数目的主要浅水分潮即可满足潮汐预报的要求。但在浅水区，常常需要考虑到六阶甚至更高阶的相互作用，才能满足潮汐预报的要求。上个世纪六十年代，一些潮汐学者试图通过扩充高频分潮的数目以使预报结果获得改进，如Zelter and Cumimngs[3]以及Rossiter and Lennon[4]曾将分潮的数目扩充到110多个，但效果并不理想。方国洪等人认为，不理想的原因在于随着频率的增加，高频分潮的数目极速的增加，不可能从中挑选出少数分潮近似代替所有分潮，难以用有限数目的浅水分潮来体现总的浅水效应。可以认为通过增加浅水分潮以改进潮汐预报，其效果可能是比较有限的。 基于以上分析，方国洪等提出了一个浅水潮汐预报的准调和分析方法[5]，可以用来推算任意时刻的潮高，也可以用来推算高、低潮，效果比传统的调和分析法有了显著的改进。但我们使用此方法对连云港的多年潮位实测数据进行分析预报时，发现最小二乘法的法方程组系数矩阵条件数很大，算法不稳定。为此，我们对浅水准调和分析模型进行了简化，简化后的模型计算更为稳定，计算时间大为减少。 2 准调和分析方法介绍 方国洪等人[5]提出的浅水预报准调和方法思路是把潮高分做两部分，一部分为低频部分，由 基金项目：水文水资源与水利工程科学国家重点实验室开放研究基金(2005407411)；中国教育部科学技术研究重点项目（104104）；中国江苏省普通高等学校高新技术产业发展项目(JH03-010) 作者简介：王如云(1963-)，教授，男，安徽芜湖人，从事计算物理学研究，E-mail:wangry@https://www.doczj.com/doc/495065127.html,
program main real*8 h(100000),h0(100000),preh(100000),prehf(100000),ph0(100000),ph(100000) !实测数据 real*8 w(1000) !分潮角速度 real*8 x(1000),y(1000),aa(1000,1000),bb(1000,1000),f1(1000),f2(1000) real*8 q(1000,1000),d(1000),v(1000,1000),aaa(1000,1000),bbb(1000,1000) integer n,ntidal,khl,kh1,n2,n12 ,np,nep,np12,nbp !n为观测数据时间长度!!ntidal为分潮数目 real*8 tdif,dt,err,pai,tbi,max,min character*15 b,char WRITE(*,111) WRITE(*,222) WRITE(*,111) 111 format(20x,'********************') 222 FORMAT(20x,'**<<潮汐调和分析>>**') pai=3.1415926 dt=1.d0 tdif=8.d0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !读取预报时段实测潮位(连云港，1979年7月) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! open(1,file='d:\Lyg1979.dat') read(1,*)char,icck,iyyk,immk,iddk,ihhk,iccj,iyyj,immj,iddj,ihhj call GDAY(iddk,immk,iyyk,icck,kd) !开始时间转换为18521015后的天数 nbp=(kd-1)*24+ihhk-tdif call GDAY(iddj,immj,iyyj,iccj,kd) !开始时间转换为18521015后的天数 nep=(kd-1)*24+ihhj-tdif np=nep-nbp+1 np12=np/12 nep=12*np12+nbp-1 np=nep-nbp+1 read(1,*) do i=1,np12 read(1,*)(ph0(j),j=(i-1)*12+1,12*i) enddo close(1) call GDAY(01,07,79,19,kd) nmlb=(kd-1)*24+00-tdif call GDAY(31,07,79,19,kd) nmle=(kd-1)*24+23-tdif nmn=nmle-nmlb+1 write(*,*)np,nmn,nmle-nmlb-1+i do i=1,nmn
中期水位资料对潮汐进行调和分析 !利用1996年7月厦门站的潮汐观测数据计算调和常数，并利用主要分潮和浅水分潮进行潮汐预报 program work implicit none character*80::a1 character(len=5),dimension(62,16)::aa integer::bb(62,12),c(62,2),caita(- 371:371),i,i1,i2,j,t1 real::N0,n(13,6),a(0:13,0:13),b(1:13,1:13),s,s0,s1,s2 ,s3,sa,hh !n代表Doodson代码；a,b为系数矩阵real::xiaoa(0:13),xiaob(13),gg1,gg2,pjchaocha,t,ma,mi !计算法方程所需的参数 real,dimension(1:13)::w,u,f,V0,f1(0:13),f2 ! f1和f2为法方程右边系数 real,dimension(13)::sita,h,g,r,h0(13),g0(13),h1(13),g 1(13) !调和常数参数 real,dimension(- 371:371)::caita1,caita3,caita4,caita8,caita9,caita5,caita 11 !主分潮、浅水分潮的潮高数值 real,dimension(:),allocatable::hightide,lowtide,chaoc ha !高低潮数值
integer,dimension(:),allocatable::hightrq,lowtrq,high t,hightt,lowt,lowtt !读取数据，把潮位数据赋值给bb，把年月份数据赋值给c open(unit=2,file='XM_July1996.dat') read(2,'(a)')a1 print*,'数据文件的第一行信息:',a1 do i=1,62 read(2,'(16a5)')aa(i,:) end do do i=1,62 read(aa(i,5:16),*)bb(i,:) read(aa(i,3:4),*)c(i,:) end do do i=1,62 c(i,2)=int(real(c(i,2))/10.0) end do close(2) !计算分潮角速率w w=(/0.002822,0.037219,0.038731,0.041781,0.163845,0.24 1534,0.078999,& &0.080511,0.083333,0.122292,0.161023,0.041553,0.08356 1/) w=360*w print* print*,'角速率w:',w !计算N0 (middle time:1996-7-16 ; data sum:744, middle number:372 ) N0=259.157-19.32818*(1996-1900)- 0.05295*(31*3+30*2+29+15+int((95.0)/4.0)) !初始升交点平 均黄经 N0=-(0.00220641*3+N0) print* !转换成格林威治时间 print*,'N0:',N0
中期水位观测资料的最小二乘分析报告摘要： 本次实验中采用了KM站（28.05N,121.17E）1997年8月的逐时潮位数据，运用中期水位观测资料的最小二乘分析方法，通过奇异值剔除、调和分析、逐时潮位回报、高低潮及余水位计算等工作，对此验潮站的数据进行了初步分析，并了解了中期水位资料分析的基本流程。 报告主要内容： （1）11 个主要分潮（MSf Q1 O1 K1 N2 M2 S2 MK3 M4 MS4）及两个随从分潮（P1 K2）的调和常数H和g （2）图像和数据文件的基本信息 （3）平均潮差和潮汐类型 （4）余水位特征分析 （5）误差分析 （6）程序的相关说明 （1）潮汐的调和常数： 利用最小二乘原理，通过引入差比关系的方法，我们可以成功得到11个主要分潮和2个随从分潮的调和常数如下： 分潮名称调和常数H调和常数g MSf 121.2222 -32.38747 Q1 62.95736 233.5120 O1 225.5294 238.7111 K1 266.1612 113.2537 N2 420.5689 167.2492 M2 1922.772 174.8581 S2 679.3940 197.3759 MK3 33.19594 252.1002 M4 32.43390 121.7806 MS4 33.60584 198.8826 M6 3.762754 94.29744 P1 73.46050 109.5160 K2 192.9479 201.4156 程序运行结果如图：
其中H关系到分潮的振幅，g关系到分潮的相位。从表中可以看出，M2 分潮的振幅最强，对当地潮位的贡献最大，这与实际情况相符，但K1分潮的调和常数H 仅有266.1412，结果偏小。 （2）图像和数据文件的基本信息： 本次报告中包含以下数据文件： 1. KM9708new.dat 数据原始文件。 2. KM9708new_02.dat 经过奇异值订正的数据文件，为方便画图时读取，没有输入数据质量信息。 3. 调和常数.txt 保存了调和常数的相关数据 4.回报值.txt 保存了利用六个主要分潮进行数据回报得到的结果，同样为方便读取，没有输入质量控制
!利用1996年7月厦门站的潮汐观测数据计算调和常数，并利用主要分潮和浅水分潮进行潮汐预报 program work implicit none character*80::a1 character(len=5),dimension(62,16)::aa integer::bb(62,12),c(62,2),caita(-371:371),i,i1,i2,j,t1 real::N0,n(13,6),a(0:13,0:13),b(1:13,1:13),s,s0,s1,s2,s3,sa,hh !n代表Doodson代码；a,b为系数矩阵 real::xiaoa(0:13),xiaob(13),gg1,gg2,pjchaocha,t,ma,mi !计算法方程所需的参数 real,dimension(1:13)::w,u,f,V0,f1(0:13),f2 !f1和f2为法方程右边系数 real,dimension(13)::sita,h,g,r,h0(13),g0(13),h1(13),g1(13) !调和常数参数 real,dimension(-371:371)::caita1,caita3,caita4,caita8,caita9,caita5,caita11 !主分潮、浅水分潮的潮高数值 real,dimension(:),allocatable::hightide,lowtide,chaocha !高低潮数值 integer,dimension(:),allocatable::hightrq,lowtrq,hight,hightt,lowt,lowtt !读取数据，把潮位数据赋值给bb，把年月份数据赋值给c