宗敬文，李厚樸，紀(jì) 兵，歐陽(yáng)永忠

1. 海軍工程大學(xué)導(dǎo)航工程教研室，湖北 武漢 430033； 2. 自然資源部海洋環(huán)境探測(cè)技術(shù)與應(yīng)用重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510300

海洋重力場(chǎng)數(shù)據(jù)的獲取手段有衛(wèi)星測(cè)高反演、船載及航空重力測(cè)量，隨著人造地球軌道衛(wèi)星技術(shù)的發(fā)展，衛(wèi)星測(cè)高技術(shù)促進(jìn)了大地測(cè)量學(xué)、地球物理學(xué)和海洋學(xué)等學(xué)科的交叉發(fā)展[1]，同時(shí)與傳統(tǒng)的海洋觀測(cè)手段相比，能夠周期性獲取除極地以外的高分辨率、全天候、長(zhǎng)時(shí)間序列的全球海洋觀測(cè)數(shù)據(jù)[2-5]。截至目前，世界上各個(gè)國(guó)家共計(jì)發(fā)射測(cè)高衛(wèi)星18顆，其中包括美國(guó)海軍的Geosat系列、美國(guó)宇航局和法國(guó)空間研究中心聯(lián)合發(fā)射的T/P系列[6]、歐空局的ERS系列[7]及中國(guó)的HY-2A[8]。衛(wèi)星測(cè)高數(shù)據(jù)遍布全球85°S—85°N的海洋區(qū)域，分辨率達(dá)1′×1′，精度達(dá)1～2 cm，利用高精度、高分辨率衛(wèi)星測(cè)高數(shù)據(jù)反演重力異常精度已達(dá)到10 mGal(1 Gal=10-2m/s2)以內(nèi)。文獻(xiàn)[9]使用逆Stokes法反演中國(guó)近海30′×30′海洋重力異常，精度達(dá)到3.5 mGal。文獻(xiàn)[10]使用Laplace方程的垂線偏差法反演全球海域的重力異常，反演的全球海域1′×1′重力場(chǎng)模型與NGDC船測(cè)數(shù)據(jù)的檢核精度達(dá)到4～8 mGal。文獻(xiàn)[11—14]在衛(wèi)星測(cè)高數(shù)據(jù)融合和精細(xì)處理方面也作出了大量的貢獻(xiàn)。

文獻(xiàn)[15]利用Stokes公式將衛(wèi)星測(cè)高數(shù)據(jù)轉(zhuǎn)換為重力異常。文獻(xiàn)[16]基于垂線偏差對(duì)地形信息敏感的特點(diǎn)，根據(jù)邊值理論由重力與地形數(shù)據(jù)確定格網(wǎng)垂線偏差模型。文獻(xiàn)[17]提出了利用垂線偏差來(lái)反演重力異常的方法，并系統(tǒng)地給出了垂線偏差、擾動(dòng)重力和重力異常之間的關(guān)系表達(dá)式，但在利用大地水準(zhǔn)面高和垂線偏差反演重力異常實(shí)際計(jì)算中，計(jì)算點(diǎn)及其附近區(qū)域到計(jì)算點(diǎn)的理論距離近似于零，會(huì)導(dǎo)致逆Stokes公式和逆Vening-Meinesz公式中的積分函數(shù)產(chǎn)生奇異。文獻(xiàn)[18—19]提出一組“非奇異變換”，系統(tǒng)地解決了物理大地測(cè)量中高程異常、垂線偏差、地形改正及重力異常梯度等泛函的中央?yún)^(qū)計(jì)算問(wèn)題。對(duì)于中央?yún)^(qū)積分問(wèn)題的研究，文獻(xiàn)[20—21]將該積分區(qū)域視為圓形域，推導(dǎo)的逆Vening-Meinesz公式被廣泛應(yīng)用于反演海洋重力異常。文獻(xiàn)[22]在利用逆Vening-Meinesz公式的FFT算法時(shí)，考慮了中央?yún)^(qū)效應(yīng)的影響，研究表明中央?yún)^(qū)效應(yīng)對(duì)重力異常的影響可以達(dá)到百微伽量級(jí)，在高精度重力場(chǎng)計(jì)算中是不能忽視的。文獻(xiàn)[23—24]分別推導(dǎo)了逆Stokes法和逆Vening-Meinesz法的衛(wèi)星測(cè)高反演重力異常解析計(jì)算公式，解決了奇異積分問(wèn)題，具有較高的計(jì)算精度，但解析計(jì)算公式系數(shù)繁多，計(jì)算量大，實(shí)際應(yīng)用中有所不便。因此，為進(jìn)一步完善衛(wèi)星測(cè)高反演重力異常的計(jì)算方法，簡(jiǎn)化計(jì)算過(guò)程，提高計(jì)算效率，本文采用數(shù)值求積公式，分別利用Simpson公式和Cotes公式對(duì)逆Stokes法和逆Vening-Meinesz法中的奇異積分問(wèn)題進(jìn)行了新的研究，系統(tǒng)地推導(dǎo)出了中央?yún)^(qū)重力異常的計(jì)算公式。此公式可直接利用格網(wǎng)節(jié)點(diǎn)處的大地水準(zhǔn)面高和垂線偏差計(jì)算重力異常值，形式簡(jiǎn)單，計(jì)算效率高，計(jì)算精度與解析法計(jì)算結(jié)果精度相當(dāng)，可以滿足實(shí)際應(yīng)用。

1 逆Stokes公式法奇異積分?jǐn)?shù)值求積公式

由文獻(xiàn)[23—24]可知，大地水準(zhǔn)面計(jì)算重力異常的逆Stokes公式平面近似形式為

(1)

(2)

假設(shè)中央?yún)^(qū)NQ可以展開(kāi)為如下Taylor級(jí)數(shù)

NQ=N(x,y)=NP+Nxx+Nyy+

(3)

由式(3)可得

(4)

設(shè)中央?yún)^(qū)積分面積為σ∈[-a≤x≤a,-b≤y≤b]，如圖1所示。為方便推導(dǎo)計(jì)算，引入新的積分變量u=x/a和v=y/b，因此中央?yún)^(qū)積分區(qū)域面積為σ′∈[-1≤u≤1,-1≤v≤1]。

(5)

圖1 中央?yún)^(qū)示意圖Fig.1 Sketch map of the innermost area

對(duì)σ1引入新的積分變量

(6)

對(duì)σ2引入新的積分變量

(7)

將式(6)和式(7)代入式(5)，可得中央?yún)^(qū)IN的表達(dá)式為

(8)

1.1 基于Simpson公式推導(dǎo)出的逆Stokes公式中央?yún)^(qū)數(shù)值求積公式

利用Simpson公式

(9)

考慮到式(4)，利用式(9)對(duì)式(8)進(jìn)行u和v方向上數(shù)值積分后可得

(10)

式中，INS表示基于Simpson公式的逆Stokes法奇異積分?jǐn)?shù)值求積公式。將N(u,v)在u,v=0,±1處的值記為Nij(i,j=0,±1)，并對(duì)式(10)在k和λ方向上繼續(xù)使用Simpson公式，其中

(11)

根據(jù)數(shù)值微分公式可由節(jié)點(diǎn)處大地水準(zhǔn)面高求出

(12)

整理后可得

(13)

將式(13)代入式(2)，最后可得中央?yún)^(qū)重力異常的數(shù)值求積公式為

Nab-4NP)

]

(14)

1.2 基于Cotes公式推導(dǎo)出的逆Stokes公式中央?yún)^(qū)數(shù)值求積公式

當(dāng)只考慮中央?yún)^(qū)范圍為σ∈[-a≤x≤a,-b≤y≤b]時(shí)，中央?yún)^(qū)的影響是不完整的，特別是在起伏變化較大的區(qū)域時(shí)，考慮將中央?yún)^(qū)的范圍擴(kuò)大為σ∈[-2a≤x≤2a,-2b≤y≤2b]，因此，本文引入了Cotes公式來(lái)計(jì)算。推導(dǎo)過(guò)程類(lèi)似于上一部分，首先利用Cotes公式

(15)

考慮到式(4)，利用式(15)對(duì)式(8)進(jìn)行u和v方向上數(shù)值積分后可得

(16)

式中，INC表示基于Cotes公式的逆Stokes法奇異積分求積公式。

如圖2所示，將N(u,v)在u,v=0,±1,±2處的值記為Nij(i,j=0,±1,±2)，對(duì)式(16)在k和λ方向上繼續(xù)使用Cotes公式，并考慮到式(11)和式(12)，整理后可得

Na-2b+N-a2b+N-a-2b-8NP)

}

(17)

將式(17)代入式(2)中，最后可得中央?yún)^(qū)重力異常的數(shù)值求積公式為

N-2a-2b+N-2a2b+N2a-2b-4NP)+

N-2a-b+Na2b+Na-2b+N-a2b+

N-a-2b-8NP)

}

(18)

圖2 Cotes公式中央?yún)^(qū)示意圖Fig.2 Sketch map of the innermost area used the Cotes formula

利用式(14)和式(18)計(jì)算中央?yún)^(qū)重力異常時(shí)，可直接使用已知網(wǎng)格節(jié)點(diǎn)處的大地水準(zhǔn)面高作為輸入值，相較于利用解析法的求解過(guò)程減少了求解多項(xiàng)式系數(shù)的步驟，大大簡(jiǎn)化了逆Stokes法反演重力異常的過(guò)程，從而提高了計(jì)算效率。

2 逆Vening-Meinesz公式法奇異積分?jǐn)?shù)值求積公式

由文獻(xiàn)[25—26]可知，垂線偏差計(jì)算重力異常的逆Vening-Meinesz公式平面近似形式為

(19)

式中，ξQ、ηQ為移動(dòng)點(diǎn)處的垂線偏差；ξP、ηP為計(jì)算點(diǎn)處的垂線偏差。為便于計(jì)算，令

(20)

將垂線偏差ξQ、ηQ在中央?yún)^(qū)展成如下Taylor級(jí)數(shù)

(21)

由式(21)得出

(22)

由于積分區(qū)域的對(duì)稱(chēng)性，并考慮到式(6)和式(7)，式(20)可改寫(xiě)為

(23)

2.1 基于Simpson公式推導(dǎo)出的逆Vening-Meinesz公式中央?yún)^(qū)數(shù)值求積公式

利用式(9)，并考慮到式(22)，對(duì)式(23)中積分部分進(jìn)行u和v方向上數(shù)值積分后可得

(24)

將ξ(u,v)、η(u,v)在u、v=0、±1處的值記為ξij、ηij(i,j=0、±1)，并考慮到根據(jù)數(shù)值微分公式可由節(jié)點(diǎn)處垂線偏差求出

(25)

最后對(duì)式(24)在k和λ方向上繼續(xù)使用Simpson公式，整理后得

(26)

式中，ΔgξS、ΔgηS表示基于Simpson公式的逆Vening-Meinesz法奇異積分?jǐn)?shù)值求積公式。因此，式(26)即為基于Simpson公式的逆Vening-Meinesz法反演中央?yún)^(qū)重力異常的數(shù)值求積公式。

2.2 基于Cotes公式推導(dǎo)出的逆Vening-Meinesz公式中央?yún)^(qū)數(shù)值求積公式

同樣，當(dāng)只考慮中央?yún)^(qū)范圍為σ∈[-a≤x≤a,-b≤y≤b]時(shí)，中央?yún)^(qū)的影響是不完整的，特別是在起伏變化較大的區(qū)域時(shí)，考慮將中央?yún)^(qū)的范圍擴(kuò)大為σ∈[-2a≤x≤2a,-2b≤y≤2b]，因此，本文引入了Cotes公式來(lái)計(jì)算。利用式(15)并考慮到式(22)，對(duì)式(23)中積分部分進(jìn)行u和v方向上數(shù)值積分后可得

(27)

如圖2所示，將ξ(u,v)、η(u,v)在u,v=0,±1,±2處的值記為ξij、ηij(i,j=0,±1,±2)。對(duì)式(27)在k和λ方向上繼續(xù)使用Cotes公式，并考慮到式(11)和式(25)，整理后可得

(28)

式中，ΔgξC、ΔgηC表示基于Cotes公式的逆Vening-Meinesz法奇異積分?jǐn)?shù)值求積公式，因此，式(28)即為基于Cotes公式的逆Vening-Meinesz法反演中央?yún)^(qū)重力異常數(shù)值求積公式。利用式(26)和式(28)計(jì)算中央?yún)^(qū)重力異常時(shí)，可直接使用已知網(wǎng)格節(jié)點(diǎn)處的垂線偏差作為輸入值，相較于利用解析法的求解過(guò)程減少了求解插值多項(xiàng)式系數(shù)的步驟，從而簡(jiǎn)化了逆Vening-Meinesz公式反演重力異常的計(jì)算過(guò)程，提高了計(jì)算效率。

3 算例分析

3.1 理論模型檢驗(yàn)

為驗(yàn)證本文推導(dǎo)的基于數(shù)值求積公式的測(cè)高反演重力異常普適計(jì)算公式的準(zhǔn)確性和可靠性，假定大地水準(zhǔn)面高和垂線偏差為理論模型值，計(jì)算了非奇異變換前、后求積公式計(jì)算結(jié)果，分析比較了非奇異變換后兩種不同普適數(shù)值求積公式的計(jì)算精度?，F(xiàn)設(shè)計(jì)如下兩個(gè)算例。

3.1.1 理論模型1

取中央?yún)^(qū)大地水準(zhǔn)面高的數(shù)學(xué)模型為

(29)

中央?yún)^(qū)垂線偏差的數(shù)學(xué)模型為

(30)

根據(jù)本文推導(dǎo)出的普適數(shù)值求積公式取a=1,b=1，因此中央?yún)^(qū)范圍為σ∈[-2≤x≤2,-2≤y≤2]，δ為平滑因子，則NP=NQ(0,0)=δ，ξP=ξQ(0,0)=0，ηP=ηQ(0,0)=0。以δ取10時(shí)為例，大地水準(zhǔn)面高如圖3所示。

圖3 大地水準(zhǔn)面高Fig.3 Geoidal height

(31)

為驗(yàn)證非奇異變換前、后積分表達(dá)式精度，分別取δ=1、5、10時(shí)，IN1和IN2在不同等分區(qū)間下的計(jì)算結(jié)果列于表1。

由表1可知，非奇異變換后IN250×50等分計(jì)算結(jié)果的準(zhǔn)確度和非奇異變換前IN1500×500等分計(jì)算結(jié)果的準(zhǔn)確度基本相同。由圖4(a)可以看出，非奇異變換前在低等分區(qū)間精度較低且收斂速度慢，非奇異變換后的收斂速度要大于變換前且計(jì)算結(jié)果更接近真值。因此，非奇異變換后的積分公式計(jì)算效率更高，收斂速度更快且精度更高。為進(jìn)一步分析比較基于Simpson公式和Cotes公式的逆Stokes法反演中央?yún)^(qū)重力異常普適數(shù)值計(jì)算公式，利用大地水準(zhǔn)面高理論模型1，取IN24000×4000等分的計(jì)算結(jié)果為準(zhǔn)確值，式(13)和式(17)的計(jì)算結(jié)果分別記為INS和INC，將δ分別取δ=3,5,8,10時(shí)的計(jì)算結(jié)果列于表2。

表1 非奇異變換前后奇異積分計(jì)算結(jié)果

圖4 δ=10時(shí)非奇異變換前后不同等分區(qū)間積分?jǐn)?shù)值Fig.4 Results of singular integral before and after the non-singular transformation

由表2可知，當(dāng)δ=3，5時(shí)，Cotes公式法的相對(duì)誤差小于Simpson公式法；當(dāng)δ=8，10時(shí)Simpson公式法的相對(duì)誤差要小于Cotes公式法。Simpson公式法數(shù)值計(jì)算公式的相對(duì)誤差波動(dòng)較大，與平滑因子δ取值有關(guān)，大地水準(zhǔn)面高變化緩和地區(qū)相對(duì)誤差較小，精度較高；Cotes公式法數(shù)值計(jì)算公式的相對(duì)誤差比較穩(wěn)定，隨著平滑因子δ取值變化，精度受影響較小，但精度稍低。相比較于Cotes公式法，Simpson公式法在平緩地區(qū)具有更好的計(jì)算精度，但其計(jì)算精度波動(dòng)較大；相比較于Simpson公式法，Cotes公式法具有更好的計(jì)算穩(wěn)定性，計(jì)算精度受地形影響較小。

(32)

為驗(yàn)證非奇異變換前、后積分表達(dá)式精度，分別取δ=1、5、10，IV1和IV2在不同等分區(qū)間下的計(jì)算結(jié)果列于表3。

表2 非奇異變換后數(shù)值解計(jì)算結(jié)果誤差

表3 非奇異變換前后奇異積分計(jì)算結(jié)果

表4 非奇異變換后數(shù)值解計(jì)算結(jié)果誤差

由表4可知，Cotes公式法數(shù)值計(jì)算公式的相對(duì)誤差明顯小于Simpson公式法數(shù)值計(jì)算公式，Simpson公式法數(shù)值計(jì)算公式的相對(duì)誤差波動(dòng)較大，與平滑因子δ取值有關(guān)；Cotes公式法數(shù)值計(jì)算公式的相對(duì)誤差比較穩(wěn)定，當(dāng)δ分別取δ=3,6,8,10時(shí)，相對(duì)誤差均小于1%，相比較于Simpson公式法數(shù)值計(jì)算公式，Cotes公式法數(shù)值計(jì)算公式具有較好的計(jì)算精度和穩(wěn)定性。

文獻(xiàn)[23]和文獻(xiàn)[24]分別推導(dǎo)了雙二次多項(xiàng)式插值表示下的大地水準(zhǔn)面高和垂線偏差非奇異積分解析計(jì)算公式。記錄其結(jié)果分別為IN3和IV3

(α02+α20))

(33)

(34)

式中,IN3為雙二次多項(xiàng)式插值表示下的大地水準(zhǔn)面高非奇異積分解析表達(dá)式；IV3為雙二次多項(xiàng)式插值表示下的垂線偏差非奇異積分解析表達(dá)式。以非奇異變換后4000×4000等分區(qū)間計(jì)算結(jié)果作為準(zhǔn)確值，其中逆Stokes法準(zhǔn)確值記為IN2，逆Vening-Meinesz法準(zhǔn)確值記為IV2，式(13)、式(17)、式(26)和式(28)的計(jì)算結(jié)果分別為INS、INC、IVS和IVC，計(jì)算結(jié)果及其相應(yīng)計(jì)算誤差見(jiàn)表5。

表5 非奇異變換后解析解與數(shù)值解的計(jì)算誤差

分析表5可知，逆Stokes法重力異常求積公式解析解與準(zhǔn)確值的相對(duì)誤差為0.64%，Simpson公式數(shù)值解與準(zhǔn)確值的相對(duì)誤差為0.15%，Simpson公式數(shù)值解精度高于解析解計(jì)算精度，二者相對(duì)誤差均小于1%，而Cotes公式數(shù)值解計(jì)算精度較低，相對(duì)誤差為1.62%；逆Vening-Meinesz法重力異常求積公式解析解與準(zhǔn)確值的相對(duì)誤差為0.62%，Cotes公式數(shù)值解與準(zhǔn)確值的相對(duì)誤差為0.46%，Simpson公式數(shù)值解與準(zhǔn)確值的相對(duì)誤差為0.50%，解析解計(jì)算精度略低于數(shù)值解計(jì)算精度，三者相對(duì)誤差均小于1%。

3.1.2 理論模型2

在理論模型1的分析基礎(chǔ)上，理論模型2僅用于分析兩種不同數(shù)值積分方法下測(cè)高反演重力異常普適計(jì)算公式的準(zhǔn)確性和可靠性。取中央?yún)^(qū)大地水準(zhǔn)面高的數(shù)學(xué)模型為

(35)

中央?yún)^(qū)垂線偏差的數(shù)學(xué)模型為

(36)

現(xiàn)取a=1,b=0.7，因此中央?yún)^(qū)范圍為σ∈[-2≤x≤2,-1.4≤y≤1.4]，δ為平滑因子，則NP=NQ(0,0)=δ，ξP=ξQ(0,0)=0，ηP=ηQ(0,0)=0。以δ取10時(shí)為例，大地水準(zhǔn)面高如圖5所示。

圖5 大地水準(zhǔn)面高Fig.5 Geoid height

(37)

(38)

為分析比較基于Simpson公式法和Cotes公式法的逆Stokes法反演中央?yún)^(qū)重力異常普適數(shù)值計(jì)算公式，利用大地水準(zhǔn)面高理論模型2，式(13)和式(17)的計(jì)算結(jié)果分別記為INS和INC，將δ分別取δ=3，5，8，10時(shí)的計(jì)算結(jié)果列于表6。

表6 非奇異變換后數(shù)值解計(jì)算結(jié)果誤差

由表6可知，當(dāng)δ=3時(shí)，Cotes公式法的相對(duì)誤差小于Simpson公式法；當(dāng)δ=5、8、10時(shí)，Simpson公式法的相對(duì)誤差小于Cotes公式法。Simpson公式法數(shù)值計(jì)算公式的相對(duì)誤差波動(dòng)較大，當(dāng)平滑因子變大時(shí)，計(jì)算精度較高，相對(duì)誤差小于0.5%；Cotes公式法數(shù)值計(jì)算公式的相對(duì)誤差比較穩(wěn)定，隨著平滑因子δ取值變化，精度受影響較小，但精度稍低。相比較于Cotes公式法，Simpson公式法在平緩地區(qū)具有更好的計(jì)算精度，但其計(jì)算精度波動(dòng)較大；相比較于Simpson公式法，Cotes公式法具有更好的計(jì)算穩(wěn)定性，計(jì)算精度受地形影響較小。

利用垂線偏差理論模型2，式(26)和式(28)的計(jì)算結(jié)果分別為IVS和IVC，將δ分別取δ=3,5,8,10時(shí)的計(jì)算結(jié)果列于表7。

表7 非奇異變換后數(shù)值解計(jì)算結(jié)果誤差

由表7可知，Cotes公式法數(shù)值計(jì)算公式的相對(duì)誤差比較穩(wěn)定，當(dāng)δ分別取δ=3,5,8,10時(shí)，計(jì)算結(jié)果的相對(duì)誤差均小于0.5%；Simpson公式法數(shù)值計(jì)算公式的相對(duì)誤差波動(dòng)較大，且當(dāng)δ=3,5時(shí)，相對(duì)誤差較大，當(dāng)δ=8,10時(shí)，相對(duì)誤差小于0.5%。對(duì)比發(fā)現(xiàn)，Cotes公式法數(shù)值計(jì)算公式相比較于Simpson公式法數(shù)值計(jì)算公式具有較好的計(jì)算精度和計(jì)算穩(wěn)定性。

將式(33)和式(34)中解析解計(jì)算結(jié)果與本文推導(dǎo)公式進(jìn)行比較，取δ=10時(shí)，以式(37)、式(38)非奇異變換后4000×4000等分區(qū)間計(jì)算結(jié)果作為準(zhǔn)確值，其中逆Stokes法準(zhǔn)確值記為IN22，逆Vening-Meinesz法準(zhǔn)確值記為IV22，IN3、INS、INC、IV3、IVS和IVC計(jì)算結(jié)果及其相應(yīng)計(jì)算誤差見(jiàn)表8。

表8 非奇異變換后解析解與數(shù)值解的計(jì)算誤差

分析表8可知，其中逆Stokes法重力異常數(shù)值求積公式解析解與準(zhǔn)確值的相對(duì)誤差為0.57%，Simpson公式法數(shù)值解與準(zhǔn)確值的相對(duì)誤差為0.15%，Simpson公式法數(shù)值解精度高于解析解計(jì)算精度，二者相對(duì)誤差均小于0.6%，而Cotes公式數(shù)值解計(jì)算精度較低，相對(duì)誤差為2.20%；逆Vening-Meinesz法重力異常數(shù)值求積公式解析解與準(zhǔn)確值的相對(duì)誤差為0.44%，Cotes公式法數(shù)值解與準(zhǔn)確值的相對(duì)誤差為0.02%，Simpson公式法數(shù)值解與準(zhǔn)確值的相對(duì)誤差為0.31%，數(shù)值解計(jì)算精度略高于解析解計(jì)算精度。

本文通過(guò)利用兩種不同理論模型，分析比較基于Simpson公式法和Cotes公式法兩種不同數(shù)值求積方法推導(dǎo)出的測(cè)高反演重力異常普適計(jì)算公式，發(fā)現(xiàn)在逆Stokes法重力異常求積公式計(jì)算中，Simpson公式法計(jì)算精度較好，但受地形影響其誤差波動(dòng)較大，適用于中央?yún)^(qū)變化比較平緩地區(qū)，Cotes公式法計(jì)算誤差穩(wěn)定性較高，但計(jì)算精度略低；在逆Vening-Meinesz法重力異常求積公式中，Cotes公式法計(jì)算精度和計(jì)算誤差穩(wěn)定性都優(yōu)于Simpson公式法，且計(jì)算精度優(yōu)于解析法計(jì)算精度，可以滿足應(yīng)用要求。

3.2 試驗(yàn)數(shù)據(jù)檢驗(yàn)

理論模型下的精度檢驗(yàn)表明，本文推導(dǎo)出的中央?yún)^(qū)奇異積分?jǐn)?shù)值求積公式有較高的計(jì)算精度且形式簡(jiǎn)單。為進(jìn)一步說(shuō)明中央?yún)^(qū)奇異積分?jǐn)?shù)值求積公式的計(jì)算精度及中央?yún)^(qū)對(duì)重力異常的影響，選定南海16°N-20°N，112°E-116°E海域作為計(jì)算區(qū)域，以EGM2008地球重力場(chǎng)模型計(jì)算得到的2′×2′分辨率的大地水準(zhǔn)面高數(shù)據(jù)作為試驗(yàn)數(shù)據(jù)，計(jì)算海域大地水準(zhǔn)面高分布圖如圖6所示，計(jì)算海域共記119×119個(gè)格網(wǎng)，中央?yún)^(qū)為8′×8′(4×4個(gè)格網(wǎng))，采用逆Stokes法利用式(33)、式(14)和式(18)實(shí)際計(jì)算了該區(qū)域的中央?yún)^(qū)重力異常，結(jié)果分別記為Δg1、Δg2、Δg3列于表9，計(jì)算結(jié)果之間的比較情況見(jiàn)表10。

由表9和表10可以看出，3種計(jì)算方法計(jì)算出的4×4個(gè)格網(wǎng)大小的中央?yún)^(qū)重力異常平均值分別為-1.088、-1.089、-1.087 mGal，三者差異較小；標(biāo)準(zhǔn)差分別為2.60、2.63、2.83 mGal；Simpson公式法與解析法計(jì)算得到的中央?yún)^(qū)重力異常差值的標(biāo)準(zhǔn)差為0.030 mGal，最大值為0.195 mGal，最小值為-0.129 mGal,因此當(dāng)考慮中央?yún)^(qū)范圍擴(kuò)大為4×4個(gè)格網(wǎng)時(shí)，Simpson公式法和解析法計(jì)算結(jié)果相當(dāng)；Cotes公式法計(jì)算出的重力異常最大值為14.39 mGal，最小值為-13.23 mGal，與Simpson公式法和解析法計(jì)算得到的中央?yún)^(qū)重力異常差值的標(biāo)準(zhǔn)差分別為0.533 mGal和0.545 mGal，最大值分別為3.806 mGal和4.254 mGal，最小值分別為-4.167 mGal和-3.879 mGal,Cotes公式法與另外兩種方法計(jì)算結(jié)果具有一定差距，這是由于Cotes公式法計(jì)算過(guò)程中計(jì)算點(diǎn)更多，反映出利用該方法計(jì)算中央?yún)^(qū)重力異常變化更為明顯,可以更好地反映出起伏較大區(qū)域中央?yún)^(qū)重力異常變化。同時(shí)，試驗(yàn)數(shù)據(jù)說(shuō)明，分別利用3種不同方法計(jì)算出的4×4個(gè)格網(wǎng)中央?yún)^(qū)重力異常范圍為-13.5 mGal～+14.5 mGal，因此在高精度重力異常反演時(shí)，中央?yún)^(qū)有不容忽視的貢獻(xiàn)。

圖6 計(jì)算海域大地水準(zhǔn)面高示意圖Fig.6 Sketch map of geoidal height in the caculating sea area

表9 逆Stokes法反演中央?yún)^(qū)重力異常統(tǒng)計(jì)情況

表10 3種方法計(jì)算結(jié)果比較

4 結(jié) 論

為完善衛(wèi)星測(cè)高反演重力異常的計(jì)算理論，本文利用Simpson公式和Cotes公式的數(shù)值求積方法，系統(tǒng)地推導(dǎo)出了逆Stokes法和逆Vening-Meinesz法反演中央?yún)^(qū)重力異常普適數(shù)值積分公式，利用理論模型和試驗(yàn)數(shù)據(jù)進(jìn)行了分析檢驗(yàn)。研究表明：

(1) 非奇異變換后公式的計(jì)算效率要高于非奇異變換前，且計(jì)算精度更高。

(2) 通過(guò)將逆Stokes法和逆Vening-Meinesz法奇異積分解析法計(jì)算結(jié)果與數(shù)值計(jì)算結(jié)果相比較，發(fā)現(xiàn)數(shù)值法計(jì)算結(jié)果精度優(yōu)于解析法計(jì)算結(jié)果，相對(duì)誤差均小于1%，完全可以滿足實(shí)際應(yīng)用。

(3) 本文推導(dǎo)出的普適數(shù)值求積公式可以直接利用網(wǎng)格節(jié)點(diǎn)處的大地水準(zhǔn)面高和垂線偏差計(jì)算重力異常，形式簡(jiǎn)單，避免了傳統(tǒng)解析表達(dá)式系數(shù)繁多、形式復(fù)雜的缺點(diǎn)，在確保了公式計(jì)算精度的同時(shí)也提高了計(jì)算效率。

(4) 試驗(yàn)數(shù)據(jù)下的分析表明，在利用逆Stokes法反演重力異常時(shí)，解析法與數(shù)值法計(jì)算結(jié)果相當(dāng)，且中央?yún)^(qū)有不容忽視的貢獻(xiàn)。