黃謨濤，鄧凱亮，吳太旗，歐陽永忠，陳欣，熊雄，劉敏，王許

1 海軍研究院，天津 300061 2 自然資源部海洋環(huán)境探測(cè)技術(shù)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室，廣州 510300 3 91001部隊(duì)，北京 100830

0 引言

外部擾動(dòng)重力計(jì)算是地球重力場(chǎng)逼近建模研究?jī)?nèi)容的重要組成部分，也是解算大地測(cè)量邊值問題的主要應(yīng)用目標(biāo)之一.外部擾動(dòng)重力信息在航空航天器飛行軌跡精密計(jì)算和空間科學(xué)技術(shù)研究中具有重要的應(yīng)用價(jià)值(Hirvonen and Moritz,1963；Moritz,1966；Heiskanen and Moritz,1967；Cruz and Laskowski,1984；黃謨濤等,2005).除代表中長(zhǎng)波長(zhǎng)的全球重力場(chǎng)球諧展開位模型外，常用于計(jì)算地球外部擾動(dòng)重力三分量的解析函數(shù)模型還包括直接積分和等效源兩大類別，前者包含Stokes積分模型、表層積分模型和向上延拓積分模型(Heiskanen and Moritz,1967；劉長(zhǎng)弘,2016)；后者包含等效重力異常模型(Bjerhammar,1964；李建成等,2003)、質(zhì)點(diǎn)模型(Sünkel,1983；吳曉平,1984；黃金水和朱灼文,1995；黃謨濤等,2005；王建強(qiáng)等,2010)和虛擬單層密度模型(許厚澤等,1984；操華勝等,1985).計(jì)算模型的完備性及決定其計(jì)算效率高低的模型結(jié)構(gòu)復(fù)雜程度，是評(píng)價(jià)模型優(yōu)劣和實(shí)用性的兩個(gè)關(guān)鍵因素.黃謨濤等(2018,2019)已對(duì)前述不同類別計(jì)算模型的技術(shù)特點(diǎn)和適用條件做了分析比較和總結(jié)，為不同應(yīng)用領(lǐng)域選用合適的計(jì)算模型提供了理論參考依據(jù).

早在1986年，許厚澤等就曾針對(duì)利用直接積分法計(jì)算高空擾動(dòng)重力的模型改化和精度評(píng)估問題做過深入研究(許厚澤等，2014).但精密計(jì)算外部擾動(dòng)重力的難點(diǎn)主要集中在接近地球表面的超低空高度段，具體體現(xiàn)為地形效應(yīng)和積分奇異性兩個(gè)方面的干擾.通過聯(lián)合采用Stokes球面積分模型和Molodensky級(jí)數(shù)改正項(xiàng)(Heiskanen and Moritz,1967)、Bjerhammar等效重力異常(Bjerhammar,1964)、Moritz解析延拓改正項(xiàng)(Moritz,1980)等方法，均可有效補(bǔ)償?shù)涂盏匦涡?yīng)對(duì)計(jì)算參量的影響；針對(duì)積分奇異性干擾，黃謨濤等(2005)推出了基于環(huán)帶扇形數(shù)據(jù)分布的非奇異積分公式.根據(jù)實(shí)用化保障需求，黃謨濤等(2018)提出了基于插值替換的非奇異分段積分改進(jìn)公式.盡管上述兩組公式都成功避開了積分奇異性問題，但前者對(duì)輸入數(shù)據(jù)分布有特殊要求，其應(yīng)用范圍受到了較大限制；后者的計(jì)算精度取決于插值高度和插值模型的選擇及計(jì)算區(qū)域重力異常場(chǎng)變化的劇烈程度，具有一定的不確定性.針對(duì)海域應(yīng)用需求，黃謨濤等(2019)推出了以重力異常和高程異常作為輸入信息、基于表層積分模型的無奇異計(jì)算公式，解決了海域外部近地空間流動(dòng)點(diǎn)擾動(dòng)重力精密快速賦值問題.實(shí)際上，該研究思路同樣適用于解決以單一重力異常作為輸入信息的Stokes球面積分模型改化問題.綜合考慮計(jì)算模型選擇的靈活性需求和數(shù)據(jù)保障種類可能面臨的局限性制約，本文沿用黃謨濤等(2019)的研究思路，通過引入基于積分恒等式變換的移去恢復(fù)技術(shù)及平面局域泰勒級(jí)數(shù)展開方法，對(duì)計(jì)算外部擾動(dòng)重力三分量的Stokes全球積分模型進(jìn)行分步改化處理，旨在消除數(shù)值積分不連續(xù)性和積分核函數(shù)奇異性的影響，提高賦值模型在超低空高度段的計(jì)算精準(zhǔn)度.

1 地球外部擾動(dòng)重力全球積分模型

由地球重力場(chǎng)位理論得知，擾動(dòng)重力可由擾動(dòng)位求偏微分得到，地球外部擾動(dòng)位的全球積分公式為(Heiskanen and Moritz,1967)：

(1)

(2)

由式(1)可推得地球外部擾動(dòng)重力三分量計(jì)算式：

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

式(3)—(5)即為計(jì)算地球外部擾動(dòng)重力矢量的全球積分模型.利用該組公式計(jì)算擾動(dòng)重力三分量的精度主要取決于模型改化誤差、輸入數(shù)據(jù)分辨率及觀測(cè)噪聲等幾個(gè)方面的因素(Heiskanen and Moritz,1967；黃謨濤等,2005).關(guān)于數(shù)據(jù)質(zhì)量影響和數(shù)據(jù)需求分析方面的內(nèi)容，黃謨濤等(2019)已經(jīng)做了比較詳細(xì)的研究和討論，這里不再重復(fù).本文主要針對(duì)式(3)—(5)實(shí)現(xiàn)過程可能引入的模型改化誤差進(jìn)行分析和探討，并提出相應(yīng)的補(bǔ)償方法.

2 計(jì)算模型分析與改化

2.1 徑向分量積分模型分析與改化

2.1.1 積分模型穩(wěn)定性分析

由式(6)得知，當(dāng)計(jì)算點(diǎn)趨近于數(shù)據(jù)網(wǎng)格點(diǎn)時(shí)，即當(dāng)r→R和ψ→0時(shí)，會(huì)出現(xiàn)分母項(xiàng)l→0，積分核函數(shù)Fr(r,ψ)發(fā)生奇異.這說明利用式(3)計(jì)算超低空外部擾動(dòng)重力徑向分量時(shí)，會(huì)出現(xiàn)由核函數(shù)奇異引起的不確定性問題，無法保證計(jì)算結(jié)果的可靠性和有效性.盡管通過從積分域中直接扣除計(jì)算點(diǎn)所在的網(wǎng)格數(shù)據(jù)塊，理論上即可消除積分奇異性問題，但因該數(shù)據(jù)塊距離計(jì)算點(diǎn)最近，對(duì)計(jì)算結(jié)果的影響也就最大，故扣除該數(shù)據(jù)塊必然會(huì)對(duì)計(jì)算結(jié)果帶來一定程度甚至是不可忽略的影響，降低計(jì)算結(jié)果的精度(劉長(zhǎng)弘,2016).另一方面，可以證明，當(dāng)r→R和ψ→0時(shí)，由式(3)確定的外部擾動(dòng)重力徑向分量δgr并不嚴(yán)格收斂于球面上的理論預(yù)估值δgR，不符合數(shù)值逼近理論最基本的連續(xù)性要求(黃謨濤等,2019).具體證明如下.

根據(jù)Heiskanen和Moritz(1967)，由式(3)定義的地球外部擾動(dòng)重力徑向分量δgr可轉(zhuǎn)換為如下的外部擾動(dòng)位和重力異常的線性組合：

(11)

式中，擾動(dòng)位Tr由式(1)確定；外部重力異常Δgr計(jì)算表達(dá)式為：

(12)

=Kr1(r,ψ)+Kr2(r,ψ),(13)

(14)

(15)

此時(shí)，式(11)可改寫為三個(gè)積分和：

(16)

下面分別推證，當(dāng)r→R時(shí)，上述三個(gè)積分項(xiàng)的逼近結(jié)果.首先，依據(jù)式(1)和式(2)，不難推得：

(17)

(18)

式中，S(ψ)為傳統(tǒng)的球面Stokes函數(shù).其次，依據(jù)式(14)和式(16)，可得：

(19)

考慮到當(dāng)r→R和ψ→0時(shí)，式(19)中的積分核函數(shù)可能出現(xiàn)分子分母同時(shí)為零的不定式(0/0)情形，故需要對(duì)其做特殊處理.把式(19)的全球積分域σ劃分為σ1和σ2兩部分，σ1代表以計(jì)算點(diǎn)為中心、s1為半徑且無限接近于計(jì)算點(diǎn)的小球冠區(qū)域，σ2代表剩余部分(σ-σ1)的區(qū)域.因當(dāng)r→R時(shí)，Kr1(r,ψ)→0，對(duì)應(yīng)于σ2部分的積分項(xiàng)為零，故我們只需要對(duì)小球冠σ1部分的積分進(jìn)行討論.在很小的σ1區(qū)塊內(nèi)，可采用極坐標(biāo)系(s,α)，對(duì)由式(14)表示的積分核函數(shù)做平面近似處理，?。?/p>

r=R+h,R2dσ≈sdsdα,

式中，h代表計(jì)算點(diǎn)的大地高.此時(shí)，由式(14)表示的積分核函數(shù)可近似表示為：

(20)

將式(20)代入式(19)，得小球冠σ1部分也就是δgr1的積分為：

(21)

考慮到σ1是一個(gè)無限小的包圍計(jì)算點(diǎn)的球冠，可將Δg看作恒等于球面計(jì)算點(diǎn)處的重力異常ΔgR的不變量，即認(rèn)為Δg=ΔgR.將其代入式(21)，完成積分后得：

(22)

由式(22)不難看出，當(dāng)h→0(即r→R)時(shí)，有：

δgr1(r→R)=-ΔgR.

(23)

(23)式說明，由式(19)定義的第二項(xiàng)積分極限值收斂于地面上已知重力異常的負(fù)值.

最后，依據(jù)式(15)和式(16)，我們不難得到式(16)右端第三項(xiàng)積分的極限值為：

(24)

因式(15)代表的是重力異常Poisson積分核函數(shù)的零階和一階項(xiàng)，故只有重力異常的零階和一階項(xiàng)為零這樣的條件得到滿足時(shí)，積分式(24)的極限值才為零.

綜合式(17)、(23)和(24)，可得到式(16)的積分極限值為：

(25)

由式(24)和(25)得知，當(dāng)r→R時(shí)，δgr(r→R)的極限值并未嚴(yán)格收斂于球面上的理論預(yù)估值：δgR=-2TR/R-ΔgR，除非滿足假設(shè)條件：重力異常不包含零階和一階項(xiàng).

實(shí)際上，式(25)的右端之所以會(huì)出現(xiàn)多余的第三項(xiàng)，是因?yàn)槲覀冊(cè)趯?duì)式(13)右端表達(dá)式做改化處理時(shí)，人為實(shí)施了“恢復(fù)-移去”兩步運(yùn)算，即首先在式(14)中恢復(fù)零階和一階項(xiàng)影響，然后在式(15)中移去相對(duì)應(yīng)的影響項(xiàng).而由Heiskanen和Moritz(1967)得知，人們?cè)谕茖?dǎo)地球外部引力位球諧函數(shù)及其徑向?qū)?shù)積分計(jì)算式時(shí)，總是習(xí)慣于事先自動(dòng)“壓制”掉計(jì)算參量的零階和一階球諧函數(shù)項(xiàng)，使得它們只適用于這樣的參考橢球：①它與大地水準(zhǔn)面具有相同的重力位；②它的質(zhì)量與地球質(zhì)量相等；③它的中心位于地球質(zhì)心.顯然，一般參考橢球不可能完全滿足這些條件，因此，恢復(fù)積分核函數(shù)的零階和一階項(xiàng)，使原積分式成為適用范圍更廣的廣義積分式才是合情合理的.正如Heiskanen和Moritz(1967)所述：“傳統(tǒng)Stokes公式給出的擾動(dòng)位Ts一般不包含零階項(xiàng)T0和一階項(xiàng)T1，將Ts加上T0和T1才能得到完整的擾動(dòng)位T”.依據(jù)這樣的思路，Pizzetti和Heiskanen曾分別推出了更為通用的Stokes積分公式(Heiskanen and Moritz,1967).

基于上述考慮，我們認(rèn)為，在前面的改化處理時(shí)，只需要對(duì)式(13)做第一步的恢復(fù)處理，不必做第二步的移去運(yùn)算.也就是說，只需考慮式(13)右端的Kr1(r,ψ)作用，不必顧及第二項(xiàng)Kr2(r,ψ)的影響.此時(shí)，Kr(r,ψ)=Kr1(r,ψ)，式(25)簡(jiǎn)化為：

(26)

式(26)推演結(jié)果說明，在由式(6)表示的核函數(shù)中恢復(fù)零階和一階項(xiàng)影響，即可確保由式(3)確定的擾動(dòng)重力徑向分量極限值δgr(r→R)收斂于球面上的理論預(yù)估值δgR.顯然，這正是我們所期待的結(jié)果，但這里需要指出的是，式(26)只是理論上的理想化推證結(jié)果，在實(shí)際應(yīng)用中，我們是無法完全按照式(20)—(23)的推證過程來實(shí)現(xiàn)擾動(dòng)重力徑向分量數(shù)值計(jì)算的.因?yàn)?，選擇一個(gè)無限小的球冠σ1既不現(xiàn)實(shí)也不符合數(shù)據(jù)實(shí)際，而當(dāng)計(jì)算點(diǎn)與數(shù)據(jù)點(diǎn)重合時(shí)，在超低空高度段，由式(14)表示的積分核函數(shù)會(huì)出現(xiàn)非常嚴(yán)重的奇異性問題，從而導(dǎo)致式(12)的數(shù)值計(jì)算結(jié)果嚴(yán)重偏離預(yù)期的理論“真值”，這些都是由地球外部擾動(dòng)位(類同于單層位)法向?qū)?shù)在邊界面存在不連續(xù)性的固有特性和積分核函數(shù)(見式(14))具有與三維空間Dirac函數(shù)相似的跳躍特性所決定的(Heiskanen and Moritz,1967).

2.1.2 積分模型精密改化

由前面的分析得知，要想獲得穩(wěn)定且精確的擾動(dòng)重力徑向分量數(shù)值計(jì)算結(jié)果，除了需要對(duì)式(13)表示的核函數(shù)做零階和一階項(xiàng)恢復(fù)處理外，還需要對(duì)積分計(jì)算式實(shí)現(xiàn)過程中的核函數(shù)奇異性影響做進(jìn)一步的改化處理.

根據(jù)式(1)—(3)、式(6)及式(11)—(15)，恢復(fù)零階和一階項(xiàng)影響后的擾動(dòng)重力徑向分量積分核函數(shù)Ftr(r,ψ)可表示為：

Ftr(r,ψ)=Fr(r,ψ)+Kr2(r,ψ)

(27)

不難看出，當(dāng)r→R和ψ→0時(shí)，式(27)中的S(r,ψ)會(huì)出現(xiàn)分母項(xiàng)l→0，Kr1(r,ψ)可能出現(xiàn)分子分母同時(shí)為零(0/0)的不確定情況.為了消除此類核函數(shù)的奇異性影響，我們可參照Heiskanen和Moritz(1967)、黃謨濤等(2019)的研究思路，對(duì)以式(27)作為新的核函數(shù)的積分計(jì)算式(3)進(jìn)行積分恒等式變換改化處理，可將其稱為另一種形式的“移去-恢復(fù)”技術(shù).首先，將式(3)改寫為：

(28)

式中，ΔgRp代表外部空間計(jì)算點(diǎn)P(r,φ,λ)在球面上的投影點(diǎn)P(R,φ,λ)處的重力異常，也就是式(22)中的ΔgR.將式(27)代入式(28)右端的第二項(xiàng)積分，不難推得：

(29)

(30)

此時(shí)，式(28)可簡(jiǎn)寫為：

(31)

不難看出，經(jīng)上述移去恢復(fù)變換后，理論上，式(31)不再存在積分奇異性問題，至少可以說它的奇異性被中和了(Heiskanen and Moritz,1967).與此同時(shí)，當(dāng)r→R時(shí)，顧及式(17)、式(23)和式(29)后，可推得由式(31)確定的徑向分量δgr將收斂于球面上的理論預(yù)估值：δgRp=-2TRp/R-ΔgRp.這個(gè)結(jié)果說明，使用式(31)計(jì)算外部擾動(dòng)重力徑向分量，不僅可以避免奇異積分的影響，同時(shí)可確保積分計(jì)算值從地球外部到球邊界面的連續(xù)性.

從前面的推導(dǎo)過程可以看出，如果不對(duì)原積分核函數(shù)式(6)做如式(27)所示那樣的零階和一階項(xiàng)恢復(fù)處理，我們也無法實(shí)現(xiàn)式(29)和式(30)那樣的第二步積分恒等式變換運(yùn)算.由此可見，對(duì)于相同類別的地球外部重力場(chǎng)參數(shù)計(jì)算模型，包括外部擾動(dòng)位、重力異常、重力梯度等積分模型，都需要在它們的原積分核函數(shù)基礎(chǔ)上，恢復(fù)其零階和一階項(xiàng)的影響.只有這樣，才能開展后續(xù)旨在消除積分奇異性和不連續(xù)性影響的改化處理.實(shí)際上，Heiskanen和Moritz(1967)利用Poisson積分公式和積分恒等式變換，對(duì)調(diào)和函數(shù)徑向?qū)?shù)計(jì)算模型所做的去積分奇異性改化，正是建立在積分核函數(shù)事先包含了零階和一階項(xiàng)基礎(chǔ)之上的.

需要指出的是，前面對(duì)式(3)所做的模型改化只是這一進(jìn)程的第一階段，受觀測(cè)數(shù)據(jù)覆蓋范圍的限制，我們?cè)趹?yīng)用式(31)時(shí)，還需要將全球積分域劃分為近區(qū)和遠(yuǎn)區(qū)處理，近區(qū)定義為以計(jì)算點(diǎn)為中心、ψ0為半徑的球冠區(qū)域σ0，以一定階次(比如N階)的重力位模型作為參考場(chǎng)，聯(lián)合采用實(shí)測(cè)重力異常數(shù)據(jù)和移去恢復(fù)技術(shù)，對(duì)近區(qū)影響進(jìn)行數(shù)值積分計(jì)算；遠(yuǎn)區(qū)影響則采用更高階次(比如L階)的重力位模型進(jìn)行補(bǔ)償.第二次引入“移去-恢復(fù)”處理模式后，還需要對(duì)積分核函數(shù)做相應(yīng)的改化處理，以滿足積分核函數(shù)與實(shí)測(cè)重力異常信息之間的頻譜匹配要求(Novák and Heck,2002；劉敏等,2016).這里統(tǒng)一使用簡(jiǎn)單實(shí)用的Wong和Gore(1969)方法對(duì)積分核函數(shù)進(jìn)行改化.經(jīng)分區(qū)處理和核函數(shù)改化后，式(31)從全球積分模型改變?yōu)榫钟蚍e分模型：

(32)

(33)

(35)

(36)

(38)

由前面的推證過程得知，式(28)能夠過渡到式(31)是因?yàn)榇嬖诤愕仁?29)和(30).值得注意的是，上述兩個(gè)恒等式成立的條件是全球積分，當(dāng)式(31)被改化為局域積分式(32)時(shí)，該條件不再滿足，兩個(gè)恒等式失效，故必須顧及非全球積分域?qū)愕仁?29)和式(30)的影響，并在改化計(jì)算式(32)中加以補(bǔ)償.不難看出，由全球積分過渡到局域積分引起的模型誤差，就是計(jì)算點(diǎn)所在數(shù)據(jù)塊重力異常ΔgRp作用于計(jì)算參量δgr的遠(yuǎn)區(qū)影響，其計(jì)算式為：

(39)

(41)

式中，δgrp(σ-σ0)代表ΔgRp在積分遠(yuǎn)區(qū)(ψ0-π)對(duì)計(jì)算參量δgr的影響.將式(6)代入式(40)，完成積分后可得：

(42)

(43)

(44)

(45)

(46)

將式(39)加入式(32)的右端，可得到計(jì)算外部擾動(dòng)重力徑向分量的嚴(yán)密改化公式：

+δgr(σ-σ0)+δgrp(σ-σ0).

(47)

在后面的數(shù)值計(jì)算試驗(yàn)中，我們將對(duì)遠(yuǎn)區(qū)修正量δgrp(σ-σ0)的大小做進(jìn)一步的分析和驗(yàn)證.

2.2 水平分量積分模型改化

實(shí)際上，擾動(dòng)重力水平分量全球積分模型同樣面臨相類似的適用性改化問題.首先，與徑向分量相同，受觀測(cè)數(shù)據(jù)覆蓋范圍限制，水平分量計(jì)算模型也需要將全球積分域劃分為近區(qū)和遠(yuǎn)區(qū)處理(黃謨濤等,2019).為了減弱遠(yuǎn)區(qū)截?cái)嗾`差的影響，同樣需要引入位模型參考場(chǎng)進(jìn)行第二步移去恢復(fù)運(yùn)算和核函數(shù)改化.略去具體推導(dǎo)過程，直接寫出式(4)和式(5)從全球積分模型過渡到局域積分模型的改化公式如下：

(48)

(49)

(50)

(51)

(52)

(53)

(54)

×Rn,m(ψ0),(55)

此外，由式(4)、式(5)和式(7)得知，當(dāng)計(jì)算點(diǎn)與數(shù)據(jù)網(wǎng)格點(diǎn)完全重合時(shí)，該網(wǎng)格數(shù)據(jù)塊對(duì)兩個(gè)水平分量都不起作用，實(shí)施計(jì)算時(shí)可將其從積分域中扣除，以避免出現(xiàn)積分奇異性問題(黃謨濤等,2019).但是，當(dāng)網(wǎng)格數(shù)據(jù)塊的面積較大且計(jì)算點(diǎn)周圍的重力異常場(chǎng)變化比較劇烈時(shí)，這種簡(jiǎn)單的處理方法可能會(huì)給計(jì)算結(jié)果帶來毫伽(mGal=10-5m·s-2)級(jí)的誤差.顯然，對(duì)于高精度要求的擾動(dòng)引力計(jì)算，這樣的影響量仍不能忽略.為此，這里參照黃謨濤等(2019)的思路，進(jìn)一步推出計(jì)算點(diǎn)所在數(shù)據(jù)塊影響的積分計(jì)算式如下.

首先考慮到在計(jì)算點(diǎn)周圍超低空高度段范圍內(nèi)，計(jì)算點(diǎn)與積分流動(dòng)點(diǎn)之間的空間距離(l)相比地球平均半徑(R)是一個(gè)很小的量，可將擾動(dòng)重力水平分量積分核函數(shù)Fψ(r,ψ)(見式(7))做簡(jiǎn)化處理，只保留其中起主導(dǎo)作用的第一項(xiàng)：

(56)

假設(shè)與計(jì)算點(diǎn)重合的網(wǎng)格數(shù)據(jù)塊半徑為ψ00，因當(dāng)前可使用的重力觀測(cè)數(shù)據(jù)分辨率已經(jīng)達(dá)到較高的水平，相對(duì)應(yīng)的數(shù)據(jù)網(wǎng)格一般可達(dá)5′×5′甚至更小，故可進(jìn)一步對(duì)式(56)表示的積分核函數(shù)做平面近似處理.這里以北向分量為例，依據(jù)式(48)，略去截?cái)嗪撕瘮?shù)式(52)右端的第二項(xiàng)影響，采用前面定義過的極坐標(biāo)系(s,α)，可將與計(jì)算點(diǎn)重合數(shù)據(jù)塊的積分式寫為：

(57)

式中，s0代表數(shù)據(jù)網(wǎng)格大小的一半，當(dāng)數(shù)據(jù)網(wǎng)格為1′×1′時(shí)，s0=0.5′.在此基礎(chǔ)上，參照Heiskanen和Moritz(1967)的思路，將重力異常Δg在空間計(jì)算點(diǎn)P的球面投影點(diǎn)Rp處展開為泰勒級(jí)數(shù)：

+y2gyy)+…

(58)

式中，x軸指向正北，y軸向東，x=scosα，y=ssinα.并且有：

將式(58)代入式(57)，不難推得：

(59)

同理可得：

(60)

假設(shè)與計(jì)算點(diǎn)重合的數(shù)據(jù)格網(wǎng)為(i,j)，則可按(61)(62)式計(jì)算水平方向一階梯度：

gx=[Δg(i+1)-Δg(i-1)]/(4s0),(61)

gy=[Δg(j+1)-Δg(j-1)]/(4s0cosφi).

(62)

將補(bǔ)償計(jì)算式(59)和式(60)分別加到式(48)和式(49)的右端，就得到擾動(dòng)重力水平分量的嚴(yán)密改化模型：

+δgφ(σ-σ0)+δgφ ψ00,(63)

+δgλ(σ-σ0)+δgλ ψ00,(64)

式中，(σ0-ψ00)代表扣除計(jì)算點(diǎn)所在網(wǎng)格數(shù)據(jù)塊后的近區(qū).

3 數(shù)值計(jì)算檢驗(yàn)與分析

3.1 數(shù)值檢驗(yàn)使用的數(shù)據(jù)及區(qū)域

為了分析驗(yàn)證前面提出的地球外部擾動(dòng)重力三分量全球積分模型的改化效果，本文采用超高階位模型EGM2008作為數(shù)值計(jì)算檢驗(yàn)的參考標(biāo)準(zhǔn)場(chǎng)(Pavlis et al.,2012)，用于模擬產(chǎn)生地球表面(近似為球面)1′×1′網(wǎng)格重力異常觀測(cè)量“真值”及地球外部不同高度面上的1′×1′網(wǎng)格擾動(dòng)重力三分量理論“真值”(這里使用1′×1′而非5′×5′網(wǎng)格數(shù)據(jù)是為了減弱積分離散化誤差的影響).由地球位模型計(jì)算外部擾動(dòng)重力三分量的公式參見式(33)、式(50)和式(51)，計(jì)算球面重力異常的公式參見黃謨濤等(2005).

為了說明檢驗(yàn)結(jié)果的代表性，這里特意選取重力異常場(chǎng)變化比較劇烈的馬里亞納海溝作為試驗(yàn)區(qū)，具體覆蓋范圍為：6°×6°(φ:10°N—16°N；λ:142°E—148°E).首先選取截?cái)嗟?60階次的位模型EGM2008作為參考場(chǎng)，即取N=360，然后選取361～2160階次的位模型EGM2008作為數(shù)值計(jì)算檢驗(yàn)的標(biāo)準(zhǔn)場(chǎng)，即取L=2160，由該標(biāo)準(zhǔn)場(chǎng)模型計(jì)算球面上的1′×1′網(wǎng)格重力異常觀測(cè)量“真值”Δgt；進(jìn)而選取ri=R+hi，R=6371 km，使用標(biāo)準(zhǔn)場(chǎng)模型EGM2008(361～2160階次)分別計(jì)算對(duì)應(yīng)于9個(gè)高度面上的1′×1′網(wǎng)格擾動(dòng)重力三分量理論“真值”δgtri、δgt φ i和δgt λ i(i=1,2,…,9),每個(gè)高度面對(duì)應(yīng)360×360=129600個(gè)網(wǎng)格點(diǎn)數(shù)據(jù)，9個(gè)高度分別取為：hi=0 km,0.1 km,0.3 km,1 km,3 km,5 km,10 km,30 km,50 km.表1列出了其中的4個(gè)高度面上的擾動(dòng)重力三分量理論“真值”和球面上的重力異常觀測(cè)量“真值”的統(tǒng)計(jì)結(jié)果，圖1—4分別給出了球面重力異常觀測(cè)量“真值”和對(duì)應(yīng)于零高度面上的擾動(dòng)重力三分量理論“真值”分布態(tài)勢(shì)圖.

表1統(tǒng)計(jì)結(jié)果和圖1—4顯示的重力參數(shù)曲線變化形態(tài)說明，盡管已經(jīng)扣除掉2～360階次頻段的位模型參考場(chǎng)，本試驗(yàn)區(qū)域重力異常場(chǎng)變化的激烈程度仍然比較顯著，可在一定程度上代表真實(shí)地球大部分局部重力場(chǎng)的變化特征.

圖1 重力異常分布圖Fig.1 Distribution of the gravity anomalies

表1 由EGM2008模型(361～2160階次)計(jì)算得到的重力異常和擾動(dòng)重力統(tǒng)計(jì)結(jié)果(單位：mGal)Table 1 Statistics of gravity anomalies and disturbing gravities obtained by the EGM2008 model of degree 361～2160(unit:mGal)

3.2 改化模型數(shù)值檢驗(yàn)方法及結(jié)果分析

為了對(duì)比分析前述不同階段擾動(dòng)重力三分量改化模型的計(jì)算效果，這里采用球面上的1′×1′網(wǎng)格重力異?！罢嬷怠宝t作為觀測(cè)量，同時(shí)使用徑向分量3種改化模型和水平分量2種改化模型，對(duì)前面選定的試驗(yàn)區(qū)對(duì)應(yīng)于9個(gè)高度面上的1′×1′網(wǎng)格擾動(dòng)重力進(jìn)行計(jì)算分析，其中，徑向分量第1模型是指直接使用式(3)作為基礎(chǔ)計(jì)算模型，并對(duì)全球積分域做了分區(qū)處理，但在實(shí)施近區(qū)計(jì)算時(shí)，扣除掉計(jì)算點(diǎn)所在的1′×1′數(shù)據(jù)塊，以避免出現(xiàn)奇異積分問題；第2模型對(duì)應(yīng)于公式(32)；第3模型對(duì)應(yīng)于公式(47)；水平兩個(gè)分量第1模型是指直接使用式(48)和式(49)作為基礎(chǔ)計(jì)算模型，但在實(shí)施近區(qū)計(jì)算時(shí)，扣除掉與計(jì)算點(diǎn)重合的1′×1′數(shù)據(jù)塊；第2模型對(duì)應(yīng)于公式(63)和式(64).將三分量對(duì)應(yīng)于前述各階段改化模型的計(jì)算值分別與相對(duì)應(yīng)的理論“真值”δgtri、δgtφi和δgtλi做比較，可獲得不同改化模型的精度評(píng)估信息，具體比對(duì)結(jié)果列于表2.這里積分半徑統(tǒng)一取為ψ0=2°，為了減小積分邊緣效應(yīng)對(duì)評(píng)估結(jié)果的影響，表2只列出中心區(qū)2°×2°方塊內(nèi)的比對(duì)結(jié)果(下同).為了定量評(píng)估由全球積分過渡到局部積分引起的徑向分量模型誤差和由計(jì)算點(diǎn)所在網(wǎng)格重力異常變化引起的水平分量模型誤差影響，表3給出了采用式(39)計(jì)算得到的兩組分別對(duì)應(yīng)于積分半徑ψ0=2°和ψ0=5°的誤差補(bǔ)償量δgrp(σ-σ0)及采用式(59)和式(60)計(jì)算得到的誤差補(bǔ)償量δgφ ψ00和δgλ ψ00的統(tǒng)計(jì)結(jié)果.

圖2 擾動(dòng)重力徑向分量分布圖Fig.2 Distribution of the radial components of disturbing gravity

圖3 擾動(dòng)重力北向分量分布圖Fig.3 Distribution of the north components of disturbing gravity

圖4 擾動(dòng)重力東向分量分布圖Fig.4 Distribution of the east components of disturbing gravity

從表2比對(duì)結(jié)果可以看出，不同階段改化模型的計(jì)算精度存在較大差異，徑向分量的差異尤為明顯.表面上看，徑向分量第1模型的誤差是由于直接扣除了計(jì)算點(diǎn)所在數(shù)據(jù)塊的影響引起，實(shí)質(zhì)上是原始積分模型在邊界面存在不連續(xù)性所致.對(duì)比表2和表1結(jié)果可以看出，第1模型在超低空高度段的誤差量值甚至超過了徑向分量自身大小，顯然，這不是忽略計(jì)算點(diǎn)所在數(shù)據(jù)塊影響所能引起的量值，而是正如式(23)所顯示的那樣，當(dāng)r→R時(shí)，第1模型的原始計(jì)算式(3)存在與邊界面重力異常ΔgR大小大致相等的數(shù)值跳躍所致，與我們前面所做的理論分析預(yù)期相吻合.這個(gè)結(jié)果說明，徑向分量原始計(jì)算模型在超低空高度段是失效的，只有在5 km以上計(jì)算高度才是可用的.第2模型是對(duì)第1模型的改化，從理論上消除了積分奇異性和數(shù)值不連續(xù)性影響，并在超低空高度段取得了比第1模型好得多的計(jì)算精度，但由于該模型的改化過程存在不可忽略的理論缺陷，在3 km以上高度，該模型的計(jì)算精度反而不及第1模型，即使到了10 km高度，該模型的比對(duì)中誤差仍然超過不可接受的2 mGal.第3模型從理論上彌補(bǔ)了第2模型的缺陷，使得該模型的計(jì)算精度得到顯著改善，在所有9個(gè)高度面，該模型計(jì)算值與比對(duì)基準(zhǔn)“真值”的最大互差均不超過1.2 mGal，互差均方根值不超過0.4 mGal.這個(gè)結(jié)果說明，我們對(duì)第2模型所做的補(bǔ)償改化處理是正確和有效的.兩個(gè)水平分量第1和第2改化模型的區(qū)別主要體現(xiàn)在，后者比前者增加了計(jì)算點(diǎn)所在數(shù)據(jù)塊的影響.表2結(jié)果顯示，相比第1模型，第2模型計(jì)算精度在300 m以下超低空高度段的提高幅度比較明顯，從毫伽級(jí)提升到1 mGal以內(nèi)，充分體現(xiàn)了該模型的改化效果.可以預(yù)見，當(dāng)采用的數(shù)據(jù)網(wǎng)格間距加大(比如從1′×1′增大到2′×2′)且計(jì)算點(diǎn)周圍的重力異常場(chǎng)變化更為劇烈時(shí)，第2模型的改化效果會(huì)更加顯現(xiàn).這里需要指出的是，本試驗(yàn)是在一個(gè)重力位標(biāo)準(zhǔn)場(chǎng)中進(jìn)行的，表2統(tǒng)計(jì)結(jié)果反映的只是單一的計(jì)算模型誤差影響，不包含數(shù)據(jù)傳播誤差，也不包含受限于數(shù)據(jù)分辨率的數(shù)據(jù)截?cái)嗾`差.另外，對(duì)比表2中的徑向分量和水平分量精密改化模型的計(jì)算效果可以看出，前者的總體精度要略優(yōu)于后者，這說明積分模型離散化誤差對(duì)后者的影響要略大于前者.

表2 由不同改化模型計(jì)算得到的9個(gè)高度面擾動(dòng)重力三分量與“真值”的比較(單位：mGal)Table 2 Differences between the three components of disturbing gravity,obtained by different modified models,and the “true values”on 9 altitude surfaces (unit:mGal)

由表3計(jì)算結(jié)果可進(jìn)一步看出，盡管徑向分量第3模型對(duì)第2模型的補(bǔ)償量均隨參考場(chǎng)階數(shù)N、積分半徑ψ0和計(jì)算高度h的增大而減小，但當(dāng)參考場(chǎng)階數(shù)取為N=360時(shí)，即使積分半徑增大到ψ0=5°，計(jì)算高度h=10 km處的誤差補(bǔ)償量均方根值仍然超過1 mGal.在超低空高度段，兩個(gè)水平分量的誤差補(bǔ)償量也都超過1 mGal.不難推斷，如果取參考場(chǎng)階數(shù)為N=180，那么，相同高度上的誤差補(bǔ)償量還會(huì)增大.這樣的結(jié)果再次說明，對(duì)于高精度要求的地球外部重力場(chǎng)賦值，對(duì)傳統(tǒng)計(jì)算模型進(jìn)行精細(xì)改化是非常必要的.

表3 模型誤差補(bǔ)償量δgrp(σ-σ0)，δgφ ψ00和δgλ ψ00計(jì)算結(jié)果統(tǒng)計(jì)(單位：mGal)Table 3 Statistics of the computational results of error compensation δgrp(σ-σ0),δgφ ψ00 and δgλ ψ00 (unit:mGal)

為了考察數(shù)據(jù)觀測(cè)誤差對(duì)改化模型解算結(jié)果的影響，在前述試驗(yàn)基礎(chǔ)上，進(jìn)一步開展有輸入數(shù)據(jù)噪聲影響條件下的數(shù)值計(jì)算檢驗(yàn).具體做法是，在前面作為觀測(cè)量的位模型剩余重力異?！罢嬷怠宝t中分別加入1 mGal和3 mGal的隨機(jī)噪聲，形成兩組新的模擬觀測(cè)量，然后按照前面相同的計(jì)算方案和流程，依次采用前述(3+2+2)種改化模型完成9個(gè)高度面上的1′×1′網(wǎng)格擾動(dòng)重力三分量計(jì)算，最后將計(jì)算結(jié)果與相對(duì)應(yīng)高度的三分量“真值”做比對(duì)評(píng)估，具體評(píng)估結(jié)果如表4所示.為節(jié)省篇幅，這里只列出其中的比對(duì)互差均方根值(RMS).

表4 誤差影響下由不同改化模型計(jì)算得到的9個(gè)高度擾動(dòng)重力三分量與“真值”的互差均方根值(單位：mGal)Table 4 RMS of the differences between the three components of disturbing gravity,obtained by different modified models with noisy data,and the “true values”on 9 altitude surfaces (unit:mGal)

從表4統(tǒng)計(jì)結(jié)果可以看出，數(shù)據(jù)噪聲對(duì)7種改化模型解算結(jié)果的影響規(guī)律是一致的，沒有因?yàn)槟Ｐ透幕问降牟煌a(chǎn)生實(shí)質(zhì)性差異.總體而言，數(shù)據(jù)誤差只對(duì)5 km高度以下的計(jì)算結(jié)果產(chǎn)生一定程度的影響.對(duì)比表4和表2計(jì)算結(jié)果不難看出，1 mGal數(shù)據(jù)噪聲對(duì)各個(gè)改化模型解算結(jié)果的影響很小，幾乎可以忽略不計(jì)；3 mGal數(shù)據(jù)噪聲對(duì)徑向分量第1和第2模型、兩個(gè)水平分量第1模型的影響相對(duì)較小，對(duì)徑向分量第3模型、兩個(gè)水平分量第2模型的影響相對(duì)較大.這是因?yàn)樵谇懊娴慕聘幕Ｐ椭?，相?duì)于數(shù)據(jù)誤差，模型誤差的影響明顯占主導(dǎo)地位；而在后面的精密改化模型中，情況正好相反，模型誤差減小后，數(shù)據(jù)誤差影響起了主導(dǎo)作用.這些結(jié)果說明，即使采用嚴(yán)密的徑向分量第3模型和水平分量第2模型進(jìn)行外部擾動(dòng)重力計(jì)算，也要盡可能將數(shù)據(jù)觀測(cè)誤差控制在較低的水平.

4 結(jié)論

為了提高地球外部重力場(chǎng)的賦值精度，本文分析研究了外部擾動(dòng)重力Stokes全球積分計(jì)算模型的技術(shù)特點(diǎn)和適用條件，指出了擾動(dòng)重力徑向分量積分模型在邊界面存在不連續(xù)性的原因，同時(shí)提出了保持其連續(xù)性的修正方法.為了消除全球積分模型向局域積分轉(zhuǎn)換中遇到的積分奇異性問題，綜合采用移去恢復(fù)運(yùn)算和積分恒等式變換技術(shù)，同時(shí)依據(jù)實(shí)測(cè)數(shù)據(jù)保障條件，分別推出了地球外部擾動(dòng)重力三分量積分模型的分步改化公式，提出了補(bǔ)償傳統(tǒng)改化模型缺陷的修正公式.采用超高階地球位模型EGM2008建立比對(duì)標(biāo)準(zhǔn)重力異常場(chǎng)，同時(shí)選擇在重力異常場(chǎng)變化比較劇烈的馬里亞納海溝區(qū)塊開展數(shù)值計(jì)算符合度檢驗(yàn)，分別對(duì)本文推出的擾動(dòng)重力徑向分量3種分步改化模型和水平分量2種分步改化模型的計(jì)算精度進(jìn)行了檢核分析和評(píng)估.試驗(yàn)結(jié)果表明，采用最終的嚴(yán)密改化模型不僅可以有效消除原計(jì)算模型固有的積分奇異性，又可顯著提高超低空擾動(dòng)重力三分量的計(jì)算精度和穩(wěn)定性.因此，新的嚴(yán)密改化模型具有較高的推廣應(yīng)用價(jià)值，可用于地球外部全高度段擾動(dòng)重力場(chǎng)的高精度賦值.