变密度球体重力位的推导及其球谐展开

范广学，钟 振

(贵州师范大学 物理与电子科学学院系，贵州 贵阳 550025)

由于天体内部密度分布不均，在天体的外部空间会引起相应的重力变化，通常称这种现象为重力异常[1-3]. 重力异常对实际应用具有重要的影响，比如对空间探测器的精度定轨，复杂的重力异常容易导致探测器的坠毁[4,5]. 但天体重力异常也有重要的应用价值，由于重力异常是天体内部密度不均的直接反映，可应用于天体内部密度的约束. 在采矿工程中，重力异常还被应用于地底剩余质量的反演，剩余质量其实就是变密度与常密度的差值[6-8]. 在地球物理和行星物理研究中，也经常使用重力异常约束天体内部结构[9,10]. 文献[10]给出了常密度下，表面地形与相应重力异常傅里叶变换的关系，该关系经严格证明具有一定的合理性，被广泛应用于重力反演[11,12]. 由于文献[10]给出的模型只适用于平面模型，对于球体模型特别是小天体(如月球和火星)将产生较大的误差，为此，文献[13]基于文献[10]，给出了球体模型下地形与重力异常两者间球谐展开的广义傅里叶关系，并被广泛地应用于地球、行星和月球物理研究中[14-16].

近年来，随着月球和火星重力数据分辨率的提高[18,19]，利用它们的表面重力异常约束其内部结构的研究越来越多[14-16]. 文献[13]给出的算法是基于常密度假设，而近年来的研究表明变密度模型更具普适性[15]. 文献[15]在研究火星物理参数时，为了估算变密度模型的重力异常，提出了一种较复杂的数值积分方法. 该方法尽管推理严谨，计算精度高，但计算速度较慢，效率低下.基于广义傅里叶变换的球谐算法，被广泛应用于行星和月球物理研究中[13-17]. 球谐算法的最大优点是运行速度明显优于普通数值积分方法[20]. 目前，尚无研究考虑变密度球体重力位球谐展开的算法，针对该问题，本文基于广义傅里叶展开的球谐算法，推导变密度与重力位两者间球谐系数的关系，以期为变密度重力正反演问题提供一定的参考.

1 变密度球体模型的重力位

变密度球体模型如图1所示，O表示球心，r表示球体的最大半径.球体内部任意一点P(r′,θ′,φ′)的密度Δρ(θ′,φ′)，仅与经度φ′和余纬度θ′有关，该点对球体外任意点Q(R,θ,φ)的重力位有贡献，将球体内所有点进行积分，可得变密度球体对其外部点Q(R,θ,φ)的重力位为

U(R,θ,φ)=

(1)

其中，γ表示两向量OP和OQ之间的夹角，d表示点P和点Q之间的距离，有如下关系：

cosγ=cosθ′cosθ+sinθ′sinθcos(φ′-φ)

(2)

(3)

根据勒让德母函数的定义[21]，可得

(4)

联合式(1)—式(4)，顾及密度在径向上无变化，并在径向进行积分，可得变密度球体对其外部任意点Q(R,θ,φ)的重力位为

(5)

图1 变密度球体模型示意图

2 重力位的球谐展开

参考文献[1-3,21]，勒让德多项式Pn(cosγ)可以表示成规格化球谐函数的线性组合，即

(6)

(7)

如图1所示，密度仅是余纬度和经度的函数，与半径无关，参考文献[1-3]，可将此类函数表示成规格化球谐函数的线性组合，即

Δρ(θ′,φ′)=

(8)

(9)

(10)

定义立体角微元dσ=sinθ′dθ′dφ′，参考文献[1-3]，规格化球谐函数具有正交性特征，即

(11)

(12)

其中，狄拉克函数δnl和δmk取值为

(13)

将式(6)—式(8)代入式(5)，并考虑到规格化球谐函数具有式(11)、式(12)的正交性特征，可得

U(R,θ,φ)=

(14)

式(14)表明，变密度异常体对其外部任意点的重力位，可以直接表示成变密度球谐展开系数的组合形式.变密度异常体对半径为R的球面的重力位U(R,θ,φ)，可以直接表示成变密度球谐展开系数的组合形式.变密度异常体对半径为R的球面的重力位，参考文献 [1-3, 10, 13]可将其表示成规格化球谐函数的线性组合，即

U(R,θ,φ)=

(15)

(16)

考虑实际重力位一般只能展开成有限项之和[1-3,10,13]，令其最大展开阶数为lmax，可得变密度球体对其外部任意点Q(R,θ,φ)的重力位为

(17)

在实际应用中，如果不使用式(17)的算法，也可以采用数值积分方法对式(1)进行计算.考虑式(2)、式(3)，可将式(1)表示成如下形式：

(18)

其中，积分函数为

f(r′,θ′,φ′)=

(19)

3 算法测试与应用

为了验证本文推导结果的有效性，假定变密度Δρ(θ′,φ′)=sinθ′ cos2φ′，利用式(9)、式(10)算出密度的球谐展开系数，然后利用式(16)、式(17)计算出重力位的球谐展开系数及重力位.利用式(9)、式(10)计算密度球谐展开系数时，需要计算积分值，如果利用普通数值计算方法会花费较多计算时间，但考虑到式(9)、式(10)的积分函数涉及球谐函数，可以利用快速傅里叶变换进行抽样[22]，因此，式(9)、式(10)在计算时效率仍然较高.在已知变密度球谐展开系数的情况下，利用式(17)，可以快速求出全球不同位置点的重力位.而使用普通数值方法的式(18)计算时，由于每个位置点均要进行积分计算而耗费大量时间，计算全球不同位置点的重力位将是一件冗长的事情.

为了方便比较，如表1所示，取引力常数G=1，r=1 m，R=2 m，经度和纬度间隔Δθ=Δφ=0.5°，即每隔0.5°计算一个点的重力位.参考文献[1-3, 11]，由于Δθ=Δφ=0.5°，可知式(17)中最大展开阶次lmax=179.本文推导的球谐算法，采用单核估算，耗时主要由式(9)、式(10)和式(17)产生，经多次试验，发现耗时1.7 s左右.采用普通数值方法估算式(18)、式(19)比较耗时，为此，采用MATLAB优化的三维数值积分函数integral3，并使用与球谐算法相同的处理器进行20核并行计算，经多次试验，发现耗时在1 120 s左右.由此可知，本文推导的球谐算法在计算效率方面远胜于普通数值积分方法，对实际问题的大规模重力正反演具有一定的参考价值.

两种算法产生的重力位如图2所示，其中图2(a)表示基于本文球谐算法的结果，图2(b)表示基于式(18)、式(19)的数值积分结果，图2(c)表示两者差值的绝对值.图中横轴表示经度的取值，纵轴表示纬度的取值(余纬度与纬度互为余角). 由图2(a)、2(b)可以看出，两种算法得到的重力位分布完全一致. 图2(c)表明两种算法重力位差值的绝对值，在赤道附近最小，而在两极附近较大，但最大值不超过0.001 5，因此，本文推导的球谐算法具有一定的合理性.

图2 不同算法的重力位及其偏差的绝对值

表1 参数取值及计算时间

应用本文推导的式(17)，参考文献[23]提供的月壳密度分布数据，图3(a)和图3(b)分别给出了月壳密度和相应的重力位分布，投影中心坐标(180°W, 0°N)，其中圆圈1表示月球南极艾德肯盆地(South Pole-Aitken)，圆圈2表示月球背面高地.在计算过程中，参考文献[23]的最新成果，取引力常数G=6.673 981011m3·kg-1·s-2.参考文献[24]，取待估重力位参考球面的半径R=1 738 km，取月壳层厚度Tc=40 km，以及月壳层外半径r2=1 737.15 km.为了计算月壳层的重力位，先计算外半径r2球体的重力位，再计算内半径r1=r2-Tc=1 697.15 km球体的重力位，两球体重力位之差即是月壳层的重力位.图3(a)给出的圆圈1显示出高密度物质，图3(b) 相应地也显示较大的重力位.图3(a)显示月球背面高地平均密度最小(如圆圈2所示)，图3(b)显示该区域的重力位也最小.除此之外，其他区域的密度分布与相应的重力位分布规律总体一致，表明本文推导的球谐算法具有一定的参考价值.

图3 月壳层密度及40 km壳层厚度的重力位分布

4 结论

因天体内部密度分布不均，天体外部重力位具有各向异性的特性.普通数值积分方法估算变密度球体重力位时，计算时间较长，效率低下，为此，本文提出基于广义傅里叶变换的球谐算法.结果表明，本文球谐算法估算的重力位分布，与普通数值积分方法所得一致，而计算时间明显低于普通数值积分方法.将本文球谐算法应用于月壳变密度问题，发现变密度月壳的重力位分布与密度分布总体一致，表明本文球谐算法具有一定的合理性和实用价值，可为变密度球体重力正反演问题提供有益的参考.