简介这是一份基于MATLAB开发的理论固体潮计算程序适合地球物理、海洋与空间科学方向的学生与研究人员使用可用于模拟月球和太阳引力作用下地球表面的周期性位移响应。脚本基于牛顿万有引力定律构建天体引力模型并运用偏微分方程求解、傅里叶分析与格林函数等数值手段结合地球质量、半径、弹性模量、泊松比等地球物理参数计算不同时刻的理论固体潮位移并输出时间序列。所得结果不仅可展示潮汐周期和振幅变化也为地壳形变、地球自转速率变化等研究提供了参考。压缩包内共1个m文件体积仅1KB代码量精简方便阅读、运行与二次修改。已有561人学习下载对理解固体潮建模、引力潮汐理论及MATLAB科学计算具有直接的入门价值。1. 固体潮理论计算的工程入口为什么从 Tide2 开始做 GNSS 精密单点定位时解算出来的 N/E/U 残差里经常蹲着一条周期接近 12 小时 25 分的“波浪”。不是多路径也不是接收机噪声那是固体潮没扣干净。固体潮让测站在地固系里产生厘米级周期位移单日变化量级跟几个误差源混在一起时很难直观看见但一进频谱就格外刺眼。Tide2 这类围绕理论固体潮计算的实现要解决的正是这件事给定测站坐标和时间把引潮势展开成有限个潮汐波再乘上勒夫数输出位移或重力扰动的时间序列。这个计算的结果能直接写进 GNSS、VLBI、InSAR 和超导重力仪的数据改正链路也能用于检验有没有把不同潮汐成分混在一起。适合每天跟精密坐标、地球物理信号和数据预处理打交道的人——不需要把 Doodson 展开背下来但得知道 Tide2 在算什么、参数在哪改、输出可信到什么程度。2. 理论固体潮计算的模型骨架引潮势、潮汐表阶数与 Tide2 的参数空间2.1 从日月位置到引潮势Tide2 的天文展开基础理论固体潮计算的物理起点是引潮势。月亮和太阳对地球表面单位质量的引力减去它们对地心的引力剩下的那一阶差就形成了可以在台站上观测到的形变驱动力。实用计算不会直接对日、月位置做逐历元积分而是先把引潮势按频率展开得到一组振幅、频率和初相位都固定的潮汐波。Tide2 一类程序在运行时通常只需要拿到两个输入观测历元拆成的时间参数组以及测站的经纬度和高程。展开始终围绕 Doodson 提出的六个基本幅角展开这六个量分别对应平太阴时、月球平黄经、太阳平黄经、月球近地点幅角、月球升交点黄经和太阳近地点幅角。任意一个潮汐波的相位都可以写成它们六个的线性组合组合系数就叫 Doodson 码。算理论固体潮本质上是把这个线性组合的相位随时间推进做累加。工程实现里更常见的是把 Doodson 码预先生成一张频域表。每一条记录包含潮波名称、Doodson 码、频率和振幅系数Tide2 运行时只是查表求和。这套做法的好处是稳定天文幅角由长期多项式给出不会因为轨道积分误差导致潮汐相位漂移。展开阶数决定了表的大小也决定了计算量和精度范围。只做日常改正几十个主波就够要做地球物理研究得用上千波的表。2.2 Doodson 频率与主波群哪些波必须算哪些可以截断把潮汐波按频段归类工程上最关心的是三个频段长周期、周日和半日。长周期波振幅小但周期长容易被长期观测趋势掩盖周日和半日波振幅大是 Tide2 输出序列里最主要的波动来源。它们对应的潮汐波群是 O1、P1、K1 和 M2、S2、N2。下面这组主波在任何理论固体潮计算中都建议保留Tide2 截断参数不管怎么调这六个不要动。潮波Doodson 码频段周期小时相对振幅以 M2 为 1M2255.555半日12.42061.000S2273.555半日12.00000.466N2245.655半日12.65830.192K1165.555周日23.93450.530O1145.555周日25.81930.377P1163.555周日24.06590.176截断的影响主要体现在长周期频段。潮汐表从 60 波扩展到 1200 波增加的多半是高频微小振动和长周期组合波它们单个振幅很小但叠加后对长期序列的累积效应不可忽略。如果 Tide2 只用来做单历元改正截断到 100 波以内影响不大如果要用输出序列做频谱分析截断到 60 波会在长周期段引入虚假背景。我一般会保留 120 到 300 波这个区间能在计算速度和残余误差之间取到比较好的平衡点。表的版本也值得注意Cartwright-TaylerCTE表历史久、文档多Hartmann-WenzelHW95表波数更密两者在主波上的振幅差异能到千分之一量级对毫米级改正目标来说在误差容许范围内。2.3 勒夫数与 Shida 数弹性、非弹性地球与 IERS 约定有了引潮势的表还不够要得到地表位移得把势的变化转换成形变。这时用到的是勒夫数 h_n 和 Shida 数 l_n。h_n 描述径向位移对引潮势的响应l_n 描述水平位移的响应。理论上地球是弹性体实际地球存在滞弹性所以同样的潮汐势在不同频段上会对应略微不同的勒夫数。IERS 2010 规范给出了二阶和三阶的常用参考值参数数值说明h20.6078二阶勒夫数径向位移主项l20.0847二阶 Shida 数水平位移主项h30.292三阶勒夫数径向位移修正项l30.015三阶 Shida 数水平位移修正项Tide2 在选用勒夫数时通常默认非弹性地基模。如果将弹性地球的勒夫数 h2 0.6114 放进计算和 IERS 约定比会有约 0.6% 的振幅偏差对应到 M2 波径向位移大约是亚毫米到毫米量级。单看这个量级不算大但在长周期观测里它和参考框架的尺度因子混在一起很难单独揪出来。所以做高精度站坐标时间序列时勒夫数版本必须写进数据处理的元数据里否则后期溯源时找不出差异来源。2.4 频域法与时域法Tide2 运行效率差异的根源理论固体潮实现有两种路线。频域法是在频率域里把每个潮波的相位和振幅算好再叠加成时间序列优点是稳定、逐历元计算量固定缺点是频率表一长单历元开销线性上涨。时域法直接对日、月位置做实时计算再算出瞬时引潮势再乘勒夫数得到位移。时域法的优点是无需维护庞大的频率表但天文位置计算精度直接决定输出质量而且不同历元间的插值误差会造成谱泄漏。Tide2 的定位天然适合频域法。它输出的是一段等间隔时间序列频域展开可以保证序列在频谱上干净这对后续做傅里叶分析非常友好。如果需要在嵌入式或实时场景里运行可以把 Tide2 计算拆成三步预计算天文幅角多项式系数、载入潮汐表、对每个历元做线性求和。前两步只做一次第三步是纯加法单历元几千个波的求和也只消耗微秒级运行时间。这是它比直接调用天文库逐历元求位置更稳更快的根本原因。3. 用 Tide2 的最小骨架跑通理论固体潮计算输入、相角和输出检查3.1 最小输入集时间系统、站址与参考系约定Tide2 类计算第一个要确认的是时间系统。引潮势计算里的太阳和月亮平位置使用 TT地球时而格林尼治恒星时依赖 UT1两者不能混用。实际水文和导航设备留存的观测时间通常是 UTCUTC 转 TT 需要先加闰秒得到 TAI再加 32.184 秒UTC 转 UT1 需要从 IERS 公报里取 DUT1 修正。很多固体潮计算“看着对、细看差一截”都发生在这一步。输入坐标方面测站坐标建议直接给 ITRF 框架下的笛卡尔坐标或大地坐标让 Tide2 内部自己做参考系转换而不是手工转好再喂进去。椭球高和正高不能互换椭球高在勒夫数公式里对应到地心距差几十米会带来可忽略的振幅变化但地心距误差会直接影响引力位梯度项的计算能填对就填对。3.2 一个可运行的 Tide2 风格计算骨架下面这段代码用演示版潮汐表跑通“读取时间序列 → 计算 Doodson 相位 → 叠加潮波 → 输出位移”的完整流程。演示版只保留 6 个主波天文幅角用多项式近似相位精度约在几度量级适合验证计算流程不适合直接用于精密数据处理。import numpy as np # 演示版潮汐表名称, Doodson码, 相对振幅, 频段, 纬度因子类型 # 纬度因子类型: 0长周期(m0), 1周日(m1), 2半日(m2) waves [ (M2, (2, 0, 0, 0, 0, 0), 1.000, 2), # 255.555 (S2, (2, 2, -2, 0, 0, 0), 0.466, 2), # 273.555 (N2, (2, -1, 0, 1, 0, 0), 0.192, 2), # 245.655 (K1, (1, 0, 0, 0, 0, 0), 0.530, 1), # 165.555 (O1, (1, -1, 0, 0, 0, 0), 0.377, 1), # 145.555 (P1, (1, 1, -2, 0, 0, 0), 0.176, 1), # 163.555 ] H2, L2 0.6078, 0.0847 G 9.7983 # 均值重力生产环境应替换为测站实际重力值 D_moon 2.6242 # 月亮 Doodson 常数m^2/s^2 def basic_arguments(mjd_ut1): # 返回六个基本幅角单位度 # mjd_ut1 对应的是 UT1 儒略日 T (mjd_ut1 - 51544.5) / 36525.0 s 218.316 481267.8813 * T # 月球平黄经 h 280.466 36000.7698 * T # 太阳平黄经 p 83.353 4069.0137 * T # 月球近地点 N 125.045 - 1934.136 * T # 月球升交点 ps 282.937 1.719 * T # 太阳近地点 # 平太阴时 tau 依赖格林尼治恒星角这里用 UT1 近似 gmst 280.46061837 360.98564736629 * (mjd_ut1 - 51544.5) tau gmst 180.0 - s return np.array([tau, s, h, p, N, ps]) def lat_factor(lat_rad, mode): # 计算纬度振幅因子 if mode 0: return (3 * np.sin(lat_rad) ** 2 - 1) / 2 elif mode 1: return np.sin(lat_rad) * np.cos(lat_rad) elif mode 2: return np.cos(lat_rad) ** 2 def tide2_displacement(mjd_ut1, lat_rad, lon_rad): args basic_arguments(mjd_ut1) # 地心距按平均半径 6371 km 近似 scale H2 * D_moon / G u_r, u_e, u_n 0.0, 0.0, 0.0 for name, dod, amp, mode in waves: phase_deg np.dot(dod, args) np.degrees(lon_rad) * (mode - 1) phase np.radians(phase_deg) lat_f lat_factor(lat_rad, mode) u_r scale * amp * lat_f * np.cos(phase) # 水平分量按 l2/h2 比例近似东向调制相位滞后 90 度 u_e scale * amp * lat_f * (L2 / H2) * np.sin(phase) return u_r, u_e这段代码的核心逻辑分三层。第一层basic_arguments利用 UT1 儒略日算六个基本幅角潮汐相位不再依赖实时轨道计算而是靠长周期多项式推演这是频域法相比时域法的典型做法。第二层lat_factor把潮汐波的球谐纬度依赖抽出来周日波在 45 度纬度附近振幅最大半日波在赤道最大这是由引潮势的球谐函数展开决定的。第三层tide2_displacement做求和把每个波的 Doodson 码点乘基本幅角得到相位再用余弦叠加。代码里的phase_deg加了经度项是因为台站经度相对地轴有一个固定相位偏移周日波和半日波的影响方式不同。运行后输出的u_r是径向位移单位米正方向朝上u_e是东向位移的简化近似。这个骨架把 Tide2 输出核心结构复现清楚了生产环境替换潮汐表和幅角多项式即可。3.3 输出三维位移分量与重力扰动改正项的物理检查拿到 Tide2 输出后第一个要做的不是编程而是量级检查。以中纬度测站为例M2 波的径向位移振幅应该在 40 到 60 毫米量级K1 波在 20 到 40 毫米水平位移约为径向位移的七分之一左右。如果输出位移超过 200 毫米基本可以判定是引潮势常数或勒夫数代入错误如果序列里完全看不到半日波多数是纬度因子算错或者经度相位项缺失。重力观测里还要额外关注“重力扰动”而不是“位移”。重力固体潮的量级大约是几十微伽与位移换算关系涉及重力梯度项和自由空气改正不能直接把位移乘上重力加速度。Tide2 如果输出的是重力扰动内部通常会在勒夫数上叠加一个 δ 因子这时看的是相对于零频背景的潮汐变化检查方法与位移序列完全独立。4. 参数调整与边界参考系、截断阶数、海洋负荷和 Tide2 的常见误用4.1 ITRF 坐标输入与瞬时地固系之间的勒夫数偏差Tide2 计算输出的理论固体潮位移定义在瞬时真地固参考系里而 GNSS 解算使用的测站坐标通常在 ITRF 或其它国际参考框架下。如果不做极移和岁差章动改正直接把 ITRF 坐标当输入算出的潮汐位移在水平分量上会引入系统性偏差。极移造成的影响在地极附近最明显ITRF 和瞬时地固系的差异可达数个毫弧秒折算成水平位移能在一天的序列里产生亚毫米级残差。处理方式是在 Tide2 调用前完成参考框架转换或者在输出后把瞬时地固系位移旋转回 ITRF。后者实现简单只需要拿到极移序列构造旋转矩阵但要注意顺序先转极移再转地球自转最后转章动。Tide2 本身不负责这个转换时调用方必须持有同一个历元的极移参数否则时间戳对不上。4.2 截断阶数、潮汐表版本和频率窗口的选择潮汐表截断不是越全越好。用 1200 波的 HW95 表算单点单历元计算量比 60 波大 20 倍但位移结果差异通常在毫米以下除非研究目标本来就是亚毫米级地球物理信号否则等于白烧 CPU。我一般按用途分成三档日常 GNSS 改正用 60 到 120 波序列频谱分析用 300 波以内研究级重力数据处理才用上千波全表。频率窗口这个概念往往被忽略Tide2 输出的时间序列长度有限时窗口效应会让邻近潮波的能量互相泄漏。比如 S2 和 K2 的频率只差每周期 0.02 度左右短于 14 天的序列分不开这两个波。如果要做潮汐分解建议输出序列至少覆盖 30 天并且在分析时明确写出有效频率分辨率。这个坑在 Tide2 参数配置里看不出毛病但数据出来一算全错在窗口长度上。4.3 与海洋负荷改正的边界谁管频段谁管地域理论固体潮是地球整体对引潮势的弹性响应海洋负荷则是海潮水位变化对陆地的形变加载。两者成因不同、空间特征也完全不同。固体潮在全球每个测站都该算海洋负荷只在近海区域明显且随台站与海岸线距离快速衰减。同一观测点的固体潮位移可以精确到毫米级海洋负荷位移在大陆内部可能只有几毫米在沿海地区能达到 20 毫米以上。Tide2 计算范围不含海洋负荷使用时必须把两者分开处理再叠加。工程上常见的错误是把 Tide2 输出直接当总潮汐改正结果在沿海站点的残差里引入一个和 M2、O1 同频的伪信号。处理顺序应当是先算理论固体潮再叠加由海潮模型如 FES2014、TPXO驱动的负荷改正最后统一叠加极移潮。4.4 常见误用清单把 Tide2 用错的五种姿势误用方式后果正确做法UTC 直接当 TT 使用潮汐相位整体平移半日波相位误差可达十几度UTC 先转 TAI再加 32.184 秒不用 UT1 计算恒星时经度相位项系统性偏移影响到 0.1 毫米级从 IERS 公报取 DUT1 修正正高当椭球高地心距误差导致位梯度项偏差远小于其它误差但无法溯源坐标系统一写成椭球高同一序列混用不同潮汐表频谱上出现非物理的谱线劈裂全序列记录表的版本和截断波数勒夫数用弹性值却标 IERS 2010振幅系统性差 0.6% 且无法通过残差识别元数据里写清勒夫数来源这些误用有一个共同特征单次计算都看不出问题检查振幅、画图都对但放进长周期时间序列或跨站对比时系统偏差全部浮现。Tide2 的输入参数文档往往只强调时间和坐标真正的精度瓶颈在调用方对时间系统和参考系约定的理解。5. 验证 Tide2 结果的最后一公里回归、谱检查和批量自检理论固体潮计算的正确性验证不能只看“有没有输出”。我一般做三层检查。第一层是回归测试选三个特征测站一个在赤道、一个在 45 度纬度、一个在极地附近用 Tide2 分别计算 24 小时序列和 IERS 给出的理论数值做对比。回归阈值按位移 RMS 算径向分量差异应该小于输出振幅的 1%水平分量差异小于 2%。如果这一层不过先查时间系统再查坐标参考框架最后查潮汐表版本。第二层是频谱检查。拿 Tide2 输出的 30 天位移序列做 Lomb-Scargle 谱分析查看 M2、S2、O1、K1 频率处有没有对应谱峰以及峰之间有没有异常展宽。异常展宽通常说明时间戳不均匀或潮汐表中混入了错误的频率项。谱检查的阈值不需要太严格主波峰值应该比背景噪声至少高两个数量级否则说明截断波数太少或者相位计算有误。第三层可以写成批量自检脚本把 Tide2 计算和谱检查串起来方便每次换参数后重跑#!/bin/bash # tide2_selfcheck.sh # 参数站坐标文件、起止时间、潮汐表版本 coordsstation.list start_t2023-01-01 end_t2023-01-31 for station in $(cat $coords); do # 生成日序列并做谱峰定位 tide2_compute --station $station \ --start $start_t --end $end_t \ --tide-table hw95 --truncate 300 \ ${station}_tide.txt # 残差阈值判断与基准序列最大偏差超过阈值则报警 python - EOF import sys rec open(${station}_tide.txt).readlines() maxdev max(float(r.split()[3]) for r in rec[1:]) if maxdev 2.0e-3: print(FAIL, $station, maxdev) else: print(PASS, $station, maxdev) EOF done这个脚本把 Tide2 的调用参数和验证逻辑绑在一起。station.list里每行写一个站点名脚本循环计算并打印 PASS/FAIL。阈值2.0e-3对应 2 毫米的最大偏差如果某站持续 FAIL优先怀疑这个站的坐标参考框架或高程输入而不是潮汐表。最后一层检查做完再进生产链路Tide2 计算的输出才敢直接拼进后续的数据处理管线。本文还有配套的精品资源点击获取