變密度球體重力位的推導(dǎo)及其球諧展開

范廣學(xué)，鐘 振

(貴州師范大學(xué) 物理與電子科學(xué)學(xué)院系，貴州 貴陽 550025)

由于天體內(nèi)部密度分布不均，在天體的外部空間會引起相應(yīng)的重力變化，通常稱這種現(xiàn)象為重力異常[1-3]. 重力異常對實(shí)際應(yīng)用具有重要的影響，比如對空間探測器的精度定軌，復(fù)雜的重力異常容易導(dǎo)致探測器的墜毀[4,5]. 但天體重力異常也有重要的應(yīng)用價(jià)值，由于重力異常是天體內(nèi)部密度不均的直接反映，可應(yīng)用于天體內(nèi)部密度的約束. 在采礦工程中，重力異常還被應(yīng)用于地底剩余質(zhì)量的反演，剩余質(zhì)量其實(shí)就是變密度與常密度的差值[6-8]. 在地球物理和行星物理研究中，也經(jīng)常使用重力異常約束天體內(nèi)部結(jié)構(gòu)[9,10]. 文獻(xiàn)[10]給出了常密度下，表面地形與相應(yīng)重力異常傅里葉變換的關(guān)系，該關(guān)系經(jīng)嚴(yán)格證明具有一定的合理性，被廣泛應(yīng)用于重力反演[11,12]. 由于文獻(xiàn)[10]給出的模型只適用于平面模型，對于球體模型特別是小天體(如月球和火星)將產(chǎn)生較大的誤差，為此，文獻(xiàn)[13]基于文獻(xiàn)[10]，給出了球體模型下地形與重力異常兩者間球諧展開的廣義傅里葉關(guān)系，并被廣泛地應(yīng)用于地球、行星和月球物理研究中[14-16].

近年來，隨著月球和火星重力數(shù)據(jù)分辨率的提高[18,19]，利用它們的表面重力異常約束其內(nèi)部結(jié)構(gòu)的研究越來越多[14-16]. 文獻(xiàn)[13]給出的算法是基于常密度假設(shè)，而近年來的研究表明變密度模型更具普適性[15]. 文獻(xiàn)[15]在研究火星物理參數(shù)時(shí)，為了估算變密度模型的重力異常，提出了一種較復(fù)雜的數(shù)值積分方法. 該方法盡管推理嚴(yán)謹(jǐn)，計(jì)算精度高，但計(jì)算速度較慢，效率低下.基于廣義傅里葉變換的球諧算法，被廣泛應(yīng)用于行星和月球物理研究中[13-17]. 球諧算法的最大優(yōu)點(diǎn)是運(yùn)行速度明顯優(yōu)于普通數(shù)值積分方法[20]. 目前，尚無研究考慮變密度球體重力位球諧展開的算法，針對該問題，本文基于廣義傅里葉展開的球諧算法，推導(dǎo)變密度與重力位兩者間球諧系數(shù)的關(guān)系，以期為變密度重力正反演問題提供一定的參考.

1 變密度球體模型的重力位

變密度球體模型如圖1所示，O表示球心，r表示球體的最大半徑.球體內(nèi)部任意一點(diǎn)P(r′,θ′,φ′)的密度Δρ(θ′,φ′)，僅與經(jīng)度φ′和余緯度θ′有關(guān)，該點(diǎn)對球體外任意點(diǎn)Q(R,θ,φ)的重力位有貢獻(xiàn)，將球體內(nèi)所有點(diǎn)進(jìn)行積分，可得變密度球體對其外部點(diǎn)Q(R,θ,φ)的重力位為

U(R,θ,φ)=

(1)

其中，γ表示兩向量OP和OQ之間的夾角，d表示點(diǎn)P和點(diǎn)Q之間的距離，有如下關(guān)系：

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

(2)

(3)

根據(jù)勒讓德母函數(shù)的定義[21]，可得

(4)

聯(lián)合式(1)—式(4)，顧及密度在徑向上無變化，并在徑向進(jìn)行積分，可得變密度球體對其外部任意點(diǎn)Q(R,θ,φ)的重力位為

(5)

圖1 變密度球體模型示意圖

2 重力位的球諧展開

參考文獻(xiàn)[1-3,21]，勒讓德多項(xiàng)式Pn(cosγ)可以表示成規(guī)格化球諧函數(shù)的線性組合，即

(6)

(7)

如圖1所示，密度僅是余緯度和經(jīng)度的函數(shù)，與半徑無關(guān)，參考文獻(xiàn)[1-3]，可將此類函數(shù)表示成規(guī)格化球諧函數(shù)的線性組合，即

Δρ(θ′,φ′)=

(8)

(9)

(10)

定義立體角微元dσ=sinθ′dθ′dφ′，參考文獻(xiàn)[1-3]，規(guī)格化球諧函數(shù)具有正交性特征，即

(11)

(12)

其中，狄拉克函數(shù)δnl和δmk取值為

(13)

將式(6)—式(8)代入式(5)，并考慮到規(guī)格化球諧函數(shù)具有式(11)、式(12)的正交性特征，可得

U(R,θ,φ)=

(14)

式(14)表明，變密度異常體對其外部任意點(diǎn)的重力位，可以直接表示成變密度球諧展開系數(shù)的組合形式.變密度異常體對半徑為R的球面的重力位U(R,θ,φ)，可以直接表示成變密度球諧展開系數(shù)的組合形式.變密度異常體對半徑為R的球面的重力位，參考文獻(xiàn) [1-3, 10, 13]可將其表示成規(guī)格化球諧函數(shù)的線性組合，即

U(R,θ,φ)=

(15)

(16)

考慮實(shí)際重力位一般只能展開成有限項(xiàng)之和[1-3,10,13]，令其最大展開階數(shù)為lmax，可得變密度球體對其外部任意點(diǎn)Q(R,θ,φ)的重力位為

(17)

在實(shí)際應(yīng)用中，如果不使用式(17)的算法，也可以采用數(shù)值積分方法對式(1)進(jìn)行計(jì)算.考慮式(2)、式(3)，可將式(1)表示成如下形式：

(18)

其中，積分函數(shù)為

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

(19)

3 算法測試與應(yīng)用

為了驗(yàn)證本文推導(dǎo)結(jié)果的有效性，假定變密度Δρ(θ′,φ′)=sinθ′ cos2φ′，利用式(9)、式(10)算出密度的球諧展開系數(shù)，然后利用式(16)、式(17)計(jì)算出重力位的球諧展開系數(shù)及重力位.利用式(9)、式(10)計(jì)算密度球諧展開系數(shù)時(shí)，需要計(jì)算積分值，如果利用普通數(shù)值計(jì)算方法會花費(fèi)較多計(jì)算時(shí)間，但考慮到式(9)、式(10)的積分函數(shù)涉及球諧函數(shù)，可以利用快速傅里葉變換進(jìn)行抽樣[22]，因此，式(9)、式(10)在計(jì)算時(shí)效率仍然較高.在已知變密度球諧展開系數(shù)的情況下，利用式(17)，可以快速求出全球不同位置點(diǎn)的重力位.而使用普通數(shù)值方法的式(18)計(jì)算時(shí)，由于每個(gè)位置點(diǎn)均要進(jìn)行積分計(jì)算而耗費(fèi)大量時(shí)間，計(jì)算全球不同位置點(diǎn)的重力位將是一件冗長的事情.

為了方便比較，如表1所示，取引力常數(shù)G=1，r=1 m，R=2 m，經(jīng)度和緯度間隔Δθ=Δφ=0.5°，即每隔0.5°計(jì)算一個(gè)點(diǎn)的重力位.參考文獻(xiàn)[1-3, 11]，由于Δθ=Δφ=0.5°，可知式(17)中最大展開階次lmax=179.本文推導(dǎo)的球諧算法，采用單核估算，耗時(shí)主要由式(9)、式(10)和式(17)產(chǎn)生，經(jīng)多次試驗(yàn)，發(fā)現(xiàn)耗時(shí)1.7 s左右.采用普通數(shù)值方法估算式(18)、式(19)比較耗時(shí)，為此，采用MATLAB優(yōu)化的三維數(shù)值積分函數(shù)integral3，并使用與球諧算法相同的處理器進(jìn)行20核并行計(jì)算，經(jīng)多次試驗(yàn)，發(fā)現(xiàn)耗時(shí)在1 120 s左右.由此可知，本文推導(dǎo)的球諧算法在計(jì)算效率方面遠(yuǎn)勝于普通數(shù)值積分方法，對實(shí)際問題的大規(guī)模重力正反演具有一定的參考價(jià)值.

兩種算法產(chǎn)生的重力位如圖2所示，其中圖2(a)表示基于本文球諧算法的結(jié)果，圖2(b)表示基于式(18)、式(19)的數(shù)值積分結(jié)果，圖2(c)表示兩者差值的絕對值.圖中橫軸表示經(jīng)度的取值，縱軸表示緯度的取值(余緯度與緯度互為余角). 由圖2(a)、2(b)可以看出，兩種算法得到的重力位分布完全一致. 圖2(c)表明兩種算法重力位差值的絕對值，在赤道附近最小，而在兩極附近較大，但最大值不超過0.001 5，因此，本文推導(dǎo)的球諧算法具有一定的合理性.

圖2 不同算法的重力位及其偏差的絕對值

表1 參數(shù)取值及計(jì)算時(shí)間

應(yīng)用本文推導(dǎo)的式(17)，參考文獻(xiàn)[23]提供的月殼密度分布數(shù)據(jù)，圖3(a)和圖3(b)分別給出了月殼密度和相應(yīng)的重力位分布，投影中心坐標(biāo)(180°W, 0°N)，其中圓圈1表示月球南極艾德肯盆地(South Pole-Aitken)，圓圈2表示月球背面高地.在計(jì)算過程中，參考文獻(xiàn)[23]的最新成果，取引力常數(shù)G=6.673 981011m3·kg-1·s-2.參考文獻(xiàn)[24]，取待估重力位參考球面的半徑R=1 738 km，取月殼層厚度Tc=40 km，以及月殼層外半徑r2=1 737.15 km.為了計(jì)算月殼層的重力位，先計(jì)算外半徑r2球體的重力位，再計(jì)算內(nèi)半徑r1=r2-Tc=1 697.15 km球體的重力位，兩球體重力位之差即是月殼層的重力位.圖3(a)給出的圓圈1顯示出高密度物質(zhì)，圖3(b) 相應(yīng)地也顯示較大的重力位.圖3(a)顯示月球背面高地平均密度最小(如圓圈2所示)，圖3(b)顯示該區(qū)域的重力位也最小.除此之外，其他區(qū)域的密度分布與相應(yīng)的重力位分布規(guī)律總體一致，表明本文推導(dǎo)的球諧算法具有一定的參考價(jià)值.

圖3 月殼層密度及40 km殼層厚度的重力位分布

4 結(jié)論

因天體內(nèi)部密度分布不均，天體外部重力位具有各向異性的特性.普通數(shù)值積分方法估算變密度球體重力位時(shí)，計(jì)算時(shí)間較長，效率低下，為此，本文提出基于廣義傅里葉變換的球諧算法.結(jié)果表明，本文球諧算法估算的重力位分布，與普通數(shù)值積分方法所得一致，而計(jì)算時(shí)間明顯低于普通數(shù)值積分方法.將本文球諧算法應(yīng)用于月殼變密度問題，發(fā)現(xiàn)變密度月殼的重力位分布與密度分布總體一致，表明本文球諧算法具有一定的合理性和實(shí)用價(jià)值，可為變密度球體重力正反演問題提供有益的參考.