重力異常向上延拓嚴(yán)密改化模型及向下延拓應(yīng)用

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

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

依據(jù)地表重力觀測研究確定地球形狀大小及其變化特性，是重力和大地測量學(xué)的核心研究任務(wù)，既包含傳統(tǒng)的大地測量邊值問題解算和局部重力場逼近，又涵蓋外部重力場逼近計(jì)算和全球重力場建模[1-5]。外部重力場逼近計(jì)算主要有兩個方面的重要應(yīng)用，一是為精密確定運(yùn)載火箭、人造衛(wèi)星、宇宙飛船、導(dǎo)彈武器、航天飛機(jī)等航天器飛行軌跡提供重力異常場補(bǔ)償[6-7]；二是為近地空間開展重力測量數(shù)據(jù)質(zhì)量評估提供比對基準(zhǔn)[1,8-9]。重力異常向上延拓是外部重力場逼近計(jì)算最重要的研究內(nèi)容之一，它既可直接用于地球外部重力異常計(jì)算，又可間接用于重力異常向下延拓迭代解算[2,10-14]。這說明互為反問題的重力異常向上和向下延拓不僅在計(jì)算地球物理學(xué)研究中發(fā)揮著不可替代的作用[15]，在大地測量邊值問題“先向下后向上”解算中也具有非常重要的應(yīng)用價(jià)值[1,2,16-17]。Poisson積分方程是計(jì)算地球外部重力異常的基礎(chǔ)數(shù)學(xué)模型，它源于大地測量底律希勒(Dirichlet)邊值問題的球面解，因其不能顧及地形高度的起伏變化，故在陸域?qū)嵤┲亓Ξ惓Ｏ蛏涎油赜?jì)算時，必須在Poisson積分模型基礎(chǔ)上補(bǔ)償?shù)匦涡?yīng)的影響。早在1966年，Moritz就推出了一組依據(jù)地面重力異常和地形高數(shù)據(jù)計(jì)算外部重力異常嚴(yán)密而復(fù)雜的積分公式[10]。但由于受數(shù)據(jù)保障條件的限制，在實(shí)際應(yīng)用中，人們更多采用的是所謂“移去-恢復(fù)”技術(shù)，即首先移去地形質(zhì)量對觀測重力異常的影響，然后使用Poisson積分完成向上延拓計(jì)算，最后在計(jì)算結(jié)果中恢復(fù)地形質(zhì)量的影響[7,17]。顯然，不管是采用顧及地形效應(yīng)的嚴(yán)密公式，還是“移去-恢復(fù)”技術(shù)，都離不開Poisson積分方程的理論支撐和模型工程化實(shí)現(xiàn)。因此，精密求解Poisson積分方程自然就成為獲取高精度重力異常向上延拓計(jì)算結(jié)果的核心關(guān)鍵。

計(jì)算地球外部重力異常的向上延拓公式是一類有代表性的全球積分模型，其積分核函數(shù)具有與三維空間Dirac函數(shù)相似的特性，當(dāng)計(jì)算點(diǎn)由外部空間趨于數(shù)據(jù)觀測界面時，積分計(jì)算會出現(xiàn)由Dirac核函數(shù)引起的數(shù)值跳躍現(xiàn)象；而當(dāng)計(jì)算點(diǎn)趨于數(shù)據(jù)點(diǎn)或與數(shù)據(jù)點(diǎn)重合時，此類積分式還會出現(xiàn)不同程度的奇異性問題，需要通過引入合適的積分恒等式變換，將原計(jì)算模型轉(zhuǎn)換為具有穩(wěn)定數(shù)值解的連續(xù)函數(shù)模型[1]。其次是，受觀測數(shù)據(jù)覆蓋范圍的限制，實(shí)施計(jì)算時需要對此類全球積分模型進(jìn)行積分域分割處理，即將全球積分劃分為近區(qū)和遠(yuǎn)區(qū)，近區(qū)采用實(shí)測重力數(shù)據(jù)進(jìn)行數(shù)值積分計(jì)算，遠(yuǎn)區(qū)則采用地球位模型進(jìn)行補(bǔ)償[4]。需要指出的是，第一階段為消除奇異積分而實(shí)施的積分恒等式變換通常是在全球積分意義下完成的，在第二階段實(shí)施全球積分域分區(qū)處理時，人們往往會忽視積分恒等式成立的全球積分條件[1,18-19]，不再關(guān)注采用局域積分對積分恒等式帶來的數(shù)值影響，從而引起不可忽略的計(jì)算模型誤差[20]。為此，本文針對地面重力異常向上延拓全球積分模型改化問題進(jìn)行分析研究，依據(jù)實(shí)際應(yīng)用中的數(shù)據(jù)保障條件和全球積分域分割處理方式，推導(dǎo)出外部重力異常全球積分模型的嚴(yán)密改化公式，同時通過數(shù)值計(jì)算檢驗(yàn)，分析驗(yàn)證采用嚴(yán)密改化模型的必要性和有效性。

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

1.1 重力異常向上延拓全球積分模型及穩(wěn)定性分析

由重力位場理論得知，依據(jù)球外Dirichlet問題球面解的Poisson積分方程，可直接推得由地面觀測重力異常計(jì)算地球外部重力異常的全球積分式為[1]

(1)

Pn(cosψ)

(2)

首先，將全球積分域σ劃分為σ1和σ2兩部分，σ1代表以計(jì)算點(diǎn)為中心、s1為半徑且無限接近于計(jì)算點(diǎn)的小球冠區(qū)域，σ2代表剩余部分(σ-σ1)的區(qū)域。因當(dāng)r→R，ψ≠0時，K(r,ψ)→0，故對應(yīng)于σ2部分的積分項(xiàng)為零?？梢?，只需分析討論對應(yīng)于小球冠σ1部分的積分項(xiàng)?？紤]到σ1是一個很小的區(qū)塊，故可在該區(qū)塊內(nèi)對由式(2)表示的積分核函數(shù)作平面近似處理，以極坐標(biāo)系(s,α)表示，可近似取

式中，h代表計(jì)算點(diǎn)的高度。此時，式(2)表示的積分核函數(shù)可近似表示為

(3)

將式(3)代入式(1)，得σ1部分的積分為

(4)

在無限小的球冠σ1內(nèi)，可將Δgq看作一個常值，恒等于計(jì)算點(diǎn)處的重力異常，即可取Δgq=ΔgRp，代入式(4)，得

(5)

由式(5)看出，當(dāng)h→0(即r→R)時，Δgpσ1→ΔgRp，即計(jì)算值收斂于地面上已知的重力異常觀測值。很顯然，這只是理論上的理想化推證結(jié)果，在實(shí)際應(yīng)用中，是無法完全按照上述推證過程來實(shí)現(xiàn)重力異常向上延拓?cái)?shù)值計(jì)算的。因?yàn)?，選擇一個無限小的球冠σ1既不現(xiàn)實(shí)也不符合數(shù)據(jù)實(shí)際，而當(dāng)計(jì)算點(diǎn)與數(shù)據(jù)點(diǎn)重合時，在超低空高度段，由式(2)表示的積分核函數(shù)會出現(xiàn)非常嚴(yán)重的奇異性問題，從而導(dǎo)致式(1)的數(shù)值計(jì)算結(jié)果嚴(yán)重偏離預(yù)期的理論“真值”，這些都是由地球外部擾動位(類同于單層位)法向?qū)?shù)在邊界面存在不連續(xù)性的固有特性和積分核函數(shù)具有與三維空間Dirac函數(shù)相似的跳躍特性所決定的[1]。盡管通過直接扣除計(jì)算點(diǎn)所在的數(shù)據(jù)塊，理論上即可消除積分奇異性問題，但因該數(shù)據(jù)塊距離計(jì)算點(diǎn)最近，對計(jì)算結(jié)果的影響也就最大，故扣除該數(shù)據(jù)塊必然會對計(jì)算結(jié)果帶來一定程度甚至是不可忽略的影響，從而顯著降低向上延拓結(jié)果的計(jì)算精度[19]。

1.2 全球積分模型嚴(yán)密改化公式

為了消除式(1)的奇異性，確保地球外部重力位場調(diào)和函數(shù)在邊界面上的連續(xù)性，Heiskanen和Moritz[1]建議采用由式(6)定義的積分恒等式對式(1)進(jìn)行改化

(6)

在式(6)兩端同乘以地面計(jì)算點(diǎn)處的已知觀測值ΔgRp，可得

(7)

將式(1)減去式(7)，并略加整理可得

(8)

式(8)即為經(jīng)去奇異性改化后的地球外部重力異常全球積分計(jì)算公式。不難看出，經(jīng)積分恒等式(6)變換后，式(8)不再存在積分奇異性問題。同時，理論上，當(dāng)r→R和ψ→0時，由式(8)確定的外部重力異常Δgp收斂于球面上已知的觀測值ΔgRp?？梢?，經(jīng)第一步改化后的計(jì)算式(8)在邊界面也不再出現(xiàn)數(shù)值跳躍現(xiàn)象。至此，就從理論上解決了式(1)存在的積分奇異和數(shù)值不連續(xù)性問題。

如前所述，受觀測數(shù)據(jù)覆蓋范圍的限制，在實(shí)際應(yīng)用式(8)時，通常需要將全球積分域劃分為近區(qū)和遠(yuǎn)區(qū)處理，近區(qū)定義為以計(jì)算點(diǎn)為中心、ψ0為半徑的球冠區(qū)域σ0，通過引入一定階次(比如N階)的重力位模型參考場和移去恢復(fù)技術(shù)，在近區(qū)采用實(shí)測數(shù)據(jù)(減去N階重力位模型)進(jìn)行計(jì)算，遠(yuǎn)區(qū)效應(yīng)則采用更高階次(比如L階)的重力位模型進(jìn)行補(bǔ)償。引入“移去-恢復(fù)”處理模式后，還需要對積分核函數(shù)作相應(yīng)的改化處理，以滿足積分核函數(shù)與觀測重力異常信息之間的頻譜匹配要求[21-22]。這里統(tǒng)一使用簡單實(shí)用的Wong和Gore[23]方法對積分核函數(shù)進(jìn)行改化。經(jīng)分區(qū)處理和核函數(shù)改化后，式(8)從全球積分模型轉(zhuǎn)換為局域積分模型，即得第2步改化公式

dσ+Δgq(σ-σ0)

(9)

(10)

(11)

(12)

(13)

(14)

需要指出的是，當(dāng)把式(9)作為計(jì)算外部重力異常的實(shí)用公式使用時，往往會忽略一個事實(shí)：式(8)是從積分恒等式(6)變換過來的，兩者都建立在全球積分基礎(chǔ)之上[1,19]。很顯然，分區(qū)改化模型計(jì)算式(9)中的局域積分并不滿足恒等式(6)的理論假設(shè)要求，因?yàn)棣q(σ-σ0)只代表對式(9)右端積分項(xiàng)Δgq在遠(yuǎn)區(qū)(ψ0-π)的補(bǔ)償，并未顧及另一積分項(xiàng)ΔgRp在遠(yuǎn)區(qū)(ψ0-π)的影響。這就意味著，一直在使用的重力異常向上延拓計(jì)算模型即式(9)是不嚴(yán)密的，存在系統(tǒng)性偏差。數(shù)值試驗(yàn)結(jié)果表明，由此引起的外部重力異常計(jì)算誤差隨觀測值ΔgRp和積分半徑ψ0取值大小變化而變化，最大可達(dá)數(shù)十mGal(mGal=10-5m/s2)。因此，要想獲得高精度的外部重力異常計(jì)算結(jié)果，必須消除該項(xiàng)誤差的影響。由前面的分析不難得知，該項(xiàng)誤差的計(jì)算式可表示為

(15)

式中，Δgp(σ-σ0)代表ΔgRp在積分遠(yuǎn)區(qū)(σ-σ0)的影響。將式(10)代入式(15)，得

(16)

完成式(16)積分后可得

(17)

(18)

(19)

(20)

由式(8)推導(dǎo)過程得知，計(jì)算點(diǎn)所在數(shù)據(jù)塊ΔgRp對計(jì)算參量Δgp的影響已經(jīng)在式(8)右端的第一項(xiàng)中得到補(bǔ)償，故在式(20)右端的積分項(xiàng)中不再體現(xiàn)數(shù)據(jù)塊ΔgRp的作用。需要指出的是，這樣的處理結(jié)果是在把數(shù)據(jù)塊ΔgRp當(dāng)成常值條件下獲得的，如果該條件不滿足，那么還應(yīng)對數(shù)據(jù)塊ΔgRp作單獨(dú)處理。假設(shè)與計(jì)算點(diǎn)重合的網(wǎng)格數(shù)據(jù)塊半徑為ψ00，考慮到該數(shù)據(jù)塊是一個很小的區(qū)域，故可采用與上一小節(jié)相同的方式對由式(2)表示的積分核函數(shù)作平面近似處理，此時，計(jì)算點(diǎn)所在數(shù)據(jù)塊的積分式可寫為

sdsdα

(21)

式中，Δgp0代表計(jì)算點(diǎn)所在數(shù)據(jù)塊對計(jì)算參量Δgp的影響；s0代表數(shù)據(jù)網(wǎng)格大小的一半，當(dāng)數(shù)據(jù)網(wǎng)格為1′×1′時，s0=0.5′。由式(21)得知，如果把數(shù)據(jù)塊ΔgRp當(dāng)成常值看待，即認(rèn)為在計(jì)算點(diǎn)所在的數(shù)據(jù)網(wǎng)格內(nèi)處處滿足Δgq=ΔgRp，則有Δgp0=0。當(dāng)計(jì)算點(diǎn)附近的重力異常場變化比較劇烈時，上述處理方法會帶來一定的計(jì)算誤差，此時，可參照文獻(xiàn)[1]的思路，將重力異常Δgq在空間計(jì)算點(diǎn)P的球面投影點(diǎn)Rp處展開為泰勒級數(shù)

(22)

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

將式(22)代入(21)，可得

(23)

假設(shè)計(jì)算點(diǎn)所在的數(shù)據(jù)格網(wǎng)為(i,j)，則有

(24)

(25)

將式(23)代入式(20)的右端，可得到計(jì)算外部重力異常的第4步也是最終的嚴(yán)密改化公式

Δgq(σ-σ0)+Δgp(σ-σ0)+Δgp0

(26)

相比傳統(tǒng)的向上延拓計(jì)算式(9)，式(26)增加了Δgp(σ-σ0)和Δgp0兩個補(bǔ)償量。

1.3 嚴(yán)密改化公式在向下延拓迭代計(jì)算中的應(yīng)用

如前所述，位場向下延拓實(shí)為同一位場向上延拓的反問題，由于向上延拓解算具有穩(wěn)定可靠的優(yōu)良特性，因此在重力場向下延拓迭代計(jì)算中具有多方面的應(yīng)用：①通過向下和向上延拓迭代計(jì)算將航空重力測量成果向下延拓到某個基準(zhǔn)面[10,22,25-26]；②將地面觀測重力異常向下延拓到地球內(nèi)部的某個球面，然后采用球面延拓公式完成向上延拓計(jì)算[10,17]；③將地面或航空重力測量成果統(tǒng)一向下延拓到大地水準(zhǔn)面，用于解算大地測量邊值問題[1,27]；④將地面或航空重力測量成果向下延拓到地球內(nèi)部更接近于位場源體的高度面，以提高位場數(shù)據(jù)解釋推斷的準(zhǔn)確性和可靠性[15,28-29]。在這些應(yīng)用中，重力異常向上延拓模型幾乎都是以同一個方式為各方所用，也就是在向下與向上的迭代計(jì)算中起到一個關(guān)鍵性的橋梁作用。以航空重力測量數(shù)據(jù)向下延拓計(jì)算為例，具體原理和流程說明如下。

假設(shè)在高度為h的等高面上完成了航空重力測量，經(jīng)處理后得到一組重力異常網(wǎng)格值Δgp，現(xiàn)需要將Δgp向下延拓到半徑為R的球面上。由式(26)得知，現(xiàn)在的問題是已知向上延拓方程左端的球面外部觀測量Δgp，需要求解方程右端第一項(xiàng)和第二積分項(xiàng)內(nèi)的球面重力異常ΔgRp

(27)

由式(27)還不能直接求得球面值ΔgRp，因?yàn)槭?27)右端仍包含未知參量ΔgRp和Δgq。為了解決此問題，一般采用迭代計(jì)算方法對式(27)進(jìn)行求解。首先取迭代計(jì)算初值

(28)

(29)

(30)

如此繼續(xù)進(jìn)行，按式(31)計(jì)算ΔgRp第k次近似值

(31)

直到前后兩次計(jì)算值的互差絕對值或互差均方根值小于某個設(shè)定的大于零的限差值ε或σ(如可取ε=0.3 mGal或σ=0.1 mGal)為止。

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

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

為了分析比較前面提出的重力異常向上延拓積分模型改化效果，本文采用超高階位模型EGM2008作為數(shù)值計(jì)算檢驗(yàn)的仿真標(biāo)準(zhǔn)場[31]，用于模擬產(chǎn)生球面及其外部設(shè)定高度的1′×1′網(wǎng)格重力異常理論“真值”(這里使用1′×1′而非5′×5′網(wǎng)格數(shù)據(jù)是為了減弱積分離散化誤差的影響)。由地球位模型計(jì)算球面及其外部設(shè)定高度面重力異常的公式為[1,4]

(32)

式中，各個符號的意義同前。

為了體現(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作為計(jì)算檢驗(yàn)的標(biāo)準(zhǔn)場，即取L=2160，進(jìn)而選取ri=R+hi，R=6371 km，使用EGM2008模型(361～2160階次)分別計(jì)算對應(yīng)于9個高度面上的1′×1′網(wǎng)格重力異常理論“真值”Δgti(i=1,2，…，9),每一個高度面對應(yīng)360×360=129 600個網(wǎng)格點(diǎn)數(shù)據(jù)，9個高度分別為：hi=0、0.1、0.3、1、3、5、10、30、50 km。表1列出了其中的5個高度面上的重力異常理論“真值”統(tǒng)計(jì)結(jié)果，圖1、圖2分別給出了對應(yīng)于0 km高度面和3 km高度面上的理論“真值”分布態(tài)勢圖。表1為由EGM2008模型(361～2160階次)計(jì)算得到的重力異常統(tǒng)計(jì)結(jié)果。

圖1 0 km高度面重力異常分布Fig.1 The gravity anomalies on 0 km altitude surface

圖2 3 km高度面重力異常分布Fig.2 The gravity anomalies on 3 km altitude surface

表1 由EGM2008模型(361～2160階次)計(jì)算得到的重力異常統(tǒng)計(jì)結(jié)果

表1統(tǒng)計(jì)結(jié)果和圖1、圖2顯示的重力異常變化形態(tài)說明，盡管已經(jīng)扣除掉2～360階次頻段的參考場，本試驗(yàn)區(qū)域重力異常場變化的激烈程度仍然非常顯著，足以代表真實(shí)地球絕大多數(shù)局部重力場的變化特征。

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

為了對比分析不同改化模型的計(jì)算效果，這里采用零高度面也就是球面上的1′×1′網(wǎng)格重力異?！罢嬷怠宝t作為觀測量，同時使用4種改化模型對前面選定的對應(yīng)于試驗(yàn)區(qū)9個高度面上的1′×1′網(wǎng)格重力異常進(jìn)行了計(jì)算分析。其中，第1模型是指直接使用式(1)作為基礎(chǔ)計(jì)算模型，并對全球積分域作了分區(qū)處理，但在實(shí)施近區(qū)計(jì)算時，直接扣除掉與計(jì)算點(diǎn)重合的1′×1′數(shù)據(jù)塊，以避免出現(xiàn)奇異積分；第2模型對應(yīng)于式(9)；第3模型對應(yīng)于式(20)；第4模型對應(yīng)于式(26)。將4種模型的計(jì)算值分別與相對應(yīng)的理論“真值”Δgti做比較，可獲得不同改化模型的計(jì)算精度評估信息，具體比對結(jié)果列于表2。這里積分半徑統(tǒng)一取為ψ0=2°，為了減小積分邊緣效應(yīng)對評估結(jié)果的影響，表2只列出試驗(yàn)中心區(qū)2°×2°方塊內(nèi)的數(shù)據(jù)比對結(jié)果(下同)。為了定量評估重力異常遠(yuǎn)區(qū)效應(yīng)大小及由全球積分過渡到局域積分引起的模型誤差影響，表3給出了采用式(11)計(jì)算得到的3組對應(yīng)于積分半徑ψ0=2°、ψ0=5°和ψ0=10°的遠(yuǎn)區(qū)效應(yīng)貢獻(xiàn)量Δgq(σ-σ0)統(tǒng)計(jì)結(jié)果，表4給出了采用式(17)計(jì)算得到的相對應(yīng)的3組誤差補(bǔ)償量Δgp(σ-σ0)統(tǒng)計(jì)結(jié)果。表5為依據(jù)式(18)計(jì)算得到的比例系數(shù)k(r,ψ0,N)隨位模型參考場階數(shù)N、積分半徑ψ0和高度h取值變化而變化的數(shù)值結(jié)果。

表2 由不同改化模型計(jì)算得到的9個高度面重力異常與“真值”比較Tab.2 Comparisons between the gravity anomalies from different modified models and the “true values” on 9 altitude surfaces mGal

表3 遠(yuǎn)區(qū)效應(yīng)貢獻(xiàn)量Δgq(σ-σ0)結(jié)果統(tǒng)計(jì)Tab.3 Statistics of the computational results of far-zone contribution Δgq(σ-σ0) mGal

表4 模型誤差補(bǔ)償量Δgp(σ-σ0)計(jì)算結(jié)果統(tǒng)計(jì)Tab.4 Statistics of the computational results of error compensationΔgp(σ-σ0) mGal

表5 比例系數(shù)k(r,ψ0,N)隨位場階數(shù)N、積分半徑ψ0和高度h變化情況Tab.5 The scale factor k(r,ψ0,N) with the changes of degree N, integral radius ψ0 and altitude h

由表2比對結(jié)果可以看出，除了第3和第4模型之間的計(jì)算結(jié)果差異較小以外，其他3個模型的計(jì)算精度差異明顯。第1模型一方面由于直接扣除了計(jì)算點(diǎn)所在數(shù)據(jù)塊的影響，另一方面在邊界面附近，受數(shù)值積分不連續(xù)性影響嚴(yán)重，故在3 km以下高度段，該模型的計(jì)算誤差起伏較大，從近30 mGal減小到不到2 mGal,變化幅度超過一個數(shù)量級，可見該模型不能直接用于這個高度段的重力異常向上延拓計(jì)算；第2模型盡管從理論上消除了核函數(shù)奇異性影響，并在超低空高度段取得了比第1模型好得多的計(jì)算精度，但由于該模型的改化過程存在不可忽略的缺陷，故即使到了10 km高度，該模型的計(jì)算誤差仍然超過不可接受的2 mGal；第3模型從理論上彌補(bǔ)了第2模型的缺陷，使得該模型的計(jì)算精度得到顯著改善，即使在零高度面，該模型計(jì)算值與比對基準(zhǔn)“真值”的最大互差也不超過1 mGal，互差均方根值不超過0.3 mGal。這個結(jié)果說明，對第2模型所作的補(bǔ)償改化處理是正確有效的，對于高精度要求的向上延拓計(jì)算應(yīng)用，更是必不可少。第4模型是對第3模型的進(jìn)一步改進(jìn)和優(yōu)化，理論上更加嚴(yán)密，其計(jì)算精度的改善幅度取決于采用的數(shù)據(jù)網(wǎng)格間距大小和計(jì)算點(diǎn)周圍重力異常場變化的劇烈程度。理論上，采用的數(shù)據(jù)網(wǎng)格間距越大，計(jì)算點(diǎn)周圍重力異常場變化越劇烈，第4模型相對第3模型的精度改善效果會越顯著。本試驗(yàn)未能體現(xiàn)該模型相對于第3模型的顯著改善效果，是因?yàn)檫@里采用的是1′×1′網(wǎng)格數(shù)據(jù)，而且是來自2160階次的位模型EGM2008仿真數(shù)據(jù)，對應(yīng)的1′×1′網(wǎng)格范圍內(nèi)的重力異常變化幅度相對較小?？梢灶A(yù)見，如果是在重力異常場變化更為劇烈的大山區(qū)且使用更大網(wǎng)格(如2′×2′或5′×5′)的實(shí)測數(shù)據(jù)，那么第4模型定會顯現(xiàn)出符合預(yù)期的優(yōu)勢。

由表3和表4計(jì)算結(jié)果可看出，重力異常遠(yuǎn)區(qū)效應(yīng)貢獻(xiàn)量和計(jì)算模型誤差補(bǔ)償量的大小均在毫伽級以上且后者大于前者，兩者的影響都不可忽略。由表4和表5結(jié)果可進(jìn)一步看出，盡管第3模型對第2模型的補(bǔ)償量均隨位場階數(shù)N、積分半徑ψ0和計(jì)算高度h的增大而減小，但當(dāng)位場階數(shù)取為N=360時，即使積分半徑增大到ψ0=10°，計(jì)算高度h=10 km處的誤差補(bǔ)償量仍然超過1 mGal；而當(dāng)積分半徑取為ψ0=2°時，即使位場階數(shù)提高到N=2160，高度h=1 km處的誤差補(bǔ)償量也可能超過1 mGal。這樣的結(jié)果再次說明，本文對傳統(tǒng)向上延拓計(jì)算模型進(jìn)行精細(xì)改化是必要且有效的。

前面完成的數(shù)值計(jì)算檢驗(yàn)都是基于輸入數(shù)據(jù)無誤差假設(shè)條件，得到的比對評估結(jié)果只是計(jì)算模型自身完備程度的反映。為了考察數(shù)據(jù)觀測誤差對向上延拓解算結(jié)果的影響，本文在前述試驗(yàn)基礎(chǔ)上，進(jìn)一步開展有輸入數(shù)據(jù)誤差影響條件下的數(shù)值計(jì)算檢驗(yàn)。具體做法是，在前面作為觀測量的位模型剩余重力異?！罢嬷怠宝t中分別加入1 mGal和3 mGal的隨機(jī)噪聲，形成兩組新的模擬觀測量，然后按照前面相同的計(jì)算方案和流程，依次采用4種改化模型完成9個高度面上的1′×1′網(wǎng)格重力異常計(jì)算，最后將計(jì)算結(jié)果與相對應(yīng)高度的重力異常“真值”做比對評估(表6、表7)。

將表6和表7的計(jì)算結(jié)果與表2作對比可知，數(shù)據(jù)誤差對4種改化模型解算結(jié)果的影響規(guī)律是一致的，沒有因?yàn)槟Ｐ透幕问降牟煌a(chǎn)生實(shí)質(zhì)性差異，總體而言，數(shù)據(jù)誤差只對5 km高度以下的計(jì)算結(jié)果產(chǎn)生有限度的影響。模型誤差和數(shù)據(jù)誤差共同作用于向上延拓解算結(jié)果，完全符合隨機(jī)誤差作用傳播規(guī)律。不同的是，模型誤差影響在第1和第2模型中占主導(dǎo)地位，相反，在第3和第4模型中，占主導(dǎo)作用的則是數(shù)據(jù)誤差的影響。同時，相比較而言，在相同誤差影響條件下，第4模型的計(jì)算效果要略優(yōu)于第3模型。這些結(jié)果說明，即使采用嚴(yán)密的第3和第4改化模型進(jìn)行重力向上延拓計(jì)算，也要盡可能將數(shù)據(jù)觀測誤差控制在較低的水平。

表6 1 mGal誤差影響下不同改化模型計(jì)算得到的9個高度面重力異常與“真值”比較Tab.6 Comparisons between the gravity anomalies from different models with 1 mGal error and the “true values” on 9 altitude surfaces mGal

表7 3 mGal誤差影響下不同改化模型計(jì)算得到的9個高度面重力異常與“真值”比較Tab.7 Comparisons between the gravity anomalies from different models with 3 mGal error and the “true values” on 9 altitude surfaces mGal

2.3 改化模型應(yīng)用于向下延拓迭代計(jì)算效果檢驗(yàn)

為了進(jìn)一步檢驗(yàn)前述嚴(yán)密改化公式應(yīng)用于向下延拓迭代計(jì)算的實(shí)際效果，這里繼續(xù)采用前面定義的4種向上延拓改化模型作為向下延拓迭代逼近的計(jì)算式，將事先計(jì)算得到的3 km高度面上的1′×1′網(wǎng)格重力異常理論“真值”向下延拓到0 km高度面，具體迭代過程見式(28)—式(31)。將由4種改化模型用于向下延拓迭代計(jì)算獲得的重力異常網(wǎng)格值分別與0 km高度面上的理論“真值”做比較，可獲得相對應(yīng)改化模型的應(yīng)用效果評估信息，具體比對評估結(jié)果列于表8。為了比較，表8同時列出了輸入數(shù)據(jù)存在觀測噪聲情況下的迭代計(jì)算比對結(jié)果。這里的積分半徑同樣統(tǒng)一取為ψ0=2°，結(jié)束迭代計(jì)算的閾值取為σ=0.1 mGal。表8所列結(jié)果只限于中心區(qū)2°×2°方塊內(nèi)的比對數(shù)據(jù)。

表8比對結(jié)果全面反映了向上延拓模型完備性對向下延拓迭代解算精度和穩(wěn)定性的影響。首先，由表2檢核評估統(tǒng)計(jì)結(jié)果得知，第1模型在3 km高度面的計(jì)算精度為1.7 mGal，第2模型為4.0 mGal，第3和第4模型均為0.1 mGal。表8比對結(jié)果說明，向上延拓計(jì)算模型誤差越大，對向下延拓解算結(jié)果的影響也越大，需要迭代計(jì)算的次數(shù)也越多，模型誤差和數(shù)據(jù)誤差共同作用的結(jié)果也符合前面所作的理論分析預(yù)期。但觀察迭代計(jì)算收斂過程可知，盡管隨著迭代次數(shù)的增加，第1和第2模型的迭代解算結(jié)果也在緩慢逼近空中的觀測重力異常，但向下延拓結(jié)果并不趨近于0 km高度面上的理論計(jì)算“真值”，表8中對應(yīng)于第1和第2模型迭代結(jié)束后的比對結(jié)果并不是它們的最好結(jié)果。實(shí)際情況是，第1模型的最好比對結(jié)果出現(xiàn)在第2次迭代結(jié)束之后，第2模型出現(xiàn)在第1次迭代之后，具體比對結(jié)果列于表9。第1模型從第3次、第2模型從第2次迭代開始，隨著迭代次數(shù)的增加，其相對應(yīng)的向下延拓解算結(jié)果實(shí)際上是越來越偏離理論計(jì)算“真值”。顯然，兩個模型出現(xiàn)最好比對結(jié)果的迭代次數(shù)是無法預(yù)估的，與計(jì)算模型誤差大小有關(guān)，具有較大的不確定性。模型誤差較小的第3和第4模型則不存在這方面的問題。這樣的結(jié)果再次說明，向上延拓計(jì)算模型完備程度直接決定向下延拓迭代計(jì)算精度水平，同時影響迭代計(jì)算結(jié)果的穩(wěn)定性。

表8 有無誤差影響下使用不同向上延拓改化模型進(jìn)行向下延拓迭代計(jì)算結(jié)果與“真值”比較

表9 有無誤差影響下使用第1和第2向上延拓改化模型進(jìn)行向下延拓迭代計(jì)算結(jié)果與“真值”比較

3 結(jié) 論

針對位場轉(zhuǎn)換中經(jīng)常遇見的重力異常向上延拓計(jì)算模型改化問題，本文分析研究了此類全球積分模型改化的技術(shù)流程和適用條件，指出了傳統(tǒng)改化方法存在的理論缺陷，同時依據(jù)實(shí)測數(shù)據(jù)保障條件和積分恒等式適用條件要求，導(dǎo)出了重力異常向上延拓積分模型的分步改化公式，提出了補(bǔ)償傳統(tǒng)改化模型缺陷的修正公式。針對應(yīng)用廣泛的位場逆轉(zhuǎn)換需求，提出了將最終的向上延拓嚴(yán)密改化模型應(yīng)用于重力異常向下延拓迭代解算的方法和步驟。采用超高階地球位模型EGM2008作為比對標(biāo)準(zhǔn)位場，同時選擇在重力異常場變化比較劇烈的馬里亞納海溝區(qū)塊開展數(shù)值計(jì)算符合度檢驗(yàn)，分別對重力異常向上延拓分步改化模型的計(jì)算精度及在向下延拓迭代解算中的應(yīng)用效果進(jìn)行了檢核分析和評估，表明采用嚴(yán)密改化模型是必要和有效的，不僅可顯著提高重力異常向上延拓的計(jì)算精度，同時有利于提高向下延拓迭代解算過程的穩(wěn)定性和可靠性。