15

我试图了解如何使用谐波常数创建潮汐预测表。

我使用了 Bridgeport 的谐波常数（http://tidesandcurrents.noaa.gov/data_menu.shtml?stn=8467150%20Bridgeport,%20CT&type=Harmonic%20Constituents）

并使用此 python 脚本对潮汐分量求和 -

import math
import time

tidalepoch = 0
epoch = time.mktime(time.gmtime()) - tidalepoch
f = open('bridgeport.txt', 'r')
M_PI = 3.14159
lines = f.readlines()
t = epoch - 24 * 3600
i = -24
while t < epoch:
  height = 0
  for line in lines:
    x = line.split()
    A = float(x[2]) # amplitude
    B = float(x[3]) # phase
    B *=  M_PI / 180.0
    C = float(x[4]) # speed
    C *= M_PI / 648000

    # h = R cost (wt - phi)
    height += A * math.cos(C * t - B)

  print str(i) + " " + str(height + 3.61999)
  i += 1
  t += 3600

对于“今天”，每小时打印一个高度。得到的高度大约在我预期的范围内，-0.5 到 7.5 英尺，但对于日期不正确。

我在正确的轨道上吗？我如何确定潮汐时代？在维基百科 ( http://en.wikipedia.org/wiki/Arthur_Thomas_Doodson ) 上的 Doodsen 示例中，他们使用 0 来获得 1991 年 9 月 1 日的结果。但是它们的谐波值与当前发布的不同，而且那个日期确实似乎不适合我。

这是我的 bridgeport.txt 文件的内容 -

1 M2                       3.251 109.6  28.9841042 
2 S2                       0.515 135.9  30.0000000 
3 N2                       0.656  87.6  28.4397295 
4 K1                       0.318 191.6  15.0410686 
5 M4                       0.039 127.4  57.9682084 
6 O1                       0.210 219.5  13.9430356 
7 M6                       0.044 353.9  86.9523127 
8 MK3                      0.023 198.8  44.0251729 
9 S4                       0.000   0.0  60.0000000 
10 MN4                      0.024  97.2  57.4238337 
11 NU2                      0.148  89.8  28.5125831 
12 S6                       0.000   0.0  90.0000000 
13 MU2                      0.000   0.0  27.9682084 
14 2N2                      0.077  65.6  27.8953548 
15 OO1                      0.017 228.7  16.1391017 
16 LAM2                     0.068 131.1  29.4556253 
17 S1                       0.031 175.5  15.0000000 
18 M1                       0.024 264.4  14.4966939 
19 J1                       0.021 237.0  15.5854433 
20 MM                       0.000   0.0   0.5443747 
21 SSA                      0.072  61.2   0.0821373 
22 SA                       0.207 132.0   0.0410686 
23 MSF                      0.000   0.0   1.0158958 
24 MF                       0.000   0.0   1.0980331 
25 RHO                      0.015 258.1  13.4715145 
26 Q1                       0.059 205.7  13.3986609 
27 T2                       0.054 106.4  29.9589333 
28 R2                       0.004 136.9  30.0410667 
29 2Q1                      0.014 238.8  12.8542862 
30 P1                       0.098 204.1  14.9589314 
31 2SM2                     0.000   0.0  31.0158958 
32 M3                       0.012 200.1  43.4761563 
33 L2                       0.162 134.1  29.5284789 
34 2MK3                     0.015 203.7  42.9271398 
35 K2                       0.150 134.7  30.0821373 
36 M8                       0.000   0.0 115.9364166 
37 MS4                      0.000   0.0  58.9841042
4

4 回答 4

4

National Tidal Datum Epoch与Unix epoch无关；它是该时期平均潮位的量级参考。Bridgeport 数据包含每个分量的幅度和相位，然后在最后一列中给出分量的频率。相位与 GMT 相关（GMT 为 0 度），因此以 GMT 为单位指定时间，或根据时区进行偏移。最后一列是组件的功能，而不是位置，因此该列对于任何位置都是相同的。所以问题是对不同频率的正弦曲线求和，比如傅里叶分析。生成的函数很可能没有24 小时周期，因为组件都有不同的（非 24 小时）周期。我用这个作为参考。

在下面的代码中，我包含了一些可选参数来添加一些偏移量修饰符（比如你的 3.16999，不确定它来自哪里）。我将读取和解析代码与计算代码分开。从返回的函数包含一个闭包get_tides中的测量数据。您不必这样做，而是可以使用类。

from __future__ import division
import math
import time

def get_tides(fn):
    """Read tide data from file, return a function that
    gives tide level for specified hour."""
    with open(fn,'r') as fh:
        lines = fh.readlines()

    measured = []
    for line in lines:
        index,cons,ampl,phase,speed = line.split()
        measured.append(tuple(float(x) for x in (ampl,phase,speed)))

    def find_level(hour,total_offset=0,time_offset=0):
        total = total_offset
        for ampl,phase,speed in measured:
            # you could break out the constituents here
            # to see the individual contributions
            cosarg = speed * (hour + time_offset) + phase
            total += ampl * math.cos(cosarg * math.pi/180) # do rad conversion here
        return total

    return find_level


bridgeport = get_tides('bridgeport.txt')

for hour in range(-24,25,1):
    print("%3sh %7.4f" % (hour,bridgeport(hour,total_offset=3.61999)))

和输出：

-24h -0.5488
-23h -1.8043
-22h -1.7085
-21h -0.3378
-20h  1.8647
-19h  4.4101
-18h  6.8374
-17h  8.5997
-16h  9.1818
-15h  8.4168
-14h  6.5658
-13h  4.1003
-12h  1.5669
-11h -0.3936
-10h -1.2009
 -9h -0.6705
 -8h  0.9032
 -7h  3.0316
 -6h  5.2485
 -5h  7.0432
 -4h  7.8633
 -3h  7.4000
 -2h  5.8028
 -1h  3.5317
  0h  1.1223
  1h -0.8425
  2h -1.7748
  3h -1.3620
  4h  0.2196
  5h  2.4871
  6h  4.9635
  7h  7.2015
  8h  8.6781
  9h  8.9524
 10h  7.9593
 11h  6.0188
 12h  3.6032
 13h  1.2440
 14h -0.4444
 15h -0.9554
 16h -0.2106
 17h  1.4306
 18h  3.4834
 19h  5.5091
 20h  7.0246
 21h  7.5394
 22h  6.8482
 23h  5.1795
 24h  2.9995

编辑： 对于特定日期或现在的时间：

import datetime

epoch = datetime.datetime(1991,9,1,0,0,0)
now = datetime.datetime.utcnow()
hourdiff = (now-epoch).days*24

for hour in range(0,25,1):
    print("%3sh %7.4f" % (hour,bridgeport(hour+hourdiff)))
于 2013-02-11T19:46:34.777 回答
4

如果您使用 C，一种方法是使用 libTCD http://www.flaterco.com/xtide/libtcd.html来存储组成数据，这为访问数据提供了一个很好的框架。

要进行计算，您需要将纪元调整为当年的开始。此代码使用来自 libTCD 的函数，并且基于 xTide 中使用的算法。

TIDE_RECORD record;
int readStatus = read_tide_record(self.stationID, &record);
int epochStartSeconds = staringSecondForYear(year);

for (unsigned constituentNumber=0; constituentNumber<constituentCount; constituentNumber++) {
    float constituentAmplitude = record.amplitude[constituentNumber];
    float nodeFactor = get_node_factor(constituentNumber, year);
    amplitudes[constituentNumber] =  constituentAmplitude * nodeFactor;

    float constituentEpoch =  deg2rad(-record.epoch[constituentNumber]);
    float equilibrium = deg2rad(get_equilibrium(constituentNumber, year));
    phases[constituentNumber] = constituentEpoch + equilibrium;
}

然后计算自年初以来偏移的潮汐：

- (float)getHeightForTimeSince1970:(long)date {
  //calculate the time to use for this calculation
  int secondsSinceEpoch = date - epochStartSeconds;

  float height = 0;
  for(int i = 0; i < constituentCount; i++) {
    float degreesPerHour = get_speed(i);
    float radiansPerSecond = degreesPerHour * M_PI / 648000.0;
    height += amplitudes[i] * cos(radiansPerSecond * secondsSinceEpoch + phases[i]);
  }

  height += record.datum_offset;
  return height;
}
于 2013-05-22T16:46:39.427 回答
1

我刚刚编写了一个名为Pytides的小型 Python 库，用于潮汐分析和预测。它仍处于早期阶段，但您可以使用它轻松实现您想要的。我已经写了一个如何做到这一点的例子。

正如 Dave 所提到的，您的方法不起作用的原因是 NOAA 阶段是相对于成分的平衡参数发布的，因此您需要一些方法来计算这些平衡参数（我建议让 Pytides 为您做这件事！）。

于 2013-11-15T19:31:08.293 回答
0

您链接上的“相位以度为单位，参考格林威治标准时间”意味着布里奇波特潮汐与格林威治子午线的理论“平衡”潮汐的“相位”度数不同步。

您还需要将平衡潮汐作为时间的函数。一种方法是从http://coastalengineeringmanual.tpub.com/Part-II-Chap5/Part-II-Chap50024.htm的 2010 列中复制常量并使用 2010-01-01T00:00:00 作为您的时间基准/纪元/零。

此外，还有一个节点因素会根据月球的节点进动来调节许多成分。

您可以编辑并运行http://adcirc.org/home/related-software/adcirc-utility-programs/中的 titan_fac.f 程序，而不是从表中复制，以在以下位置生成一组节点因子和平衡潮汐成分您自己选择的时间基准/纪元/零。

以下是修改后的tide_fac.f 程序的输出，用于生成 2013-01-01T00:00:00 的节点和平衡潮汐相位参数：

TIDAL FACTORS STARTING:  HR- 0.00,  DAY-  1,  MONTH-  1  YEAR- 2013

FOR A RUN LASTING   365.00 DAYS


CONST   NODE     EQ ARG (ref GM)
NAME   FACTOR    (DEG) 


M2    1.02718     270.21
S2    1.00000       0.00
N2    1.02718      16.12
K1    0.92336      17.70
M4    1.05510     180.42
O1    0.87509     248.94
M6    1.08377      90.63
MK3   0.94846     287.91
S4    1.00000       0.00
MN4   1.05510     286.33
NU2   1.02718      73.02
S6    1.00000       0.00
MU2   1.02718     178.93
2N2   1.02718     122.03
OO1   0.63334     333.60
LAMB  1.02718     287.40
S1    1.00000     180.00
M1    0.00000     159.37
J1    0.89282     275.36
MM    1.09369     254.09
SSA   1.00000     201.62
SA    1.00000     280.81
MSF   1.02718      91.28   
MF    0.74476     312.33
RHO1  0.87509      51.75   
Q1    0.87509     354.85
T2    1.00000       2.35
R2    1.00000     177.65
2Q1   0.87509     100.76
P1    1.00000     349.19
2SM2  1.02718      89.79
M3    1.04114     225.31
L2    0.00000     348.15
2MK3  0.97424     162.72
K2    0.81662     214.58
M8    1.11323       0.84
MS4   1.02718     270.21
于 2013-05-22T04:55:25.663 回答