一種球面地形重力位的計算方法

朱化強

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

在月球和類地行星表面地形重力異常的估算中，早期研究采用低階次近似，常常導致較大的誤差[1].在直角坐標系中采用高階近似，可以得到較好的精度[2]，被廣泛地應用于重力異常的正反演問題[3-5].這種高階近似算法，主要依據快速傅里葉變換，通過地形頻譜即可以迅速求解相應的重力異常，再經過傅里葉逆變換即可求得研究區域表面地形引起的重力異常.這種方法不但能保持快速計算，還便于重力異常信息的濾波處理，具有較大的實用價值，在地球物理研究中被廣泛地應用[3-5].隨著月球探測數據的增加，發現諸如月球這類小天體，其曲率產生誤差將不能被忽略[6-7].文獻[2]的模型對月球和火星這類小天體重力異常的正反演，會產生不可剔除的模型誤差，導致結果解釋不合理等現象[6-7].近年來，隨著火星和月球探測活動的開展，相關重力數據越來越豐富[8-9]，有必要構建適用于這類小天體的球坐標重力正反演模型.為此，本文提出一種球面地形重力位的估算方法，以期為小天體重力正反演提供一定的參考.

1 直角坐標平面模型的重力位

參考文獻[2]，假設直角坐標系下平面模型的重力位為Δg，產生該重力位的地形高程為h(r)，對兩者進行傅里葉展開，可得重力位的傅里葉變換F[Δg]與n階ρ地形傅里葉變換F[hn(r)]的關系為

其中，r表示坐標原點至地形參考點的位置矢量，n表示傅里葉展開的階數，z0表示參考平面的位置高階，k表示波矢量，G表示萬有引力常量，ρ表示地形密度.上式對于大天體而言，界面地形的重力位因界面弧度的影響較小，但對小天體而言，地形界面因弧度對重力位的貢獻必須考慮進去.

2 球面地形的重力位

如圖1所示，假定天體的球心為O，平均參考半徑為R，相對該參考球面距離為h(θ,φ)的地形點B(半徑為r).對于半徑為r0的參考球面(r0>R)的點A，點B對點A的重力位為(如圖1虛線所示)

(1)

其中，ρ(r,θ,φ)表示天體表面地形的密度函數，G為萬有引力常量，dV表示天體表面地形的體積微元；r0和r分別表示點A和點B的徑向矢量，兩矢量的夾角為γ，兩矢量差的模為|r-r0|.

圖1 界面地形重力位估算示意圖

假定表面地形的密度為常數，即ρ(r,θ,φ)=ρ0，并考慮到勒讓德多項式母函數的關系[10-11]，以及積分與求和進行互換，對式(1)在徑向進行積分得

UA(r0,θ0,φ0)=

Rn+3]Pn(cosγ)sinθdθdφ

(2)

上式中，Pn(cosγ)表示n階勒讓德多項式.利用二項式公式，可以將式(2)進一步化簡成

UA(r0,θ0,φ0)=

(3)

(4)

為了進一步化簡式(3)，參考文獻[10]，將n階勒讓德多項式Pn(cosγ)表示成n階m次正則化球諧函數Yinm的關系，即

(5)

其中

(6)

上式中，Pnm(cosθ)表示正則化的n階m次連帶勒讓德函數，cosmφ和sinmφ表示余弦和正弦函數，θ和φ分別表示余緯度和經度.參考文獻[12]，正則化的球諧函數Yinm滿足如下正交性原則

(7)

上式中狄拉克函數δnl和δmk滿足如下關系

(8)

如圖1所示，參考文獻[3-5，10]可將地形高程h(θ,φ)展開l階p次成球諧函數的形式，假定地形高程球諧展開系數為hilp，即

(9)

(10)

(11)

將式(11)代入式(3)，顧及式(5)—(8)，可得

(12)

(13)

其中，M表示天體在參考半徑R內，平均密度為ρ0的總質量，而且必須滿足r0>r.

3 討論與應用

有關地形引起的重力位，其位系數式(13)涉及求和運算，而且隨著階次n的增加，計算量不斷增大，但這實際上沒必要考慮太多高階項.式(13)來源于式(4)，其中h/R表示表面地形高程與天體平均半徑的比值，此項的值一般較小，特別是對于高階項.以火星表面最高的奧林匹斯山(Olympus)為例，其高程接近22 km[12-14]，而火星的平均半徑為3389.5 km[15]，h/R的高階項趨近于零，一般計算時不超過7階.圖2表示求和項式(4)中每一項的值隨k和n的分布，由圖可知，當k超過5時，式(4)的高階項幾乎沒有變化.圖3表示式(4)中所有項之和，由圖可知，對于同一個球諧階數n，k取k=3和k=7時并沒有太大的區別，因此，實際應用中，式(4)的求和項數最多取前7項即可滿足精度要求.

圖2 重力位求和式(4)中每一項隨k和n變化圖

圖3 重力位求和式(4)隨k和n變化的總和分布

為了應用本文算法，參考文獻[16，17]，取火星行星殼密度ρc=2900 kg·m-3，取火星重力位參考球面半徑r0=3396 km[17]，地形數據取自文獻[16].在計算過程中取式(13)最大求和項數為5項，火星表面地形產生的重力位分布如圖4所示，其投影中心坐標為(-90°W，0°E)，左邊是火星表面著名的奧林匹斯山脈(如圖4白色區域所示).圖4表明火星北半球的重力位明顯低于其南半徑，這主要是由于北半球因低洼的盆地產生負的重力位，而南半球因高原產生正重力位所致，這與火星表面如文獻[16]給出的地形特征一致，表明本文構建的算法具有一定的合理性.

圖4 火星表面地形取火星殼密度2900 kg·m-3時的重力位分布

為了比較直角坐標模型與本文球諧模型的差異，圖5給出了本文算法的重力位UA與直角坐標平面模型重力位Δg的差值UA-Δg.火星南半球為高地，北半球為平原，赤道附近聚積了較多火山.圖5表明兩種算法對兩極，以及赤道附近火山區域重力位的計算存在一定的偏差.文獻[18]的研究表明，火星扁率對其重力位的影響較大.圖5也顯示兩種算法的重力位在兩極存在一定的差異，表明火星表面曲率對重力位的影響不能忽略.因此，本文基于球坐標系推導的球諧模型更適合于球狀天體表面地形重力位的估算.

圖5 直角坐標模型與本文球諧模型求解的重力位偏差圖(火星殼密度取2900 kg·m-3)

4 結論

由于月球和火星此類小天體曲率特征的影響，目前地學界常用的直角坐標模型不適合于這類小天體的重力位估算.特別是近年來火星和月球探測數據越來越豐富，有必要構建適用于這類小天體的球坐標重力正反演模型.本文基于牛頓引力位理論，應用球諧變換對重力位進行展開，通過特定的數值推理，發現表面地形引起的重力位可以用地形的高階展開系數直接計算，進而簡化了傳統數值積分復雜的求積過程.同時，研究發現高階地形產生的影響不超過7階，在實際應用過程中取5階以內的高階地形就能滿足精度要求.