简介:基于Python的pygtide模块源码包,面向地球物理与海洋学的科研人员与高年级学生,依托成熟的ETERNA PREDICT算法完成高精度地球引力潮计算,可支撑地震预测、海平面变化分析、潮汐建模等研究课题。压缩包总计38个文件,约18MB,主要包含6个Python源码、13个dat数据文件、XML配置、pyd预编译模块、PDF使用手册、示例图片和说明文本,层次清晰,便于直接导入工程。目前已有166人浏览学习,说明其在相关领域具备一定参考热度。资源不仅提供完整的核心算法代码与历史潮汐数据,还附带用户指南、README及测试脚本,用户可直接调用或在此基础上调整参数,省去重复编译与数据整理环节,适合作为高校科研项目或课程设计的实现基础。
1. 处理重力固体潮数据时,你首先扣掉的这段信号就是 pygtide 算出来的
做过重力观测数据处理的人都有体会:仪器记录的原始序列里,能量最大的那一项不是地壳运动,也不是仪器漂移,而是日月相对位置变化引起的潮汐形变,幅度能到200微伽量级。做地震前兆、地球自由振荡或海平面变化研究之前,第一步就是把这项理论潮汐从观测值里扣干净。pygtide 这个基于 Python 的模块,就是把经典 ETERNA PREDICT 算法封装成标准 Python 接口,配好潮波展开系数、极移、闰秒等数据文件之后,几行代码就能得到某个台站的理论引力潮时间序列。源码包里同时给了纯 Python 实现和编译好的计算扩展,适合地球物理、海洋学、大地测量领域的研究者直接拿来做数据预处理,也适合想研究潮汐预测算法工程实现的人拆开看。
2. ETERNA PREDICT 的算法骨架与 pygtide 源码工程结构
2.1 潮汐预测的物理模型:引潮位展开、Tamura 潮波表与勒夫数
地球引力潮计算的本质是求解引潮位。日月在地球表面任一点产生的引潮位可以表示成一系列频率固定的谐波之和,每个分波对应一个天文频率,这就是潮波展开。实际展开中,Tamura(1987)发表的展开表包含1200多个分波,是固体潮预测领域使用最广的展开体系之一;pygtide 源码里的tamura.py实现的正是这套展开系数表,它把每个潮波的天文幅角与振幅系数整理成程序可直接引用的数据结构。
有了引潮位还不够,地球不是刚体,地表任一点对引潮力的响应需要用勒夫数(Love numbers)修正。不同频率分波对应的勒夫数取值不同,尤其在近周日频段,近周日共振效应会让某些分波的响应与弹性地球模型有明显偏差。ETERNA PREDICT 的处理方式是把勒夫数展开成频率的函数,在预测阶段对每个分波分别应用对应的响应系数,再合成地表理论值。理解这层逻辑对后面调参很重要:如果你拿 pygtide 的计算结果去对比其他软件,差异主要来自潮波表截断阶数和勒夫数模型的选取。
2.2 源码包里的 39 个文件分别承担什么职责
拿到上传的源码包后,第一件事是先对照文件清单搞清楚模块边界。这个工程的文件虽然多,但层次非常清楚,核心计算、数据文件、构建配置和用户文档分得很开。
| 路径 | 类型 | 作用 |
|---|---|---|
earthtides/pygtide.py | Python 源码 | 主模块,封装潮汐预测主流程,提供计算入口 |
earthtides/tamura.py | Python 源码 | Tamura 潮波展开系数实现,供预测模块调用 |
earthtides/etpred.pyd | 编译扩展 | ETERNA PREDICT 核心计算引擎的 Windows 64 位版本 |
earthtides/commdat/*.dat | 数据文件 | EOP 地球自转参数、闰秒、极移、潮波参数等 |
earthtides/commdat/pygtide.waves.ini | 配置文件 | 定义参与计算的潮波目录 |
pygtide_update_data.py | 脚本 | 拉取并更新 EOP 与闰秒数据文件 |
test.py | 脚本 | 冒烟测试,验证模块计算链路是否正常 |
PREDICT/PyGTide_user-guide.pdf | 文档 | 官方用户指南 |
setup.py/MANIFEST.in | 构建配置 | 声明打包与安装规则 |
重点留意etpred.pyd的文件名后缀:cp36-win_amd64明确标记了这是 CPython 3.6、Windows 64 位环境下的编译产物。如果你当前环境的 Python 版本或操作系统不匹配,这个扩展是直接 import 不进来的。此时不必慌,pygtide.py加上tamura.py的组合提供了完整的纯 Python 计算路径,只是速度慢一些,工程上完全可以接受。
2.3 pygtide.waves.ini 潮波目录的配置机制
commdat/目录下最能体现这个模块设计思路的文件是pygtide.waves.ini。这个 INI 文件不是随便写几个分波名,而是整个预测程序的分波字典:程序启动时读取它,获取参与计算的分波频率、振幅修正标志和相位基准,然后与tamura.py里的展开表对应起来。
实际工程里,一般不会直接把所有 1200 个分波全部压进预测计算,那样计算开销不小,而且部分高频分波对重力固体潮的影响远低于观测噪声。常见做法是在waves.ini里维护一个较全的目录,而在计算时通过参数控制是否启用某些频段。修改这个文件时要注意频率参数是天文常数,不能随意改动;可以调整的是某个分波组是否参与合成,而不是去改频率值本身,否则理论潮汐序列的相位会对不上真实天文周期。
3. 环境匹配、安装导入与第一次理论潮汐计算
3.1 解包后先做版本和环境匹配检查
拿到源码包后不要急着 import。先把upload.zip解压,确认目录结构里是否有earthtides这个包目录,然后检查你的 Python 版本。如果你还没装 Python,先按当前操作系统的 python 安装教程把环境配好,推荐直接用官方安装包,避免用系统自带的老版本。
打开终端或命令行,进入项目根目录,先做一次环境探测:
python -c "import sys; print(sys.version)" python -c "import struct; print(struct.calcsize('P') * 8, 'bit')"第一行确认 Python 主版本号,第二行确认解释器位数。etpred.pyd是为 64 位 Python 3.6 编译的,所以如果你这边是 3.8、3.10 或更高的 Python 版本,这个扩展是加载不了的。探测结果会直接决定后面走哪条计算路径。
3.2 安装方式与扩展加载失败的降级策略
正常安装直接用setup.py即可:
python setup.py install或使用开发模式,方便后续改动源码直接生效:
pip install -e .安装完成后做一个快速验证:
import pygtide print(pygtide.__file__)3.2.1 验证 etpred 编译扩展是否可用
单独验证编译扩展能否加载:
import importlib ext = importlib.import_module('etpred') print("编译扩展可用:", ext)如果这一步抛出ImportError,并且错误信息里带DLL load failed或者No module named 'etpred',基本就是 Python 版本或平台位数不匹配。看到cp36-win_amd64这种文件名,可以直接判断这个编译模块只在 Python 3.6 的 Windows 64 位环境里生效。
3.2.2 纯 Python 回退路径
扩展不可用时,pygtide.py会自动回退到纯 Python 计算路径,依赖tamura.py提供的展开系数。使用上没有差别,只是大规模批量计算时耗时会明显增加。测试阶段完全够用,等到正式批处理大批台站时再考虑统一到 Python 3.6 环境还是重写扩展的逻辑。
3.3 最小可运行脚本:计算一个台站的理论固体潮序列
下面是一个最小可运行脚本,计算北京某个点位 2024 年 1 月 1 日 0 时起的理论引力潮垂直分量:
import sys sys.path.append('./earthtides') from pygtide import PyGTide tide = PyGTide() tide.set_tide_input( lat=39.9042, # 纬度,北纬为正,单位:度 lon=116.4074, # 经度,东经为正,单位:度 alt=44.0, # 台站海拔,单位:米 date='20240101', # 起始日期,格式 YYYYMMDD time='000000', # 起始时刻,格式 HHMMSS unit='ugal' # 输出单位:微伽(1 uGal = 10 nm/s^2) ) result = tide.predict( steps=1440, # 输出点数 interval=60 # 采样间隔,单位:秒 ) print(result[:3])这段代码的逻辑是:先实例化PyGTide对象,通过set_tide_input把台站坐标、起算时间和输出单位传入,最后调用predict生成理论潮汐序列。steps与interval相乘就是总时长,这里 1440 点乘以 60 秒正好是一天的逐分钟序列。返回的result是带表头的文本形式的序列,前三条打印出来可以直接看到时间戳和幅度值。
参数上提三个关键点。第一,lat和lon建议保留至少 3 位小数,经纬度 0.001 度的偏差在高频分波上会造成相位偏移;第二,alt对重力固体潮的影响比很多人想象的小,但会参与自由空气梯度修正,几十米的误差可以忽略,上千米的台站必须填准;第三,unit选ugal时结果直接对应重力仪观测单位,如果后续要和 GNSS 位移序列对比,你可能需要换算成 nm/s^2,两者差一个 10 的因子。
4. EOP 数据与闰秒更新对预测精度的影响及校验
4.1 哪些数据文件会随时间失效
潮汐预测看似只依赖日月历表,实际上还依赖地球自转参数。commdat/下那批.dat文件里,[raw]_finals2000A.dat、[raw]_eopc04_IAU2000.dat存放的是 IERS 发布的地球自转参数,包括 UT1-UTC 差值、极移分量;[raw]_Leap_Second_History.dat是闰秒历史表。这三个文件直接参与天文幅角的计算,因为潮汐分波的相位基准是基于 UT1 而不是 UTC 的。
这类文件是时效性数据,时间越长误差积累越明显。UT1-UTC 的日变化虽然只有几毫秒到几十毫秒,但换算到高频潮波相位上,会让部分分波的预测值产生微伽级的偏差。做高精度台站数据处理时,这个量级不能被忽略,因此更新数据文件和重新计算是配套动作。
4.2 用项目自带的更新脚本刷新数据
pygtide 考虑到了这个问题,源码包根目录下有pygtide_update_data.py,用来拉取最新数据替换本地文件:
python pygtide_update_data.py脚本会访问 IERS 之类的数据服务,把新的finals2000A.dat、Leap_Second_History.dat等内容下载下来,直接覆盖commdat/下的同名文件。离线环境下没法跑这个脚本,可以手动从数据源下载对应文件,替换时注意保持原始文件命名和格式不变,最好先备份旧文件再覆盖,方便回滚对比。
这个更新动作到底改变了什么,用脚本验证最直观:更新前后分别计算同一个台站同一时段的理论潮汐序列,然后做差,观察差异随时间的漂移。如果日期离数据文件末尾越远,差值越大,说明数据时效性对你的应用场景确实有影响。研究近海海平面变化这种需要长期稳定参考的场景,保证 EOP 数据新鲜是基本功。
4.3 更新数据的差分校验脚本
写一个简单的更新前后对比:假设result_old是更新前计算的序列,result_new是更新后的序列,解析数值列后做差:
import numpy as np def parse_ugal(seq): vals = [] for line in seq.strip().splitlines()[1:]: parts = line.split() if len(parts) >= 3: vals.append(float(parts[-1])) return np.array(vals) old = parse_ugal(result_old) new = parse_ugal(result_new) diff = old - new print("最大差异: %.4f uGal" % np.max(np.abs(diff))) print("均方根差异: %.4f uGal" % np.sqrt(np.mean(diff ** 2)))这段代码先从模块输出的文本序列里提取最后一列数值,再对比两个版本结果的差异。如果当前日期距离数据文件末尾时间不久,差分结果应该非常小;如果差异显著,优先确认是不是 EOP 文件长期未更新导致的。顺便说一句,用均方根而不是平均差来评估,是因为潮汐序列是振荡信号,正负偏差会抵消,均方根才能反映真实的偏离幅度。
5. 批量计算多个台站与避开常见坑位
5.1 批量生成整年理论潮汐序列
实际科研场景里,单个台站的计算很少,通常是几十个台站、一整年逐分钟的数据。批量处理时复用同一个PyGTide实例,循环里只更新坐标参数,性能会好很多:
stations = [ {"name": "WUHN", "lat": 30.53, "lon": 114.35, "alt": 30.0}, {"name": "KUNM", "lat": 25.03, "lon": 102.80, "alt": 1990.0}, ] for st in stations: tide = PyGTide() tide.set_tide_input( lat=st["lat"], lon=st["lon"], alt=st["alt"], date='20240101', time='000000', unit='ugal' ) series = tide.predict(steps=366 * 24, interval=3600) with open(st["name"] + "_tide_2024.txt", "w") as f: f.write(series) print(st["name"], "计算完成")循环里每个台站单独 new 一个实例,避免共享内部状态导致的参数污染;采样间隔设 3600 秒、点数设 366 乘 24,正好覆盖闰年整年小时值。输出到带台站名的文本文件,后续做残差分析时直接读取即可。
5.2 四个最容易踩的坑
- 时间基准混淆:
date和time传入的是 UTC 时刻,内部计算用的是 UT1。如果你输入的是北京时间,记得先减 8 小时再传入,否则整条序列会偏一个时区。 - 扩展模块版本匹配:
etpred.pyd绑定 CPython 3.6 和 Windows 64 位,换版本后 import 报错是预期行为,直接走纯 Python 路径,不要因此怀疑源码有问题。 - 单位选择不一致:和重力仪原始观测值对比时,仪器输出可能是 mV 或 10^-8 m/s^2,需要先统一单位再求残差,否则误差项里会混入一个常数倍率。
- 坐标精度不够:纬度、经度只填到小数点后两位,重力潮汐在纬度方向上梯度明显,高频分波对坐标误差的敏感度远高于低频分波,至少保留 3 位小数。
5.3 一个有效的交叉验证技巧
没有实测数据时,如何确认计算结果可信?找一个公开的在线固体潮预测服务,选同一个台站、同一时段、同一单位做交叉对比。两者在低频分波上应该高度一致,差异主要来自潮波展开的截断阶数。把差值序列画出来,如果呈现明显的周期性,说明两个算法在某个分波的勒夫数取值或频率基准上有差异;如果只是整体常偏,大概率是基准重力值或单位换算问题。这个技巧在项目验收和论文复核阶段非常实用。
本文还有配套的精品资源,点击获取