梁偉，李建成,b，徐新禹,b,*，張勝軍，趙永奇

a School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China

b Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China

c School of Resources and Civil Engineering, Northeastern University, Shenyang 110004, China

1. 引言

高分辨率地球重力場(chǎng)模型可用于高精度大地水準(zhǔn)面確定、全球高程基準(zhǔn)統(tǒng)一、動(dòng)態(tài)海面地形確定和地球內(nèi)部結(jié)構(gòu)探測(cè)等科學(xué)和工程應(yīng)用研究[1-5]。隨著新一代衛(wèi)星重力計(jì)劃（Challenging Minisatellite Payload,CHAMP）的出現(xiàn)[6]以及重力場(chǎng)與氣候?qū)嶒?yàn)衛(wèi)星（Gravity Recovery and Climate Experiment, GRACE）[7]和重力場(chǎng)與海洋環(huán)流探測(cè)衛(wèi)星（Gravity field and Steady-State Ocean Circulation Explorer, GOCE）[8]的成功發(fā)射，地球重力場(chǎng)模型的中長(zhǎng)波精度得到了顯著提高[9-16]；同時(shí)，地面重力數(shù)據(jù)、衛(wèi)星測(cè)高數(shù)據(jù)、航空重力數(shù)據(jù)和地形數(shù)據(jù)包含了高精度的短波信號(hào)[17]，因此可以聯(lián)合衛(wèi)星、測(cè)高和地面等多源重力數(shù)據(jù)解算高精度、高分辨率的地球重力場(chǎng)模型。

國際地球重力場(chǎng)模型中心（International Center for Global Earth Models, ICGEM; http://icgem.gfz-potsdam.de/tom_longtime）發(fā)布的高分辨率地球重力場(chǎng)模型有：Earth Gravitational Model (EGM) 2008 [18]、European improved gravity model of the earth by new techniques (EIGEN)-6C4 [19]、GOCE and EGM2008 Combined model（GECO） [20]、gravity observation combination (GOCO)模型GOCO05c [21]、Experimental Gravity Field Model(XGM) 2016 [22]和SGG-UGM-1 [23]，表1給出了這些模型的屬性。由于這些模型使用了新一代重力衛(wèi)星的數(shù)據(jù)，模型的中長(zhǎng)波精度得到了顯著提高。EGM2008是當(dāng)前使用最多的重力場(chǎng)模型，它的計(jì)算中使用了幾乎精度最高的5′×5′分辨率全球重力異常和ITG-GRACE03S模型的法方程，但是并沒有用到GOCE衛(wèi)星的數(shù)據(jù)；其中全球重力異常數(shù)據(jù)融合了地面重力數(shù)據(jù)、衛(wèi)星測(cè)高數(shù)據(jù)和地形數(shù)據(jù)。EIGEN-6C4是EIGEN系列模型的代表，相較于EGM2008，它使用了GOCE衛(wèi)星數(shù)據(jù)和Laser Geodynamics Satellite (LAGEOS)衛(wèi)星數(shù)據(jù)，不過地面重力異常數(shù)據(jù)由EGM2008模型計(jì)算。EGM2008和EIGEN-6C4的高階次系數(shù)都是采用塊對(duì)角最小二乘法解算得到，而720階次GOCO05c和719階次的XGM2016是聯(lián)合衛(wèi)星和重力異常法方程，采用嚴(yán)格最小二乘方法估計(jì)的。GECO和SGG-UGM-1都是通過使用GOCE數(shù)據(jù)精化EGM2008模型得到[22-24]。文獻(xiàn)[23]的檢核結(jié)果表明EGM2008、EIGEN-6C4和SGG-UGM-1在美國的精度接近，而包含GOCE數(shù)據(jù)的EIGEN-6C4和SGG-UGM-1在中國的精度優(yōu)于EGM2008。

本文將介紹的SGG-UGM-2模型與SGG-UGM-1同屬一個(gè)系列的高分辨率重力場(chǎng)模型，它是SGG-UGM-1模型的進(jìn)一步發(fā)展，采用橢球諧分析方法計(jì)算，并使用新的海域重力異常數(shù)據(jù)，以及融入GRACE衛(wèi)星觀測(cè)法方程，以此提升模型的精度。

由于地球更接近于橢球，球近似的誤差要大于橢球近似，所以在全球重力場(chǎng)模型計(jì)算中往往將地面重力數(shù)據(jù)延拓到參考橢球面上，如Ohio States University (OSU)91 [25]、EGM96 [26]和EGM2008 [18]，因此橢球諧分析比球諧分析更適合于橢球面數(shù)據(jù)的處理[27]。Hotine[28]和Jekeli [29,30]對(duì)第二類勒讓德函數(shù)進(jìn)行了規(guī)格化，并推導(dǎo)了球諧系數(shù)與橢球諧系數(shù)之間的轉(zhuǎn)換公式（下文稱為“Jekeli轉(zhuǎn)換”或“Jekeli系數(shù)轉(zhuǎn)換”）。Gleason [31]通過數(shù)值試驗(yàn)驗(yàn)證了Jekeli轉(zhuǎn)換的精度。Sebera等[32]將勒讓德函數(shù)的計(jì)算擴(kuò)展至二階導(dǎo)數(shù)，并減少了計(jì)算超幾何變換需要的迭代次數(shù)。

在基于橢球面重力異常構(gòu)建全球重力場(chǎng)模型時(shí)，一般首先對(duì)橢球面的重力異常數(shù)據(jù)進(jìn)行橢球諧分析得到橢球諧系數(shù)，然后將橢球諧系數(shù)轉(zhuǎn)換為擾動(dòng)位的球諧系數(shù)（Csnm），本文將這種計(jì)算重力場(chǎng)模型的方法簡(jiǎn)稱為EHACT (Ellipsoidal Harmonic Analysis and Coefficients Transformation) 方法。EHA-CT方法可能有不同具體實(shí)現(xiàn)策略，例如，可以通過最小二乘方法進(jìn)行橢球諧分析，也可以使用數(shù)值積分方法實(shí)現(xiàn)橢球諧分析。Rapp和Pavlis[33]推導(dǎo)了使用橢球面平均重力異常數(shù)據(jù)計(jì)算系數(shù)Csnm的計(jì)算公式，該公式用在了OSU89 [33]、OSU91 [25]、EGM96 [26]、IGG97LB [34]和MOD99 [35]等模型的計(jì)算中。前述ICGEM發(fā)布的高分辨率重力場(chǎng)模型中，只有EGM2008是使用橢球諧分析計(jì)算的，在EGM2008的求解中，首先使用了最小二乘方法進(jìn)行橢球諧分析得到橢球諧系數(shù)（），然后使用Jekeli轉(zhuǎn)換將轉(zhuǎn)換為。本文將首先介紹EHA-CT方法，并推導(dǎo)相關(guān)的嚴(yán)密公式，及引入Driscoll和Healy [36]的采樣理論，得到適用于離散點(diǎn)值的嚴(yán)密積分公式。需要指出的是，本文推導(dǎo)的適用于格網(wǎng)均值的積分公式和文獻(xiàn)[33]推導(dǎo)的公式略有不同。此外，需要說明的是：由橢球面重力異常數(shù)據(jù)計(jì)算重力場(chǎng)模型屬于橢球邊界面大地測(cè)量邊值問題的范疇，除EHA-CT方法之外，文獻(xiàn)[36-43]也給出了其他解算方法。

海域面積占地球表面的71%，高精度的海域重力異常是構(gòu)建高分辨率重力場(chǎng)模型的關(guān)鍵，越來越多的衛(wèi)星測(cè)高數(shù)據(jù)可以用于海域重力異常的確定，如Geosat GM /ERM (17 d), ERS-1/GM (168 d), ERS/ERM (35 d), T/P/T/P Tandem (10 d), Jason-1/ERM (10 d), Envisat (35d/30 d),Jason-2/ERM (10 d), Jason-1/GM (406 d), CryoSat-2 (369 d),SARAL/AltiKa ERM(35 d), HY-2A (14 d), Jason-2/GM 和SARAL/AltiKa GM。由這些多星測(cè)高數(shù)據(jù)可以使用數(shù)值分析方法[44,45]或者最小二乘配置方法[46,47]計(jì)算得到緯度-80.738°～80.738°范圍內(nèi)1′×1′空間分辨率的重力異常數(shù)據(jù)，用于高分辨率地球重力場(chǎng)模型的解算。本文使用衛(wèi)星測(cè)高數(shù)據(jù)計(jì)算了海洋區(qū)域1′×1′的重力異常數(shù)據(jù)，并與EGM2008計(jì)算的陸地重力異常數(shù)據(jù)融合，組成了全球重力異常格網(wǎng)數(shù)據(jù)集。除此之外，GRACE衛(wèi)星提供了長(zhǎng)達(dá)15年的衛(wèi)星重力觀測(cè)數(shù)據(jù)，使用這些數(shù)據(jù)可以計(jì)算高精度的重力場(chǎng)長(zhǎng)波信號(hào)。ITSG-Grace2018 [48]是ITSG(Institute of Theoretical Geodesy and Satellite Geodesy)系列GRACE重力場(chǎng)模型的最新模型，模型的作者將模型的法方程與模型同時(shí)發(fā)布，本文將利用該法方程與其他類型的重力數(shù)據(jù)一起用于高分辨率重力場(chǎng)模型的解算。SGGUGM-1中使用的海域重力異常來自EGM2008，且沒有用到GRACE數(shù)據(jù)，SGG-UGM-2的解算中將使用新的海域重力異常數(shù)據(jù)和ITSG-Grace2018模型法方程。

表1 ICGEM發(fā)布的高分辨率重力場(chǎng)模型的主要屬性

本文分為6個(gè)部分，首先在第2部分介紹EHA-CT方法的基本原理，并推導(dǎo)相應(yīng)的計(jì)算公式；在第3部分使用數(shù)值試驗(yàn)對(duì)本文推導(dǎo)的公式進(jìn)行驗(yàn)證；第4部分介紹SGG-UGM-2模型的具體計(jì)算過程；第5部分給出模型的檢核結(jié)果；最后在第6部分進(jìn)行總結(jié)。

2. 基本原理與方法

2.1. 由橢球面重力異常數(shù)據(jù)求解重力場(chǎng)模型系數(shù)的EHA-CT方法

任意滿足拉普拉斯方程的調(diào)和函數(shù)f在其調(diào)和區(qū)域內(nèi)可用橢球諧函數(shù)級(jí)數(shù)表示為[30,32]：

式中，(u, δ, λ)為橢球坐標(biāo)[1]；E為參考橢球的線性偏心率；b為參考橢球的短半軸；Ynm(δ, λ)為完全正規(guī)化的面球諧函數(shù)：

調(diào)和函數(shù)rΔg的橢球諧展開式為[31,33]：

式中，rΔg為r和Δg的積，r為地心向徑，Δg為重力異常；為rΔg的橢球諧展開系數(shù)；a為參考橢球的長(zhǎng)半徑。

當(dāng)重力異常數(shù)據(jù)位于參考橢球面上時(shí)，式（2）可簡(jiǎn)化為：

式中，[]E表示數(shù)據(jù)位于參考橢球表面。

式（4）所示的系數(shù)轉(zhuǎn)換方法為，首先使用Jekeli系數(shù)轉(zhuǎn)換將系數(shù)轉(zhuǎn)為系數(shù)g的球諧展開式系數(shù)，然后再將系數(shù)轉(zhuǎn)換為待求的擾動(dòng)位的球諧系數(shù)，下節(jié)將給出這些轉(zhuǎn)換的計(jì)算公式。

在重力場(chǎng)模型求解中EHA-CT方法會(huì)有不同的實(shí)現(xiàn)形式，如果利用數(shù)值積分方法進(jìn)行橢球諧分析，則該實(shí)現(xiàn)形式稱為EHA-CT方法的積分形式；如果使用最小二乘方法進(jìn)行橢球諧分析，則稱為EHA-CT方法的最小二乘形式。而且模型求解時(shí)使用的離散的重力異常觀測(cè)數(shù)據(jù)可以是格網(wǎng)點(diǎn)值數(shù)據(jù)或者格網(wǎng)平均值數(shù)據(jù)，針對(duì)不同的離散重力數(shù)據(jù)，EHA-CT方法也有不同的計(jì)算公式，2.1.1節(jié)和2.1.2節(jié)將給出這些計(jì)算公式。

2.1.1. 積分形式的EHA-CT方法

根據(jù)面球諧函數(shù)的正交性，由式（3）可以推導(dǎo)出使用橢球面重力異常數(shù)據(jù)計(jì)算系數(shù)的公式：

式中，文獻(xiàn)[30]中有Λnmk的計(jì)算方法。

式中，GM為地球引力常數(shù)。

將式（5）和式（6）代入式（7）可得：

式中，γ=GM/a2。

雖然式（8）是嚴(yán)密的，但該式中的數(shù)據(jù)為參考橢球面連續(xù)的重力異常數(shù)據(jù)，而重力場(chǎng)模型解算可使用的數(shù)據(jù)一般為離散重力數(shù)據(jù)，所以該式不能直接用于重力場(chǎng)模型求解，需要把它離散化。針對(duì)不同形式的離散數(shù)據(jù)，式（8）有不同的離散化形式，有等角格網(wǎng)離散點(diǎn)值數(shù)據(jù)的離散化形式，也有等角格網(wǎng)離散平均值數(shù)據(jù)的形式。

如果離散數(shù)據(jù)為等角格網(wǎng)的重力異常平均值，式（8）的離散化形式為：

式中，i和j表示數(shù)據(jù)位于橢球面上的第i (i = 0, …, N-1)個(gè)緯度帶和第j (j = 0, …, 2N-1)個(gè)子午帶，這樣的全球等角格網(wǎng)共有N行和2N列，i_max=N-1, j_max=2N-1；[ri]E為橢球面上第i個(gè)緯度帶上格網(wǎng)中心的地心向徑；為格網(wǎng)平均重力異常值是格網(wǎng)內(nèi)[rΔg]E的平均值；是勒讓德函數(shù)的積分和的計(jì)算公式為Appendix A中有該式更為詳細(xì)的推導(dǎo)。

如果離散數(shù)據(jù)為等角格網(wǎng)的重力異常點(diǎn)值數(shù)據(jù)，式（8）的離散化形式為：

盡管式（11）中引入了平滑因子，但由于勒讓德函數(shù)的正交性問題并未完全解決，計(jì)算得到的系數(shù)仍然有離散化誤差，當(dāng)前沒有閱讀到相關(guān)文獻(xiàn)中有這個(gè)問題的解決辦法，本文將在第3節(jié)中分析該離散化誤差的大小。

然而，針對(duì)格網(wǎng)點(diǎn)值，可采用特殊的采樣和賦權(quán)方法解決數(shù)據(jù)離散化中的勒讓德函數(shù)正交性問題，如Gauss-Legendre (GL) [50]積分方法和Driscoll/Healy (DH)[36]積分方法等，本文將該方法引入到基于離散重力異常點(diǎn)值估計(jì)模型位系數(shù)。根據(jù)DH方法，可將橢球面劃分為[N×N] (nmax= N/2-1)的格網(wǎng)，其中，nmax為待求系數(shù)的最高階，緯度方向的格網(wǎng)間隔為Δδ = 180°/N，經(jīng)度方向的格網(wǎng)間隔為Δλ = 360°/N；也可劃分為Δδ = Δλ =180°/N的等角格網(wǎng)。DH方法定義數(shù)據(jù)的權(quán)為：

將式（12）中的權(quán)因子代入到式（10），可得橢球面點(diǎn)值數(shù)據(jù)計(jì)算系數(shù)的公式為：

與DH方法類似，GL方法同樣通過定義特殊的格網(wǎng)劃分方式和權(quán)因子來解決勒讓德函數(shù)的正交性問題[50]。GL方法中定義的格網(wǎng)是不規(guī)則的格網(wǎng)，緯度方向上格網(wǎng)的采樣間隔相等，經(jīng)度方向上的采樣點(diǎn)為Pn|m|函數(shù)的零值點(diǎn)。Hirt等[51]使用GL方法對(duì)卷積積分進(jìn)行數(shù)值分析，文獻(xiàn)[50,51]中有GL方法的格網(wǎng)劃分和賦權(quán)方法的詳細(xì)介紹。

使用式（11）和式（13），可以分別由橢球面格網(wǎng)重力異常均值和點(diǎn)值計(jì)算重力場(chǎng)模型系數(shù)。Rapp和Pavlis [33]中的式（20）也是使用橢球面重力異常均值數(shù)據(jù)計(jì)算系數(shù)的公式，其推導(dǎo)方法與本文的式（11）相同，然而式（11）和文獻(xiàn)[33]中的式（20）略有不同，其差異體現(xiàn)在本文推導(dǎo)式（11）和文獻(xiàn)[33]中的一項(xiàng)因子分別為1/(n-1)和1/(n-2k-1)。從上文的推導(dǎo)可知，該差異項(xiàng)表征了系數(shù)和系數(shù)之間的關(guān)系，而從式（7）可知，這兩個(gè)系數(shù)只與n有關(guān)。該項(xiàng)的差異會(huì)對(duì)系數(shù)求解產(chǎn)生影響，本文將在第3節(jié)分析該影響的大小。

2.1.2. EHA-CT方法的最小二乘計(jì)算公式

其中

根據(jù)式（14），在Gauss-Markov模型下，使用橢球面數(shù)據(jù)的函數(shù)模型和隨機(jī)模型為：

同時(shí)，將式（6）和式（7）代入式（14）中，可以得到橢球面重力異常均值和點(diǎn)值數(shù)據(jù)與系數(shù)Cnsm的關(guān)系：

根據(jù)式（17）可以使用系數(shù)計(jì)算參考橢球面的重力異常均值和點(diǎn)值數(shù)據(jù)，同樣可以使用該式由橢球面重力異常數(shù)據(jù)建立求解系數(shù)的觀測(cè)方程。

EGM2008的計(jì)算中也使用了最小二乘形式的EHACT方法[18]，計(jì)算策略如下：

對(duì)比式（4）和式（18）可知，該策略與本文使用最小二乘方法的策略不同，如式（18）所示，Pavlis等[18]使用最小二乘方法首先求得系數(shù)，然后再使用Jekeli系數(shù)轉(zhuǎn)換將系數(shù)轉(zhuǎn)換為系數(shù)。

2.2. 重力異常數(shù)據(jù)和衛(wèi)星重力數(shù)據(jù)的聯(lián)合求解方法

本文使用最小二乘方法聯(lián)合多類的重力數(shù)據(jù)，如重力異常數(shù)據(jù)、GOCE衛(wèi)星數(shù)據(jù)和GRACE衛(wèi)星數(shù)據(jù)求解重力場(chǎng)模型。在最小二乘方法中，當(dāng)多類數(shù)據(jù)不相關(guān)時(shí)，使用最小二乘方法聯(lián)合多類數(shù)據(jù)得到的聯(lián)合解為：

式中，Ai為根據(jù)第i類數(shù)據(jù)建立的設(shè)計(jì)矩陣，它表征了待求參數(shù)與第i類數(shù)據(jù)之間的關(guān)系；Pi為數(shù)據(jù)i的權(quán)陣；li是觀測(cè)數(shù)據(jù)向量；是數(shù)據(jù)i的初始方差分量；Ni和Ui是數(shù)據(jù)的法方程矩陣和法方程向量；是數(shù)據(jù)i的初始相對(duì)權(quán)。

由于聯(lián)合不同法方程時(shí)使用的初始方差分量一般不準(zhǔn)確，本文在聯(lián)合多個(gè)法方程時(shí)使用方差分量估計(jì)方法（variance component estimation method, VCE）[53]通過迭代估計(jì)數(shù)據(jù)的方差分量，并在迭代的過程中求解得到多個(gè)法方程的聯(lián)合解。在VCE方法中，方差分量的計(jì)算公式為[53]：

式中，k表示第k次迭代；vk,i為殘差向量；rk,i為冗余度；文獻(xiàn)[53]中有更為詳盡的計(jì)算公式。

然而，在使用VCE方法聯(lián)合多個(gè)法方程時(shí)會(huì)遇到一些問題。一方面，如果使用的法方程非自主構(gòu)建時(shí)，法方程的提供者有時(shí)不能同時(shí)提供殘差的加權(quán)和，這時(shí)就不能直接使用式（20）計(jì)算方差分量；另一方面，如果這些法方程的誤差隨機(jī)模型不同且不能準(zhǔn)確建模時(shí)，VCE方法不能得到最優(yōu)解[54]。Jean等[55]等提出了改進(jìn)的VCE方法，在使用VCE方法時(shí)規(guī)避了直接使用觀測(cè)值殘差計(jì)算方差分量，而是根據(jù)每個(gè)法方程的解相對(duì)聯(lián)合解的“誤差”計(jì)算每個(gè)法方程的方差分量，計(jì)算公式如下[55]：

這樣在每次迭代過程中，可以根據(jù)此次迭代的聯(lián)合解，解算每個(gè)法方程的方差分量，并用于下一次迭代，直到迭代收斂并解得最終的聯(lián)合解。

解算SGG-UGM-1模型時(shí)，由于重力異常數(shù)據(jù)法方程和衛(wèi)星重力數(shù)據(jù)法方程的待求參數(shù)均為重力場(chǎng)模型的球諧系數(shù)，將這些法方程加權(quán)求和即可得到不同數(shù)據(jù)的聯(lián)合觀測(cè)法方程。然而求解SGG-UGM-2模型時(shí)，重力異常數(shù)據(jù)法方程的待求參數(shù)為橢球諧系數(shù)，衛(wèi)星重力數(shù)據(jù)法方程的參數(shù)為球諧系數(shù)在聯(lián)合這些法方程之前首先應(yīng)把法方程的參數(shù)統(tǒng)一。本文在法方程聯(lián)合之前，將參數(shù)為球諧系數(shù)的衛(wèi)星觀測(cè)法方程轉(zhuǎn)換為參數(shù)為橢球諧系數(shù)的法方程。

式中，下標(biāo)“s”表示相應(yīng)的參量為衛(wèi)星觀測(cè)數(shù)據(jù)的參量；x^s為待求系數(shù)向量。

式中，Tse為根據(jù)式（6）和式（7）計(jì)算的系數(shù)轉(zhuǎn)換矩陣；x^e為系數(shù)向量。

將式（23）代入式（22）中，可得：

在最小二乘準(zhǔn)則下，可得式（24）的法方程為：

分析式（25）可知，該式給出了不同參數(shù)形式的法方程之間的關(guān)系，因此可以使用該式將參數(shù)為球諧系數(shù)的衛(wèi)星觀測(cè)法方程轉(zhuǎn)換為參數(shù)為橢球諧系數(shù)的法方程，這樣將法方程的參數(shù)形式統(tǒng)一后，再將這些法方程加權(quán)求和即可得到聯(lián)合觀測(cè)法方程。

3. EHA-CT方法計(jì)算公式的檢驗(yàn)

本節(jié)的主要目的是通過數(shù)值試驗(yàn)檢驗(yàn)第2節(jié)推導(dǎo)的計(jì)算公式。一方面，這可以保證用于SGG-UGM-2模型求解的計(jì)算公式是嚴(yán)密的；另一方面，2.1.1節(jié)的分析表明本文推導(dǎo)的由橢球面重力異常均值計(jì)算系數(shù)的公式與Rapp和Pavlis [33]的相應(yīng)公式存在差異，通過數(shù)值試驗(yàn)可以驗(yàn)證公式的嚴(yán)密性，并分析該差異對(duì)求解系數(shù)的影響，為研究人員選擇合適的公式計(jì)算重力場(chǎng)模型。

首先，使用重力異常均值數(shù)據(jù)通過式（11）計(jì)算得到一個(gè)重力場(chǎng)模型，稱為Model1。圖1中的紅色曲線展示了Model1相對(duì)EGM2008模型的系數(shù)階誤差。2.1.1節(jié)的分析表明，盡管式（11）中引入了平滑因子，該式計(jì)算的系數(shù)仍有離散化誤差，圖1中Model1的誤差反映了離散化誤差的大小，該項(xiàng)誤差對(duì)模型系數(shù)的影響最大可達(dá)10-11量級(jí)，這個(gè)誤差在模型解算中不可忽略。除特殊說明外，本文中模型的系數(shù)誤差均為絕對(duì)誤差。

然后，使用重力異常點(diǎn)值數(shù)據(jù)分別通過式（13）和式（10）計(jì)算了模型Model2和Model3，式（13）和式（10）的區(qū)別在于是否使用DH權(quán)因子。模型Model2和Model3的系數(shù)階誤差分別為圖1中的藍(lán)色和深灰色曲線。由圖1可知，模型Model2全頻段的系數(shù)階誤差均小于10-17量級(jí)，可以認(rèn)為該誤差是計(jì)算機(jī)截?cái)嗾`差引起的。模型Model3的系數(shù)階誤差相對(duì)較大，而且在奇數(shù)階和偶數(shù)階之間有很大的波動(dòng)。模型Model3的誤差遠(yuǎn)大于模型Model1和Model2。結(jié)果說明式（13）是嚴(yán)密的，可用于模型反演，同時(shí)也說明，本文采用DH權(quán)因子有效解決了積分離散化后勒讓德函數(shù)的正交性問題。

同時(shí)，使用模擬的重力異常點(diǎn)值數(shù)據(jù)，通過2.1.2節(jié)中的最小二乘方法計(jì)算了模型Model4，具體步驟為：首先，建立了式（14）所示的觀測(cè)方程，然后求解該觀測(cè)方程得到了系數(shù)；然后使用式（6）和式（7）將系數(shù)轉(zhuǎn)換為系數(shù)。求解系數(shù)時(shí)使用了BDLS方法和OpenMP [57]并行計(jì)算技術(shù)提高計(jì)算效率。模型Model4的系數(shù)階誤差RMS為圖1中的品紅色曲線，從圖中可以看出Model4模型的誤差很小，可以忽略，這充分表明EHA-CT方法最小二乘形式的公式是嚴(yán)密的，可以用于SGGG-UGM-2模型的求解。而且，相對(duì)于積分方法，最小二乘方法更方便地聯(lián)合解算多類觀測(cè)數(shù)據(jù)，因此本文將使用最小二乘方法求解SGG-UGM-2。此外，我們也按照計(jì)算Model4的方法，使用重力異常均值數(shù)據(jù)計(jì)算了相應(yīng)的模型，得到了與點(diǎn)值數(shù)據(jù)一致的結(jié)果，在本文中未做展示。

由于本文推導(dǎo)的式（11）和文獻(xiàn)[33]推導(dǎo)的式（20）均存在離散化誤差，其對(duì)解算系數(shù)的影響和公式差異對(duì)求解系數(shù)的影響會(huì)混疊在一起，在數(shù)值模擬試驗(yàn)中很難分辨由式（11）估計(jì)系數(shù)的誤差是來源于系數(shù)公式嚴(yán)密性的影響，還是來源于離散化的影響，因此很難以此檢驗(yàn)公式的嚴(yán)密性?？紤]到式（13）與式（11）類似，都是使用離散重力異常計(jì)算系數(shù)的公式，它們的區(qū)別在于式（13）中的數(shù)據(jù)為重力異常點(diǎn)值數(shù)據(jù)，它們有相同的因子項(xiàng)1/(n-1)。類似地，如果文獻(xiàn)[33]中的式（20）是正確的，那么將式（13）中的1/(n-1)替換為1/(n-2k-1)，就可以得到與該式相對(duì)應(yīng)的且采用了DH權(quán)的，由重力異常點(diǎn)值數(shù)據(jù)計(jì)算模型系數(shù)的公式：

式（16）的嚴(yán)密性可以反映其因子項(xiàng)1/(n-2k-1)的正確性，進(jìn)而反映文獻(xiàn)[33]中公式（20）中1/(n-2k-1)的正確性，這是本文嚴(yán)格檢驗(yàn)該公式與本文公式嚴(yán)密性的思路。為此，我們使用式（26），由格網(wǎng)重力異常點(diǎn)值數(shù)據(jù)計(jì)算了模型Model5，Model5的階誤差RMS為圖1中的綠色曲線。從圖中可知，在低階次部分，Model5的系數(shù)誤差較大，在5階附近可達(dá)10-9量級(jí)。低階部分（小于50階）的誤差大于10-11，這個(gè)誤差不可以忽略。因此可以推斷文獻(xiàn)[33]中式（20）的1/(n-2k-1)是不準(zhǔn)確的。

為了進(jìn)一步分析文獻(xiàn)[33]中的式（20）的嚴(yán)密性，分別使用模型Model2和Model5計(jì)算了參考橢球面上的重力異常和大地水準(zhǔn)面誤差，表2給出了誤差的統(tǒng)計(jì)結(jié)果，圖2給出了其空間分布。由表2可知，模型Model5在2160階的重力異常和大地水準(zhǔn)面的誤差分別為0.18 mGal (1 mGal = 1×10-5m·s-2)和11 cm，這遠(yuǎn)遠(yuǎn)大于Model2的誤差。Model5在階次100～2160頻段的大地水準(zhǔn)面誤差的最大值為3.5 cm，在階次200～2160頻段的大地水準(zhǔn)面誤差最大值為1.9 cm；在階次100～2160頻段的重力異常誤差的最大值為1.4 mGal，在階次200～2160頻段的重力異常誤差最大值為1.2 mGal。由圖2可知，Model5有較大且為系統(tǒng)性的大地水準(zhǔn)面誤差，而Model2的大地水準(zhǔn)面誤差遠(yuǎn)小于Model5。這些結(jié)果反映了式（26）中1/(n-2k-1)的影響不可忽略。

4. SGG-UGM-2高分辨率重力場(chǎng)模型的解算

本文聯(lián)合衛(wèi)星測(cè)高數(shù)據(jù)、衛(wèi)星重力數(shù)據(jù)和重力異常數(shù)據(jù)計(jì)算了2190階2160次的高分辨率重力場(chǎng)模型SGGUGM-2。本節(jié)將首先簡(jiǎn)要介紹不同重力數(shù)據(jù)（衛(wèi)星重力和衛(wèi)星測(cè)高）的處理策略，然后介紹多個(gè)法方程的聯(lián)合解算策略。

4.1. 建立GOCE和GRACE衛(wèi)星觀測(cè)法方程

圖1. 模型Model 1（紅色）、Model 2（藍(lán)色）、Model 3（深灰）、Model 4（品紅）和Model 5（綠色）相對(duì)EGM2008的系數(shù)階誤差RMS。（a）從2階到2160階；（b）2階到200階。黑線是EGM2008模型系數(shù)的大小。

表2 模型Model2和Model5相對(duì)EGM2008的重力異常和大地水準(zhǔn)面誤差的統(tǒng)計(jì)結(jié)果

圖2. 模型Model 2（a）和Model 5（b）的大地水準(zhǔn)面誤差的空間分布（2160階）。

本文使用的GOCE衛(wèi)星重力觀測(cè)數(shù)據(jù)為2009年11月1日至2012年5月31日的EGG_NOM_2 (GGT, PRD,IAQ, PRM)數(shù)據(jù)和2009年11月1日至2010年7月5日的SST_PSO_2 (PKI, PCV, CCD)數(shù)據(jù)，數(shù)據(jù)的采樣間隔為1 s [58]。EGG_NOM_2數(shù)據(jù)主要包含梯度儀坐標(biāo)系（gradiometer reference frame, GRF）下的重力梯度張量數(shù)據(jù)（gravity gradient tensor, GGT），用于慣性參考系（inertial reference frame, IRF）和GRF坐標(biāo)系轉(zhuǎn)換的四元數(shù)數(shù)據(jù)EGG_IAQ_2，以及加速度數(shù)據(jù)EGG_CCD_2C。SST_PSO_2數(shù)據(jù)包含衛(wèi)星運(yùn)動(dòng)學(xué)軌道SST_PKI_2（PKI軌道）、精密PKI軌道的方差-協(xié)方差信息SST_PCV_2、簡(jiǎn)化動(dòng)力學(xué)軌道SST_PRD_2及地球固定參考系（earth-fixed reference frame, EFRF）與IRF之間轉(zhuǎn)換的四元數(shù)SST_PRM_2。

模型反演中僅使用了GGT數(shù)據(jù)中精度較高的對(duì)角線分量（Vxx, Vyy, Vzz）用于法方程的構(gòu)建[24]。使用衛(wèi)星重力梯度數(shù)據(jù)（satellite gravity gradient, SGG）和高低衛(wèi)星跟蹤數(shù)據(jù)（satellite-to-satellite in high-low mode, SST-hl）分別構(gòu)建了法方程，兩個(gè)法方程參數(shù)的最高階分別為220和130。建立GOCE衛(wèi)星觀測(cè)法方程的主要解算策略如下。

（1）首先對(duì)SGG數(shù)據(jù)和SST-hl數(shù)據(jù)進(jìn)行數(shù)據(jù)內(nèi)插、粗差探測(cè)等預(yù)處理。

（2）使用直接法由SGG數(shù)據(jù)建立法方程[24]，其中在建立觀測(cè)方程時(shí)進(jìn)行了通帶為5～41 mHz的ARMA(auto regressive moving-average)濾波[59]，以解決SGG數(shù)據(jù)的有色噪聲問題。根據(jù)公式fmax=Nmax/Tr計(jì)算的通帶近似最高頻率為fmax= 41 mHz，其中Tr=5383 s，是衛(wèi)星繞地球旋轉(zhuǎn)一周的所用時(shí)間）[60]，Nmax= 220，是求解模型的最高階。

（3）采用點(diǎn)域加速度方法建立GOCE衛(wèi)星SST-hl法方程，求解SST法方程得到觀測(cè)值殘差[60-62]。衛(wèi)星的加速度是根據(jù)衛(wèi)星運(yùn)動(dòng)學(xué)軌道使用EDF5 (extended differentiation fi lter)方法計(jì)算的，其中采樣間隔Δt = 5 s [62]。

（4）根據(jù)SGG觀測(cè)法方程和SST-hl觀測(cè)法方程的方差分量將兩個(gè)法方程聯(lián)合，文獻(xiàn)[24]中有法方程聯(lián)合的具體實(shí)現(xiàn)方法。

這是本文建立GOCE衛(wèi)星觀測(cè)法方程的簡(jiǎn)要介紹，該法方程也是純GOCE衛(wèi)星模型GOSG01S和高分辨率模型SGG-UGM-1計(jì)算中使用的法方程。文獻(xiàn)[24]中有建立該法方程的數(shù)據(jù)處理方法和計(jì)算流程的詳細(xì)介紹。

本文沒有直接使用GRACE衛(wèi)星的原始觀測(cè)數(shù)據(jù)構(gòu)建法方程，而是使用GRACE靜態(tài)重力場(chǎng)模型ITSGGrace2018 [48]的法方程作為GRACE衛(wèi)星觀測(cè)法方程。ITSG-Grace2018模型的研制團(tuán)隊(duì)發(fā)布了SINEX [63]格式的法方程，本文將SINEX格式的法方程轉(zhuǎn)換格式后，將其用于SGG-UGM-2的聯(lián)合解算。

4.2. 全球海域重力異常的確定

本文使用Geosat、ERS-1、Envisat、T/P、Jason-1、CryoSat-2和SARAL/AltiKa等多顆衛(wèi)星的測(cè)高數(shù)據(jù)恢復(fù)了全球海域重力異常[45]，表3給出了粗差剔除、低通濾波等預(yù)處理前后的衛(wèi)星測(cè)高觀測(cè)值數(shù)目。從衛(wèi)星測(cè)高觀測(cè)值中獲取的大地水準(zhǔn)面高和垂線偏差是計(jì)算海洋重力異常的重要起算數(shù)據(jù)，其中相鄰觀測(cè)值之間差分計(jì)算垂線偏差的過程可以有效抑制軌道徑向誤差及其他長(zhǎng)波誤差項(xiàng)。因此，無需進(jìn)行交叉點(diǎn)平差，之前的數(shù)值測(cè)試表明該方法與大地水準(zhǔn)面解算方法解算精度相當(dāng) [64]。

我們對(duì)表3所述的多顆衛(wèi)星測(cè)高數(shù)據(jù)進(jìn)行聯(lián)合處理，獲取垂線偏差信息，進(jìn)而計(jì)算海域重力異常[45]。這些處理包括：首先采用兩步擬合法對(duì)多個(gè)測(cè)高衛(wèi)星的原始波形進(jìn)行重跟蹤處理[45]，并將沿軌觀測(cè)值統(tǒng)一重采樣至5 Hz，提高觀測(cè)值的精度和有效觀測(cè)值數(shù)量；然后顧及數(shù)據(jù)產(chǎn)品自帶的路徑延遲及地球物理環(huán)境改正項(xiàng)信息，獲取沿軌海表面高度和梯度信息，基于EGM2008模型梯度和3倍標(biāo)準(zhǔn)差準(zhǔn)則剔除粗差；顧及沿軌差分過程中會(huì)放大高頻噪聲，采用Parks-McClellan低通濾波器濾波得到梯度觀測(cè)值；分別使用DOT2008A模型和EGM2008模型考慮海面地形效應(yīng)和大地水準(zhǔn)面長(zhǎng)波信號(hào)，結(jié)合測(cè)高衛(wèi)星地面軌跡處經(jīng)向、緯向速度的近似計(jì)算公式，得到“移去”參考場(chǎng)貢獻(xiàn)的沿軌殘余垂線偏差信息，并計(jì)算格網(wǎng)點(diǎn)處的殘余垂線偏差方向分量；根據(jù)拉普拉斯方程導(dǎo)出的重力異常與垂線偏差的關(guān)系，進(jìn)一步計(jì)算得到格網(wǎng)化殘余重力異常；最后，“恢復(fù)”EGM2008參考場(chǎng)的長(zhǎng)波信號(hào)后就得到了全球1′×1′分辨率的海洋重力異常數(shù)據(jù)。

為了驗(yàn)證解算結(jié)果可靠性，采用美國國家地球物理數(shù)據(jù)中心（National Geophysical Data Center, NGDC）提供的船測(cè)重力數(shù)據(jù)進(jìn)行精度分析，并基于水深信息將檢核數(shù)據(jù)劃分為淺水區(qū)、非淺水區(qū)和深水區(qū)三類，分別與DTU10、DTU13和SS V23.1進(jìn)行對(duì)比，如表4所示。結(jié)果表明，本文計(jì)算的海洋重力異常數(shù)據(jù)精度是可靠的，在非淺水區(qū)和開闊海域精度與同期的DTU13 and SS V23.1模型相當(dāng)，優(yōu)于早期的DTU10模型。

4.3. GOCE和GRACE衛(wèi)星觀測(cè)法方程與重力異常觀測(cè)法方程的聯(lián)合解算

衛(wèi)星重力觀測(cè)數(shù)據(jù)和重力異常數(shù)據(jù)對(duì)重力場(chǎng)不同頻段信號(hào)的敏感度不同，因此不同的數(shù)據(jù)適于解算重力場(chǎng)的頻段也不相同，充分利用各類觀測(cè)數(shù)據(jù)包含的重力場(chǎng)信號(hào)，將多類觀測(cè)數(shù)據(jù)進(jìn)行最優(yōu)聯(lián)合解算是高精度高分辨率重力場(chǎng)模型解算的關(guān)鍵。通常認(rèn)為衛(wèi)星重力對(duì)重力場(chǎng)中長(zhǎng)波信號(hào)敏感，重力異常對(duì)高頻及超高頻信號(hào)敏感。本文采用法方程的聯(lián)合解算實(shí)現(xiàn)多類觀測(cè)數(shù)據(jù)的聯(lián)合求解，解算中認(rèn)為GRACE衛(wèi)星重力數(shù)據(jù)、GOCE衛(wèi)星觀測(cè)數(shù)據(jù)和重力異常數(shù)據(jù)三類數(shù)據(jù)之間不相關(guān)。圖3給出了衛(wèi)星觀測(cè)法方程和重力異常法方程的聯(lián)合求解策略，該策略參考了EGM96 [26]和EIGEN系列[19]模型解算的法方程聯(lián)合解算方法。

表3 恢復(fù)海域重力異常使用的衛(wèi)星測(cè)高數(shù)據(jù)

表4 在典型區(qū)域使用NGDC船測(cè)重力數(shù)據(jù)對(duì)模型進(jìn)行檢核的結(jié)果（單位：mGal）

GRACE觀測(cè)法方程、GOCE觀測(cè)法方程和重力異常法方程中位系數(shù)參數(shù)的最高階分別為200、220和2159。為提高計(jì)算效率，本文構(gòu)建重力異常法方程時(shí)使用了BDLS方法，因此重力異常法方程是塊對(duì)角形式的法方程。如圖3所示，這些法方程的聯(lián)合可以分為高階次部分和低階次部分，在高階次部分，模型的251～2159階ˉ系數(shù)是求解重力異常塊對(duì)角形式法方程得到的；而在低階次部分，模型的2～250階ˉ系數(shù)是通過聯(lián)合全階次的衛(wèi)星法方程和塊對(duì)角形式的重力異常法方程得到的。250階稱為過渡階，一般大于衛(wèi)星法方程的最高階，這樣可以使系數(shù)及其誤差在此處的變換較為平滑。假設(shè)模型的兩部分的系數(shù)是不相關(guān)的，才能將模型分為兩部分單獨(dú)解算，當(dāng)然，這樣的假設(shè)會(huì)帶來一些近似誤差。選取合適的過渡階會(huì)盡量減小該近似誤差，一方面選取的階次越高，通過多類觀測(cè)數(shù)據(jù)聯(lián)合求解得到的系數(shù)越多，但是低階次部分的聯(lián)合法方程的規(guī)模也會(huì)變大，相應(yīng)的計(jì)算量也會(huì)增大；另一方面，過渡階提高到一定的階，階次增大對(duì)模型的結(jié)果產(chǎn)生影響很小?？紤]本研究使用的計(jì)算平臺(tái)和計(jì)算精度的需要，本文選取的過渡階為250階。

圖3. GRACE、GOCE和重力異常法方程的聯(lián)合解算方法。

考慮到衛(wèi)星觀測(cè)數(shù)據(jù)求解的重力場(chǎng)低階部分系數(shù)的精度高于重力異常數(shù)據(jù)求解的系數(shù)，低階的聯(lián)合觀測(cè)法方程101階以下部分僅使用了衛(wèi)星法方程，沒有使用重力異常法方程。這個(gè)階次的選擇是通過分析衛(wèi)星模型ITSG-Grace2018、GOSG01S以及EGM2008模型的大地水準(zhǔn)面誤差得到的，EGM2008反映了由重力異常數(shù)據(jù)解算的系數(shù)精度。如圖4所示，EGM2008的大地水準(zhǔn)面誤差在100階以下的誤差均大于GOSG01S模型，在20～100階的誤差大于ITSG-Grace2018，從100階開始它們誤差的差異開始減小，所以認(rèn)為EGM2008在100階以內(nèi)的表現(xiàn)較差。由于本文使用的陸地重力異常數(shù)據(jù)是由EGM2008計(jì)算得到的，因此認(rèn)為重力異常求解的系數(shù)精度與EGM2008的特性一致，所以重力異常法方程沒有用于構(gòu)建聯(lián)合法方程的2～100階部分，而是用于構(gòu)建聯(lián)合法方程的101～250階部分。

4.4. SGG-UGM-2模型的解算流程

圖5給出了SGG-UGM-2模型的具體解算流程，模型解算中每一步的計(jì)算過程如下。

（1）根據(jù)4.1節(jié)的處理策略，使用GOCE衛(wèi)星的SGG和SST-hl觀測(cè)數(shù)據(jù)建立了220階次的GOCE衛(wèi)星觀測(cè)法方程，法方程的參數(shù)為系數(shù)。同時(shí)，解算該法方程得到純GOCE模型GOSG01S。ITSG-Grace2018模型法方程的原始格式為SINEX格式，對(duì)其進(jìn)行格式轉(zhuǎn)換，轉(zhuǎn)換后的格式與GOCE法方程一致。

（2）使用2.2節(jié)介紹的法方程參數(shù)轉(zhuǎn)換方法，將步驟（1）中構(gòu)建的兩個(gè)法方程轉(zhuǎn)換為參數(shù)是橢球諧系數(shù)的法方程。

（3）根據(jù)4.2節(jié)中的海洋重力異常數(shù)據(jù)求解方法，使用多顆測(cè)高衛(wèi)星的數(shù)據(jù)解算了5′×5′分辨率的海洋重力異常數(shù)據(jù)。

（4）根據(jù)式（27），使用EGM2008模型計(jì)算了參考橢球面上5′×5′分辨率的陸地重力異常點(diǎn)值數(shù)據(jù)，其中，GM和a參數(shù)分別 為3.986 004 415×1014m3·s-2和6 378 136.3 m，這與衛(wèi)星法方程中相應(yīng)參數(shù)一致。將步驟（3）計(jì)算的海洋重力異常數(shù)據(jù)和EGM2008計(jì)算的陸地重力異常數(shù)據(jù)融合得到全球重力異常數(shù)據(jù)，命名為ΔgF[65]。

圖4. EGM2008、GOSG01S和ITSG-Grace2018模型的大地水準(zhǔn)面階誤差。

（5）由數(shù)據(jù)ΔgF使用塊對(duì)角最小二乘方法構(gòu)建2～2159階系數(shù)塊對(duì)角形式的法方程，在計(jì)算中引入OpenMP并行計(jì)算技術(shù)提高計(jì)算效率，同時(shí)求解該塊對(duì)角法方程得到一組2～2159階的系數(shù)。為滿足塊對(duì)角最小二乘方法的要求，用單位陣作為數(shù)據(jù)的權(quán)陣；根據(jù)EGM2008模型在2190階的大地水準(zhǔn)面累積誤差，選取(5 mGal)2為法方程的初始方差分量。

（8）按照步驟（7）所述方法求得低階次部分的系數(shù)后，使用2.2節(jié)中的VCE方法更新各個(gè)法方程的相對(duì)權(quán)。然后使用新的相對(duì)權(quán)對(duì)步驟（7）和步驟（8）進(jìn)行迭代計(jì)算，直到法方程的相對(duì)權(quán)收斂。在SGG-UGM-2的計(jì)算中，最終進(jìn)行了4次迭代，最終得到的重力異常數(shù)據(jù)、GOCE數(shù)據(jù)和GRACE數(shù)據(jù)法方程的相對(duì)權(quán)分別為1.0、1.9和1.3。

（10）使用式（6）和式（7）將步驟（9）中得到的2～2159階橢球諧系數(shù)轉(zhuǎn)換為2～2190階0～2159次的球諧系數(shù)系數(shù)為擾動(dòng)位的球諧系數(shù)，加上GRS80橢球的正常位系數(shù)后就可以得到最終解算的SGG-UGM-2模型系數(shù)。

圖5. SGG-UGM-2模型的計(jì)算流程。

5. SGG-UGM-2模型的精度分析

5.1. 與EGM2008在頻域和空域的對(duì)比

為了分析SGG-UGM-2模型系數(shù)的精度，首先計(jì)算了它與EGM2008的系數(shù)差異，圖6中的紅色曲線為差異的階RMS，圖6也給出了另外兩個(gè)高分辨率模型EIGEN-6C4和GECO與EGM2008的差異，圖7給出了模型SGG-UGM-2、EIGEN-6C4和GECO與EGM2008系數(shù)階差異的譜圖。由圖6和圖7可知，在160階以下，這3個(gè)模型與EGM2008的差異非常接近，尤其是SGG-UGM-2和EIGEN-6C4，這是由于這3個(gè)模型都使用了GOCE衛(wèi)星數(shù)據(jù)。從160階開始，這些模型的差異開始變大，這是由于這些模型使用了不同的重力異常數(shù)據(jù)和不同的數(shù)據(jù)聯(lián)合解算方法。本文使用的海域重力異常數(shù)據(jù)是自主計(jì)算的，而GECO使用的是EGM2008格網(wǎng)重力異常，EIGEN-6C4使用的重力異常也與EGM2008重力異常接近，因此相較SGG-UGM-2，EIGEN-6C4在360階以上的系數(shù)與EGM2008更為接近，GECO模型360階以上的系數(shù)與EGM2008幾乎一致。

圖6. 模型SGG-UGM-2、EIGEN-6C4和GECO與EGM2008的系數(shù)階差異RMS。

圖7. 模型SGG-UGM-2 (a)、EIGEN-6C4 (b)和GECO (c)與EGM2008的系數(shù)階差異譜圖。

此外，使用SGG-UGM-2計(jì)算了全球重力異常數(shù)據(jù)，并計(jì)算了它與EGM2008模型重力異常和EIGEN-6C4重力異常的差異，差異的空間分布如圖8所示，可知SGGUGM-2與EGM2008差異主要體現(xiàn)在青藏高原、南美洲、非洲中部等區(qū)域，這些區(qū)域的重力異常數(shù)據(jù)都很稀少，SGG-UGM-2與EIGEN-6C4在這些區(qū)域的差異相對(duì)較小，這反映了GOCE數(shù)據(jù)對(duì)SGG-UGM-2與EIGEN-6C4的貢獻(xiàn)，GOCE數(shù)據(jù)改善了這兩個(gè)模型在重力稀疏區(qū)域的精度。除這些區(qū)域之外，SGG-UGM-2與其他兩個(gè)模型在沿海區(qū)域的差異也較大，造成該差異的原因是SGG-UGM-2解算中使用了與其他兩個(gè)模型不同的海域重力異常數(shù)據(jù)。

5.2. 中國和美國的GPS/水準(zhǔn)數(shù)據(jù)的檢核結(jié)果

首先使用了中國大陸區(qū)域的649個(gè)GPS/水準(zhǔn)點(diǎn)[66]和美國區(qū)域的6169個(gè)水準(zhǔn)點(diǎn)[67]檢驗(yàn)了SGG-UGM-2模型的精度，需要說明的是，中國區(qū)域GPS/水準(zhǔn)數(shù)據(jù)的高程基準(zhǔn)是正常高，美國區(qū)域GPS/水準(zhǔn)數(shù)據(jù)的高程系統(tǒng)為正高。使用SGG-UGM-2、SGG-UGM-1、EGM2008、GECO和EIGEN-6C4模型計(jì)算了中國區(qū)域的似大地水準(zhǔn)面高和美國區(qū)域的大地水準(zhǔn)面高，通過與GPS/水準(zhǔn)數(shù)據(jù)對(duì)比，計(jì)算模型的（似）大地水準(zhǔn)面誤差，計(jì)算中使用的潮汐均為無潮汐系統(tǒng)[67-69]。He等[69]的分析表明美國區(qū)域的GPS/水準(zhǔn)數(shù)據(jù)在東西方向有約70 cm的傾斜，而中國區(qū)域數(shù)據(jù)在東西方向的傾斜為9 cm。表5和表6給出了模型的（似）大地水準(zhǔn)面高在美國區(qū)域和中國區(qū)域誤差的統(tǒng)計(jì)結(jié)果，在美國區(qū)域和中國區(qū)域誤差的直方圖分別為圖9和圖10，這些差異沒有進(jìn)行任何的改正。

根據(jù)表5、表6、圖9和圖10，通過分析模型誤差的STD可知，在美國，EGM2008、SGG-UGM-1、SGG-UGM-2、GECO和EIGEN-6C4等模型的精度非常接近，模型誤差STD之間的差異小于7 mm，這些模型誤差的直方圖的分布形態(tài)也很接近，展示了較大誤差離散度，這可能是因?yàn)槊绹腉PS/水準(zhǔn)數(shù)據(jù)在不同區(qū)域的質(zhì)量不一致引起的。SGG-UGM-2在美國的精度最好，誤差的STD為0.277 m，而EGM2008的精度最差，誤差的STD為0.284 m。然而，這些模型在中國區(qū)域的檢核結(jié)果差異較大，模型誤差STD不盡相同，STD為0.157～0.240 m。使用了GOCE數(shù)據(jù)的SGG-UGM-1、SGG-UGM-2、GECO和EIGEN-6C4模型的精度較為接近，EIGEN-6C4的檢核結(jié)果最好，它的誤差STD為0.157 m，而EGM2008的檢核結(jié)果最差，它的誤差STD為0.240 m。從圖10可以看出，使用了GOCE數(shù)據(jù)的模型在中國的誤差直方圖也相對(duì)接近，EGM2008誤差直方圖的離散度最高，它的特性與其他幾個(gè)模型明顯不同，這可能是因?yàn)镋GM2008構(gòu)建中使用GOCE觀測(cè)數(shù)據(jù)，且使用中國區(qū)域的重力數(shù)據(jù)較為稀疏，而其他幾個(gè)模型通過使用GOCE衛(wèi)星數(shù)據(jù)提高了模型長(zhǎng)波系數(shù)的精度，改善了模型在中國的精度，所以這些模型在中國區(qū)域的精度優(yōu)于EGM2008。另外，SGG-UGM-1和SGG-UGM-2模型在中國和美國的檢核結(jié)果都較為理想。同時(shí)，對(duì)比這兩個(gè)模型可知，由于SGG-UGM-2中使用了最新的GRACE衛(wèi)星法方程和海域重力異常數(shù)據(jù)，它在美國和中國區(qū)域的精度都優(yōu)于SGG-UGM-1。EIGEN-6C4和SGG-UGM-2使用的觀測(cè)數(shù)據(jù)和模型解算方法比較接近，所以它們?cè)趦蓚€(gè)區(qū)域的精度也比較接近。兩者的不同之處在于，EIGEN-6C4使用了更多的衛(wèi)星重力數(shù)據(jù)，如LAGEOS衛(wèi)星數(shù)據(jù)，以及采用了不同于本文的海洋重力異常數(shù)據(jù)；SGG-UGM-2解算中各類數(shù)據(jù)的相對(duì)權(quán)是使用VCE方法計(jì)算的，而EIGEN-6C4使用的相對(duì)權(quán)是通過經(jīng)驗(yàn)方法確定的[12,19]。

圖8. SGG-UGM-2模型重力異常與EGM2008 (a)和EIGEN-6C4 (b)重力異常差異的空間分布。

表5 模型大地水準(zhǔn)面與美國GPS/水準(zhǔn)數(shù)據(jù)差異的統(tǒng)計(jì)結(jié)果（6169個(gè)點(diǎn)）（單位：m）

表6 模型似大地水準(zhǔn)面與中國GPS/水準(zhǔn)數(shù)據(jù)差異的統(tǒng)計(jì)結(jié)果（649個(gè)點(diǎn)）（單位：m）

圖9. EGM2008 (a)、 EIGEN-6C4 (b)、SGG-UGM-1 (c)、SGG-UGM-2 (d)和GECO (e)模型與美國區(qū)域GPS/水準(zhǔn)數(shù)據(jù)差異的直方圖。

同時(shí)，本文也使用了這些模型大地水準(zhǔn)面誤差的變異函數(shù)檢核模型精度，模型誤差的變異函數(shù)反映了模型誤差方差與GPS/水準(zhǔn)點(diǎn)間距離的關(guān)系。本文使用了文獻(xiàn)[70]中經(jīng)驗(yàn)變異函數(shù)的計(jì)算方法，通過文獻(xiàn)[71,72]可知，“gammavariance” (GVR)代表了在給定的距離上模型誤差的方差。圖11給出了模型在中國和美國區(qū)域誤差的經(jīng)驗(yàn)變異函數(shù)。在中國區(qū)域，EGM2008和GECO的經(jīng)驗(yàn)變異函數(shù)與其余3個(gè)模型有明顯的差異，其在中國的GVR都較大。SGG-UGM-1、SGG-UGM-2和EIGEN-6C4在兩個(gè)區(qū)域的GVR都接近。此外，在3000 km距離上，EIGEN-6C4模型在美國和中國的GVR均最小，這在一定程度上說明它在該波段的系數(shù)精度最優(yōu)。

5.3. 青島和臺(tái)灣的GPS/水準(zhǔn)數(shù)據(jù)的檢驗(yàn)結(jié)果

為了檢驗(yàn)?zāi)Ｐ驮谘睾^(qū)域和海島的精度，使用了中國青島和中國臺(tái)灣兩個(gè)區(qū)域的GPS/水準(zhǔn)數(shù)據(jù)對(duì)模型進(jìn)行檢驗(yàn)，這兩個(gè)區(qū)域都是沿海區(qū)域。青島和臺(tái)灣的GPS/水準(zhǔn)數(shù)據(jù)中分別有152個(gè)點(diǎn)和88個(gè)點(diǎn)。表7和表8給出了模型在青島和臺(tái)灣的GPS/水準(zhǔn)檢驗(yàn)的結(jié)果，圖12和圖13是在青島和臺(tái)灣區(qū)域模型大地水準(zhǔn)面誤差的直方圖，圖14是模型在這兩個(gè)區(qū)域誤差的經(jīng)驗(yàn)變異函數(shù)。青島和臺(tái)灣的GPS/水準(zhǔn)數(shù)據(jù)的高程系統(tǒng)分別為正常高系統(tǒng)和正高系統(tǒng)，所以分別計(jì)算了模型在兩個(gè)區(qū)域的似大地水準(zhǔn)面高和大地水準(zhǔn)面高。

由表7和表8可知，在青島區(qū)域SGG-UGM-2模型誤差的STD小于SGG-UGM-1，但在臺(tái)灣區(qū)域卻大于SGGUGM-1。一方面，青島區(qū)域的結(jié)果表明SGG-UGM-2中使用了新的海洋重力異常數(shù)據(jù)提升了其在沿海區(qū)域的精度；但臺(tái)灣區(qū)域的精度卻沒有提高，這可能是由于SGUGM-1和SGG-UGM-2使用的陸地重力異常是EGM2008計(jì)算的，而EGM2008模型構(gòu)建時(shí)使用了臺(tái)灣區(qū)域的精度很好的地面重力數(shù)據(jù)，新的衛(wèi)星重力數(shù)據(jù)和海域重力數(shù)據(jù)不能進(jìn)一步提高其在臺(tái)灣的精度，這也造成了SGG-UGM-1、EIGEN-6C4和GECO在臺(tái)灣的精度均不如EGM2008。在青島，這些模型誤差的直方圖比較接近，這與表7中的統(tǒng)計(jì)結(jié)果一致。而在臺(tái)灣這些模型的直方圖有明顯的不同，這是由于臺(tái)灣的地形較為復(fù)雜造成的。與表7和表8類似，在圖14中，EGM2008的經(jīng)驗(yàn)變異函數(shù)最小，尤其是在80～140 km的空間分辨率上，近似對(duì)應(yīng)于模型的140～250階次，從采用了GOCE數(shù)據(jù)的超高階重力場(chǎng)模型在海洋和美國本土的區(qū)域檢驗(yàn)可以說明，這部分頻段的主要貢獻(xiàn)來自于地面重力數(shù)據(jù)，因此GOCE衛(wèi)星重力和新導(dǎo)出的海洋重力異常對(duì)模型在臺(tái)灣區(qū)域精度的提升無實(shí)質(zhì)性貢獻(xiàn)。

圖10. EGM2008 (a)、 EIGEN-6C4 (b)、SGG-UGM-1 (c)、SGG-UGM-2 (d)和GECO (e)模型與中國區(qū)域GPS/水準(zhǔn)數(shù)據(jù)差異的直方圖。

圖11. EGM2008、EIGEN-6C4、GECO、SGG-UGM-1和SGG-UGM-2模型與中國（a）和美國（b）區(qū)域GPS/水準(zhǔn)差異的經(jīng)驗(yàn)變異函數(shù)。

表7 模型似大地水準(zhǔn)面與青島GPS/水準(zhǔn)數(shù)據(jù)差異的統(tǒng)計(jì)結(jié)果（152個(gè)點(diǎn)）（單位：m）

表8 模型大地水準(zhǔn)面與臺(tái)灣GPS/水準(zhǔn)數(shù)據(jù)差異的統(tǒng)計(jì)結(jié)果（88個(gè)點(diǎn)）（單位：m）

圖12. EGM2008 (a)、 EIGEN-6C4 (b)、SGG-UGM-1 (c)、SGG-UGM-2 (d)和GECO (e)模型與青島GPS/水準(zhǔn)數(shù)據(jù)差異的直方圖。

圖13. EGM2008 (a)、 EIGEN-6C4 (b)、SGG-UGM-1 (c)、SGG-UGM-2 (d)和GECO (e)模型與臺(tái)灣GPS/水準(zhǔn)數(shù)據(jù)差異的直方圖。

圖14. EGM2008、EIGEN-6C4、GECO、SGG-UGM-1和SGG-UGM-2模型與青島（a）和臺(tái)灣（b）GPS/水準(zhǔn)差異的經(jīng)驗(yàn)變異函數(shù)。

6. 結(jié)論

本文深入研究了EHA-CT方法的基本原理，并給出了實(shí)用化解算策略；推導(dǎo)了根據(jù)EHA-CT方法使用橢球面重力異常均值數(shù)據(jù)和點(diǎn)值數(shù)據(jù)反演重力場(chǎng)模型的嚴(yán)密計(jì)算公式，在由橢球面點(diǎn)值計(jì)算重力場(chǎng)模型的公式中引入了DH采樣和賦權(quán)理論[36]，數(shù)值試驗(yàn)結(jié)果表明，本文推導(dǎo)公式是嚴(yán)密的，可用于重力場(chǎng)模型的解算；對(duì)Rapp和Pavlis [33]實(shí)現(xiàn)EHA-CT方法的積分公式進(jìn)行了驗(yàn)證，數(shù)值試驗(yàn)表明該式不正確，對(duì)重力場(chǎng)模型的中長(zhǎng)波分量影響較大。

同時(shí)，聯(lián)合GOCE衛(wèi)星的SGG和SST-hl觀測(cè)數(shù)據(jù)、ITSG-Grace2018模型的法方程、多代測(cè)高衛(wèi)星數(shù)據(jù)和EGM2008重力異常數(shù)據(jù)解算了一個(gè)新的5′×5′分辨率的地球重力場(chǎng)模型SGG-UGM-2，模型完全展開到2190階2159次。在頻域和空域內(nèi)與EGM2008的對(duì)比分析說明，SGG-UGM-2模型的精度可靠。中國和美國GPS/水準(zhǔn)數(shù)據(jù)檢核的模型大地水準(zhǔn)面誤差的統(tǒng)計(jì)結(jié)果、直方圖和經(jīng)驗(yàn)變異函數(shù)表明，在中國和美國區(qū)域SGG-UGM-2優(yōu)于SGG-UGM-1、GECO和EGM2008模型，在美國區(qū)域SGG-UGM-2模型的精度最優(yōu)。由于SGG-UGM-2模型解算中使用了新的GRACE法方程和新反演的海域重力異常數(shù)據(jù)，因此SGG-UGM-2模型在中國大陸和美國的精度均略優(yōu)于SGG-UGM-1。這些結(jié)果均說明SGG-UGM-2模型精度是可靠的，可用于相關(guān)的地球科學(xué)研究和工程應(yīng)用；同時(shí)，研究結(jié)果也表明解算SGG-UGM-2模型使用的方法是嚴(yán)密的，是后續(xù)研制和發(fā)布SGG-UGM系列高分辨率重力場(chǎng)模型的重要基礎(chǔ)。

致謝

感謝Mayer-Gürr和Kvas提供ITSG-Grace2018模型的法方程。本研究獲得國家自然科學(xué)基金創(chuàng)新研究群體項(xiàng)目（41721003）和面上項(xiàng)目（41574019、41774020）資助，并獲德國DAAD主題網(wǎng)絡(luò)項(xiàng)目（57421148）、高分辨率對(duì)地觀測(cè)系統(tǒng)重大專項(xiàng)、中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金（N170103009）資助。我們也感謝編輯和匿名審稿人的建設(shè)性意見，幫助我們提高了稿件的質(zhì)量。

Compliance with ethics guidelines

Wei Liang, Jiancheng Li, Xinyu Xu, Shengjun Zhang,and Yongqi Zhao declare that they have no con fl ict of interest or fi nancial con fl icts to disclose.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2020.05.008.