馬志偉，陸　洋，涂　弋，朱傳東，郗　慧

1. 中國科學(xué)院測量與地球物理研究所大地測量與地球動力學(xué)國家重點實驗室，湖北 武漢 430077； 2. 中國科學(xué)院大學(xué)，北京 100049

?

利用Abel-Poisson徑向基函數(shù)模型化局部重力場

馬志偉1，2，陸洋1，2，涂弋1，2，朱傳東1，2，郗慧1，2

1. 中國科學(xué)院測量與地球物理研究所大地測量與地球動力學(xué)國家重點實驗室，湖北 武漢 430077； 2. 中國科學(xué)院大學(xué)，北京 100049

多種類型高分辨率重力場數(shù)據(jù)的不斷增加，使得在局部范圍內(nèi)精化重力場模型成為了可能。本文采用Abel-Poisson核將重力場量表示成有限個徑向基函數(shù)線性求和的形式，對局部區(qū)域的多種重力場數(shù)據(jù)進(jìn)行聯(lián)合建模。為了提高運算速度，運用了基于自適應(yīng)精化格網(wǎng)算法的最小均方根誤差準(zhǔn)則(RMS)來求解徑向基函數(shù)平均帶寬。以南海核心地區(qū)為例，聯(lián)合兩種不同類型、不同分辨率的重力場資料(大地水準(zhǔn)面起伏6′×6′、重力異常2′×2′)，構(gòu)建了局部區(qū)域高分辨率的重力場模型。所建模型表示的重力場參量達(dá)到了2′×2′的分辨率，對原始的重力異常數(shù)據(jù)(2′×2′)擬合的符合程度達(dá)到±0.8×10-5m/s2。結(jié)果表明，利用徑向基函數(shù)方法進(jìn)行局部重力場建模，避免了球諧函數(shù)建模收斂慢的問題，有效提高了模型表示重力場的分辨率。

局部重力場建模;重力異常;徑向基函數(shù);數(shù)據(jù)自適應(yīng);最小均方根誤差準(zhǔn)則

作為大地測量學(xué)一個重要的物理量，地球重力場對研究地球形狀及其內(nèi)部構(gòu)造具有重要的意義，特別是最近幾十年，隨著專題重力衛(wèi)星CHAMP、GRACE和GOCE的相繼投入使用和衛(wèi)星測高、重力儀觀測技術(shù)的不斷改善，人們對于求解局部高精度、高分辨率的重力場模型的需求日益迫切。因此，如何充分利用觀測數(shù)據(jù)進(jìn)行有效地反演地球重力場，成為近年來研究的一個重要課題：如文獻(xiàn)[1]研究了衛(wèi)星測高誤差對重力場反演的影響；文獻(xiàn)[2]利用兩種衛(wèi)星數(shù)據(jù)從譜組合法的角度對其展開了討論等。通常，對于重力場的表示是利用球諧函數(shù)將其展開至某些特定階次，例如EGM08展開至2190階，EIGEN-GL04S展開至150階。這種傳統(tǒng)的重力場表示方法有著應(yīng)用上的便利性，但由于其全局緊支撐特性，任何一個球諧系數(shù)的變化都會引起整個重力場的改變；另外，它也不能顧及局部數(shù)據(jù)分布的不均勻性，不適用于局部重力場的精化。由于其本身的局部緊支撐特性，徑向基函數(shù)能夠捕捉到細(xì)節(jié)的重力場信號，可應(yīng)用于局部重力場建模。

近些年來，眾多學(xué)者對于徑向基函數(shù)表示局部重力場作出了嘗試，并在基函數(shù)構(gòu)建以及基函數(shù)格網(wǎng)設(shè)計算法方面提出了各自的理論?；瘮?shù)的構(gòu)建大致可分為以下幾類：點質(zhì)量基函數(shù)法[3]、多級小波基函數(shù)法[4-5]、Blackman 小波基函數(shù)法[6-8]、球諧樣條核函數(shù)法[9]。點質(zhì)量基函數(shù)法在一定程度上可以減少基函數(shù)個數(shù)，但是該方法實施難度較大，且計算復(fù)雜度較其他方法并無改善；小波基函數(shù)和球諧樣條核函數(shù)法由于采用固定類型的基函數(shù)格網(wǎng)設(shè)計方案，為了捕捉細(xì)節(jié)信息，需要不斷細(xì)化基函數(shù)格網(wǎng)，進(jìn)而產(chǎn)生大量的基函數(shù)格網(wǎng)點數(shù)，這一方面對計算機(jī)性能提出了嚴(yán)峻的挑戰(zhàn)，另一方面不可避免地對某些區(qū)域造成過度參數(shù)化，從而加劇了法方程矩陣的病態(tài)性，最終達(dá)不到較好的建模效果。

在基函數(shù)格網(wǎng)設(shè)計算法上，由于牽涉基函數(shù)位置、數(shù)量、帶寬等諸多因素，目前還沒有一個統(tǒng)一的標(biāo)準(zhǔn)來模型化局部重力場。文獻(xiàn)[9]利用多種準(zhǔn)則(最小距離、均勻性、靈活性等)分析了多種格網(wǎng)的優(yōu)劣，結(jié)果表明，Reuter格網(wǎng)和三角頂點格網(wǎng)最適合應(yīng)用于模型化局部重力場。文獻(xiàn)[10]研究了不同的徑向基函數(shù)在不同類型的格網(wǎng)下表示重力場的差異和精度, 所得結(jié)論是：Abel-Poisson核、Blackman核和三次多項式核在空域內(nèi)具有較小的震蕩，更適合于局部重力場建模。文獻(xiàn)[11]指出徑向基函數(shù)的數(shù)量與研究區(qū)域內(nèi)數(shù)據(jù)的分辨率及其區(qū)域大小有關(guān)，并且觀測值的數(shù)量并不直接影響基函數(shù)的數(shù)量。文獻(xiàn)[8—9]將基函數(shù)格網(wǎng)與已知數(shù)據(jù)的球諧階次聯(lián)系起來，認(rèn)為所選取的基函數(shù)格網(wǎng)點數(shù)要大于同等階次的球諧函數(shù)系數(shù)個數(shù)。文獻(xiàn)[5]利用多級小波族形成的框架聯(lián)合模型化局部重力場，并且根據(jù)冗余度值(小波個數(shù)與產(chǎn)生的Hilbert空間的維度比值)來選定基函數(shù)的格網(wǎng)密度參數(shù)，這種方法可以提供一個相對穩(wěn)定的重力場表達(dá)，但是所需要的小波基函數(shù)數(shù)量仍然很大。文獻(xiàn)[6]提出了數(shù)據(jù)自適應(yīng)精化格網(wǎng)算法，在顯著減少基函數(shù)個數(shù)的同時，還可以靈活地選取與觀測信號相適應(yīng)的基函數(shù)位置和帶寬。

本文在文獻(xiàn)[6]自適應(yīng)算法的基礎(chǔ)上，運用最小均方根誤差準(zhǔn)則(RMS準(zhǔn)則)來求解平均帶寬，并采用方差分量估計法，對各個方差因子作出估計?；诖?，聯(lián)合南海地區(qū)兩種不同的重力場數(shù)據(jù)，構(gòu)建該區(qū)域高分辨率的徑向基函數(shù)地球重力場模型。

1　球面徑向基函數(shù)建模理論及其求解

球面徑向基函數(shù)(spherical radial basis functions，SRBF)，簡稱徑向基函數(shù)或基函數(shù)，是地球向徑方向上的局部化對稱函數(shù)，其函數(shù)值大小只與兩點間的球面角距有關(guān)，它的特點是具有很好的空間局部化特性，即大部分能量集中于基函數(shù)點附近的局部區(qū)域。圖1繪出了一個空間域正則化Abel-Poisson徑向基函數(shù)的基本形狀。

球面上任意一個位置向量x可以表示為

x=r[cosφcosλcosφsinλsinφ]T=rr

(1)

式中，φ為緯度；λ為經(jīng)度；r為向徑；r為x方向上的單位向量。

考慮兩個位置向量x、xi，其中x位于半徑為R(Bjerhammar球半徑)的球面σR上或外部，xi位于球面上,則極點位于xi處的外部徑向基函數(shù)可表示為式(2)

(2)

式中，Pn為n階勒讓德多項式Pn0；r、ri分別為x、xi方向上的單位向量；kn為基函數(shù)核，不同的基函數(shù)核會得到不同的徑向基函數(shù)，筆者取核函數(shù)為Abel-Poisson核[10,15]，其具體表達(dá)形式為

kn=e-np

(3)

p為Abel-Poisson核的帶寬參數(shù)，不同的帶寬參數(shù)值對應(yīng)不同的頻譜，進(jìn)而導(dǎo)致徑向基函數(shù)(SRBF)在空域內(nèi)不同的空間形狀。圖2和圖3給出了Abel-Poisson核在不同帶寬參數(shù)下的頻譜覆蓋情況及SRBF的空間形狀變化。從圖2 可以看出，帶寬參數(shù)p越大，頻譜覆蓋范圍越?。籶越小，頻譜覆蓋范圍越大。圖3中，帶寬參數(shù)p越大，SRBF越寬闊，能量越分散；p越小，SRBF越狹窄，能量越集中。

半徑為R的球面外調(diào)和函數(shù)T可以表示為N個徑向基函數(shù)的和的形式

(4)

式中，αi為徑向基函數(shù)系數(shù)；T(x)為擾動位；GM為地球引力系數(shù)；RE為地球平均半徑。

地球上同一區(qū)域可能有多種類型的重力觀測數(shù)據(jù)，比如重力異常Δg、大地水準(zhǔn)面起伏N(或者高程異常)，它們都可以表示為擾動位的線性函數(shù)[12,16]

(5)

(6)

與擾動位類似，其他重力觀測量同樣也可以表示為有限個徑向基函數(shù)疊加求和的形式，下面直接給出用徑向基函數(shù)表示的重力異常Δg、大地水準(zhǔn)面起伏N公式

(7)

(8)

式中，Γi(x,xi)、Θi(x,xi)為徑向基函數(shù)，它們與Ψi(x,xi)的線性關(guān)系可由式(5)和式(6)導(dǎo)出，其他符號與前面公式相同。

式(4)、式(7)和式(8)都具有相同的未知基函數(shù)系數(shù)αi(i=1,2,…，N)，因此，可以利用不同類型的觀測數(shù)據(jù)，共同求解未知系數(shù)。

不同類型的重力觀測數(shù)據(jù)可組成如式(9)的線性方程

(9)

(10)

2　數(shù)據(jù)自適應(yīng)精化格網(wǎng)算法(data-adaptive network design，DAND)

徑向基函數(shù)建模，首先需要設(shè)置基函數(shù)格網(wǎng)，基函數(shù)格網(wǎng)有多種選擇，如等經(jīng)緯度格網(wǎng)、等三角形格網(wǎng)或Reuter格網(wǎng)等。本文主要采用Reuter格網(wǎng)[18]。Reuter格網(wǎng)在全球范圍可以實現(xiàn)良好的均勻分布，并且格網(wǎng)點的生成比較靈活，Reuter格網(wǎng)的密度由分辨率水平L決定， 圖4(a)、(b)繪出了L=20和L=30時的全球Reuter格網(wǎng)分布情況。

從圖4可以看出，Reuter格網(wǎng)分布呈從赤道到兩級遞減的趨勢，不會因為靠近極點而逐漸變密；L越大，格網(wǎng)點間距越小，分布越密。但是，單一的基函數(shù)格網(wǎng)并不能取得良好的建模效果，當(dāng)觀測值分布較不均勻時，基函數(shù)格網(wǎng)與觀測值之間便不能很好地匹配，往往會導(dǎo)致一些地區(qū)過參數(shù)化，而另一些地區(qū)欠參數(shù)化情況的出現(xiàn)；另外，基函數(shù)的分布也必須依據(jù)重力信號變化的強(qiáng)弱作出調(diào)整，單一的Reuter格網(wǎng)并不能滿足這一要求。因此，有必要尋求與數(shù)據(jù)分布、信號變化相適應(yīng)的基函數(shù)格網(wǎng)設(shè)計算法。

針對上述問題，文獻(xiàn)[6]提出了依賴于觀測數(shù)據(jù)的自適應(yīng)精化格網(wǎng)算法，它可以明顯地改善以往單一格網(wǎng)算法欠參數(shù)化的情況，能夠與觀測信號之間進(jìn)行良好地匹配，并且可以聯(lián)合多種數(shù)據(jù)共同建模，在保證計算精度的基礎(chǔ)上，減少了過度擬合現(xiàn)象的發(fā)生。另外，為了提高計算速度，本文采用最小均方根誤差準(zhǔn)則(RMS)替代Klees的廣義交叉驗證法(GCV)來求解基函數(shù)平均帶寬，具體步驟如下。

(1) 定義球面粗格網(wǎng)，在格網(wǎng)點設(shè)置徑向基函數(shù)，用阻尼最小二乘法求解基函數(shù)系數(shù)，用最小均方根誤差準(zhǔn)則(RMS)求解基函數(shù)平均帶寬

(11)

RMS準(zhǔn)則的判斷標(biāo)準(zhǔn)如下。

(a) 根據(jù)數(shù)據(jù)頻譜信息，估計帶寬p的變化范圍，并指定適當(dāng)?shù)牟骄嘧兓?/p>

(c) 每一個帶寬對應(yīng)一個剩余值，其中最小RMS剩余所對應(yīng)的帶寬即為初始格網(wǎng)最佳帶寬。

在數(shù)據(jù)量較大的情況下，RMS準(zhǔn)則相對GCV明顯提高了收斂速度。表1給出了由EGM2008模型計算的南海地區(qū)不同階次的(300、400、500)重力異常數(shù)據(jù)利用這兩種準(zhǔn)則判斷基函數(shù)格網(wǎng)最佳帶寬的實例，結(jié)果顯示：兩者得到的最佳帶寬基本一致，但是在運算時間上RMS卻明顯少于GCV。需要說明的是，這里的改進(jìn)是指在判斷精化基函數(shù)是否應(yīng)該存在的準(zhǔn)則上，因為減少了計算的循環(huán)次數(shù)，在計算速度上有很大提高，但在精度方面，由于兩者篩選結(jié)果基本一致，所得精度也極為相近。

表1RMS和GCV確定最佳帶寬的運算時間和最佳p值的比較

Tab.1Comparisons of the time and optimum bandwidth of RMS and GCV

D/ON.obsRMSGCVOpt.ptime/sOpt.ptime/s300609 0.040～0.060690.040～0.055 176840018620.017～0.0191370.018～0.020395650036570.015～0.0189940.016～0.01813485

注：D/O代表數(shù)據(jù)階次；N.obs為觀測值個數(shù)；Opt.p為最佳帶寬。

(2) 求解(1)的最小二乘剩余，從最大到最小，依次判定是否其可被選為精化基函數(shù)，每選定一個基函數(shù)，在局部區(qū)域內(nèi)更新一次“觀測值剩余”，直至選不出滿足條件的徑向基函數(shù)為止。

(3) 將(1)中的粗格網(wǎng)基函數(shù)和(2)中的精化基函數(shù)結(jié)合起來，作為對該局部重力場的完全參數(shù)化，最后利用所有的觀測數(shù)據(jù)，用方差分量估計法(VCE)共同求解基函數(shù)系數(shù)。

精化基函數(shù)算法輸入?yún)?shù)對建模結(jié)果的好壞起決定性作用，判斷準(zhǔn)則主要有以下幾點。

(1) 待定SRBF的剩余值應(yīng)大于一定的閾值τ1；

(2) 待定SRBF的周圍(半徑為ψc的球蓋)應(yīng)至少存在q個足夠大的剩余觀測值，足夠大意味著q個點的平均值大于另一閾值τ2；

(3) 新增加的基函數(shù)點與已有的基函數(shù)點之間應(yīng)大于一定的距離(ψmin)。

τ1、τ2、ψc、q、ψmin通常在試驗中嘗試確定，一般的，τ1、τ2要足夠小，依據(jù)q和數(shù)據(jù)分辨率可確定ψc，ψmin不大于初始粗格網(wǎng)間距。

3　實例計算

選取南海地區(qū)約10°×8°的研究區(qū)域作為重力場建模試驗對象，具體范圍為北緯8°N—18°N，東經(jīng)108°E—116°E。輸入數(shù)據(jù)有兩種，一種是丹麥空間研究所提供的2′×2′的自由空氣重力異常數(shù)據(jù)(DTU13)，總共72 000個數(shù)據(jù)；另一種是由EIGEN-GL04S重力場模型計算得到的大地水準(zhǔn)面起伏數(shù)據(jù)。EIGEN-GL04S由GRACE衛(wèi)星數(shù)據(jù)和LAGEOS數(shù)據(jù)聯(lián)合解算得到，但是其在高階精度稍差，因此在建模之前，有必要對其可靠階次進(jìn)行分析。圖5繪出了EIGEN-GL04S模型大地水準(zhǔn)面階方差頻譜。在2～36階范圍內(nèi)，大地水準(zhǔn)面誤差階方差呈先下降后上升的趨勢(藍(lán)線)，累計誤差階方差(黑線)緩慢上升，前36階對應(yīng)的累計誤差方差為0.006 m；在36 階以后，誤差階方差和累計誤差階方差都迅速增大；在125階左右，累計誤差階方差與信號階方差(紅線)相當(dāng)，達(dá)到0.08 m；此后，累計階方差開始大于信號階方差。因此，為了更好地利用EIGEN-GL04S模型的中低階信息，本文采用“移去-恢復(fù)法”，以EIGEN-GL04S模型的前36階作為參考，而高階截斷至120階(累計大地水準(zhǔn)面誤差為0.06 m)?；诖擞嬎懔搜芯繀^(qū)域內(nèi)6′×6′的高程異常數(shù)據(jù)，數(shù)據(jù)總量為8181個。

圖1　空間域正則化球面徑向基函數(shù)Fig.1　A regularized spherical radial basis function in the spatial domain

圖2　Abel-Poisson核函數(shù)不同帶寬參數(shù)下的頻譜Fig.2　Normalized spectrum of Abel-Poisson kernel of different bandwidth

圖3　Abel-Poisson正則化基函數(shù)不同帶寬參數(shù)下的空間形狀Fig.3　The space shape of the Abel-Poisson base function under different bandwidth

圖4　全球Reuter格網(wǎng)Fig.4　Global Reuter grids

圖5　EIGEN-GL04S大地水準(zhǔn)面階方差頻譜Fig.5　EIGEN-GL04S geoid degree variance spectrum

EIGEN-GL04S在中低階波段信號精度較高，但分辨率較低； DTU13重力異常主要由衛(wèi)星測高數(shù)據(jù)得到，分辨率較高，包含了豐富的高頻信號。因此，聯(lián)合這兩種數(shù)據(jù)，彌補(bǔ)彼此之間的不足，便可能得到一個高質(zhì)量、高分辨率的重力場模型。圖6繪制了移去低階重力場貢獻(xiàn)后的兩種數(shù)據(jù)的分布情況。

考慮到徑向基函數(shù)建模普遍出現(xiàn)的“邊緣效應(yīng)”，將基函數(shù)格網(wǎng)建模區(qū)域四周各擴(kuò)展1°，即基函數(shù)格網(wǎng)覆蓋區(qū)域為7°N—19°N，107°E—117°E。初始粗格網(wǎng)采用Reuter格網(wǎng)[18]，取格網(wǎng)密度參數(shù)L=1800，格網(wǎng)總數(shù)為11 768。接下來便是精化格網(wǎng)算法輸入?yún)?shù)的確定，經(jīng)過多次嘗試，當(dāng)取τ1=1、τ2=0.5時，此時得到的重力場反演結(jié)果與觀測值最為相近；另外，為了防止孤立的大剩余觀測值的出現(xiàn)(可能不是真實的重力場信號)，待定基函數(shù)格網(wǎng)點附近應(yīng)存在若干較大的剩余觀測值，本例中取q=3；依據(jù)數(shù)據(jù)分辨率和q的大小，取精化格網(wǎng)點選擇半徑ψc=0.04°；為了防止SRBF過度集中，任意兩個SRBF的球面距應(yīng)大于一定閾值ψmin，設(shè)置為粗格網(wǎng)間距的一半，即ψmin=0.05°。依據(jù)上述篩選原則，得到精化基函數(shù)點數(shù)為2719。最終利用所有選定的基函數(shù)(粗+精)，共同構(gòu)建局部重力場。

表2方差分量估計法正則化條件數(shù)和方差因子

Tab.2Condition number and variance factors of variance component estimate

方法條件數(shù)σ21(m2/s4)σ22/m2σ2μ正則化前5.85×1012111正則化后3.73×1031.6940.1347.42×10-2

本文主要采用Abel-Poisson徑向基函數(shù)建模(Abel-Poisson basis function modeling，APBF)。圖7繪制出了建模后恢復(fù)的重力異常數(shù)據(jù)的誤差分布情況，圖8為該地區(qū)對應(yīng)的SRTM_PLUS30海底地形。

從圖7可以看到，絕大部分地區(qū)誤差值較小，且分布均勻，較大的誤差主要集中于區(qū)域東南角的南沙群島、西部的長山山脈和北部的西沙群島，這主要歸因于這些區(qū)域復(fù)雜的海底地形。南沙群島由于其島礁眾多、分布較廣，且位于海盆的邊緣地帶，地形最為復(fù)雜，建模后擬合誤差最大值也出現(xiàn)在這一區(qū)域，分別為6.80×10-5m/s2、-7.54×10-5m/s2(表3)；區(qū)域東部狹長的長山山脈，地質(zhì)構(gòu)造復(fù)雜，高峰多達(dá)1500～2000 m以上，地勢崎嶇，誤差達(dá)到5.47×10-5m/s2；西沙群島和中沙群島之間東北向的巨型海槽和群島東南側(cè)的中央深海海盆，使得該地區(qū)的重力異常變化比周圍大多數(shù)地方都大，因此該地區(qū)也有較大的誤差；另外，在中沙群島與南沙群島的中間地帶，海底地形起伏雖小于前3個地區(qū)，但仍可以看到明顯的剩余誤差。除上述幾個地區(qū)外，絕大部分區(qū)域建模誤差均小于±1.5×10-5m/s2，區(qū)域整體擬合誤差RMS為±0.80×10-5m/s2(表3)，驗證了應(yīng)用Abel-Poisson徑向基函數(shù)和數(shù)據(jù)自適應(yīng)算法局部構(gòu)建重力場模型的有效性。

EGM08模型是目前已經(jīng)得到公認(rèn)的高階次重力場模型，可以作為檢驗新模型質(zhì)量的參考。但EGM08模型只能展開至2190階(約5′×5′),而本文采用的DTU13數(shù)據(jù)以及APBF模型的分辨率為2′×2′。由于這些數(shù)據(jù)的空間分辨率不同，因此必須轉(zhuǎn)化為相同的空間尺度，再進(jìn)行相互比較。表3列出了EGM08模型與APBF模型以及DTU13數(shù)據(jù)的偏差統(tǒng)計結(jié)果。表3中，APBF模型與EGM08模型的偏差絕對值不超過9.56×10-5m/s2，平均值為-0.076×10-5m/s2，RMS為0.75×10-5m/s2，兩者之間有一定的偏差，但相差不大，進(jìn)一步證明了APBF 模型的有效性和可靠性。

表3APBF、EGM08和DTU13三者之間重力異常偏差比較

Tab. 3Statistics of gravity anomalies bias computed from APBF, EGM08 and DTU13

10-5m/s2

由于兼顧了大地水準(zhǔn)面和重力異常不同頻譜的影響，建模解得的徑向基函數(shù)系數(shù)α不再與先前任何單一數(shù)據(jù)的頻譜相對應(yīng)，而是融合了兩種數(shù)據(jù)的共同特性。因此,依據(jù)這些系數(shù)，不僅可以得到與EIGEN-GL04S相對應(yīng)的低階(120)大地水準(zhǔn)面起伏數(shù)據(jù)，還可以得到與DTU13重力異常頻譜相適應(yīng)的更高階次的其他重力場相關(guān)量。圖9繪制了徑向基函數(shù)模型APBF計算的剩余大地水準(zhǔn)面起伏數(shù)據(jù)圖，相關(guān)的誤差統(tǒng)計見表4。

圖6　南海部分地區(qū)的重力場剩余輸入數(shù)據(jù)Fig.6　The input residual gravity data of parts of the South China Sea

圖7　APBF模型重力異常數(shù)據(jù)恢復(fù)誤差Fig.7　Gravity anomaly recovery error of APBF model

圖8　南海部分地區(qū)海底地形(SRTM_PLUS30)Fig.8Submarine topography of the South China Sea(SRTM_PLUS30)

圖9　APBF模型計算的剩余大地水準(zhǔn)面起伏Fig.9　Geoid undulation residuals recovered by APBF model

圖9(a)為徑向基函數(shù)APBF模型恢復(fù)的120階剩余大地水準(zhǔn)面起伏，通過與EIGEN-GL04S輸入大地水準(zhǔn)面起伏數(shù)據(jù)(圖6(b))比較，兩者符合程度較好，偏差最大為0.453m，最小為-0.471m，RMS為±0.066m(表4)，均保持在較低的水平。

圖9(b)為APBF模型得到的高分辨率大地水準(zhǔn)面起伏(約3.6km的空間分辨率)，為了評估APBF模型的高階大地水準(zhǔn)面起伏數(shù)據(jù)的質(zhì)量，選取EGM08模型作為參考對象，兩個模型同時計算至2190(分辨率約9km)，統(tǒng)計結(jié)果列于表4。

表4APBF與EIGEN-GL04S和EGM08大地水準(zhǔn)面起伏偏差統(tǒng)計

Tab.4Statistics of geoid undulation bias computed from APBF ,EIGEN-GL04S and EGM08

m

由表4可以看到， APBF模型與EGM08模型的偏差不大，從另一角度表明，APBF模型能夠有效地表示高階、高分辨率的大地水準(zhǔn)面。

4　結(jié)　論

本文利用Abel-Poisson 徑向基函數(shù)和改進(jìn)的數(shù)據(jù)自適應(yīng)精化格網(wǎng)算法，聯(lián)合兩種不同類型、不同分辨率的重力場數(shù)據(jù)，對南海部分地區(qū)進(jìn)行建模。徑向基函數(shù)模型APBF對原始重力異常數(shù)據(jù)擬合的符合程度達(dá)到了±0.80×10-5m/s2(分辨率為2′×2′)，絕大部分地區(qū)擬合誤差值小于±1.5×10-5m/s2，分布均勻；由于APBF模型兼顧了兩種數(shù)據(jù)的不同頻譜特性，不僅能較好地恢復(fù)該地區(qū)的低階大地水準(zhǔn)面數(shù)據(jù)，而且還可以得出高階、高分辨率的大地水準(zhǔn)面信息，從而達(dá)到了高精度、高分辨率的局部重力場建模的目的。在數(shù)據(jù)處理中，采用RMS準(zhǔn)則相對GCV可以明顯提高收斂速度，減少計算量。

與其他建模方法一樣，徑向基函數(shù)建模在一定程度上不可避免地受到復(fù)雜地形的影響。但是，應(yīng)用數(shù)據(jù)自適應(yīng)算法，絕大部分誤差均可以控制在較小的范圍內(nèi)，從而保證局部重力場建模的有效性。

基函數(shù)模型可以融合不同來源、不同分布的重力數(shù)據(jù)，如GRACE、衛(wèi)星測高等，從而得到高分辨率、高精度的局部重力場，可以用來精化局部大地水準(zhǔn)面；利用徑向基函數(shù)可以方便地開展重力場多尺度分解方面的研究，結(jié)合其他地球物理資料，有助于更好地了解地球的內(nèi)部構(gòu)造和動力學(xué)變化機(jī)制。

致謝：感謝丹麥技術(shù)大學(xué)空間研究所提供的DTU自由空氣重力異常數(shù)據(jù)。

[1]彭富清, 陳雙軍, 金群峰. 衛(wèi)星測高誤差對海洋重力場反演的影響[J]. 測繪學(xué)報, 2014, 43(4): 337-340. DOI: 10.13485/j.cnki.11-2089.2014.0050.PENG Fuqing, CHEN Shuangjun, JIN Qunfeng. Influence of Altimetry Errors on Marine Geopotential Recovery[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(4): 337-340. DOI: 10.13485/j.cnki.11-2089.2014.0050.

[2]鐘波, 羅志才, 李建成, 等. 聯(lián)合高低衛(wèi)-衛(wèi)跟蹤和衛(wèi)星重力梯度數(shù)據(jù)恢復(fù)地球重力場的譜組合法[J]. 測繪學(xué)報, 2012, 41(5): 735-742.

ZHONG Bo, LUO Zhicai, LI Jiancheng, et al. Spectral Combination Method for Recovering the Earth’s Gravity Field from High-low SST and SGG Data[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 735-742.

[3]HEIKKINEN M. Solving the Shape of the Earth by Using Digital Density Models[R]. Helsinki: Finnish Geodetic Institute, 1981.

[4]HOLSCHNEIDER M, CHAMBODUT A, MANDEA M. From Global to Regional Analysis of the Magnetic Field on the Sphere Using Wavelet Frames[J]. Physics of the Earth and Planetary Interiors, 2003, 135(2-3): 107-124.

[5]CHAMBODUT A, PANET I, MANDEA M, et al. Wavelet Frames: An Alternative to Spherical Harmonic Representation of Potential Fields[J]. Geophysical Journal International, 2005, 163(3): 875-899.

[6]SCHMIDT M, FABERT O, HAN S C, et al. Gravity Field Determination Using Multiresolution Techniques[C]∥Proceedings of the 2nd International GOCE User Workshop. Frascati: GOCE, 2004.

[7]SCHMIDT M, FABERT O, SHUM C K. On the Estimation of a Multi-resolution Representation of the Gravity Field Based on Spherical Harmonics and Wavelets[J]. Journal of Geodynamics, 2005, 39(5): 512-526.

[8]SCHMIDT M, FENGLER M, MAYER-GüRR T, et al. Regional Gravity Modeling in Terms of Spherical Base Functions[J]. Journal of Geodesy, 2007, 81(1): 17-38.

[9]EICKER A. Gravity Field Refinement by Radial Basis Functions from In-situ Satellite Data[D]. Bonn: Universit?t Bonn, 2008.

[10]BENTEL K, SCHMIDT M, GERLACH C. Different Radial Basis Functions and Their Applicability for Regional Gravity Field Representation on the Sphere[J]. GEM-International Journal on Geomathematics, 2013, 4(1): 67-96.

[11]WITTWER T B. Regional Gravity Field Modelling with Radial Basis Functions[D]. Delft: Delft University of Technology, 2009.

[12]KLEES R, TENZER R, PRUTKIN I, et al. A Data-driven Approach to Local Gravity Field Modelling Using Spherical Radial Basis Functions[J]. Journal of Geodesy, 2008, 82(8): 457-471.

[13]PANET I, KUROISHI Y, HOLSCHNEIDER M. Wavelet Modelling of the Gravity Field by Domain Decomposition Methods: An Example Over Japan[J]. Geophysical Journal International, 2011, 184(1): 203-219.

[14]FENGLER M J, FREEDEN W, KOHLHAAS A, et al. Wavelet Modeling of Regional and Temporal Variations of the Earth’s Gravitational Potential Observed by GRACE[J]. Journal of Geodesy, 2007, 81(1): 5-15.

[15]GRAFAREND E W, KRUMM F. The Abel-Poisson Kernel and the Abel-Poisson Integral in a Moving Tangent Space[J]. Journal of Geodesy, 1998, 72(7-8): 404-410.

[16]TENZER R, KLEES R. The Choice of the Spherical Radial Basis Functions in Local Gravity Field Modeling[J]. Studia Geophysica et Geodaetica, 2008, 52(3): 287-304.

[17]KOCH K R, KUSCHE J. Regularization of Geopotential Determination from Satellite Data by Variance Components[J]. Journal of Geodesy, 2002, 76(5): 259-268.

[18]REUTER R. über Integralformeln der Einheitssph?re und Harmonische Splinefunktionen[M]. Rheinisch: Hochschule Aachen, 1982.

[19]SCHMIDT M, HAN S C, KUSCHE J, et al. Regional High-resolution Spatiotemporal Gravity Modeling from GRACE Data Using Spherical Wavelets[J]. Geophysical Research Letters, 2006, 33(8). DOI: 10.1029/2005GL025509.

[20]KOCH K R. Parameter Estimation and Hypothesis Testing in Linear Models[M]. Berlin Heidelberg: Springer, 2013.

[21]PANET I, KUROISHI Y, HOLSCHNEIDER M. Flexible Dataset Combination and Modelling by Domain Decomposition Approaches[M]∥SNEEUW N, NOVK P, CRESPI M, et al. VII Hotine-Marussi Symposium on Mathematical Geodesy. Berlin: Springer, 2012: 67-73.

[22]汪海洪. 小波多尺度分析在地球重力場中的應(yīng)用研究[D]. 武漢: 武漢大學(xué), 2005.

WANG Haihong. Research on Applications of Wavelet Multiscale Analysis in the Earth’s Gravity Field[D]. Wuhan: Wuhan University, 2005.

[23]王虎彪, 王勇, 陸洋, 等. 用衛(wèi)星測高和船測重力資料聯(lián)合反演海洋重力異常[J]. 大地測量與地球動力學(xué), 2005, 25(1): 81-85.

WANG Hubiao, WANG Yong, LU Yang, et al. Inversion of Marine Gravity Anomalies by Combining Multi-altimeter Data and Shipborne Gravimetric Data[J]. Journal of Geodesy and Geodynamics, 2005, 25(1): 81-85.

(責(zé)任編輯：陳品馨)

Regional Gravity Field Modeling with Abel-Poisson Radial Basis Functions

MA Zhiwei1,2，LU Yang1,2，TU Yi1,2，ZHU Chuandong1,2，XI Hui1,2

1. State Key laboratory of Geodesy and Earth’s Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China; 2. University of Chinese Academy of Sciences, Beijing 100049, China

With the increasing number of various types of high-resolution gravity observations, earth gravity models can be regionally refined. We use Abel-Poisson kernel to represent the gravity as the linear summation of finite radial basis functions and combine the multiple gravity data to build a regional gravity model with high resolution. The minimum root mean square criterion based on the data adaptive algorithm is proposed to calculate the base function, which promote the speed of computation significantly. Taking the central South China Sea as an example, two different types of gravity data, namely geoid undulations with resolution of 6′×6′ and gravity anomaly with resolution of 2′×2′, are used to construct the high-resolution regional gravity model. The model has a resolution of 2′×2′, and has a great agreement with original gravity anomaly, reaching to ±0.8×10-5m/s2.Our results show that using radial basis functions to construct the regional gravity field can avoid the problem of slow convergence of spherical harmonic functions, and can improve the resolution remarkably.

local gravity field modeling; gravity anomaly; radial basis function; data adaptive; root mean square error criterion

The National Basic Research Program of China (973) (No.2013CB733301); The National Natural Science Foundation of China (Nos. 41274025; 41174064)

馬志偉，陸洋，涂弋，等.利用Abel-Poisson徑向基函數(shù)模型化局部重力場[J].測繪學(xué)報，2016,45(9)：1019-1027.

10.11947/j.AGCS.2016.20150519.

MA Zhiwei, LU Yang, TU Yi, et al.Regional Gravity Field Modeling with Abel-Poisson Radial Basis Functions[J]. Acta Geodaetica et Cartographica Sinica,2016,45(9):1019-1027. DOI:10.11947/j.AGCS.2016.20150519.

P223

A

1001-1595(2016)09-1019-09

國家973計劃(2013CB733301)；國家自然科學(xué)基金(41274025;41174064)

2015-10-14

馬志偉(1986—)，男，博士生，主要從事多源重力數(shù)據(jù)融合研究工作。First author： MA Zhiwei(1986—)，male, PhD candidate, majors in multi-source gravity data fusion

E-mail： jzmazhiwei@163.com

修回日期： 2016-06-27