朱化強(qiáng)

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

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

1 直角坐標(biāo)平面模型的重力位

參考文獻(xiàn)[2]，假設(shè)直角坐標(biāo)系下平面模型的重力位為Δg，產(chǎn)生該重力位的地形高程為h(r)，對(duì)兩者進(jìn)行傅里葉展開(kāi)，可得重力位的傅里葉變換F[Δg]與n階ρ地形傅里葉變換F[hn(r)]的關(guān)系為

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

2 球面地形的重力位

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

(1)

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

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

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

UA(r0,θ0,φ0)=

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

(2)

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

UA(r0,θ0,φ0)=

(3)

(4)

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

(5)

其中

(6)

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

(7)

上式中狄拉克函數(shù)δnl和δmk滿足如下關(guān)系

(8)

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

(9)

(10)

(11)

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

(12)

(13)

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

3 討論與應(yīng)用

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

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

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

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

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

為了比較直角坐標(biāo)模型與本文球諧模型的差異，圖5給出了本文算法的重力位UA與直角坐標(biāo)平面模型重力位Δg的差值UA-Δg.火星南半球?yàn)楦叩兀卑肭驗(yàn)槠皆?，赤道附近聚積了較多火山.圖5表明兩種算法對(duì)兩極，以及赤道附近火山區(qū)域重力位的計(jì)算存在一定的偏差.文獻(xiàn)[18]的研究表明，火星扁率對(duì)其重力位的影響較大.圖5也顯示兩種算法的重力位在兩極存在一定的差異，表明火星表面曲率對(duì)重力位的影響不能忽略.因此，本文基于球坐標(biāo)系推導(dǎo)的球諧模型更適合于球狀天體表面地形重力位的估算.

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

4 結(jié)論

由于月球和火星此類(lèi)小天體曲率特征的影響，目前地學(xué)界常用的直角坐標(biāo)模型不適合于這類(lèi)小天體的重力位估算.特別是近年來(lái)火星和月球探測(cè)數(shù)據(jù)越來(lái)越豐富，有必要構(gòu)建適用于這類(lèi)小天體的球坐標(biāo)重力正反演模型.本文基于牛頓引力位理論，應(yīng)用球諧變換對(duì)重力位進(jìn)行展開(kāi)，通過(guò)特定的數(shù)值推理，發(fā)現(xiàn)表面地形引起的重力位可以用地形的高階展開(kāi)系數(shù)直接計(jì)算，進(jìn)而簡(jiǎn)化了傳統(tǒng)數(shù)值積分復(fù)雜的求積過(guò)程.同時(shí)，研究發(fā)現(xiàn)高階地形產(chǎn)生的影響不超過(guò)7階，在實(shí)際應(yīng)用過(guò)程中取5階以內(nèi)的高階地形就能滿足精度要求.