0 引言

精密引潮力的计算是地球物理、大地测量、天体测量、空间科学等相关学科研究的基本理论问题(张捍卫等, 2005; 郗钦文, 2007; 张晶等, 2007; 徐建桥等, 2014).张捍卫等(2004)根据弹性力学的基本理论，从地球的微小弹性运动方程出发，给出了引潮力的严格数学定义，抛弃了前人对引潮力定义所给定的所有假设条件，回避了引潮力不同定义容易引起的歧义性.Hartmann和Wenzel (1995a, b)、Kudryavtsev (2004)基于引潮位理论公式和不同的天体历表，将月球、太阳、近地行星的引潮位分别展开至6阶、3阶、2阶，计算得到引潮位及引潮力的时间序列.地球扁率对引潮力的计算存在影响(Bartels, 1957; Wilhelm, 1983; Dahlen, 1993)，郗钦文(1999)和张捍卫等(2010)分别指出，地球扁率对引潮力的影响不超过5×10^-11 ms^-2(10^-11 ms^-2=1 nGal)，这在高精度的潮汐分析、引潮力计算中是必须要考虑的；Hartmann和Wenzel (1995a, b)通过对公式中参数的数量级估计，表明当引潮力计算截断阈值为0.001×10^-11 ms^-2时，需要考虑地球扁率对月球和太阳引潮力的影响，对其他行星的影响则不需考虑.通常情况下，引潮力被认为是引潮位的径向梯度值，而重力仪的测量数据是基于铅垂线方向，因此需将径向元素转换为铅垂线方向的元素，Wenzel (1974)指出此处可忽略垂线偏差的影响，只需将基于球面的径向元素转换为基于椭球面的法向元素即可.

在引潮力的计算过程中，天球参考系转换、勒让德函数计算是两大核心工作，在以往的计算中(Hartmann和Wenzel, 1995a, b; Roosbeek, 1996; Kudryavtsev, 2004)，对于前者采用经典的基于春分点的岁差章动转换，对于后者采用手工推导出各阶次的具体表达式；这种做法存在两个不足：一是缺乏检核条件，无法快速确定计算过程与结果的准确性；二是勒让德函数表达式推导的难度随着阶次升高将会倍增.

在本文中，对于天球参考系转换，分别采用基于春分点的岁差章动转换、基于CIO的无旋转原点转换两种方法，转换参数分别采用IERS 2003、2010规范中推荐的4组参数模型；对于勒让德函数计算，分别采用勒让德函数的递推算法、完全规格化缔合勒让德函数的递推算法两种数学模型.通过以上方法确保了引潮力计算过程与结果的准确性.在天球参考系转换过程中，涉及到历表更新、参考框架偏差、岁差章动模型更新、极移、时间尺度转换等内容，都会对引潮力的计算带来影响.另外，地球扁率、径向与法向的转换也会对引潮力的计算产生影响，本文对相关因素对引潮力计算的影响进行了计算与分析，并求得了具体的影响数值.

1 理论公式

在图 1所示的地心天球中，某历元天体对地球上某测站点的引潮位V^〈1〉为(Melchior, 1983)：

(1)

(1)式中，GM_J为万有引力常数与天体J的质量之积，R_J、r分别表示天体、测站点的地心距，Z_J为天体与测站之间的地心天顶距，P_n(x)为n阶勒让德函数.

由球面三角学中边的余弦公式得：

(2)

(2)式中，(α_J, δ_J)、(λ, φ)分别表示天体、测站点在ITRS (International Terrestrial Reference System，国际地球参考系)中的地心经度、地心纬度，并有(Hofmann-Wellenhof and Moritz, 2006)：

(3)

(3)式中P_n^m(x)为n阶m次完全规格化缔合勒让德函数.

则(1)式可表达为

(4)

(4)式中a为地球参考椭球长半径.

地球扁率对引潮位的影响V^〈2〉为(郗钦文, 1999; 张捍卫等, 2010)：

(5)

(5)式中，J_2⊕为地球重力场二阶带球谐系数.

则天体对测站点总的引潮位V为V^〈1〉和V^〈2〉之和，即(4)式和(5)式之和.

记：

(6)

则g_{r, S}、g_{λ, S}、g_{φ, S}分别表示测站点在球坐标系下的径向引潮力(径向重力固体潮)、径向经度固体潮、径向纬度固体潮.

为了确保引潮位及引潮力计算过程的准确性，本文同时还采用(7)式对V进行检核计算，公式为

(7)

(7)式中的cosZ_J由(2)式给出，此时g_{r, S}、g_{λ, S}、g_{φ, S}的计算公式可由(6~7)式得到.

记Δφ为测站点处的大地纬度B与球心纬度φ之差(即径向与法向的夹角)，则测站点在椭球坐标系下的引潮力(法向重力固体潮)g_{r, E}、法向经度固体潮g_{L, E}、法向纬度固体潮g_{B, E}为(Wenzel, 1974)：

(8)

由上可知，天体在ITRS中的地心坐标、勒让德函数及其导数、完全规格化缔合勒让德函数及其导数的计算是引潮力计算工作的基础，下面将对它们分别进行介绍.

2 天球参考系转换

天体坐标基于美国JPL (Jet Propulsion Laboratory，喷气推进实验室)发布的DE (Development Ephemeris，展开历表)历表计算得到的，以地球为中心星、各天体为目标星的历表计算结果是天体在GCRS (Geocentric Celestial Reference System，地心天球参考系)中的坐标，还需要将其转换为ITRS中的坐标.

DE历表中的时间尺度为TDB (Barycentric Dynamical Time，质心力学时)，ITRS、岁差章动模型中的时间尺度均为TT (Terrestrial Time，地球时)，TDB与TT之间的差值不超过±1.7ms (millisecond，毫秒)，此偏差在历表计算中可忽略不计(Soffel and Langhans, 2013; 中国科学院紫金台天文台, 2014)，在本文计算中，时间尺度均为TT.

2.1 天球参考系转换模型

GCRS与ITRS间有两种转换方法：基于春分点的岁差章动转换、基于CIO (Celestial Intermediate Origin，天球中间零点)的无旋转原点转换，IERS 2003、2010规范针对两种转换模型分别推荐了相应的转换参数，两种转换方法相应的转换过程如图 2所示：

在图 2的两种转换过程中，涉及到三个中间参考系：瞬时真赤道参考系、CIRS (Celestial Intermediate Reference System，天球中间参考系)、TIRS (Terrestrial Intermediate Reference System，地球中间参考系).其中瞬时真赤道参考系的Z轴指向CIP (Celestial Intermediate Pole，天球中间极)，X轴指向瞬时真春分点；CIRS的Z轴指向CIP，X轴指向CIO；TIRS的Z轴指向CIP，X轴指向TIO (Terrestrial Intermediate Origin，地球中间零点).

2.1.1 基于春分点的岁差章动转换模型

在t时刻，基于春分点的岁差章动转换的数学模型为(McCarthy and Petit, 2004; Petit and Luzum, 2010)：

(9)

(9)式中t为从J2000.0起算的TT儒略世纪数(以下相同)，[ITRS]、[GCRS]分别表示相应参考系中的三维直角坐标向量，、、分别称为框架偏差-岁差-章动矩阵、地球自转矩阵、极移矩阵，三个矩阵的计算模型分别为

(10)

在中，为参考框架偏差矩阵(简称框架偏差)，为岁差矩阵，为章动矩阵，它们的计算公式分别为

(11)

在(11)式中，ξ₀、η₀表示CIP在GCRS中J2000.0时刻的天极偏差，dα₀为瞬时平春分点的赤经偏差值；ψ_A、ω_A为基本赤道岁差参数，χ_A为黄道岁差导出量，ε_A为瞬时黄赤交角；Δψ为黄经章动，Δε为交角章动.

在中，GST (Greenwich Sidereal Time，格林尼治真恒星时)可根据GMST (Greenwich Mean Sidereal Time，格林尼治平恒星时)、EE (Equation of the Equinoxes，赤经章动)以及EECT (Complementary Terms of the Equation of the Equinoxes，赤经章动补充项)计算公式为

(12)

在中，x_p、y_p是CIP在ITRS中极移的两个分量，s′为TIO定位角，表示TIO在CIP赤道上的位置.计算公式分别为

(13)

(14)

(13)式中等号右端三项分别表示IERS (International Earth Rotation and Reference Systems Service，国际地球自转与参考系服务)数据中心发布的极移数据、由海洋潮汐引起的极移量、由章动引起的极移量.

2.1.2 基于CIO的无旋转原点转换模型

在t时刻，基于CIO的无旋转原点转换的数学模型为(McCarthy and Petit, 2004; Petit and Luzum, 2010)：

(15)

(15)式中、也分别称为框架偏差-岁差-章动矩阵、地球自转矩阵，计算模型分别为

(16)

(17)

(16)式中，s为CIO定位角，表示CIO在CIP赤道上的位置，ERA (Earth Rotation Angle)为地球自转角.

极移矩阵的计算模型同上.

2.2 转换参数及其计算

在两种天球参考系转换模型中，涉及到众多转换参数，IERS 2003、2010规范对这些转换参数分别推荐了相应的计算模型，现概述如下，分别列于表 1、表 2中，具体内容可查阅IERS规范(McCarthy and Petit, 2004; Petit and Luzum, 2010)，相关数据可从ftp://tai.bipm.org/iers或ftp://toshi.nofs.navy.mil下载获取.

2.3 计算过程中的注意事项

在两种参考系转换模型中，均涉及到了ERA的计算，两个规范中的计算公式均为

(18)

(18)式中T_u是从J2000.0起算的儒略日数，时间系统为UT1(Universal Time, 世界时)，这就涉及到TT与UT1之间的转换，记二者的差值为ΔT=TT-UT1.

ΔT采用雷伟伟等(2015)给出的经验公式计算得到，该经验公式是基于最小二乘法拟合得到一组5次多项式，不仅适用于公元纪年、JD、MJD以及从J2000.0起算的儒略世纪数等四种时间格式，便于程序设计，而且精度远高于目前国际上通用的Morrison & Stephenson经验公式(李广宇, 2015)的结果.

两种坐标转换模型实现的关键在于转换参数的计算，参数计算形式主要有多项式、周期项、泊松项三类，而周期项又可视为泊松项的特例，因此泊松项的计算就是参数计算的关键所在.在程序设计时，对IERS两个规范中的章动模型、EECT、X、Y、s、(x_p, y_p)_{ocean tides}、(x_p, y_p)_nutation中的展开项部分进行了数据格式的预处理，将其中的周期项、泊松项统一表达为的形式，s_{i, j}、c_{i, j}为相应的系数，并注意统一系数的单位，arg为基本天文辐角量的线性组合.

3 勒让德函数及其导数的递推算法

在P_n(sinθ)、P_n^m(sinθ)中，θ为图 1中的地心纬度，θ∈[-90°, 90°](即对于天体θ=δ、对于测站θ=φ).在计算引潮力的过程中，还需要计算P_n(sinθ)对sinθ的一阶导数P′_n(sinθ)，以及P_n^m(sinθ)对θ的一阶导数P_n^m(1)(sinθ).

3.1 P_n(sinθ)、P′_n(sinθ)的递推算法

P_n(sinθ)存在如下递推关系(Hofmann-Wellenhof and Moritz, 2006)，公式为

(19)

(19)式中n≥2.

P′_n(sinθ)的递推公式可由(19)式获得.

3.2 P_n^m(sinθ)、P_n^m(1)(sinθ)的递推算法

P_n^m(sinθ)、P_n^m(1)(sinθ)的递推算法主要有标准向前列递推算法和标准向前行递推算法，其中前者在适用范围、计算速度等方面均优于后者(Holmes and Featherstone, 2002; Fantino and Casotto, 2009; 雷伟伟等, 2016)，在本文计算中，采用标准向前列递推算法，算法公式为

(20)

(21)

(22)

(23)

(20)式中n≥2且n>m；(21)式中m≥2；(22)式中n>m；(23)式中m≥0.

递推初值为

(24)

4 影响因素计算与分析

Hartmann and Wenzel (1995a, b)首先以德国Black Forest Observatory (BFO) Schiltach观测站为例进行引潮位展开研究与计算，随后Roosbeek (1996)、Kudryavtsev (2004)均以此测站为例进行研究与计算，为了与国际研究保持一致，本文也沿用此测站位置坐标进行计算.BFO测站点的球心坐标为λ=8.3300°E、φ=48.139419°N、r=6366836.969 m，并有Δφ=B-φ=0.191181°.

各天体在GCRS中的坐标由DE431历表计算得到，然后根据天球参考系转换模型，得到天体在ITRS中坐标(α_J, δ_J, R_J)，再基于勒让德函数及其导数、完全规格化缔合勒让德函数及其导数的递推公式，分别基于(4~6)式、(6~7)式计算得到V、g_{r, S}、g_{λ, S}、g_{φ, S}，再由(8)式转换为法向结果，得到g_{r, E}、g_{L, E}、g_{B, E}的具体数值.参数a、J_2⊕采用IERS 2003、2010规范中的推荐值，各天体GM_J值采用DE431历表“头文件”中给出的天文常数数值.

在GCRS与ITRS的两种坐标转换模型中，均需要知道极移参数x_p、y_p的具体数值，但由于地球自转的不规则性和不可预知性，目前无法准确预测二者的具体数值，只能根据IERS公布的EOP模型中已发布的数据进行计算.目前最新的EOP 08 C04模型中，发布了从1962年1月1日至今每日一组的数据，且每天进行更新，由于是后处理数据，故在时间上延迟一个月.在本文计算中，时间范围均从1962年1月1日至2015年8月1日，该范围与当前EOP 08 C04模型发布数据的时间范围一致，计算时间间隔均为1小时，分别计算得到V、g_{r, E}、g_{L, E}、g_{B, E}序列.

4.1 各天体展开阶次及天体取舍

在进行精密引潮力的计算工作之前，首先需要确定参与计算的天体及其展开阶数，基于现阶段重力仪测量精度水平，并兼顾重力仪未来的发展，以≥0.001×10^-11 ms^-2 (0.001 nGal)为计算截断阈值，由于：

(25)

取r_max=a，对于x∈[-1, 1]，由(19)式可得≤1，对于天体J而言，其第n阶参与计算的条件为

(26)

基于DE431历表，计算得到1000AD~3000AD共计2000年内时间间隔为1天的各天体地心距序列，统计得到最小地心距min (R_J)列于表 3中.

根据表 3及DE431历表“头文件”中给出的GM_J数值，即可得到各天体在各阶次引潮力最大值.

由表 4可得，在精密引潮力计算时，月球需计算至6阶，太阳需计算至3阶，水星、金星、火星、木星、土星需计算至2阶，而天王星、海王星、冥王星则不参与计算.

在以下计算过程中，理论上引潮力的计算至0.001×10^-11 ms^-2即可，但为了分析各类影响因素对引潮力计算影响的具体量值，实际计算至0.000001×10^-11 ms^-2.

4.2 计算过程准确性的验证

为确保计算过程与结果的准确性，对天球参考系转换，分别采用基于春分点的岁差章动转换、基于CIO的无旋转原点转换两种方法，转换模型中的参数又分别基于IERS 2003、2010规范推荐的四组参数计算模型；对引潮力计算公式，分别采用基于勒让德函数及其导数递推算法、完全规格化缔合勒让德函数及其导数递推算法两种数学模型得到的(4~6)式、(6~7)式计算得到.

在下列各项影响因素计算过程中，均采用(4~6)式、(6~7)式分别独立计算得到V、g_{r, S}、g_{λ, S}、g_{φ, S}序列，然后由(8)式转换为V、g_{r, E}、g_{L, E}、g_{B, E}序列.计算结果表明：对所有序列而言，两组公式结果间的差值ΔV、Δg_{r, E}、Δg_{L, E}、Δg_{B, E}序列均为零，这就从勒让德函数递推计算方面确保了计算过程与计算结果的准确性.

在天球参考系转换方面，下文(见本文4.4节)计算结果表明：两个规范中转换参数模型间的差异、两种坐标转换方法间的差异对引潮力计算的影响非常微小，可完全忽略不计，这就从天球参考系转换方面确保了计算过程与计算结果的准确性.

4.3 地球扁率的影响

引潮位V^〈1〉是将地球作为均质球体得到的，而真实的地球则是一个两级略扁的旋转椭球体，因此地球扁率就对引潮位的计算带来影响.基于(5~6)、(6~7)式，计算得到由地球扁率而引起的各天体对测站引潮位与引潮力影响的具体数值.

由表 5可得，由地球扁率引起的月球对测站引潮位与引潮力的影响分别在0.00013 m²s^-2、1.890×10^-11 ms^-2、1.994×10^-11 ms^-2、1.909×10^-11 ms^-2以内；由地球扁率引起的太阳对测站引潮位与引潮力的影响分别在0.00000011 m²s^-2、0.002×10^-11 ms^-2、0.002×10^-11 ms^-2、0.002×10^-11 ms^-2以内；其他天体的影响均为零.

因此，在精密引潮力的计算过程中，只需考虑地球扁率对月球、太阳的影响即可.

4.4 岁差章动模型更新及两种天球参考系转换方法的影响

前已述及，针对两种坐标转换方法，IERS 2003、2010规范分别推荐了两套4组不同的转换参数，为明确岁差章动模型更新、不同坐标转换方法间的差异对引潮力计算的影响，分别基于不同的岁差章动模型和坐标转换方法，计算得到相应序列.为便于表达，本文对坐标转换模型进行如下命名：

分别基于2010A、2010B、2003A、2003B四种坐标转换模型，计算得到V、g_{r, E}、g_{L, E}、g_{B, E}各自的四组序列值，差值序列即反映了岁差章动模型更新、不同坐标转换方法间差异对引潮位与引潮力计算的影响，结果列于表 7.

由表 7可得，岁差章动模型更新(2010A-2003A、2010B-2003B)对引潮位与引潮力的影响分别在2×10^-8 m²s^-2、0.00061×10^-11 ms^-2、0.00047×10^-11 ms^-2、0.00040×10^-11 ms^-2以内；不同坐标转换方法(2010A-2010B、2003A-2003B)的影响量值分别在2×10^-10 m²s^-2、0.000005×10^-11 ms^-2、0.000004×10^-11 ms^-2、0.000003×10^-11 ms^-2以内.

这说明岁差章动模型更新、两种坐标转换方法间的差异对引潮位与引潮力计算的影响甚微，实际计算时可完全忽略不计.在引潮位与引潮力的计算中，可采用任意岁差章动模型及坐标转换方法.

4.5 径向法向转换的影响

由(4~6)式、(6~7)式计算得到的g_{r, S}、g_{λ, S}、g_{φ, S}是径向梯度结果，需根据(8)式将其转换为法向结果.若忽略径向到法向的转换，就对引潮位与引潮力的计算带来影响.

由表 8可得，若忽略径向到法向的转换，对V和g_{L, E}的数值没有影响，对g_{r, E}、g_{B, E}影响的最大值分别达到476×10^-11 ms^-2、500×10^-11 ms^-2，该影响量值远远超出允许的精度范围，即计算结果错误.

因此，在引潮力计算中，一定不能忽略将径向结果转换为法向结果，否则结果错误.

4.6 参考框架偏差的影响

参考框架偏差是由于在标准历元J2000.0时刻岁差章动模型给出的天极、分点与GCRS的天极、经度零点不重合造成的(Hilton and Hohenkerk, 2004).IERS 2003、2010规范中考虑了参考框架偏差对天球参考系转换的影响，但在IERS 1996规范中并没有考虑该影响.

以2010A模型为例，计算了在基于春分点的岁差章动转换过程中，当忽略参考框架偏差时，对引潮位与引潮力计算影响的具体数值.

由表 9可得，在天球参考系转换过程中如果忽略参考框架偏差的影响，对引潮位与引潮力的影响量值分别在1×10^-6m²s^-2、0.032×10^-11 ms^-2、0.032×10^-11 ms^-2、0.021×10^-11 ms^-2以内.

4.7 地球时与世界时转换的影响

在计算ERA时，时间尺度为世界时UT1，如果误将时间尺度用作TT，即忽略了TT与UT1的差值ΔT，就会对引潮位与引潮力计算结果带来影响.

由表 10可得，如果忽略地球时与世界时间的转换，对引潮位及引潮力的影响量值分别达到0.03 m²s^-2、916×10^-11 ms^-2、1144×10^-11 ms^-2、696×10^-11 ms^-2，该影响量值远远超出允许的精度范围，即计算结果错误.

因此，在引潮力计算中，一定不能忽略地球时与世界时转换的影响，否则结果错误.

4.8 极移的影响

由于极移量x_p、y_p的具体数值目前无法准确预测，只能根据EOP 08 C04模型中已发布的相关数据进行计算.中国天文年历(中国科学院紫金台天文台, 2014)在进行天球参考系转换时，并没有考虑极移的影响，实际上是将天体坐标转换至TIRS.由于在两种坐标转换方法中，极移影响与岁差章动模型无关，因此极移对坐标转换的影响量对两种坐标转换方法、两个规范而言均是相同的.本文计算了当忽略极移时，对引潮位与引潮力计算影响的具体数值.

由表 11可得，如果忽略极移的影响，对引潮位与引潮力的影响量值分别在2×10^-5m²s^-2、0.586×10^-11 ms^-2、0.543×10^-11 ms^-2、0.341×10^-11 ms^-2以内.

4.9 历表更新的影响

随着天体测量理论与技术的进步，尤其是精度更高的天体测量数据的积累，JPL不断的对DE历表进行更新，发布了一系列版本的历表，其中比较著名的有DE200、DE403、DE405、DE406、DE413、DE421、DE423等，目前最新的是2014年发布的DE432历表.在IERS 2003、2010规范中，分别将DE405、DE421历表作为ICRS的动力学实现.

在后续的引潮位频谱法展开研究中，采用DE431历表，其覆盖时间范围为30000年(13200BC~17191AD).为明确历表更新对引潮力的影响，分别基于DE405、DE421、DE431历表，计算了相应的序列，并统计出各序列间的差异.

由表 12可得，历表更新对引潮位与引潮力的影响量值分别在4×10^-7m²s^-2、0.012×10^-11 ms^-2、0.013×10^-11 ms^-2、0.013×10^-11 ms^-2以内.

历表更新对引潮力的影响有两方面原因：一是由于历表自身精度的提高；二是不同历表采用的天文常数系统不同，在DE405、DE421、DE431三个历表中，各天体GM值以及定义距离的AU (Astronomical Unit，天文单位)都是不同的.

5 引潮位与引潮力的具体数值

在确保计算过程的准确性及分析各类影响因素对引潮力计算影响量值的基础上，以德国BFO为例，采用(4~6)式，基于DE431历表计算得到各个天体坐标，并通过基于春分点的岁差章动转换方法进行转换，参考框架偏差、岁差章动模型采用IERS 2010规范推荐值，极移数据基于EOP 08 C04模型，各天体GM采用DE431历表“头文件”中的数值，最终计算得到1962年1月1日至2015年8月1日之间每隔1小时的V、g_{r, E}、g_{L, E}、g_{B, E}序列，各序列的统计结果见表 13.

为能更清晰直观的看出V、g_{r, E}、g_{L, E}、g_{B, E}数值及其随时间的分布情况，现将四个序列分别绘于图 3~6中.

6 结束语 6.1

天球参考系转换、勒让德函数的计算是引潮力计算的两个核心工作.为确保引潮力计算过程与结果的准确性，对于天球参考系转换，分别采用基于春分点的岁差章动转换和基于CIO的无旋转原点转换两种方法，转换参数又分别采用IERS 2003、2010规范推荐的4组参数模型；对于勒让德计算，分别采用勒让德函数及其导数递推算法、完全规格化缔合勒让德及其导数的标准向前列递推算法两种函数计算模型.

6.2

在1962年1月1日至2015年5月1日间的54年内，位于地球中纬度地区的德国BFO测站，V、g_{r, E}、g_{L, E}、g_{B, E}的量值分别在4.7578173471 m²s^-2、149459.841×10^-11 ms^-2、142582.390×10^-11 ms^-2、142717.163×10^-11 ms^-2以内.

6.3

岁差章动模型更新、两种坐标转换方法间的差异对引潮力计算结果的影响可完全忽略不计；地球扁率、极移、参考框架偏差、历表更新对引潮力的影响分别在1.891×10^-11 ms^-2、0.586×10^-11 ms^-2、0.032×10^-11 ms^-2、0.012×10^-11 ms^-2以内，这些因素在精密引潮力的计算中不能忽略；径向法向转换、地球时与世界时转换对引潮力的影响分别达到476×10^-11 ms^-2、916×10^-11 ms^-2，如果忽略这两个因素，计算结果错误.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！