應(yīng)用格林積分直接以地面邊值確定外部擾動(dòng)重力場(chǎng)

田家磊，吳曉平，李?yuàn)檴?/p>

信息工程大學(xué)地理空間信息學(xué)院，河南 鄭州450001

1 擾動(dòng)重力位的地面邊值問題

確定地球外部重力場(chǎng)和大地水準(zhǔn)面是大地測(cè)量學(xué)的主要任務(wù)之一［1－5］。確定地球外部重力場(chǎng)和大地水準(zhǔn)面的斯托克斯理論需要將地面觀測(cè)的重力異常歸算至大地水準(zhǔn)面［6］，再采用調(diào)和函數(shù)球面邊值的解式（如Stokes公式）求得大地水準(zhǔn)面及其外部的擾動(dòng)重力位［7－8］。歸算將涉及對(duì)大地水準(zhǔn)面至地面的質(zhì)量遷移，對(duì)大地水準(zhǔn)面產(chǎn)生間接影響，而且由于歸算對(duì)質(zhì)量進(jìn)行了調(diào)整，改變了外部擾動(dòng)重力場(chǎng)，因此歸算到大地水準(zhǔn)面上的重力異常用以確定外部擾動(dòng)重力位會(huì)導(dǎo)致結(jié)果的歪曲［9－10］。

直接以地面重力異常為邊值的Molodensky問題從理論上避免了歸算的困難，成為近代外部重力場(chǎng)研究的理論基石［11］。然而，由于地球表面的復(fù)雜性，給這一問題的求解帶來極大難度。Molodensky基于基本積分方程的小參數(shù)解法得到地面擾動(dòng)位的級(jí)數(shù)解式［12－13］。文獻(xiàn)［10］提出將地面重力異常解析地延拓到一點(diǎn)的水準(zhǔn)面上，再采用球面的Stokes積分得到地面擾動(dòng)位，其結(jié)果同樣是級(jí)數(shù)的形式。文獻(xiàn)［15］也研究得到類似的級(jí)數(shù)解。文獻(xiàn)［16］則提出將地面重力異常調(diào)和地延拓到一個(gè)內(nèi)部球面上，再由球面邊值問題解逼近外部擾動(dòng)位，其調(diào)和延拓需要求解Poisson積分方程。盡管這些理論解的途徑有所不同，但在一定前提下它們是等價(jià)的。雖然經(jīng)過線性化和地球橢球作球近似后的所謂簡(jiǎn)單Molodensky問題的研究已得到幾近完美的理論結(jié)果，但它們的實(shí)現(xiàn)仍具有相當(dāng)大的困難［17］。由于受到數(shù)據(jù)和高階項(xiàng)計(jì)算穩(wěn)定性的限制，目前在實(shí)際上通常只能考慮到一階項(xiàng)［10］。

對(duì)于確定地球外部擾動(dòng)重力場(chǎng)問題，上述解在應(yīng)用上受到一定的限制。像Molodenky解通常應(yīng)用于地面，Moritz的解析延拓解和Bjeharmmar解雖然可以拓展到外部空中，但邊值的延拓仍是一個(gè)較復(fù)雜的過程［10］。本文側(cè)重于應(yīng)用的需要，討論直接由地面邊值確定外部擾動(dòng)位的方法。

2 地面邊值問題的格林公式表示

確定地球外部擾動(dòng)重力位T歸結(jié)為下面的邊值問題［18］

由位理論知，T作為調(diào)和函數(shù)可以由格林第三恒等式表示為［19］

式中，l是計(jì)算點(diǎn)p至地面Σ上面元dΣ的空間距離；n是相對(duì)于調(diào)和空間的邊界面外法線方向。取局部北東天坐標(biāo)系，求法線方向?qū)?shù)得［20］根據(jù)位理論為擾動(dòng)重力矢量為法線的方向余弦，因此可知即為擾動(dòng)重力在法線方向的分量，如圖1。

對(duì)于所謂的簡(jiǎn)單Molodensky問題，即將地球橢球近似為球面時(shí)，上述各元素的幾何關(guān)系見圖2。由距離公式

式中，l為P點(diǎn)與dΣ單元處的距離；rp為P點(diǎn)的球心距離；r為dΣ單元處的球心距離；ψ為P點(diǎn)到dΣ單元處的極距，求導(dǎo)可得cos（l，n），因此可得

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

應(yīng)用格林公式可以在實(shí)際地球表面上計(jì)算外部擾動(dòng)位。其條件是需要同時(shí)具有地面上的擾動(dòng)位和擾動(dòng)重力矢量的觀測(cè)值。這在實(shí)際應(yīng)用中是有困難的。一方面，所需的邊值條件很難滿足。另一方面地球表面非常復(fù)雜，這就使得在地表起伏較大地區(qū)該式中的法線方向變化劇烈，其計(jì)算相當(dāng)困難。盡管如此，格林公式提供了不需作任何邊值的歸算或延拓而以地球自然表面上的邊值條件確定外部擾動(dòng)重力場(chǎng)的唯一可能的解析形式。

3 基于格網(wǎng)數(shù)字地形面的擾動(dòng)重力位邊值問題解式

對(duì)于四周的側(cè)面積

在實(shí)際應(yīng)用中，地面通常是采用格網(wǎng)化的數(shù)字形式表示的。在這些格網(wǎng)單元中，其高度取作一個(gè)常數(shù)，因此這些格網(wǎng)顯示的只是一個(gè)個(gè)平面單元。而在重力場(chǎng)實(shí)際數(shù)據(jù)處理過程中，由于觀測(cè)數(shù)據(jù)是離散的，不可能達(dá)到地面邊值的連續(xù)分布，通常將離散的觀測(cè)數(shù)據(jù)（如重力異常）處理成不同規(guī)格的格網(wǎng)平均值，它們可視為相應(yīng)于平面格網(wǎng)地形面上的值。另外，在重力計(jì)算時(shí)，也是根據(jù)精度要求選取最優(yōu)的積分半徑和格網(wǎng)大小，而非采用越小越好的格網(wǎng)選擇策略。所以在實(shí)際上面臨的是由不同平面格網(wǎng)單元組合的地球表面的邊值問題。這種數(shù)字格網(wǎng)地形面是與實(shí)際地球表面最接近且為應(yīng)用上必然的形式。在此情況下，格林公式形式的簡(jiǎn)化有助于其應(yīng)用的實(shí)現(xiàn)。

設(shè)地形表面由經(jīng)緯度方向微小格網(wǎng)構(gòu)成，每一個(gè)格網(wǎng)為高程是常數(shù)的球面體單元（如圖3所示），該格網(wǎng)中的邊值元素為一定值，得到離散的格林積分為

式中，δgx、δgy是該格網(wǎng)點(diǎn)擾動(dòng)重力的兩個(gè)水平分量；A是該格網(wǎng)點(diǎn)至計(jì)算點(diǎn)P的方位角。

圖3 格網(wǎng)數(shù)字地形面示意圖

在每一格網(wǎng)單元中，對(duì)于朝上的格網(wǎng)頂面積，被積元素為

注意到，法線方向n都是指向網(wǎng)格內(nèi)的（如圖3所示），當(dāng)該格網(wǎng)點(diǎn)處于地形高處的情況下，相對(duì)的兩個(gè)側(cè)面積擾動(dòng)重力分量取相反的符號(hào)，其積分是抵消的，對(duì)于在地形持續(xù)升高或降低的地區(qū)，格網(wǎng)通常有1—2個(gè)單向的側(cè)面積需要計(jì)算。考慮到在重力計(jì)算情況下，格網(wǎng)間的起伏（平均坡度）通常小于5%，即使重力水平分量與垂直分量數(shù)值大小相當(dāng)（一般后者要大于前者），其積分貢獻(xiàn)也較小，故可忽略側(cè)面積的擾動(dòng)重力水平分量，得到

根據(jù)重力異常與擾動(dòng)重力的關(guān)系，式（12）也可表示成重力異常的形式

由式（4）的關(guān)系，式（11）還可表示成

在式（14）中，需要有地面擾動(dòng)位和重力異常（或擾動(dòng)重力垂直分量）兩個(gè)邊值條件。它是由一個(gè)重力異常（或擾動(dòng)重力垂直分量）的單層位和擾動(dòng)位的Poisson核的位兩部分組成?？梢宰C明，在球面情況下，它們分別是這兩種邊值問題的解，即

可見，格林積分在球面情況下即退變成兩個(gè)單邊值問題解之和。

目前的邊值問題理論多是研究以重力異常為單一邊值條件，地面擾動(dòng)位通常是作為要取得的結(jié)果。但對(duì)于確定外部擾動(dòng)重力場(chǎng)而言，為了能應(yīng)用格林公式直接由地面的邊值求取外部擾動(dòng)重力場(chǎng)，需要有地面擾動(dòng)位作為邊值條件，為此要用到已有的地面擾動(dòng)位的結(jié)果。

因?yàn)閿_動(dòng)位與高程異常有以下關(guān)系

式中，γ是正常重力；ζ是高程異常。故可將式（14）中的擾動(dòng)位替代為高程異常

式（17）可由地面重力異常和高程異常直接確定外部擾動(dòng)位，而不需要經(jīng)過延拓或歸算。注意到，根據(jù)格林公式的性質(zhì)，該公式適用于地面外部的空間，如果計(jì)算點(diǎn)位于地面Σ上時(shí)，其公式應(yīng)表示為

格林公式在從空間到層面上是不連續(xù)的，另外，應(yīng)用式（16）求取外部擾動(dòng)位，要用到地面高程異常。目前我國(guó)地面高程異常亦取得較好的精度（東部地區(qū)優(yōu)于10 cm，西部地區(qū)優(yōu)于30 cm），可以作為格林公式中的已知地面元素。但對(duì)式（17）分析可見，含高程異常項(xiàng)的是Poisson積分核函數(shù)，它對(duì)于計(jì)算點(diǎn)附近的高程異常的精度要求較高。假如平均重力異常的精度為3 mGal（1 mGal＝10－5m／s2），對(duì)于1 km 高度的計(jì)算點(diǎn)，則要求積分范圍10 km以內(nèi)的高程異常要優(yōu)于30 cm的精度才能與平均異常精度匹配。當(dāng)高度增加或減小時(shí)，這個(gè)范圍要求作相應(yīng)擴(kuò)大或減小。為了降低對(duì)高程異常的精度要求，需要對(duì)格林積分式作以下改進(jìn)。

由位理論知，當(dāng)p點(diǎn)位于地面外部時(shí)有

則式（2）可表示成

式中，p0點(diǎn)是計(jì)算點(diǎn)沿向徑至地面上的點(diǎn)?？梢宰C明，該式在空間和地面上都是適用的，據(jù)此，式（14）變?yōu)?/p>

上述改進(jìn)公式的作用在于：因?yàn)椴捎昧讼鄬?duì)于p0點(diǎn)的高程異常差作為輸入數(shù)據(jù)，它消除了高程異常的系統(tǒng)性誤差（包含高程基準(zhǔn)的系統(tǒng)差）和長(zhǎng)波長(zhǎng)的誤差。特別對(duì)于距離p0較近的點(diǎn)，高程異常差可以消除許多公共性誤差，這些誤差在高程異常誤差中占很大比重，從而可以保證計(jì)算對(duì)高程異常數(shù)據(jù)的精度要求。另外，p0點(diǎn)附近高程異常在原格林積分中貢獻(xiàn)占優(yōu)，特別是計(jì)算點(diǎn)接近地面時(shí)積分強(qiáng)奇異。而在改進(jìn)的格林積分中，由于p0點(diǎn)附近高程異常差接近于零，降低了積分的奇異性，以精度更高的重力異常起主要作用。另外注意到，經(jīng)過改化的格林公式不僅適用于外部空間，同時(shí)也適用于地面，即從地面到外部是連續(xù)的。而原始格林公式在外部至地面的極限是不同的（相差1／2，這是由于單層位的導(dǎo)數(shù)在層面上的極限不連續(xù)）。

4 由地面邊值確定擾動(dòng)重力場(chǎng)元

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

利用式（29），即可由地面上的重力異常和相對(duì)p0的高程異常差計(jì)算地面或空間一點(diǎn)的擾動(dòng)重力矢量的3個(gè)分量。

5 試驗(yàn)計(jì)算分析

下面進(jìn)行兩個(gè)模擬計(jì)算。試驗(yàn)1驗(yàn)證了邊界面為球面時(shí)改進(jìn)的格林公式和Stokes公式的等效性，試驗(yàn)2驗(yàn)證了在實(shí)際地形面為邊界面的情況下格林積分公式優(yōu)于Stokes公式。

5．1 改進(jìn)的格林公式與Stokes公式的等效性驗(yàn)證

采用2160階EGM2008地球重力場(chǎng)模型模擬區(qū)域21°N—31°N，109°N—119°N，大地高為0的分辨率為2′×2′的重力異常和高程異常網(wǎng)格數(shù)據(jù)?？罩袛?shù)據(jù)用EGM2008模型模擬區(qū)域23．5°N—28．5°N，111．5°E—116．5°E，高度為3000 m，分辨率為2′×2′的150×150格網(wǎng)擾動(dòng)重力三分量，將其作為“真值”檢驗(yàn)格林積分計(jì)算外部擾動(dòng)重力的質(zhì)量（見表1）。利用格林公式選取積分半徑為2．5°采用移去恢復(fù)法計(jì)算外部擾動(dòng)重力（見表2）。

可以看出利用改進(jìn)的格林公式大大地提高了徑向分量的精度。利用Stokes公式計(jì)算相同區(qū)域相同高度的擾動(dòng)重力三分量，結(jié)果如表3所示。

采用同樣的地面重力數(shù)據(jù)分別利用Stokes公式和改進(jìn)的格林公式計(jì)算同樣范圍的1500 m高度的擾動(dòng)重力三分量得到的結(jié)果如表4、表5所示。

表1 未改進(jìn)的格林公式計(jì)算的高度3000 m擾動(dòng)重力差值Tab．1 Calculation of disturbing gravity difference at the altitude of 3000 m by using unimproved Green formula mGal

表2 改進(jìn)的格林公式計(jì)算的高度3000 m擾動(dòng)重力差值Tab．2 Calculation of disturbing gravity difference at the altitude of 3000 m by using improved Green formula mGal

表3 Stokes公式計(jì)算的高度3000 m擾動(dòng)重力差值Tab．3 Calculation of disturbing gravity difference at the altitude of 3000 m by using Stokes formula mGal

表4 改進(jìn)的格林公式計(jì)算的高度1500 m擾動(dòng)重力差值Tab．4 Calculation of disturbing gravity difference at the altitude of 1500 m by using improved Green formula mGal

表5 Stokes公式計(jì)算的高度1500 m擾動(dòng)重力差值Tab．5 Calculation of disturbing gravity difference at the altitude of 1500 m by using Stokes formula mGal

5．2 實(shí)際地形面為邊界面下格林積分公式優(yōu)越性驗(yàn)證

為了檢驗(yàn)格林公式在實(shí)際地形面上的計(jì)算效果，選取了采用2160階EGM2008地球重力場(chǎng)模型模擬區(qū)域32°N—36°N，108°E—114°E，分辨率為2′×2′的重力異常和高程異常網(wǎng)格數(shù)據(jù)。高程信息采用實(shí)際的高程，分辨率為2′×2′。試驗(yàn)區(qū)高程的統(tǒng)計(jì)信息：最大值2 463．5 m、最小值55 m、平均值686 m、標(biāo)準(zhǔn)差459．6 m。同樣采用2160階EGM2008模型計(jì)算的距離地面1500 m高度的擾動(dòng)重力三分量作為“真值”檢驗(yàn)改進(jìn)格林公式以及Stokes公式的計(jì)算精度（見表6、表7）。

表6 積分半徑10′時(shí)Stokes公式計(jì)算的高度1500 m擾動(dòng)重力差值Tab．6 Calculation of disturbing gravity difference at the altitude of 1500 m by using Stokes formula with the integral radius of 10′ mGal

表7 積分半徑10′時(shí)改進(jìn)的格林公式計(jì)算的高度1500 m擾動(dòng)重力差值Tab．7 Calculation of disturbing gravity difference at the altitude of 1500 m by using improved Green formula with the integral radius of 10′ mGal

可以看出在實(shí)際的地形區(qū)，在徑向方向改進(jìn)的格林公式計(jì)算的精度要大大優(yōu)于Stokes公式，水平方向分量精度與Stokes公式相當(dāng)。結(jié)果表明以實(shí)際地形面為邊界面求解外部擾動(dòng)重力，改進(jìn)的格林公式相比Stokes公式結(jié)果較好。

為了檢驗(yàn)不同積分半徑、不同計(jì)算高度對(duì)于改進(jìn)的格林積分公式求解外部擾動(dòng)重力三分量結(jié)果的影響，分別以1°、30′、20′為積分半徑計(jì)算3000 m高度擾動(dòng)重力三分量，結(jié)果統(tǒng)計(jì)如表8—表10所示。

根據(jù)以上數(shù)據(jù)可以看出應(yīng)用改進(jìn)的格林公式計(jì)算擾動(dòng)重力三分量可以取得比較理想的結(jié)果，在1°的積分半徑內(nèi)，隨著積分半徑的增大，改進(jìn)的格林公式計(jì)算擾動(dòng)重力徑向分量的精度得到了提高。這是由于在1°的范圍內(nèi)，地面點(diǎn)與空中點(diǎn)具有一定的相關(guān)性，因此隨著半徑的增大，精度也就相應(yīng)地提高。

表9 積分半徑30′時(shí)改進(jìn)的格林公式計(jì)算的高度3000 m擾動(dòng)重力差值Tab．9 Calculation of disturbing gravity difference at the altitude of 3000 m by using improved Green formula with the integral radius of 30′ mGal

表10 積分半徑1°時(shí)改進(jìn)的格林公式計(jì)算的高度3000 m擾動(dòng)重力差值Tab．10 Calculation of disturbing gravity difference at the altitude of 3000 m by using improved Green formula with the integral radius of 1° mGal

6 結(jié) 論

格林公式提供了以地面作為邊界確定外部擾動(dòng)位的解析表達(dá)，它要求有地面擾動(dòng)位及擾動(dòng)重力矢量邊值條件。本文應(yīng)用格林公式，以格網(wǎng)數(shù)值地形面為邊界面，忽略較小的垂直面積上的擾動(dòng)重力水平分量積分，推導(dǎo)出以擾動(dòng)重力垂直分量（或重力異常）和高程異常差直接求取擾動(dòng)重力場(chǎng)元的公式，該結(jié)果規(guī)避了Molodenky問題經(jīng)典解法需要進(jìn)行地面數(shù)據(jù)的延拓和級(jí)數(shù)解的形式。這不僅減少了延拓過程的巨大工作量，也避免了延拓產(chǎn)生的不適定性誤差。核函數(shù)簡(jiǎn)單，易于計(jì)算，為確定地球外部重力場(chǎng)提供了一種思路。

［1］ SJ?BERG L E．On the Discrete Boundary Value Problem of Physical Geodesy with Harmonic Reduction to an Internal Sphere［R］．Stockholm：Royal Institute of Technology，1975．

［2］ WANG Zhengtao，DANG Yamin，CHAO Dingbo．Theory and Methodology of Ultra－h(huán)igh－degree Geopotential Model Determination［M］．Beijing：Surveying and Mapping Press，2011．（王正濤，黨亞民，晁定波．超高階地球重力位模型的確定的理論與方法［M］．北京：測(cè)繪出版社，2011．）

［3］ CHENG Luying．General Kernel Functions Based on Integral Transformation among Different Disturbing Potential Elements［J］．Acta Geodaetica et Cartographica Sinica，2013，42（2）：203－210．（程蘆穎．不同擾動(dòng)位泛函間的積分變換廣義核函數(shù) ［J］．測(cè)繪學(xué)報(bào)，2013，42（2）：203－210．）

［4］ SEEBER G．Satellite Geodesy［M］．2nd ed．Berlin：Walter de Gruyter，2003．

［5］ LIN Miao，DENKER H，MüLLER J．Regional Gravity Field Modeling Using Free－positional Point Masses［J］．Studia Geophysica et Geodaetica，2014，58（2）：206－227．

［6］ XU Houze，LU Zhonglian．Earth’s Gravity Feld and Geoid in China［M］．Beijing：The People’s Liberation Army Press，1997．（許厚澤，陸仲連．中國(guó)地球重力場(chǎng)與大地水準(zhǔn)面［M］．北京：解放軍出版社，1997．）

［7］ STOKES G G．On the Variation of Gravity on the Surface of the Earth ［J］． Transactions of the Cambridge Philosophical Society，1849，8（5）：672－695．

［8］ DENKER H，Regional Gravity Field Modeling：Theory and Practical Results［M］∥XU Guochang．Sciences of Geodesy－Ⅱ． Berlin Heidelbreg： Springer Verlag， 2013：185－291．

［9］ HEISKANEN W A，MORITZ H．Physical Geodesy［M］．LU Fukang，HU Guoli trans．Beijing：Surveying and Mapping Press，1979．（海斯卡涅 W A，莫里斯H．物理大地測(cè)量學(xué)［M］．盧?？?，胡國(guó)理，譯．北京：測(cè)繪出版社，1979．）

［10］ MORITZ H．Advance Physical Geodesy［M］．London：Abacus Press，1980．

［11］ HOFM ANN－WELLEN HOF B，MORIZE H．Physical Geodesy［M］．Heidelberg：Springer，2005．

［12］ MOLODENSKY M S，EREMEEV V F，YOURKINA M I．Methods for Study of the External Gravitational Field and Figure of the Earth．Works of Central Research Institute of Geodesy．Aerial Photography and Cartography［M］．Moscow：Geodetic Literature Press，1960．

［13］ MOLODENSKY M S，EREMEEV V F，YOURKINA M I．An Evaluation of Accuracy of Stokes’Series and of Some Attempts to Improve His Theory［J］．Bulletin Géodésique，1962，63（1）：19－37．

［14］ SU NKEL H．Point Mass Models and the Anomalous Gravitational Field of the Earth：Report No．328［R］．Columbus：Ohio State University，1981．

［15］ KATSAMBALOS K E．Simulation Studies on the Computation of the Gravity Vector in Space from Surface Data Considering the Topography of the Earth：Report No．323［R］．Columbus：Ohio State University，1981．

［16］ BJERHAMMAR A．A New Theroy of Geodetic Gravity［M］．Stockholm：Tekniska Hogskolan，1964．

［17］ FORSBERG R．A Study of Terrain Reductions，Density Anomalies，and Geophysical Inversion Methods in Gravity Field Modelling：Report No．355［R］．Columbus：Ohio State University，1984．

［18］ JIANG Tao，LI Jiancheng，DANG Yamin，et al．Regional Gravity Field Modeling Based on Rectangular Harmonic Analysis［J］．Science China：Earth Sciences，2014，57（7）：1637－1644．

［19］ LU Zhonglian．Theory and Method of Earth’s Gravity［M］．Beijing：The People’s Liberation Army Press，1996．（陸仲連．地球重力場(chǎng)理論與方法［M］．北京：解放軍出版社，1996．）

［20］ LU Zhonglian，WU Xiaoping，DING Xingbin，et al．Ballistic Missile Gravity［M］．Beijing：The Ba Yi Press，1993．（陸仲連，吳曉平，丁行斌，等．彈道導(dǎo)彈重力學(xué)［M］．北京：八一出版社，1993．）