高精度重力测量的潮汐改正软件及其评估分析

2023-02-12朱传东明晓冉庞柳青田兴宇苏建文

大地测量与地球动力学 2023年1期

朱传东陈铭明晓冉庞柳青田兴宇苏建文

1 中国地震局第一监测中心，天津市耐火路7号，300180

潮汐是对日月等天体引潮力的响应，包括固体潮、海潮和大气潮，潮汐改正是重力观测数据处理的基础。重力固体潮与海洋潮汐负荷是影响重力观测的主要干扰因素，二者最大变化量分别可达300 μGal与10 μGal。通常使用引潮位对重力固体潮进行分析计算，主要计算方法有天顶距直接法和分波法[1-2]。由于潮汐效应具有相同的频率，常规的滤波方式无法分离海洋潮汐负荷，因此需要利用高精度海潮模型进行改正[3-6]。

本文基于MATLAB GUI界面自主编写一套高精度重力测量的潮汐改正软件GTIDE，并结合国际主流软件对其可靠性与适用性进行分析，以期为重力观测的数据处理与分析提供灵活可靠的软件平台。

1 数据处理方法

1.1 重力固体潮理论值

根据潮汐力学理论可知，重力固体潮理论值是将月亮与太阳的引潮位分别取至三阶与二阶，再由天顶距直接计算求得(向上为正)[7]。GTIDE软件计算给出的重力固体潮理论值为扣除永久潮汐引力直接部分后的零潮汐改正值(单位μGal)，软件可以输入由统计经验模型获得的平均潮汐因子(默认1.16)或者由固体潮观测台站数据计算得到的实测潮汐因子[7-8]。月亮与太阳的杜德森常数以及平均地球半径都与采用的参考椭球和天文参数有关，GTIDE软件数据取自1967年国际大地测量协会[7]。其他的主要参数与月球、太阳的6个轨道运动天文参数有关，在J2000.0系统下可表示为儒略世纪数(以力学时为参改)的线性函数。此外，GTIDE软件可以手动输入或者批量导入计算过程中需要用到的测站地理经纬度坐标以及观测时间等数据。

1.2 海洋潮汐负荷

海洋潮汐负荷改正可以由海潮模型与重力格林函数的褶积积分求得(向上为正)[9]：

L(φ,λ,t)=

(1)

式中，φ和λ分别为重力测站的经度和纬度，φ′和λ′分别为负荷点的经度和纬度，ρw为海水密度，r为地球半径，H(φ′,λ′)为全球或近海区域的瞬时潮高，G(ψ)是负荷点到测站角距为ψ的重力负荷格林函数。

从式(1)可以看出，海潮模型与重力格林函数是影响海潮负荷计算精度的主要因素。GTIDE软件可以选择性输入NAO99b、TPXO、FES、EOT、DTU和CSR等一系列高精度全球海潮模型或NAO99Jb近海海潮模型。选中近海海潮模型后，软件可将其与全球海潮模型相应近海区域值进行替换。GTIDE软件统一将瞬时潮高表示为8个主要调和谐波分量(M2、S2、N2、K2、K1、O1、P1和Q1)的总和，同时为保证每个分潮波均满足海水质量守恒定律，需分别扣除等厚度的水质量层以修正全球海潮模型[10]。对于重力格林函数，GTIDE软件可以选择性输入特定坐标系下的G-B、1066A、PREM、PREMsoft、PREMhard、iasp91和ak135[11]。在实际计算过程中，式(1)中的离散褶积积分区域共划分为4个部分，GTIDE软件可手动输入每个积分区域的边界范围及其离散网格大小，也可以在近区范围内选择性输入SRTM DEM数据以改善海陆边界的识别精度。为改善重力格林函数直接引力项的计算精度，采用海水质量引力效应对测站垂向分量进行计算[12]。

基于上述设置可计算出各分潮波的负荷振幅Lp与负荷相位αp，通过叠加计算可获得海潮负荷重力效应随时间变化的特征：

Gp(φ,λ,t)=

(2)

式中，fp与up分别为交点因子和交点订正角；t为观测时间，可在软件界面中手动输入或者批量导入，由此得到任意时刻的负荷值。各分潮波的负荷振幅、负荷相位以及负荷值时间序列可在选定的输出文件中保存。

2 结果与分析

2.1 重力固体潮理论值

为验证GTIDE软件计算得到的理论重力固体潮精度，将其与ETERNA软件包中PREDICT程序计算得到的结果进行比较。PREDICT程序采用分波法计算得到理论重力固体潮，在实际计算过程中可提供7种引潮位展开表。本文选用Tamura的引潮位展开表，其计算精度可达μGal级[13]。

以分布于不同纬度的4个重力测站为例(永兴岛站、厦门站、佘山站、蓟县站)，对上述2种软件计算得到的理论重力固体潮结果进行比较分析。计算过程中，将采样时间跨度设置为2021年全年，采样时间间隔设置为1 min，采样参考时间选取世界时(UTC)，平均潮汐因子设置为1.16，同时保留永久潮汐引力直接部分，并从时间序列中扣除该时间段内的平均值。表1为2种软件得到的理论重力固体潮结果，图1为8月蓟县站2种软件的理论重力固体潮及差值。从全年采样数据来看，理论重力固体潮在各个测站上的峰对峰值最大可达326.9 μGal，说明重力固体潮会对重力观测造成较大影响。2种软件在各个测站上的结果差值均在0.7 μGal以内，相关系数均为1.0，呈现出较为一致的变化特征，说明本文GTIDE软件的计算结果具有可靠性。2种软件计算结果的差异可能与选用的参考椭球、天文参数、引潮位截断以及引潮位展开表有关[14]。

表1 GTIDE与PREDICT理论重力固体潮计算结果

图1 GTIDE与PREDICT的理论重力固体潮及其差值Fig.1 The theoretical gravity solid tide and its difference between GTIDE and PREDICT

2.2 海洋潮汐负荷

为验证GTIDE软件计算得到的海洋潮汐重力负荷的精度，将其与Bos-Scherneck网站计算得到的结果进行比较分析。Bos-Scherneck网站的海洋潮汐重力负荷值同样由海潮模型与重力格林函数的褶积积分计算得到，在实际计算过程中可提供G-B、STW105重力格林函数与多种全球海潮模型的褶积计算结果。本文仍以上述4个重力测站为例(永兴岛站、厦门站、佘山站、蓟县站距海岸线分别约1 km、25 km、70 km和100 km)，对2种软件计算得到的海洋潮汐重力负荷进行比较分析。实际计算过程中，默认重力测站的高程为0 m，将重力格林函数设置为G-B模型，将海潮模型设置为FES2004模型。表2为2种软件得到的海洋潮汐重力负荷结果。可以看出，各分潮波的重力负荷最大可达3.5 μGal，说明海洋潮汐会对重力观测造成一定的影响，且这种影响会随台站与海洋之间距离的增大而减小。从各分潮波振幅结果来看，2种软件在各个测站上的差值均在0.07 μGal以内，除K2与N2分潮波外，其余占比均在自身量级的10%以内；从各分潮波的相位结果来看，除K2与N2分潮波外，其余差值基本在6.5°以内。2种软件计算结果的差异可能与海岸线数据以及褶积积分方法有关[9]。值得注意的是，K2与N2分潮波的相位差值整体较大，可能是因为二者的振幅量级相对较小，对上述影响更为敏感，但其对重力负荷计算结果的影响非常有限。

表2 GTIDE与Bos-Scherneck海洋潮汐重力负荷计算结果

图2为2021-08四个测站2种软件的海洋潮汐重力负荷时间序列，其中采样时间间隔设置为1 min，采样参考时间选取世界时(UTC)。2种软件在永兴岛站、厦门站、佘山站和蓟县站的序列差值分别在0.25 μGal、0.19 μGal、0.29 μGal和0.14 μGal以内，相关系数均为1.0左右，呈现出较为一致的变化特征，说明GTIDE软件的计算结果具有可靠性。

图2 GTIDE与Bos-Scherneck的海洋潮汐重力负荷时间序列Fig.2 The time series of ocean tide loading of GTIDE and Bos-Scherneck

基于GTIDE软件，本文进一步分析海洋潮汐模型与重力格林函数对海洋潮汐重力负荷的影响。限于篇幅，本文仅选择佘山站和厦门站的3种海洋潮汐模型与重力格林函数模型进行比较分析。表3为FES2004、FES2014b和TPXO9_atlas海洋潮汐模型的重力负荷计算结果，重力格林函数均设置为G-B模型。其中，FES2004、FES2014b与TPXO9_atlas模型的空间分辨率分别为1/8°、1/16°和1/30°。从表3可以看出，不同海洋潮汐模型的重力负荷计算结果差异较为明显。相比于FES2004，FES2014b与TPXO9_atlas在2个测站上的最大振幅差值可达0.80 μGal，最大相位差值可达49.4°(对应振幅为0.46 μGal)，上述差异在沿海且海岸线较为复杂的厦门站更加明显。综上可知，海洋潮汐重力负荷计算结果对海潮模型的精度较为敏感，因此需要采用更高精度的区域海潮模型、验潮站数据以及海岸线数据进行模型构建。

表3 不同海潮模型的海洋潮汐重力负荷计算结果

表4对不同格林函数的重力负荷计算结果进行比较分析，其中格林函数的地球模型分别为G-B、PREM和PREMhard，海潮模型均设置为FES2004。3种地球模型的区别主要体现在地壳上地幔区域，不同于G-B、PREM这2种全球平均地球模型，PREMhard是以Crust2.0固结沉积层为最外层的地球模型，能够有效改善PREM纵向地壳弹性结构参数的精度。由表可见，相比于G-B 模型，PREM与PREMhard在2个测站上的最大振幅差值可达0.14 μGal，最大相位差值可达12.6°(对应振幅为0.25 μGal)。由此可知，虽然不同地球模型得到的重力格林函数差异明显(地壳上地幔参数的改进显著改善了高阶负荷勒夫数和测站近区的格林函数)，但通过褶积积分计算可知，其对海洋潮汐重力负荷计算结果的影响有限。