以擾動重力為邊值條件確定外部擾動重力場

田家磊，李新星，劉曉剛，馮進(jìn)凱，范 雕

（1. 信息工程大學(xué) 地理空間信息學(xué)院，鄭州 450001；2. 63850部隊，吉林 白城 137000；3. 武漢大學(xué) 測繪學(xué)院，武漢 430079；4. 西安測繪研究所，西安 710054；5. 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室，西安 710054；6. 大地測量與地球動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，武漢 430077）

地球重力場是地球的基本物理場，同時也是大地測量學(xué)、海洋學(xué)、地震學(xué)、國防科學(xué)、空間科學(xué)、地球物理學(xué)等學(xué)科重要的研究對象[1-2]。地球重力場及其時變反映了地球內(nèi)部質(zhì)量分布情況和運(yùn)動狀態(tài)，同時決定了大地水準(zhǔn)面形狀與隨時間的變化[3-4]。物理大地測量學(xué)的兩大基本任務(wù)就是確定地球的形狀和地球外部擾動重力場[5]。

導(dǎo)彈、飛機(jī)、火箭、航天飛機(jī)、衛(wèi)星、宇宙探測器等飛行器在飛行過程中都要受到地球外部重力場的影響，所以它們的發(fā)射、制導(dǎo)與控制等過程需要精確的重力場信息來提供這些運(yùn)動物體力學(xué)的重力場約束[6]。隨著慣性器件的精度越來越高，擾動重力場對慣性導(dǎo)航系統(tǒng)（INS）的影響越發(fā)地突出，想要進(jìn)一步的提高導(dǎo)航精度就需要在導(dǎo)航的過程中對擾動重力進(jìn)行補(bǔ)償[7]，由此也加大了對高精度、高分辨率的擾動重力場模型的需求[8]。重力輔助導(dǎo)航也是近些年研究的熱點(diǎn)之一[9]。

目前地球外部擾動重力場的模型構(gòu)建和等效逼近理論主要包括兩類：模型逼近和算法逼近。模型逼近主要有Sttokes理論和Molodensky理論等。算法逼近主要有貝亞哈馬解、虛擬點(diǎn)質(zhì)量解[10]、虛擬單層密度解、最小二乘配置解及球諧函數(shù)解等。雖然上述方法可以用來確定地球外部擾動重力場，但在應(yīng)用上均受到一定的限制。

傳統(tǒng)的重力測量及數(shù)據(jù)處理得到的是重力異常，經(jīng)典的物理大地測量理論也都是以重力異常為邊值條件。但是隨著 GNSS、航空重力等技術(shù)的發(fā)展，高精度、高分辨率擾動重力數(shù)據(jù)的獲取變得經(jīng)濟(jì)、便捷與高效。目前利用擾動重力確定區(qū)域大地水準(zhǔn)面的研究比較熱門，然而對于利用擾動重力構(gòu)建擾動重力場的研究比較少。下面結(jié)合實(shí)際應(yīng)用的需要，討論利用擾動重力確定外部擾動重力三分量的方法。

1 格林積分法

球近似下地球各元素之間的幾何關(guān)系如圖1所示。

由位理論可得，擾動位T是調(diào)和函數(shù)，根據(jù)調(diào)和函數(shù)的性質(zhì)，擾動位可由格林第三式表示為：

圖1 球近似下各元素之間的幾何關(guān)系Fig.1 Geometric relations among the elements under the condition of the ball approximation

式中，∑表示地面，d∑為地面上的面元，l是計算點(diǎn)p與面元的距離，n是相對于調(diào)和空間的邊界面外法線方向。ψ為計算點(diǎn)p與地面∑上的面元d∑的球心角，?T?n和T分別表示地面的擾動重力和擾動位。式(1)就是應(yīng)用格林公式以實(shí)際地球表面為邊界面，以擾動位和擾動重力矢量為邊值條件，計算外部擾動位的公式。它需要具有二者的觀測值才能進(jìn)行解算，地形復(fù)雜時法線方向的計算是比較困難，這兩個因素限制了格林公式的實(shí)際應(yīng)用。盡管如此，格林公式提供了一種以地球自然表面為邊界面（自然表面上的擾動位與擾動重力矢量為邊值），而不需對任何邊值進(jìn)行延拓或歸算，確定地球外部擾動重力場的解析形式[11]。

實(shí)際應(yīng)用中的地形與重力數(shù)據(jù)都是格網(wǎng)化的形式，通過分析與簡化可以得到基于格網(wǎng)數(shù)字地形面的格林公式：

式中，δgi為積分點(diǎn)的擾動重力，lpi表示計算點(diǎn)與地面點(diǎn)之間的距離，Ti為地面點(diǎn)的擾動位，rrp和ri分別表示計算點(diǎn)和地面點(diǎn)的地心向徑，Δσi=cosφi? Δφi?Δλi為格網(wǎng)的單位球面積，λ、φ分別為經(jīng)度和緯度。

式(2)表明利用格林公式求解地球外部擾動位同時需要有地面的擾動位以及擾動重力的垂直分量這兩個邊值條件。分析式(2)的形式可以看出它是由擾動重力垂直分量的單層位和擾動位的泊松核的位兩部分組成。在球近似的情況下它們分別是這兩種邊值問題的解。對式(2)做一定的改進(jìn)可得：

式中，p0是計算點(diǎn)沿著向徑至地面上的點(diǎn)，Tp0表示p0點(diǎn)的擾動位，γ是正常重力，ζ是高程異常。

式(3)采用了相對于p0點(diǎn)的高程異常差作為輸入數(shù)據(jù)，這樣做的好處有：它消除了高程異常的系統(tǒng)性誤差和長波長的誤差，p0點(diǎn)附近高程異常差接近于零，降低了積分的奇異性，精度更高的擾動重力可以起到更大的作用。

地球外部擾動重力矢量通常采用沿球心向徑、緯度和經(jīng)度三個方向的分量表示。它們與擾動位有以下關(guān)系：

其中，Api是計算點(diǎn)p與流動面元之間的方位角。將式(3)帶入式(4)得到：

利用式(5)，即可由地面上的擾動重力和相對0p的高程異常差計算出地面或空間一點(diǎn)的擾動重力矢量的三個分量。

2 豪汀積分法

第二邊值問題是指根據(jù)邊界面上調(diào)和函數(shù)的法線導(dǎo)數(shù)推求界面內(nèi)或者界面外的調(diào)和函數(shù)的問題，亦稱Neuman邊值問題，球外部第二邊值問題的解稱為Hotine積分。

式中，δg為擾動重力，H(rp,ψpi)為Hotine核，具體形式為：離散化可得：

由式(4)可得：

根據(jù)圖1所示的球近似下各元素之間的幾何關(guān)系可得：

將式(9)帶入式(8)可得由豪汀公式計算擾動重力三分量的公式：

3 點(diǎn)質(zhì)量模型法

瑞典的地球物理學(xué)家貝亞哈馬于 1964年提出的一種理論，其核心思想就是把地球?qū)ν獠奎c(diǎn)的擾動位當(dāng)成完全處于地球內(nèi)部一個虛擬球（貝亞哈馬球，見圖2）對它的引力位。基于貝亞哈馬理論，后來有學(xué)者提出用位于貝亞哈馬球上面的虛擬單層位模型等效逼近地球外部重力場。單層位模型的積分核函數(shù)是空間距離的倒數(shù)，便于計算。若將虛擬單層質(zhì)量離散化就可以得到虛擬點(diǎn)質(zhì)量模型。

現(xiàn)假定點(diǎn)質(zhì)量模型由分布于地球內(nèi)部的k個擾動質(zhì)量為Mj(j=1,2,…,k)的離散擾動質(zhì)點(diǎn)組成。根據(jù)牛頓萬有引力定律可以算出地球表面外任意一點(diǎn)P的擾動位，即

式中，G是引力常數(shù)。擾動重力場可以由埋藏在地球

圖2 貝亞哈馬球Fig.2 Bjerhammar sphere

式中，BR為貝亞哈馬球的半徑，δig為第i個擾動重力觀測量。擾動質(zhì)點(diǎn)的質(zhì)量可以由地面上已知的擾動重力來推求。由式(12)可以推導(dǎo)出點(diǎn)質(zhì)量模型計算擾動重力三分量的式：

4 擾動重力的向上延拓

擾動重力和擾動重力的徑向分量數(shù)值相等方向相反，因此擾動重力的向上延拓是求擾動重力徑向分量的實(shí)用方法。向上延拓屬于第一邊值問題，泊松積分是第一邊值問題的解。對于調(diào)和函數(shù)的向上延拓的實(shí)用公式為：

rδg是調(diào)和函數(shù)，根據(jù)調(diào)和函數(shù)的性質(zhì)得到關(guān)于擾動重力的向上延拓式：

式中，P是計算點(diǎn)，P0是P沿向徑方向在球面上的投影點(diǎn)，Q為流動點(diǎn)。

5 實(shí)驗(yàn)與分析

下面采用2160階EGM2008地球重力場模型模擬計算北緯 23.5°～28.5°，東經(jīng) 111.5°～116.5°區(qū)域，大地高為 0、分辨率為2′×2′、大小為150×150的格網(wǎng)擾動重力和高程異常網(wǎng)格數(shù)據(jù)?？罩袛?shù)據(jù)采用2160階EGM2008模型模擬計算北緯 24°～28°，東經(jīng) 112°～116°區(qū)域，高度為3000 m和1500 m、分辨率為2′×2′、大小為120×120的格網(wǎng)擾動重力三分量，將其作為“真值”，檢驗(yàn)各種方法計算外部擾動重力的精度。

采用同樣的地面重力數(shù)據(jù)分別利用三種方法計算同樣范圍的擾動重力矢量然后和“真值”做差，結(jié)果統(tǒng)計如表1～6所示。

可以看出，隨著計算高度的增加，三種方法的結(jié)果都變好，1500 m高度計算的結(jié)果點(diǎn)質(zhì)量模型的精度是最好的，豪汀積分法次之但是其計算的徑向分量結(jié)果不好，格林積分法的結(jié)果稍差。3000 m高度的水平分量計算結(jié)果豪汀積分法最好，格林積分法次之，點(diǎn)質(zhì)量模型稍差，但是結(jié)果都比較好。相比另外兩種方法點(diǎn)質(zhì)量模型在兩個高度的徑向分量的計算都取得了較好的結(jié)果。

表1 格林積分法計算的3000 m高度擾動重力三分量差值Tab.1 Differences between true values and three components of disturbance gravity calculated by Green integral method at height of 3000 m

表2 豪汀積分法計算的3000 m高度擾動重力三分量差值Tab.2 Differences between true values and three components of disturbance gravity calculated by Hotine integral method at height of 3000 m

表3 點(diǎn)質(zhì)量模型法計算的3000 m高度擾動重力三分量差值Tab.3 Differences between true values and three components of disturbance gravity calculated by Point mass method at height of 3000 m

表4 格林積分法計算的1500 m高度擾動重力三分量差值Tab.4 Differences between true values and three components of disturbance gravity calculated by Green integral method at height of 1500 m

表5 豪汀積分法計算的1500 m高度擾動重力三分量差值Tab.5 Differences between true values and three components of disturbance gravity calculated by Hotine integral method at height of 1500 m

表6 點(diǎn)質(zhì)量模型法計算的1500 m高度擾動重力三分量差值Tab.6 Differences between true values and three components of disturbance gravity calculated by Point mass method at height of 1500 m

為了研究擾動重力誤差對計算擾動重力三分量的影響，并且實(shí)測的擾動重力數(shù)據(jù)由于觀測條件的限制必然會有一定的誤差，在原始數(shù)據(jù)加入3 mGal的高斯白噪聲，然后計算擾動三分量差值，結(jié)果統(tǒng)計如表7～10所示。

從表 7～10可以看出，原始地面加入誤差后，三種方法計算擾動重力三分量仍然取得了較好的效果，說明三種方法都有一定的抗差能力。同時利用向上延拓算法也可以很好地完成擾動重力徑向分量的計算。

表7 格林積分法加入3 mGal誤差計算的1500 m擾動重力三分量差值Tab.7 Differences between true values and three components of disturbance gravity calculated by Green integral method with 3mGal error added at the height of 1500 m

表8 豪汀積分法加入3 mGal誤差計算的1500 m擾動重力三分量差值Tab.8 Differences between true values and three components of disturbance gravity calculated by Hotine integral method with 3mGal error added at the height of 1500 m

表9 點(diǎn)質(zhì)量模型法加入3 mGal誤差計算的1500 m擾動重力三分量差值Tab.9 Differences between true values and three components of disturbance gravity calculated by Point mass method with 3 mGal error added at the height of 1500 m

表10 向上延拓法計算擾動重力差值Tab.10 Differences between disturbance gravity calculated by upward continuation method

6 結(jié) 論

精確地確定地球的外部擾動重力場對提高慣性導(dǎo)航的精度愈發(fā)重要。隨著GNSS技術(shù)的發(fā)展，重力測量以及數(shù)據(jù)處理中越來越廣泛地采用大地高系統(tǒng)，因而獲得了越來越多的擾動重力數(shù)據(jù)。本文研究了幾種以擾動重力為邊值條件確定地球外部擾動重力場的方法。經(jīng)過實(shí)驗(yàn)說明三種方法都可以有效地計算擾動重力三分量。在低空領(lǐng)域點(diǎn)質(zhì)量模型法與豪汀積分法（計算水平分量）加擾動重力的向上延拓（計算徑向分量）的組合是比較好的方法。雖然在低空領(lǐng)域格林積分方法的結(jié)果稍差，但是其計算過程中不需要數(shù)據(jù)的歸算，并且其核函數(shù)簡單便于計算，在滿足其條件的情況下也可以有效地確定地球的外部擾動重力場。