利用重力场恢复与气候试验卫星(Gravity Recovery and Climate Experiment, GRACE)及其后续实验卫星(GRACE Follow-on)^[1-2]，可以全天候探测地球重力场中长波分量及其随时间的变化。其时变重力场可以用于反演全球季节性和年际气候变化伴随的地表多相态水体(冰川、陆地水、非潮汐海水质量)储量的变化^[3-4]，推动水文大地测量学的发展；也可用于固体地球过程的反演(例如大地震位错变形^[5]和冰川均衡调整^[6])，推动了地球动力学的研究。

近20 a卫星重力技术的发展为全球地表质量变化的研究提供了机会。GRACE双星采用近极轨道，轨道高度485 km，双星相距约220 km，用微波干涉仪测量双星间的距离及其速度变化，用高轨全球导航卫星系统(Global Navigation Satellite System, GNSS)卫星确定轨道位置和轨道速度，由加速度计测量非保守力影响。在GRACE双星的运行期间(自2002年3月至2017年7月)，有效累积超过14 a的月重力场时间序列。根据GRACE观测解算的时变重力场反演的多相态水体储量变化的空间分辨率达300 km，探测精度达到10 cm等效水柱或1 cm/a趋势速率。GRACE Follow-on卫星已发射，设计运行寿命5 a，目标大于10 a，双星轨道高更低(250 km)，双星距离更短(约50 km)，且使用激光干涉测距仪，能够较精确地测量双星间距离和星间速度，地球重力场观测精度比GRACE至少提高10倍，反演的多相态水体储量变化的空间分辨率达100 km，精度将大有提高。Swarm是为观测地球磁场设计的组合星座，由3颗位于不同高度的极轨卫星组成(一颗轨道高530 km，两颗450 km)，该星座2013年11月发射，迄今状态良好。由于携带加速度计、激光测距仪和GNSS接收机，能够恢复地球时变重力场，可弥补GRACE卫星及Follow-on卫星计划之间的重力观测空白。

GRACE时变重力已较成功应用于反演地表多相态水体储量变化和固体地球动力学过程，在陆地水或地下水储量变化方面，研究了印度北部^[7]、华北^[8]、北美、北欧等冰川均衡调整(glacial isostatic adjustment，GIA)近场区^[9-10]、青藏高原及其周边地区^[11]、三峡地区^[12]、加州山谷地区^[13]、全球^[14]、亚马逊流域和赞比西河流域^[15]等的陆地水储量变化；在冰川/冰盖消融量方面，研究了南极冰盖^{[16, 4]}和高亚洲冰川^[17]的质量变化；在固体地球动力学过程方面，主要有GIA^[18-20]、大地震同震变形^[5]的研究；此外还有质量海平面变化方面的研究^[21]等。

一般用地表质量异常表示多相态水体储量变化。必须注意，由于卫星重力观测采用地球系统质心(center of mass, CM)坐标，因此其观测的一阶斯托克斯系数为零，但这并不意味着地表质量异常没有一阶项，因为地表水文、海洋观测和模型是相对地球外表面的构形中心(center of figure, CF)坐标框架下的结果。此外，二阶零次斯托克斯系数可能受潮汐混频影响，具有很大的不确定性^[22]。所以，实践中使用其他方法独立估计一阶、二阶斯托克斯系数，用以补充和替换GRACE时变重力模型的结果^[23-24]。

忽略一阶项贡献对高纬度质量异常变化和大尺度流域-海洋质量交换的估计具有较大影响^[25]。例如，单独基于GRACE与基于测高和比容估计的平均海平面变化的季节项振幅相差15%，而补充GRACE的一阶斯托克斯系数后，其差异减少到1%^[26]；若忽略其质量异常一阶项趋势影响，可使基于GRACE的南北极和全球山地冰盖/冰川消融引起的海平面上升低估30%^[27]。

最近，独立估计时变重力场一阶、二阶斯托克斯系数越来越受到关注。继Swenson等^[23]提出GRACE-OBP法估计一阶斯托克斯系数、Cheng等^[24]用卫星激光测距(satellite laser ranging, SLR)资料估计二阶斯托克斯系数后，Sun等^[28-30]考虑陆地交换水自吸引和旋转反馈效应，利用时变重力模型的协方差信息，进一步完善一阶、二阶斯托克斯系数的估计，给出了另一独立的结果。此外，Chanard等^[31]利用GRACE二阶及以上斯托克斯系数，结合Swenson模型一阶斯托克斯系数结果^[23]，在全球689个连续GNSS观测站点上，计算由全球地表多相态水体储量变化(其中也包含非潮汐大气)引起的径向和水平位移，发现水平位移与GNSS观测的振幅和相位结果有较大差异，这在一定程度说明Swenson模型可能有较大误差。

因此，本文重新审视时变重力场一阶、二阶斯托克斯系数的估计方法，特别是对比不同模型及其转换的质量异常结果，为当前GRACE和下一代卫星重力的成功应用，具体如何进行一阶、二阶斯托克斯系数补充及替换，提供参考。

1 一阶、二阶斯托克斯系数的反演方法

1.1   时变重力场的球谐展开式

地球重力场通常由大地水准面的形状描述，而大地水准面对应全球海洋平均海平面的等位面，其形状可用球谐展开表示^[32]。本文所关注的时变重力场可以理解为不同时间大地水准面的变化，或同一时间相对平均大地水准面的变化。因此，在地表及其外部任意点(θ, φ)，地球时变重力场由大地水准面形状变化表示，其球谐展开式为：

式中，(ΔC′_lm, ΔS′_lm)为由卫星重力观测的l阶m次斯托克斯系数的变化(简称斯托克斯系数)；${{\tilde P}_{lm}}\left( {\cos \theta } \right)$为l阶m次的连带勒让德多项式；a=6 378.136 3 km，为赤道半径。由于卫星重力观测采用CM坐标，故一阶斯托克斯系数为零。此外，由于地球系统质量守恒，零阶斯托克斯系数也为零。受分辨率限制，卫星重力只能提供有限阶次的斯托克斯系数的结果，例如，当前德克萨斯大学空间研究中心(the Center of Space Research, CSR)重力模型最大可达96阶。

采用CF为原点的笛卡尔坐标框架，则式(1)的一阶斯托克斯系数不为零，依此非零系数可求得CM相对于CF沿3个坐标轴向的位移^[33]：

根据时变重力场二阶、零次斯托克斯系数，可以简单求得时变引力位二阶带谐系数：

式(3)转换便于与SLR观测结果进行对比。

1.2   全球地表质量异常的球谐展开式

首先，假设时变重力场的一阶和二阶零次斯托克斯系数已经补充和替换，其他斯托克斯系数由GRACE等卫星时变重力模型(GSM+GAC)给出；其次，所有斯托克斯系数已经过固体地球过程影响的改正，即去除GIA、大地震位错变形的影响，则改正后的斯托克斯系数不仅反映了地表多相态水体(也包含非潮汐大气)储量变化相关的地表质量异常的直接牛顿引力贡献, 也反映了负荷形变激发的地球内部质量重新分布贡献。地表质量异常用地表面密度表示，根据文献[34]，地表面密度球谐展开形式为：

式中，地表面密度球谐系数由重力场斯托克斯系数转换求得：

式中，k_l为l阶与位扰动相关的弹性负荷勒夫数^[35-36]；地球平均密度ρ_ave= 5.517 g/cm³。地表面密度除以水密度(1.0 g/cm³)即是等效水柱高。

值得说明的是，当前该负荷勒夫数一般采用根据参考地球模型PREM(preliminary reference earth model) ^[37]计算的结果(如Wang等^[36]的结果)，但对其一阶负荷勒夫数则有特殊的考虑，即取CF坐标框架下的结果：

式中，Blewitt^[38]采用Farrell^[39]根据古登堡-布伦地球模型计算的地球中心(center of earth, CE)为原点坐标框架下的勒夫数h₁=-0.290和l₁=0.113，计算出CF框架下k₁=0.021；根据PREM^[37]的CE框架下的h₁=-0.285 7和l₁=0.103 6^[36]，得出CF框架下k₁=0.026 2，依此并根据式(5)计算的一阶项地表质量异常比k=0.021时小0.5%。然而，当前的相关研究文献和发布的一阶斯托克斯系数模型一般建议使用k₁=0.021，但本文采用相对较新的PREM地球模型，且为了与计算二阶及以上阶地表质量异常广泛采用的PREM地球模型保持一致，因此采用k₁=0.026 2，提供了一个与PREM模型一致的可替代的结果。

1.3   C₁₀、C₁₁、S₁₁和C₂₀的反演方法

为了方便，省去时变重力场斯托克斯系数和地表面密度的球谐系数中的Δ。人们知道，GRACE提供了3类时变重力场斯托克斯系数模型：①GSM模型主要反映陆地多相态水储量变化(不含非潮汐大气质量信号)，即反映土壤湿度、地表水、地下水、冰川、积雪、冻土变化，此外包含GIA效应、海洋的大陆径流补给的淡水质量重新分布(陆地交换水)；②GAC模型，模拟的全球非潮汐大气、海洋质量贡献；③GAD模型，是GAC去除陆地大气影响的结果，也即指模拟的非潮汐海洋质量和非潮汐海洋大气质量的贡献。

在大地水准面公式(1)中，使用GRACE的时变重力场斯托克斯系数GSM模型，用式(4)表达全球地表面密度。取全球地表面密度海洋中的结果，即取反映陆地交换水的部分，取陆地结果为零，定义陆地交换水的球谐系数：

其中，dΩ=sinθdθdφ。

假设将陆地水总储量的亏损量补给到海洋，GRACE-OBP算法在海洋各处等量补给，而改进的GRACE-OBP算法则考虑自吸引效应和旋转反馈效应。因此，已知海洋的陆地交换水的球谐系数C₁₀^ocean、C₁₁^ocean、S₁₁^ocean和C₂₀^ocean，则根据式(7)和式(8)列出已知量和待求量(C₁₀、C₁₁、S₁₁和C₂₀)的关系^{[23, 28]}：

用矩阵方程表示为：

解方程(13)，得地表面密度的低阶球谐系数(C₁₀、C₁₁、S₁₁和C₂₀)，进一步通过式(14)转换得到地球时变重力场的GSM模型的相应斯托克斯系数：

如果时变重力场的斯托克斯系数需要包含全球非潮汐大气、海洋质量贡献，则在式(14)左边的结果中加上GRACE GAC模型的相应阶次斯托克斯系数。在后面的讨论中，忽略式(14)左边的“′”上标。

2 不同模型的C₁₀、C₁₁、S₁₁和C₂₀比较

在Chambers Model中，Swenson和Chambers用GRACE-OBP法^[23]，结合2002年7月至2016年12月二阶以上GRACE的CSR 60阶的斯托克斯系数月解，其中C₂₀由Cheng等^[24]给出的SLR结果替换，计算了一阶斯托克斯系数(C₁₀、C₁₁、S₁₁)月解。在Sun Model中，Sun等^[28-30]用改进的GRACE-OBP法, 结合二阶以上(C₂₀除外)GRACE的CSR 60阶同期斯托克斯系数月解、CSR 96阶斯托克斯系数月解、德国地学研究中心(German Research Center for Geoscien-ce, GFZ)的90阶斯托克斯系数月解和喷气推进实验室(Jet Propulsion Laboratory, JPL)的90阶斯托克斯系数月解, 分别计算了一阶斯托克斯系数(C₁₀、C₁₁、S₁₁)和C₂₀的月解。分别对所求月解的时间序列进行最小二乘回归分析，得到相应趋势速率、周年的振幅和相位及其各自回归的方差结果。

表 1给出了不同模型C₁₀的结果。对比Chambers Model和Sun Model/CSR/60(表示模型主要开发者为Sun、采用CSR 2~60阶的GRACE时变重力模型，下同)的结果，Chambers Model减少趋势速率比Sun Model增加27%，而周年振幅则减少36%，周年相位相差较小(即0.005 a)；对比Sun Model/CSR/60和Sun Model/CSR/96的结果，GRACE截断阶数增加36，对趋势信号和周年信号影响非常小；对比Sun Model/CSR/96、Sun Model/GFZ/90和Sun Model/JPL/90的结果，CSR与JPL趋势速率相差最大(达17%)，周年信号相差很小。

表  1  不同模型C₁₀的趋势速率、周年振幅和周年相位

Table  1.  Trend Rate, Annual Amplitude and Annual Phase of Different Models C₁₀

模型	趋势速率/10-10·a-1	周年振幅/10-9	周年相位/a
Chambers Model	-0.189 58±0.011 90	0.168 11±0.067 52	0.009±0.006
Sun Model/CSR/60	-0.150 36±0.016 22	0.254 78±0.009 11	0.014±0.006
Sun Model/CSR/96	-0.145 59±0.017 18	0.256 93±0.009 43	0.017±0.006
Sun Model/GFZ/90	-0.159 92±0.019 32	0.264 20±0.010 85	0.017±0.007
Sun Model/JPL/90	-0.182 92±0.016 38	0.257 33±0.009 21	0.004±0.006

下载: 导出CSV 

| 显示表格

表 2给出了C₁₁的结果。类似的比较表明，Chambers Model较Sun Model的减少趋势增加266%，而振幅则减少21%，周年相位相差较小(即0.002 a)；GRACE截断阶数增加36，对趋势信号影响为19%，周年信号影响非常小；CSR、GFZ与JPL趋势速率相差较大，CSR与JPL相差最大(达63%)，但周年信号相差都较小。

表  2  不同模型C₁₁的趋势速率、周年振幅和周年相位

Table  2.  Trend Rate, Annual Amplitude and Annual Phase of Different Models C₁₁

模型	趋势速率/10-11·a-1	周年振幅/10-9	周年相位/a
Chambers Model	-1.105 01±0.095 17	0.115 08±0.005 39	-0.012±0.008
Sun Model/CSR/60	-0.325 10±0.137 73	0.143 68±0.007 74	0.009±0.009
Sun Model/CSR/96	-0.269 83±0.141 31	0.141 60±0.007 75	0.011±0.009
Sun Model/GFZ/90	-0.501 81±0.153 78	0.145 68±0.008 64	0.005±0.010
Sun Model/JPL/90	-0.707 83 ±0.144 66	0.137 76±0.008 13	-0.001±0.010

下载: 导出CSV 

| 显示表格

表 3是S₁₁的结果。Chambers Model的减少趋势和周年振幅远远小于Sun Model的结果，周年相位则相差较小；GRACE截断阶数增加36, 对趋势信号影响为10%，周年信号影响则非常小；CSR、GFZ的结果非常相近，但它们与JPL的趋势速率相差较大(32%)，但周年信号相差很小。

表  3  不同模型S₁₁的趋势速率、周年振幅和周年相位

Table  3.  Trend Rate, Annual Amplitude and Annual Phase of Different Models S₁₁

模型	趋势速率/10-11·a-1	周年振幅/10-9	周年相位/a
Chambers Model	-0.088 99±0.108 34	0.128 47±0.006 14	0.477±0.008
Sun Model/CSR/60	0.557 83±0.126 42	0.149 95±0.007 11	0.431±0.008
Sun Model/CSR/96	0.614 25±0.137 43	0.153 61±0.007 55	0.424±0.008
Sun Model/GFZ/90	0.607 76±0.141 33	0.149 27±0.007 95	0.427±0.009
Sun Model/JPL/90	0.880 19±0.150 84	0.140 74±0.008 50	0.425±0.010

下载: 导出CSV 

| 显示表格

表 4是C₂₀的结果。类似的比较表明, Cheng Model的趋势速率比Sun Model/CSR/60小44%，周年振幅小20%，周年相位相差0.03 a；GRACE截断阶数增加36, 对趋势速率和周年信号影响均很小；CSR、GFZ和JPL的趋势和周年信号相差都很小。

表  4  不同模型C₂₀的趋势速率、周年振幅和周年相位

Table  4.  Trend Rate, Annual Amplitude and Annual Phase of Different Models C₂₀

模型	趋势速率/10-11·a-1	周年振幅/10-10	周年相位/a
Cheng Model	-0.106 67± 0.072 37	0.671 25± 0.041 82	0.025±0.010
Sun Model/CSR/60	-0.187 19±0.008 11	0.846 46± 0.045 62	0.054±0.009
Sun Model/CSR/96	-0.188 47±0.008 14	0.863 53± 0.044 71	0.052±0.008
Sun Model/GFZ/90	-0.190 11±0.008 30	0.847 95± 0.046 70	0.059±0.009
Sun Model/JPL/90	-0.186 42±0.081 91	0.858 61± 0.046 11	0.051±0.009

下载: 导出CSV 

| 显示表格

3 不同模型的C₁₀、C₁₁、S₁₁和C₂₀质量异常的比较

利用上述6个不同模型给出的C₁₀、C₁₁、S₁₁和C₂₀月解，通过式(4)分别计算全球地表质量异常一阶项和二阶项的时间序列，通过最小二乘回归分析，得到全球地表质量异常一阶项、二阶项的趋势速率、周年振幅和周年相位，一阶项的结果见图 1至图 3，二阶项的结果见图 4。

图  1  不同模型的地表质量异常一阶项趋势速率

Figure  1.  Trend Rates of Degree One Surface Mass Anomaly from Different Models

下载: 全尺寸图片 幻灯片

图  2  不同模型的地表质量异常一阶项周年振幅

Figure  2.  Annual Amplitudes of Degree One Surface Mass Anomaly from Different Models

下载: 全尺寸图片 幻灯片

图  3  不同模型的地表质量异常一阶项周年相位

Figure  3.  Annual Phases of Degree One Surface Mass Anomaly from Different Models

下载: 全尺寸图片 幻灯片

图  4  不同模型C₂₀的质量异常的趋势速率、周年振幅和周年相位

Figure  4.  Trend Rates, Annual Amplitudes and Annual Phases of Surface Mass Anomaly from Different Models C₂₀

下载: 全尺寸图片 幻灯片

图 1显示了5个模型的一阶项质量异常速率。可以看出，Chambers Model的一阶项质量异常速率基本绕零度经线对称分布，其正负极值位于(180°E，60°S)和(0°N，60°N)，而4个Sun Model则基本具有一致的一阶项质量异常速率分布形态，其正负极值一致位于(120°E，75°S)和(75°N，60°W)；Sun Model(图 1(b))与Chambers Model(图 1(a))的结果差异见图 1(f)，最大差异位于以红海为中心的欧亚非大陆和太平洋，幅值达0.6 mm/a；图 1(c)与图 1(b)的比较说明GRACE截断阶数增加36使极值范围有所缩小；图 1(c)、1(d)、1(e)的比较说明，CSR与GFZ极值相当(1.1 mm/a)，但JPL的极值增加到1.3 mm/a。子图顶部显示所用模型和高原内的数值范围。

图 2给出了5个模型的一阶项质量异常周年振幅。可以发现，Chambers Model的一阶项质量异常周年振幅具有与其他4个Sun Model完全一致的分布形态和最大值位置；Chambers Model周年振幅极值(图 2(a))比Sun Model(图 2(b))小26%，Sun Model与Chambers Model的结果差异见图 2(f)，最大差异位于两极地区，幅值达5.5 mm；图 2(c)与图 2(b)的比较显示GRACE截断阶数增加36不影响周年振幅；图 2(c)、2(d)、2(e)的比较说明CSR、GFZ与JPL的周年振幅相当(19.2 mm)。

图 3给出了5个模型的一阶项质量异常周年相位。可以看出，Chambers Model周年相位(图 3(a))与Sun Model(图 3(b))分布形态相似，Sun Model与Chambers Model的结果差异见图 3(f)，除非洲南端-南极-太平洋的窄条带外，其他地区差异一般小于0.02 a；GRACE截断阶数增加36不影响周年相位；CSR、GFZ与JPL的周年相位基本相同。

图 4给出了5个模型的二阶项质量异常的趋势速率、周年振幅和周年相位。可以看出，4个Sun Model的结果十分接近，完全无法分辨它们之间的差异；二阶项质量异常的趋势速率、周年振幅在35.26°N和35.26°S取零值，而相位则发生相差0.5 a的跳变；趋势速率(图 4(a))在赤道和两极取正负极值，Cheng Model的正负极值幅度分别比4个Sun Model小44%和43%；周年振幅(图 4(b))在赤道和两极取极值，Cheng Model的正负极值幅度分别比4个Sun Model低21%和21%；图 4(c)中，Cheng Model比4个Sun Model的相位延迟0.029 a。

4 结语

本文简略地介绍了GRACE-OBP算法及其改进的算法，有助于估计GRACE时变重力场的3个一阶斯托克斯系数(C₁₀、C₁₁、S₁₁)和C₂₀系数；并指出在全球地表质量异常的球谐表达式中，各项勒夫数需采用一致的地球模型，因一般较高阶项采用PREM的结果，故在CF框架下一阶勒夫数应取PREM模型的结果(0.026 2)。

本文运用基于GRACE-OBP算法的Chambers Model给出的一阶斯托克斯系数月解、基于SLR观测的Cheng Model给出的C₂₀月解和基于改进的GRACE-OBP算法的4个Sun Model给出的一阶斯托克斯系数月解和C₂₀月解，分别计算了趋势项、周年项，考察了传统方法较改进方法结果的差异、GRACE截断阶数增加36对结果的影响以及运用3个机构(CSR/GFZ/JPL)时变重力模型的结果差异。不同模型间存在显著差异：① Chambers Model C₁₀的减少趋势速率比Sun Model增加27%，周年振幅减少36%，CSR、GFZ的趋势速率与JPL的趋势速率相差达17%；② Chambers Model C₁₁的减少趋势速率比Sun Model大26%，振幅小21%，截断阶数增加36对趋势信号影响为19%，CSR、GFZ与JPL的趋势速率相差较大，最大达63%；③ Chambers Model S₁₁的减少趋势速率和周年振幅均远远小于Sun Model，截断阶数增加36对趋势信号影响为10%，CSR、GFZ与JPL的速率相差达32%；④ Cheng Model C₂₀的趋势速率比Sun Model小44%，周年振幅小20%。

运用上述6个不同模型给出的C₁₀、C₁₁、S₁₁和C₂₀月解，计算了全球地表质量异常一阶项和二阶项的趋势速率、周年振幅和周年相位，不同模型结果主要差异如下:①对一阶项质量异常趋势速率，Chambers Model与Sun Model具有不同的分布，二者之间沿纬线存在60°的平移，最大差异在欧亚非大陆和太平洋；②对一阶项质量异常周年振幅，Chambers Model与Sun Model最大差异位于两极地区，达5.5 mm，最大值小26%；③对二阶项质量异常趋势速率，Cheng Model的正负极值幅度比4个Sun Model小44%和43%；④对二阶项质量异常周年振幅，Cheng Model的正负极值幅度比4个Sun Model均低21%。

由上可知，基于GRACE-OBP算法的C₁₀、C₁₁、S₁₁和SLR观测的C₂₀月解所直接计算的趋势项和周年项，或计算的地表质量异常的一阶项和二阶项的趋势项和周年项，比用改进的GRACE-OBP算法的结果在趋势项上有很大差异，在周年项上差异相对较小。因此，在计算陆地水储量变化时，一定要合理选择一阶斯托克斯系数和C₂₀，本文建议使用基于改进的GRACE-OBP算法的Sun Model系列的结果。此外，运用GRACE模型的不同截断阶数和不同的模型发布机构对趋势项和周年项都有一定影响，但趋势项影响更大，所以必须选择最好的GRACE模型且具有较大的阶次。