重力異常垂向梯度嚴(yán)密改化模型及應(yīng)用

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

海軍研究院，天津 300061

0 引言

重力異常場是地球內(nèi)部質(zhì)量分布不均勻性及外部地形變化形態(tài)的綜合反映（Dehlinger,1978；Torge,1989）.利用地表重力觀測研究確定地球形狀大小及變化特性，是重力與大地測量學(xué)的核心任務(wù),其研究內(nèi)容包含大地測量邊值問題解算、局部和外部重力場逼近（Heiskanen and Moritz,1967；李建成等,2003；Sansò and Sideris,2013）；依據(jù)地表重力觀測研究確定地球內(nèi)部物質(zhì)結(jié)構(gòu)、形態(tài)及物理特性，則是地球物理學(xué)的重要研究主題之一，其研究內(nèi)容涉及利用數(shù)值計(jì)算方法解決地球重力場的正演和反演問題（Moritz,1990；馬在田等,1997）.由此可見，地球重力場研究不僅跨越了地球物理學(xué)和大地測量學(xué)兩個(gè)學(xué)科的研究領(lǐng)域，同時(shí)跨越了地球內(nèi)部、地表及外部空間三個(gè)不同特征的研究空域，其研究內(nèi)容除了涉及不同空域的多維度重力場建模技術(shù)外，還包含地表重力觀測數(shù)據(jù)的向上和向下延拓歸算技術(shù)（Moritz,1980；Cruz and Laskowski,1984；梁錦文,1989；王興濤等,2004；劉敏等,2018）.實(shí)際上，大地測量邊值問題解算和地球外部重力場逼近計(jì)算都可以歸結(jié)為廣義的重力位場向上和向下延拓問題（Moritz,1966；Heiskanen and Moritz,1967；Moritz,1980；于錦海等,2001；Sansò and Sideris,2017）.可見，重力向上和向下延拓技術(shù)在地球重力場研究中具有非常重要的應(yīng)用價(jià)值.實(shí)施向上和向下延拓解算除了可以依托傳統(tǒng)的Poisson積分方程外（Heiskanen and Moritz,1967），還可以采用計(jì)算穩(wěn)定性更好的Taylor級數(shù)展開模型（Peters,1949）.精確求取位場各階垂向?qū)?shù)（也稱垂向梯度）是運(yùn)用Taylor級數(shù)展開延拓法的關(guān)鍵，國內(nèi)外諸多學(xué)者為此提出了許多解決該問題的方法（Moritz,1980；劉保華等,1995；姚長利等,1997；Fedi and Florio,2002；Wei,2014；張沖等,2017；黃謨濤等,2018；Tran and Nguyen，2020）.實(shí)際上，除了向上和向下延拓解算，位場垂向梯度還可以直接應(yīng)用于地質(zhì)結(jié)構(gòu)反演和地球資源勘探，因?yàn)橹亓μ荻葦?shù)據(jù)更能精細(xì)反映地球內(nèi)部淺層地質(zhì)結(jié)構(gòu)密度的變化形態(tài)（王虎彪等,2009；劉金釗等,2013,2020）.

在前述各種應(yīng)用中，重力異常一階垂向梯度是應(yīng)用面最廣泛的重力場特征參數(shù)之一，在各類重力歸算和大地測量邊值解算中發(fā)揮著關(guān)鍵性作用（Heiskanen and Moritz,1967；Moritz,1980；Wang et al.,2012）.求解觀測重力異常的全球積分是獲取重力異常垂向梯度的主要手段，但實(shí)施此類積分模型數(shù)值計(jì)算時(shí)均涉及全球積分模型改化問題，一方面需要對全球積分計(jì)算式進(jìn)行積分域分割處理，通常做法是將全球積分劃分為近區(qū)和遠(yuǎn)區(qū)，近區(qū)采用實(shí)測數(shù)據(jù)進(jìn)行數(shù)值積分計(jì)算，遠(yuǎn)區(qū)則采用地球位模型進(jìn)行補(bǔ)償（黃謨濤等,2020）；另一方面需要對積分核函數(shù)進(jìn)行去奇異性處理，傳統(tǒng)做法是通過引入合適的積分恒等式變換，將原積分計(jì)算模型轉(zhuǎn)換為具有穩(wěn)定數(shù)值解的連續(xù)函數(shù)模型（Heiskanen and Moritz,1967）.需要指出的是，在實(shí)際應(yīng)用中，人們在實(shí)施全球積分域分割處理時(shí)，往往會忽視積分恒等式成立的全球積分條件，不再關(guān)注采用局域積分對積分恒等式帶來的數(shù)值影響（Heiskanen and Moritz,1967；Moritz,1980；歐陽明達(dá)等,2014；翟振和等,2015；劉長弘,2016），從而引起不可忽略的計(jì)算誤差（黃謨濤等，2019）.考慮到重力異常一階垂向?qū)?shù)是計(jì)算二階及其更高階導(dǎo)數(shù)的基礎(chǔ)，本文主要針對球面及外部重力異常一階垂向?qū)?shù)全球積分模型改化問題進(jìn)行分析研究和試驗(yàn)驗(yàn)證，依據(jù)實(shí)測數(shù)據(jù)保障條件和積分域分割處理方式，分別推出兩類參數(shù)全球積分模型的嚴(yán)密改化公式，同時(shí)通過數(shù)值計(jì)算和向下延拓應(yīng)用對比分析，進(jìn)一步驗(yàn)證采用嚴(yán)密改化模型的可行性和有效性.

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

1.1 地球外部重力異常垂向梯度嚴(yán)密改化模型

由Heiskanen和Moritz（1967）得知，地球外部重力異常全球積分計(jì)算式可表示為

（1）

（2）

由式（1）不難看出，當(dāng)計(jì)算點(diǎn)趨近于數(shù)據(jù)點(diǎn)時(shí)，即當(dāng)r→R和ψ→0時(shí)，會出現(xiàn)分子分母項(xiàng)同時(shí)為零的不定式（0/0）情況，積分核函數(shù)K（r,ψ）發(fā)生奇異.為了消除積分式（1）的奇異性，通常做法是從積分域中直接扣除計(jì)算點(diǎn)所在的網(wǎng)格數(shù)據(jù)塊，從而避免出現(xiàn)積分核函數(shù)分母項(xiàng)l→0現(xiàn)象.這種處理方法看似比較簡單明了，但往往會給計(jì)算結(jié)果帶來不可忽略的誤差.為此，Heiskanen和Moritz（1967）建議采用式（3）表示的積分恒等式對式（1）進(jìn)行改化：

（3）

略去推導(dǎo)過程，直接寫出改化公式如下：

（4）

式中，ΔgRp代表空間計(jì)算點(diǎn)P（r,φ,λ）在地面投影點(diǎn)P0（R,φ,λ）處的已知觀測值.式（4）即為人們常用的經(jīng)去奇異性改化后的地球外部重力異常全球積分計(jì)算公式.從形式上不難看出，引入積分恒等式（3）變換后，計(jì)算式（4）不再存在積分奇異性問題，同時(shí)確保當(dāng)r→R時(shí)，積分計(jì)算值Δgp收斂于球面觀測量ΔgRp.Heiskanen和Moritz（1967）認(rèn)為，經(jīng)式（3）變換后，至少可以說式（4）核函數(shù)的奇異性被中和了.于錦海等（2001）從理論上證明了式（4）右端奇異積分項(xiàng)在Chauchy主值意義下的存在性，從而為實(shí)施式（4）及其微分算子的數(shù)值計(jì)算提供了理論依據(jù).

地球外部重力異常徑向偏導(dǎo)數(shù)計(jì)算模型可直接由式（4）求微分得到

（5）

依據(jù)式（5）可推出重力異常一階徑向偏導(dǎo)數(shù)Δg′rp計(jì)算模型為

（6）

=K1（r,ψ）.

（7）

如前所述，受觀測數(shù)據(jù)覆蓋范圍的限制，在實(shí)際使用式（6）計(jì)算地球外部重力異常一階徑向偏導(dǎo)數(shù)時(shí)，通常需要將全球積分域劃分為近區(qū)和遠(yuǎn)區(qū)處理，近區(qū)定義為以計(jì)算點(diǎn)為中心、ψ0為半徑的球冠區(qū)域σ0，剩下的部分稱為遠(yuǎn)區(qū)（σ-σ0）.一般以一定階次（比如N階）的重力位模型作為參考場，聯(lián)合采用實(shí)測重力數(shù)據(jù)和移去恢復(fù)技術(shù)，對近區(qū)數(shù)據(jù)影響進(jìn)行數(shù)值積分計(jì)算；遠(yuǎn)區(qū)效應(yīng)則采用更高階次（比如L階）的重力位模型進(jìn)行補(bǔ)償（黃謨濤等，2020）.引入基于參考場的“移去-恢復(fù)”處理模式后，還需要對積分核函數(shù)作相應(yīng)的改化處理，以滿足積分核函數(shù)與觀測重力異常信息之間的頻譜匹配要求（Novák and Heck,2002；劉敏等，2016）.這里統(tǒng)一使用簡單實(shí)用的Wong和Gore（1969）方法對積分核函數(shù)進(jìn)行改化，即從原核函數(shù)中截?cái)嗟襞c位模型參考場相同階次的球諧展開式.經(jīng)分區(qū)處理和核函數(shù)改化后，計(jì)算式（6）從全球積分模型轉(zhuǎn)換為局域積分模型，也就是第一步改化模型：

（8）

×Pn（cosψ），

（9）

（10）

（11）

（12）

（13）

需要指出的是，式（8）并不是嚴(yán)密的改化模型，還不能作為最終的實(shí)用化公式使用.這是因?yàn)?，式?）是由全球積分式（6）改化為局域積分得來的，式（6）又是從積分恒等式（3）變換過來的，而恒等式（3）成立的前提條件是積分域覆蓋全球，由分區(qū)改化得來的計(jì)算式（8）對應(yīng)的是局域積分，顯然不滿足積分恒等式（3）的假設(shè)條件要求.盡管在式（8）的右端已經(jīng)通過超高階位模型計(jì)算值Δg′rq（σ -σ0）顧及了遠(yuǎn)區(qū)效應(yīng)的影響，但Δg′rq（σ -σ0）只代表式（8）右端積分項(xiàng)Δgq在遠(yuǎn)區(qū)（σ-σ0）的補(bǔ)償，并未顧及另一積分項(xiàng)ΔgRp在遠(yuǎn)區(qū)（σ-σ0）的影響.這就意味著，當(dāng)我們采用局部區(qū)域觀測數(shù)據(jù)完成式（8）計(jì)算時(shí)，其計(jì)算結(jié)果必然存在一定大小的系統(tǒng)性模型偏差.數(shù)值試驗(yàn)結(jié)果表明，要想獲得高精度的外部重力異常垂向梯度計(jì)算結(jié)果，必須消除該項(xiàng)誤差的影響.

從前面的分析得知，由全球積分過渡到局域積分引起計(jì)算式（8）的模型誤差可用公式表示為

（14）

（15）

（16）

式中，Δg′rp（σ-σ0）代表ΔgRp在積分遠(yuǎn)區(qū)（σ-σ0）對計(jì)算參量Δg′rp的影響.將式（7）代入式（15），并進(jìn)行積分運(yùn)算，可推得

+rRcosψ0（r2-5R2）]，

（17）

（18）

在式（8）右端加入模型誤差修正項(xiàng)Δg′rp（σ -σ0），可得到計(jì)算外部重力異常一階徑向偏導(dǎo)數(shù)的第二步改化模型：

+Δg′rq（σ -σ0）+Δg′rp（σ -σ0）.

（19）

按照同樣的思路可推得外部重力異常二階及以上高階偏導(dǎo)數(shù)相對應(yīng)的改化公式，但考慮到由式（5）定義的高階偏導(dǎo)數(shù)涉及復(fù)雜的觀測函數(shù)連續(xù)性和核函數(shù)強(qiáng)奇異性問題，我們擬另文作專題研討，這里不再展開討論.需要指出的是，在重力異常場變化比較劇烈的區(qū)域，使用式（19）計(jì)算重力異常垂向梯度還會帶來一定的模型誤差.這是因?yàn)椋谑剑?9）右端的積分項(xiàng)中，我們是把計(jì)算點(diǎn)所在的網(wǎng)格數(shù)據(jù)塊重力異常當(dāng)作常值ΔgRp處理的，該數(shù)據(jù)塊對計(jì)算參量Δg′rp的影響已經(jīng)在式（19）右端的第一項(xiàng)和最后一項(xiàng)中得到補(bǔ)償，在積分項(xiàng)中不再體現(xiàn)該數(shù)據(jù)塊的影響.但當(dāng)計(jì)算點(diǎn)附近的重力異常場變化比較劇烈時(shí)，再將計(jì)算點(diǎn)所在數(shù)據(jù)塊當(dāng)成常值處理可能帶來不可忽略的誤差，必須對其作相應(yīng)的補(bǔ)償.假設(shè)與計(jì)算點(diǎn)重合的網(wǎng)格數(shù)據(jù)塊半徑為ψ00，考慮到該數(shù)據(jù)塊是一個(gè)很小的區(qū)域，故可采用極坐標(biāo)系（s,α）對積分核函數(shù)作平面近似處理.取

略去（h/R）2及以上高階項(xiàng)影響，可將由式（7）表示的積分核函數(shù)近似表示為

（20）

此時(shí)，計(jì)算點(diǎn)所在數(shù)據(jù)塊的積分式可寫為

（21）

式中，Δg′rp0代表計(jì)算點(diǎn)所在數(shù)據(jù)塊重力變化特征對計(jì)算參量Δg′rp的影響；s0代表數(shù)據(jù)網(wǎng)格大小的一半，當(dāng)數(shù)據(jù)網(wǎng)格為1′×1′時(shí)，s0=0.5′.由式（21）得知，如果把中心數(shù)據(jù)塊的重力異常當(dāng)成常值ΔgRp看待，即認(rèn)為在計(jì)算點(diǎn)所在的數(shù)據(jù)網(wǎng)格內(nèi)處處滿足Δgq=ΔgRp，則有Δg′rp0=0.當(dāng)計(jì)算點(diǎn)附近的重力異常場變化比較劇烈時(shí)，可參照Heiskanen和Moritz（1967）的思路，將重力異常Δgq在空間計(jì)算點(diǎn)P的球面投影點(diǎn)Rp處展開為泰勒（Taylor）級數(shù)：

+…，

（22）

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

將式（22）代入式（21），同時(shí)考慮到修正項(xiàng)Δg′rp0本身的量值一般比較小，Taylor級數(shù)展開式（22）中的三階及以上各項(xiàng)的綜合影響可以忽略不計(jì)（Heiskanen and Moritz，1967），由此不難推得

（23）

假設(shè)與計(jì)算點(diǎn)重合的數(shù)據(jù)格網(wǎng)為（i,j），可按式（24）、（25）計(jì)算gxx和gyy：

（24）

（25）

將式（23）加入式（19）的右端，就得到計(jì)算外部重力異常垂向梯度的嚴(yán)密改化模型：

+Δg′rq（σ-σ0）+Δg′rp（σ-σ0）+Δg′rp0.

（26）

1.2 地面重力異常垂向梯度嚴(yán)密改化模型

相比地球外部重力異常垂向梯度，在地球重力場逼近計(jì)算實(shí)際應(yīng)用中，人們更關(guān)注的是地球表面或重力觀測面（近似為球面）上的垂向梯度.為此，這里專門給出球面重力異常一階徑向偏導(dǎo)數(shù)嚴(yán)密改化公式.

實(shí)際上，球面重力異常垂向梯度只是外部重力異常垂向梯度的一個(gè)特例.在式（6）和式（7）中，令r=R，可推得球面上的重力異常一階徑向偏導(dǎo)數(shù)Δg′Rp計(jì)算模型：

（27）

（28）

l0=2Rsin（ψ/2），

（29）

將式（28）代入式（27）得

（30）

式（30）就是重力歸算應(yīng)用中最常見的球面重力異常一階徑向偏導(dǎo)數(shù)計(jì)算模型（Heiskanen and Moritz,1967）.

顯然，在實(shí)際應(yīng)用中，同樣需要對式（30）右端的全球積分做分區(qū)改化處理，同時(shí)需要引入位模型參考場對已知參量和待求參量進(jìn)行移去恢復(fù)計(jì)算和積分核函數(shù)截?cái)嗵幚?基于與式（6）同樣的改化流程，首先將式（30）右端的全球積分改化為近區(qū)積分和遠(yuǎn)區(qū)球諧函數(shù)展開兩部分：

（31）

（32）

（33）

式中，Δg′Rq（σ-σ0）代表球面重力異常一階徑向偏導(dǎo)數(shù)遠(yuǎn)區(qū)效應(yīng)位模型計(jì)算值；Qn（Δg′R）為相對應(yīng)的Poisson積分核截?cái)嘞禂?shù)；其他符號意義同前.很顯然，式（31）—（33）可以直接由式（8）—（11）令r=R得到.類似于式（8），經(jīng)第一步改化后的式（31）同樣存在由全球積分過渡到局域積分引起的模型誤差，該誤差大小可由式（34）確定：

（34）

l0（ψ0）=2Rsin（ψ0/2），

（35）

式中，Δg′Rp（σ-σ0）代表ΔgRp在積分遠(yuǎn)區(qū)（σ-σ0）對計(jì)算參量Δg′Rp的影響.同樣可以證明，式（34）可直接由式（14）令r=R得到.

在式（31）的右端加入模型誤差修正項(xiàng)Δg′Rp（σ -σ0），可得到計(jì)算球面重力異常一階偏導(dǎo)數(shù)的第二步改化模型：

顯然，計(jì)算球面重力異常一階偏導(dǎo)數(shù)的嚴(yán)密改化公式還應(yīng)當(dāng)包含與式（23）相對應(yīng)的計(jì)算點(diǎn)所在數(shù)據(jù)塊的影響（Heiskanen and Moritz,1967）.在式（23）中令h=0（即r=R），可直接求得

（37）

式中，Δg′Rp0代表計(jì)算點(diǎn)所在數(shù)據(jù)塊重力變化特征對計(jì)算參量Δg′Rp的影響.式（37）與Heiskanen和Moritz（1967）的推導(dǎo)結(jié)果完全一致，將其加入式（36）的右端，就得到計(jì)算球面重力異常垂向梯度的嚴(yán)密改化模型：

+Δg′Rq（σ-σ0）+Δg′Rp（σ-σ0）+Δg′Rp0，

（38）

式中，Δg′Rq（σ-σ0）和Δg′Rp（σ-σ0）分別代表Δgq和ΔgRp在積分遠(yuǎn)區(qū)（σ-σ0）對計(jì)算參量Δg′Rp的影響，分別由式（32）和式（34）計(jì)算；Δg′Rp0代表計(jì)算點(diǎn)所在數(shù)據(jù)塊對計(jì)算參量Δg′Rp的影響，由式（37）計(jì)算.后面的數(shù)值計(jì)算檢驗(yàn)，將進(jìn)一步驗(yàn)證第一步改化公式（31）增加模型誤差兩個(gè)修正項(xiàng)Δg′Rp（σ-σ0）和Δg′Rp0的必要性及有效性.

1.3 重力異常垂向梯度向下延拓應(yīng)用

如前所述，向下延拓是重力異常垂向梯度最重要的應(yīng)用方向之一.除了我們常見的航空重力向地面延拓計(jì)算外（王興濤等,2004；黃謨濤等,2018），重力異常垂向梯度最具標(biāo)志性的應(yīng)用場景是，將地面重力異常延拓歸算到海平面或過計(jì)算點(diǎn)的水準(zhǔn)面，進(jìn)而用于大地測量邊值問題解算（Moritz,1980；于錦海等,2001）.實(shí)際上，基于現(xiàn)代邊值問題理論的Molodensky零階加一階項(xiàng)級數(shù)解可解釋為，首先將地面重力異常向下解析延拓到海平面，用Stokes積分求得海平面上的高程異常，再將該結(jié)果向上延拓到地面（Heiskanen and Moritz,1967；李建成等，2003）.其中，地面空間重力異常Δg向海平面延拓的計(jì)算公式可表達(dá)為

（39）

式中，Δg*代表海平面上的重力異常；h代表地面點(diǎn)的正常高；（?Δg/?h）為地面重力異常的垂向梯度，通常使用前面介紹的一階徑向?qū)?shù)（?Δg/?r）來替代.美國從20世紀(jì)90年代開始，每隔3至5年就會對作為國家高程基準(zhǔn)的大地水準(zhǔn)面模型進(jìn)行更新?lián)Q代，新發(fā)布的USGG2009模型在構(gòu)建過程中，就使用了式（39）作為地面重力異常的歸算模型（Wang et al.,2012）.顯然，與完整的向下解析延拓模型相比較（Moritz,1980），式（39）已經(jīng)事先省略了二次及更高次項(xiàng)的影響，關(guān)于這些高次項(xiàng)影響的討論已經(jīng)超出本文的研究范圍，故不再做更多的評述，這里主要就一階徑向?qū)?shù)（?Δg/?r）計(jì)算模型的完備性對重力異常向下延拓精度的影響進(jìn)行分析和驗(yàn)證，具體見后面的數(shù)值計(jì)算檢驗(yàn)環(huán)節(jié).

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

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

為了分析比較前述不同階段改化模型的計(jì)算效果，本文采用超高階位模型EGM2008作為數(shù)值計(jì)算檢驗(yàn)的標(biāo)準(zhǔn)場（Pavlis et al.,2012），用于模擬產(chǎn)生球面1′×1′網(wǎng)格重力異常觀測量“真值”（這里使用1′×1′而非5′×5′網(wǎng)格數(shù)據(jù)是為了減弱積分離散化誤差的影響），同時(shí)產(chǎn)生球面及外部設(shè)定高度的重力異常垂向梯度理論“真值”.由重力位模型計(jì)算地球外部重力異常及一階徑向偏導(dǎo)數(shù)的公式為（Heiskanen and Moritz,1967；黃謨濤等，2005）

（40）

（41）

式中各個(gè)符號意義同前.在式（40）和（41）中令r=R，可得到計(jì)算地球表面（球面）重力異常及其一階徑向偏導(dǎo)數(shù)的公式.

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

表1統(tǒng)計(jì)結(jié)果和圖1、圖2顯示的重力異常及一階徑向偏導(dǎo)數(shù)變化形態(tài)說明，盡管已經(jīng)扣除掉2～360階次頻段的參考場，本試驗(yàn)區(qū)域標(biāo)準(zhǔn)場重力變化的激烈程度仍然非常顯著，可在一定程度上代表真實(shí)地球大部分局部重力場的變化特征.

表1 由EGM2008模型（361～2160階次）計(jì)算得到的重力異常及其一階導(dǎo)數(shù)統(tǒng)計(jì)結(jié)果Table 1 Gravity anomalies and their first order derivatives obtained by the EGM2008 model of degree 361～2160

2.2 垂向梯度改化模型檢驗(yàn)結(jié)果分析

為了對比分析不同改化模型的計(jì)算效果，這里首先采用零高度面上的1′×1′網(wǎng)格重力異?！罢嬷怠宝t0作為觀測量，同時(shí)使用前述4種地球外部重力異常垂向梯度（Δg′rp）改化模型，對前面選定的試驗(yàn)區(qū)對應(yīng)于7個(gè)高度面上的1′×1′網(wǎng)格重力異常一階偏導(dǎo)數(shù)進(jìn)行數(shù)值計(jì)算檢驗(yàn)和分析，其中，第1模型是指直接對式（1）求徑向偏導(dǎo)數(shù)作為基礎(chǔ)計(jì)算模型，并對全球積分域作了分區(qū)處理，但在實(shí)施近區(qū)計(jì)算時(shí)，扣除掉與計(jì)算點(diǎn)重合的那個(gè)1′×1′數(shù)據(jù)塊，以避免出現(xiàn)奇異積分；第2模型對應(yīng)于公式（8）；第3模型對應(yīng)于公式（19）；第4模型對應(yīng)于公式（26）.將4種改化模型的計(jì)算值分別與相對應(yīng)高度面的理論“真值”Δg′ti作比較，可獲得不同改化模型的精度評估信息，具體比對結(jié)果列于表2.這里積分半徑統(tǒng)一取為ψ0=2°，為了減小積分邊緣效應(yīng)對評估結(jié)果的影響，表2只列出中心區(qū)2°×2°方塊內(nèi)的比對結(jié)果（下同）.為了定量評估由全球積分過渡到局域積分引起的模型誤差影響，這里同時(shí)給出了采用式（14）計(jì)算得到的兩組分別對應(yīng)于ψ0=2°和ψ0=5°的誤差補(bǔ)償量Δg′rp（σ-σ0）統(tǒng)計(jì)結(jié)果，具體見表3，該補(bǔ)償量是本文推出的嚴(yán)密改化公式中最重要的修正項(xiàng).

圖1 球面重力異常分布Fig.1 The gravity anomalies on the sphere

圖2 球面重力異常一階徑向偏導(dǎo)數(shù)分布Fig.2 The first order radial partial derivatives of the gravity anomaly on the sphere

在前述試驗(yàn)基礎(chǔ)上，我們進(jìn)一步采用同一高度面上的1′×1′網(wǎng)格重力異?！罢嬷怠宝ti作為觀測量，同時(shí)使用前述4種地面重力異常垂向梯度（Δg′Rp）改化模型，對相對應(yīng)7個(gè)高度面上的1′×1′網(wǎng)格重力異常一階偏導(dǎo)數(shù)進(jìn)行數(shù)值計(jì)算檢驗(yàn)和分析，其中，第1模型是指直接對式（1）求徑向偏導(dǎo)數(shù)并令r=R作為基礎(chǔ)計(jì)算模型，且對全球積分域作了分區(qū)處理，但在實(shí)施近區(qū)計(jì)算時(shí)，扣除掉與計(jì)算點(diǎn)重合的1′×1′數(shù)據(jù)塊，以避免出現(xiàn)奇異積分；第2模型對應(yīng)于公式（31）；第3模型對應(yīng)于公式（36）；第4模型對應(yīng)于公式（38）.將4種改化模型的計(jì)算值分別與相對應(yīng)高度面的理論“真值”Δg′ti作比較，可獲得不同改化模型的精度評估信息，具體比對結(jié)果列于表4.

首先，從表2互比結(jié)果可以看出，我們對地球外部重力異常垂向梯度積分模型所作的分階段改化處理，取得了符合預(yù)期的解算效果.第1模型的誤差看似主要源于直接扣除了計(jì)算點(diǎn)所在數(shù)據(jù)塊的影響，實(shí)質(zhì)上是由于該積分模型在邊界面存在不連續(xù)性所致.對比表2和表1結(jié)果可以看出，第1模型在1 km以下超低空高度段的誤差量值遠(yuǎn)遠(yuǎn)超過了垂向梯度自身大小，顯然，這不是忽略計(jì)算點(diǎn)所在數(shù)據(jù)塊影響所能引起的量值，而是第1模型原始計(jì)算式（直接對式（1）求徑向偏導(dǎo)數(shù)得到）在邊界面存在比較顯著的類似于質(zhì)面和質(zhì)體位那樣的數(shù)值跳躍所致，這是由地球重力位在邊界面存在不連續(xù)特性所決定的（Heiskanen and Moritz,1967）.這個(gè)結(jié)果說明，重力異常垂向梯度原始計(jì)算模型在超低空高度段是失效的，只有在5 km以上計(jì)算高度才是可用的.第2模型從理論上消除了第1模型的積分奇異性和數(shù)值不連續(xù)性影響，計(jì)算精度得到顯著提升，其相對檢核精度（指互差均方根/垂向梯度自身）都控制在20%以內(nèi)，但由于該模型的改化過程存在不可忽略的理論缺陷，在10 km以上高度，該模型的計(jì)算精度反而不及第1模型.第3模型從理論上彌補(bǔ)了第2模型的缺陷，使得該模型的計(jì)算精度得到進(jìn)一步改善，在超低空高度段，該模型的相對計(jì)算精度不低于7%，在1 km以上高度段，相對精度優(yōu)于3%.這個(gè)結(jié)果說明，我們對第2模型所作的修正和補(bǔ)償處理是正確且有效的.第4模型是在第3模型基礎(chǔ)上，增加了計(jì)算點(diǎn)所在數(shù)據(jù)塊重力變化特征對計(jì)算參量的影響，表2結(jié)果顯示，相比第3模型，第4模型計(jì)算精度在300 m以下超低空高度段又得到一定程度的提升，在零高度面，其相對計(jì)算精度從6.7%提升到4.0%，提升幅度超過40%，充分體現(xiàn)了該模型的改化效果.可以預(yù)見，當(dāng)采用的數(shù)據(jù)網(wǎng)格間距加大（比如從1′×1′增大到2′×2′）且計(jì)算點(diǎn)周圍的重力異常場變化更為劇烈時(shí)，第4模型的改化效果會更加顯現(xiàn).

表2 由不同外部改化模型計(jì)算得到的7個(gè)高度面重力異常一階導(dǎo)數(shù)與“真值”比較（單位：mGal·km-1）Table 2 Comparisons between the first order radial partial derivatives of gravity anomalies,obtained by different modified models outside the earth,and the “true values”on 7 altitude surfaces （unit:mGal·km-1）

表3 模型誤差補(bǔ)償量Δg′rp（σ-σ0）計(jì)算結(jié)果統(tǒng)計(jì)（單位：mGal·km-1）Table 3 Statistics of the computational results of model error compensation Δg′rp（σ-σ0） （unit:mGal·km-1）

表4 由不同地面改化模型計(jì)算得到的7個(gè)高度面重力異常一階導(dǎo)數(shù)與“真值”比較（單位：mGal·km-1）Table 4 Comparisons between the first order radial partial derivatives of gravity anomalies,obtained by different modified models on the surface,and the “true values”on 7 altitude surfaces （unit:mGal·km-1）

由表3計(jì)算結(jié)果可進(jìn)一步看出，盡管第3模型對第2模型的補(bǔ)償量均隨參考場階數(shù)N、積分半徑ψ0和計(jì)算高度h的增大而減小，但當(dāng)參考場階數(shù)取N=360時(shí)，即使積分半徑增大到ψ0=5°，7個(gè)高度面的誤差補(bǔ)償量均方根值仍然接近甚至超過垂向梯度自身大小的10%.這樣的結(jié)果再次說明，對于高精度要求的地球重力場逼近計(jì)算，對重力異常垂向梯度傳統(tǒng)積分模型進(jìn)行精細(xì)改化和校正是非常必要的.

對比表4和表2計(jì)算結(jié)果可以看出，在相同的數(shù)據(jù)分辨率和精度保障條件下，利用某一高度面的重力異常觀測數(shù)據(jù)計(jì)算該高度面外部的重力異常垂向梯度，其精度都要比利用本高度面觀測數(shù)據(jù)計(jì)算本高度面的垂向梯度精度高.這個(gè)結(jié)果顯然跟重力異常垂向梯度積分計(jì)算模型誤差隨計(jì)算高度升高而衰減有關(guān).相比較而言，因第1模型在邊界面存在較大的數(shù)值跳躍問題，故利用該模型和本高度面觀測數(shù)據(jù)計(jì)算本高度面垂向梯度的結(jié)果偏離理論“真值”的幅度最大，完全失去了其使用價(jià)值.第2至第4模型兩種方式計(jì)算得到的垂向梯度精度差異相對較小，但在條件允許情況下，仍應(yīng)優(yōu)先采用前一種脫離邊界面的方式進(jìn)行垂向梯度計(jì)算.

2.3 改化模型向下延拓應(yīng)用效果分析

為了考察垂向梯度計(jì)算模型改化誤差對重力異常向下延拓解算結(jié)果的影響，這里特別設(shè)計(jì)如下試驗(yàn)流程：

步驟一：使用地球外部6個(gè)高度面上的重力異常一階徑向偏導(dǎo)數(shù)理論“真值”Δg′ti，依據(jù)公式（39）將對應(yīng)于6個(gè)高度面上的重力異常Δgti向下延拓到零高度面，分別將各個(gè)高度面的延拓計(jì)算值與零高度面的理論“真值”Δgt0作比較，計(jì)算互差統(tǒng)計(jì)結(jié)果.

步驟二：將式（39）中的垂向梯度替換為與表2統(tǒng)計(jì)結(jié)果相對應(yīng)的外部重力異常一階徑向偏導(dǎo)數(shù)（Δg′rp）4種改化模型的計(jì)算結(jié)果，重復(fù)步驟一的試驗(yàn).

步驟三：將式（39）中的垂向梯度替換為與表4統(tǒng)計(jì)結(jié)果相對應(yīng)的地面重力異常一階徑向偏導(dǎo)數(shù)（Δg′Rp）4種改化模型的計(jì)算結(jié)果，重復(fù)步驟一的試驗(yàn).

前述三步驟計(jì)算統(tǒng)計(jì)結(jié)果列于表5，為節(jié)省篇幅，表中只列出其中的對比互差均方根值（RMS）.

從表5統(tǒng)計(jì)結(jié)果可以看出，垂向梯度計(jì)算模型誤差直接影響重力異常向下延拓的解算精度，對積

表5 不同改化模型計(jì)算一階導(dǎo)數(shù)用于重力異常向下延拓與“真值”比較（互差均方根值）（單位：mGal）

分計(jì)算模型的修正和改化效果已經(jīng)在向下延拓的解算結(jié)果中得到充分體現(xiàn).不難看出，使用球面外部第3和第4模型計(jì)算得到的垂向梯度（Δg′rp）實(shí)施重力異常向下延拓解算的效果，已經(jīng)完全等同于使用垂向梯度理論“真值”（Δg′ti）獲得的解算效果，兩個(gè)模型計(jì)算值與使用“真值”計(jì)算結(jié)果的差異不超過0.1 mGal；而使用第1模型時(shí)，兩者的差異超過10 mGal；使用第2模型時(shí)，兩者的差異也超過1 mGal,這些結(jié)果再次驗(yàn)證了嚴(yán)密改化模型的有效性.需要指出的是，重力異常向下延拓的解算精度除了與垂向梯度計(jì)算精度水平密切相關(guān)外，還取決于向下延拓模型自身的完備性、延拓高度大小及計(jì)算區(qū)域重力場變化的劇烈程度.就本試驗(yàn)而言，由表5可以看出，將第1模型排除在外，使用只顧及到一階項(xiàng)的向下延拓模型即公式（39）進(jìn)行重力異常歸算，要想得到優(yōu)于1 mGal的解算精度，必須將向下延拓高度控制在1 km以內(nèi)，否則，需要將向下延拓模型拓展到更高階次（黃謨濤等,2018）.這方面的內(nèi)容已經(jīng)超出本文的研究范圍，這里不再做深入討論.

3 結(jié)論

將全球積分模型改化為局域模型是實(shí)現(xiàn)重力異常垂向梯度精密計(jì)算的前提條件.本文分析研究了重力異常垂向梯度全球積分計(jì)算模型的技術(shù)特點(diǎn)和適用條件，指出了開展積分模型精密改化的必要性和可行性.針對全球積分模型向局域積分轉(zhuǎn)換中遇到的積分奇異性和不連續(xù)性問題，綜合采用積分恒等式變換和移去恢復(fù)換算技術(shù)，同時(shí)依據(jù)實(shí)測數(shù)據(jù)保障條件，分別推出了計(jì)算地球外部及地面重力異常垂向梯度全球積分模型的分步改化公式，提出了補(bǔ)償傳統(tǒng)改化模型理論缺陷的修正公式.采用超高階地球位模型EGM2008作為模型比對標(biāo)準(zhǔn)重力異常場，同時(shí)選擇在重力異常場變化比較劇烈的馬里亞納海溝區(qū)塊開展數(shù)值計(jì)算試驗(yàn)，分別對本文推出的重力異常垂向梯度兩類8種分步改化模型的計(jì)算精度及向下延拓應(yīng)用效果進(jìn)行了檢核分析和評估.試驗(yàn)結(jié)果表明，采用最終的嚴(yán)密改化模型不僅可以有效消除原計(jì)算模型固有的積分奇異性和數(shù)值跳躍問題，又可顯著提高超低空重力異常垂向梯度的計(jì)算精度和穩(wěn)定性，有效提升重力異常向下延拓的解算精度水平.因此，新的嚴(yán)密改化模型具有較高的推廣應(yīng)用價(jià)值，可用于地球表面及外部重力場的高精度逼近計(jì)算.