不同地球模型对中国陆区大气负荷位移建模的影响

李长海， 钟萍， 姜中山， 汤苗， 杨兴海， 由晓文

西南交通大学地球科学与环境工程学院， 成都 611756

0 引言

环境负载如大气、海洋、陆地水等的变化会使固体地球表面产生明显的位移变化，并同步引起重力场变化.环境负载引起的地表形变的准确计算有利于相关地球物理信号的研究(Van Dam et al., 1994; 王敏等，2005; 袁林果等，2008).相比于其他环境负载，地表气压的观测数据、数值模式和同化方案都较为完善(赵天保等, 2010)，并且当前大气负荷效应的研究也日趋成熟.相关研究表明，大气负荷引起的垂向地表形变具有明显的季节性特征且在中国地区的变化幅度最大超过20 mm(张诗玉等, 2006; 周江存和孙和平, 2007; Yue et al., 2020)，对GPS站点位移进行大气负荷修正可以明显减小数据残差(Van Dam et al., 1994; Tregoning and Watson, 2009; 姜卫平等, 2016).因此，对于地壳构造运动的研究及大地参考框架的建立，精确改正大气负荷效应是十分必要的.

负荷位移建模方法主要有格林函数法和球谐函数法，两者在数学上是等价的，其精度在误差水平上也是一致的(沈迎春等, 2017).格林函数法是应用卷积计算全球或区域格点负荷对测站的影响，在负荷位移建模中被广泛采用(Farrell, 1972; van Dam et al., 1994; Dill et al., 2015; Li et al., 2020).早期研究采用一维模型来表征地球内部结构及密度和速度等特征(Dziewonski and Anderson, 1981)，当前进行负荷位移建模也普遍采用一维地球模型.但是随着理论研究的深入和观测精度的不断提高，建立更加完整和精确的负荷模型愈发重要.冯锐(1985)的研究表明中国陆区地壳厚度分布不均，东部变化平缓，西部变化复杂.基于对各种地球物理数据的统计分析，Laske等(2013)提出了三维地球模型CRUST1.0，该模型最深可达上地幔，考虑了地壳的区域差异.Wang等(2012)计算了不同地球模型(PREM、iasp91、ak135+crust2.0)的负荷勒夫数，这些勒夫数在超过200阶后差异较大.他们认为这可能是由于模型参数在220 km深度处不连续和不同模型在地壳中弹性参数差异的影响所致，并建议当使用格林函数角距Θ<1°时，用CRUST2.0的地壳结构代替PREM或ak135的上层.贾路路等(2014)的研究表明，全球平均的地球模型不能很好地反映中国陆区特别是青藏高原地区的地壳结构.Dill等(2015)采用局部负荷格林函数来计算负荷形变，结果表明，区域地壳结构对负荷位移建模的影响不容忽视.

中国陆区的地壳结构复杂，是研究地球模型差异对负荷建模影响的天然试验场.为了研究地球模型差异对利用地表气压数据估算中国陆区地表垂向形变的影响，本文基于不同地球模型计算了对应的负荷格林函数，并与全球地表气压再分析数据估算中国陆区大气负荷垂向形变，对比分析不同地球模型的负荷位移建模差异.结合研究区域内GPS参考站的去线性时间序列数据，根据大气负荷修正后GPS残差WRMS的改正情况来评估不同地球模型的建模差异.

1 数据

1.1 地球模型

本文讨论了四种地球模型，包括一维地球模型PREM、STW105、AK135和三维地球模型CRUST1.0.以上地球模型的模型参数(不同地表深度下的纵波波速VP、横波波速VS和岩石密度ρ)由地震联合研究机构(Incorporated Research Institutions for Seismology，IRIS)提供.其中CRUST1.0模型由冰层、水层、沉积层和地壳层组成.在计算时，本文将CRUST1.0模型地幔层及以下的结构用PREM模型中的对应部分补充.不同地球模型的模型参数如图1所示，当深度超过1000 km后，不同地球模型的模型参数基本一致，为了便于展示，我们仅绘制了1000 km以上的部分.

图1 不同地球模型的模型参数纵坐标表示地球内部一点至地表的深度.横坐标表示该深度下的模型参数值.红、绿、蓝三种实线分别表示AK135、PREM和STW105模型的模型参数(包括纵波波速VP，横波波速VS和密度ρ)，黄色虚线表示CRUST1.0模型在中国陆区3°×3°各格网点的模型参数.Fig.1 Model parameters of different earth models The Y-axis represents the depth from a point inside the earth to the surface. The X-axis represents the model parameter value at that depth. The three solid lines of red, green and blue represent the model parameters of the AK135, PREM and STW105 models (including longitudinal wave velocity, shear wave velocity and density). The yellow dotted line represents the model parameters of the CRUST1.0 model at 3°×3° grid points in Chinese continent.

1.2 GPS残差时间序列

为了正确分析大气负荷效应对中国陆区GPS测站的影响，本文采用国家地震科学数据中心提供的单天解去线性项的精密坐标时间序列(http:∥www.eqdsc.com/data/pgv-sjxl.htm)，该数据使用麻省理工学院提供的GAMIT/GLOBK软件进行解算，采用IGS提供的精密星历信息，采用GMF模型改正先验对流层延迟，使用IERS03模型改正固体潮，使用otl_FES2004全球海潮模型改正海潮负荷效应.数据时间跨度为2011—2017年.由于本文获取的GPS数据已经扣除了部分地球物理信号(固体潮、海潮负荷)，因此文中将其称为GPS残差.为了获得更可靠的统计结果，本文剔除了GPS残差时间序列中超过三倍中误差的历元数据，筛选出了数据完整率大于70%的243个测站.

1.3 地表气压数据

本文采用美国国家环境预报中心(National Centers for Environmental Prediction, NCEP)和美国国家大气研究中心(National Center for Atmospheric Research, NCAR)提供的全球地表气压再分析数据NCEP R-1(https:∥psl.noaa.gov/data/gridded/data.ncep.reanalysis.surface.html)计算大气负荷引起的地表形变，数据时间跨度为2011—2017年.该数据的空间分辨率为2.5°×2.5°，时间分辨率为6 h，有3～5天的数据延迟.由于原GPS数据未扣除大气潮汐效应，为了保持一致，本文在计算气压变化量时未扣除地表气压数据中的潮汐项(如S1、S2).

2 计算方法

2.1 负荷格林函数计算

本文通过负荷格林函数与全球格网点负荷卷积计算得到大气负荷引起的地表形变.负荷格林函数表示单位点负荷变化引起的地表形变，其中垂向位移格林函数可以通过下式计算：

(1)

其中，Gr表示垂向位移格林函数，Θ表示负荷源与形变点之间的角距离，RE为地球半径，me为地球质量，n为计算阶数，hn为负荷勒夫数，Pn为勒让德多项式.以地表至6371 km深度的地球模型参数作为输入数据，本文使用Martens等(2019)提供的LoadDef程序中的run_ln计算了PREM、STW105、AK135和CRUST1.0模型10000阶次的负荷勒夫数，并将该结果作为输入数据，使用run_gf计算了相应的负荷格林函数.对于CRUST1.0地球模型，本文用CRUST1.0模型每个2.5°×2.5°网格点的地壳层参数代替PREM模型的地壳层参数，计算了全球陆地2.5°×2.5°的结果，即每个格网点的地球模型所对应的负荷勒夫数和负荷格林函数.

2.2 大气负荷位移计算

本文采用格林函数方法(Farrell, 1972; van Dam et al., 1994; Dill et al., 2015; Li et al., 2020)计算地表气压变化引起的垂向地表形变，公式如下：

(2)

其中，u表示形变点r′处的垂向位移，ΔP表示负荷源r处的气压变化量，g表示重力加速度，Ω表示积分区域.考虑到海洋的逆气压效应，需要对海洋区域的气压进行反变气压改正，本文采用Spiridonov和Vinogradova(2019)、Li等(2020)和van Dam and Ray(2010)的改正策略，即仅对陆地区域进行积分计算，并使用van Dam提供的海陆掩膜数据.在计算气压变化时，本文取20年(1991—2010年)的气压均值作为参考气压(van Dam and Ray, 2010).需要说明的是，一维地球模型的负荷格林函数仅与负荷点和形变点之间的角距有关，而用CRUST1.0模型替换PREM模型的地壳结构后，该地球模型的负荷格林函数还与负荷源所处的位置有关，本文用CRUST1.0地球模型的负荷格林函数进行负荷位移建模时，参照了Dill等(2015)的计算方法，公式如下：

(3)

与(2)式不同的是，(3)式中的Gr为负荷源r所处格网的局部负荷格林函数.在实际计算时，两种方法的原理是相同的.

2.3 负荷修正效果评估

本文采用WRMS分析法(Jiang et al., 2013)，通过扣除GPS时间序列中的大气负荷效应，统计扣除前后GPS残差的加权均方根值(weighted root mean square, WRMS)的相对变化量来评估大气负荷效应修正对测站坐标时间序列的修正效果.令ω(i)=1/sigu(i)2，WRMS相对变化量可通过下式计算：

WRMS(diff percentage)%=

(4)

WRMS(var)=

(5)

式中var(i)和sigu(i)分别表示i时刻的变量及其不确定度；Dgps表示GPS垂向残差；Dload表示大气负荷引起的垂向位移；ndat表示测站的观测值数量.计算不同地球模型的大气负荷效应改正所对应的WRMS值时只需替换对应的Dload分量即可.由式(4)—(5)可知，WRMS相对变化量大于0表示大气负荷校正可以降低GPS时间序列的RMS.

2.4 周年振幅估计

本文根据

y(ti)=a+bti+csin(2πti)+dcos(2πti)+vi,

(6)

采用最小二乘方法对周年振幅进行估计.式中，y表示垂向位移，ti为以年为单位的时间历元，a为常数项，b为线性速度项，c和d为周年运动的系数项，vi代表残差序列，A为拟合得到的周年振幅.

2.5 相关系数计算

本文采用Pearson相关系数(Correlation Coefficient,CC)来衡量不同数据之间的一致性，计算方法如下：

(7)

3 结果和讨论

3.1 计算结果验证

为了验证本文大气负荷位移计算结果的正确性，根据全球地球物理流体中心(Global geophysical fluids center, GGFC)提供的大气负荷位移格网数据(van Dam and Ray, 2010)，本文计算了与之对应的2.5°×2.5°中国陆区2011年大气负荷垂向位移的时间序列并将二者进行比较.由于GGFC计算大气负荷位移时用低通滤波扣除了地表气压数据的高频项，因此本文比较两种结果时使用的是日均值.由图2可以看出中国陆区大部分地区本文的计算结果和GGFC提供的产品精度相当.由于计算时采取的积分间隔不同，因此东部沿海地区的负荷位移建模结果存在略微差异，但并不影响不同地球模型间的比较.

图2 本文计算结果与GGFC提供的产品对比(a) 本文用PREM模型计算的与GGFC提供的2011年中国陆区2.5°×2.5°各格网点大气负荷位移时间序列差值的均方根; (b) 三个具有代表性的格网点的时间序列结果(包括本文计算结果、GGFC提供的产品及两者的差值)及所有格网点的RMS统计结果.Fig.2 Comparison between our results and GGFC′s products(a) The RMS of the difference between the 2.5° × 2.5° atmospheric load displacement time series over the year 2011 at each grid point in Chinese continent calculated by the PREM model in this paper and the GGFC; (b) shows the time series of three representative ones and the RMS statistical results of all grid points. The time series results of the grid points including our results, the products provided by GGFC and the difference.

3.2 大气负荷位移特征分析

为了对中国陆区大气负荷垂向位移的空间分布有一个直观的了解，本文首先用PREM模型计算了1°×1°中国陆区2011—2017年大气负荷垂向位移的时间序列，并根据公式(7)采用最小二乘估计方法拟合得到了周年振幅.计算结果如图3所示.从图中可以发现大气负荷效应最大的是我国的华中北部地区、华东西北部地区及华北东南部地区，周年振幅可达5.0 mm以上.其次是华北中部地区、华中南部地区、西北北部地区及华东中部地区，周年振幅超过4.0 mm.影响最小的是青藏高原地区，周年振幅小于2.0 mm.该结果与相关研究(王敏等, 2005; Yue et al., 2020)的结论基本一致.大气负荷效应呈现复杂多样的分布情况主要与中国陆区地形变化、跨越多个纬度带以及季风气候有关.一方面，气压的季节性变化受气压带和季风环流的影响.另一方面，受纬度、海拔等因素的影响，气压变化会随着纬度的升高而增大(姜卫平等，2016)，而随着海拔升高、大气密度减小、气压变化减小，负荷形变也随之减小.

图3 中国陆区大气负荷垂向位移周年振幅图中标注了4个典型测站(SDJX、YNLJ、FJXP、XJJJ).Fig.3 Annual amplitude of PREM-derived ATML vertical displacement in China 4 typical stations (SDJX，YNLJ，FJXP，XJJJ) are marked in the figure.

为了分析大气负荷效应对中国陆区不同区域GPS测站的影响，本文用PREM地球模型计算了中国陆区243个GPS测站的大气负荷垂向位移，通过WRMS分析法来评估大气负荷效应修正效果，并用Pearson相关系数衡量GPS垂向残差与大气负荷垂向形变之间的一致性.两个数据集的相关性及WRMS相对修正量如图4所示.统计结果表明，对于中国陆区的243个GPS测站，大气负荷垂向位移与GPS垂向残差的平均相关系数为0.44，其中强相关(0.6<|CC|<0.8)的测站有58个，主要分布在华北地区及西北北部.中等程度相关(0.4<|CC|<0.6)的测站有74个，主要分布在华中、华北和西北地区.弱相关(0.2<|CC|<0.4)的测站有91个，主要分布在西南地区、西北南部以及华东地区.极弱相关(|CC|<0.2)的测站有20个，主要分布在新疆西南部、西藏西部以及华南地区.在加入大气负荷效应改正后，WRMS的平均相对改正量为10.59%.其中，222个测站在改正后GPS残差的RMS得到不同程度的降低.效果最好的是NMEL站(二连浩特)，其WRMS相对改正量为39.59%.另外，也存在部分测站在大气负荷修正后WRMS反而变大.当数据一致性越高时，大气负荷效应在GPS残差中的比例越高，WRMS的修正效果也越好.

图4 大气负荷垂向位移与GPS垂向残差的关系(a) 大气负荷垂向形变与GPS垂向残差的相关性系数； (b) 大气负荷效应改正后各测站WRMS的修正百分比.Fig.4 The relationship between ATML displacement and GPS vertical residual(a) Correlation coefficient of ATML displacement and GPS vertical residual； (b) Corrected percentage of WRMS of each station after atmospheric loading effect correction.

本文在图3中标注了4个典型测站SDJX、YNLJ、XJJJ和FJXP，下面通过位移时间序列进一步分析GPS残差中大气负荷效应的影响，结果如图5所示.其中，嘉祥站(SDJX)位于华北地区，该测站的大气负荷垂向位移与GPS垂向残差具有很好的一致性，说明在扣除固体潮和海潮等影响后该区域的非构造形变以大气负荷为主.另外，从图中还可以看出，嘉祥站的GPS残差时间序列在2014年1月和2015年1月前后存在明显的跳动，而大气负荷位移时间序列的变化基本一致.这说明在其他信号影响小或改正较好的情况下，GPS观测对气压变化响应十分灵敏，同时大气负荷效应也可以用来解释GPS监测到的地表垂向形变.丽江站(YNLJ)位于西南地区，受热带海洋气团和极地大陆气团的交替控制，属于亚热带季风气候区，兼具海洋性气候和大陆性气候特征，干湿季节分明，其GPS残差时间序列具有明显的季节性且变化幅度远大于大气负荷形变，该地区非构造形变的成因除了气压变化以外，受地表水负荷、土壤水迁移等其他环境因素的影响较大，因此GPS残差与大气负荷形变的一致性较低.霞浦站(FJXP)位于东南沿海地区，受大气负荷的影响较小，并且由于其他环境负荷及海潮模型不准确的影响，其GPS残差时间序列除了表现出明显的季节性以外还存在较强的高频噪声，因此大气负荷效应修正对WRMS的改正较小.芨芨站(XJJJ)位于新疆东部，与前者相比，该区域年降水量少且远离海洋，受陆地水负荷及海潮负荷的影响小.从时间序列来看，两者在峰谷值所处时刻基本一致，大气负荷效应可以解释GPS观测结果的主体部分.但是该站的GPS残差数据中存在大量噪声，因此与嘉祥站相比，芨芨站的数据一致性及WRMS改正量均不高，分别为0.57和7.83%.

图5 部分典型测站的位移时间序列，其中红线表示GPS垂向残差，蓝线表示由PREM模型计算的大气负荷垂向形变Fig.5 Time series of some stations with PREM-derived ATML displacement (blue curve) together with its GPS residuals (red curve) for the Up component

3.3 不同地球模型的建模差异

为了研究地球模型对大气负荷位移建模的影响，本文分别用三种地球模型(AK135、STW105和CRUST1.0)计算了2011—2017年中国陆区1°×1°的大气负荷垂向位移，接着与PREM模型的计算结果进行对比.图6分别为2011年不同模型的负荷位移的最大差值.其中，PREM与AK135模型的负荷位移建模结果存在区域性差异，西北、东北、华中和华东地区的建模差异最大，超过0.2 mm，青藏高原地区和华南地区的建模差异较小，仅0.1 mm左右.PREM和STW105的负荷位移建模的差异分布与前者基本一致，但是后者的计算结果更接近，最大建模差异小于0.05 mm.PREM与CRUST1.0模型的建模结果也存在区域性差异，但与前两者的分布不同，其中青藏高原地区的差异最大，超过0.4 mm，西北地区次之，在0.2 mm左右，东部地区最小，低于0.05 mm.

图6 不同模型计算的大气负荷位移在2011年的最大差值的绝对值(a、b、c) 分别为STW105、AK135和CRUST1.0与PREM地球模型的最大负荷位移差.Fig.6 The absolute value of the maximum differences between the ATML displacement using different earth models over 2011(a, b，c) are the largest load displacement difference between STW105, AK135, CRUST1.0 and the PREM earth model, respectively.

使用一维地球模型进行负荷位移建模时无法考虑浅层地壳结构的区域差异，单一的结构使不同一维地球模型的建模差异具有相同的分布特征即主要集中在负荷较大的区域.青藏高原地区的大气负荷较小，因此一维地球模型在该区域的负荷建模结果基本一致.CRUST1.0模型考虑了地壳结构的区域差异，表现出与前者不同的差异分布，尤其是青藏高原地区.为了比较不同地球模型在该区域的负荷位移建模精度，本文用AK135、STW105和CRUST1.0模型计算了青藏高原地区44个GPS测站2011—2017年的大气负荷垂向位移及WRMS的相对改正量，并与PREM模型的改正效果进行对比，如图7所示.图中的值为

图7 WRMS相对改正量差值(a、b、c) 分别为STW105、AK135和CRUST1.0与PREM模型对青藏高原地区不同测站WRMS相对改正量之差.Fig.7 Relative correction difference of WRMS(a, b, c) are the differences of relative WRMS correction amount of different stations in the Qinghai-Tibet Plateau region between STW105, AK135, CRUST1.0 and the PREM earth model, respectively.

(8)

其中，每个地球模型的WRMS(diff percentage)通过公式(4)计算得到.结果为正表示相比于PREM模型，使用其他地球模型计算的大气负荷位移可以更好地降低GPS残差的WRMS，结果为负则表示改正效果不如PREM模型.

由图7可知，STW105和PREM模型对WRMS的改正效果相当.与PREM相比，改正最好的是炉霍站(SCLH)，WRMS减小了0.07%，最差的是同仁站(QHTR)，WRMS增大了0.01%，平均降低了0.02%.AK135模型的改正效果比PREM模型差，其中改正最好的是祁连站(QHQL)，WRMS减小了0.20%，最差的是拉萨站(LHAZ)，WRMS增大了0.40%，平均升高了0.16%.CRUST1.0模型的改正效果明显优于PREM，改正最好的是拉萨站(LHAZ)，WRMS减小了0.93%，最差的是祁连站(QHQL)，WRMS增大了0.30%，平均降低了0.30%.由前面的分析可知，青藏高原地区的垂向地表形变受大气负荷效应的影响较小，受其他环境负荷的影响更大，因此不同模型对残差WRMS的改正较小.但从该结果可以看出，考虑地壳浅层局部差异的CRUST1.0模型与一维地球模型的负荷位移建模结果存在明显差异.但是对于不到1%的改正量来讲，其效果虽然有，但相比于其他的误差和影响来说，三维模型并不占优.

3.4 负荷格林函数差异

负荷位移建模中不同地球模型的差异主要体现为负荷格林函数的差异，而由负荷格林函数差异引起的建模差异可以简化为下式：

(9)

其中,Δu为位移差，S为负荷作用区域的面积，ΔG为不同地球模型的负荷格林函数差值.

由该式可知，不同地球模型的建模差异主要与负荷作用面积、负荷格林函数差值以及气压变化量有关.图8为不同地球模型的标准化的垂向负荷格林函数，其中CRUST1.0模型表示的是中国陆区2.5°×2.5°的结果.可以发现，不同地球模型的负荷格林函数在远场基本一致，差异集中在近场(θ<1°).因此，在比较不同地球模型的大气负荷位移建模结果时，主要考虑近场的气压变化量及负荷格林函数差异.当近场气压变化量固定时，建模差异与负荷格林函数之差成正比.当负荷格林函数之差固定时，建模差异与近场气压变化量成正比.对于STW105模型和AK135模型，由于东部地区和西北地区的气压变化量大，青藏高原地区的气压变化量小，因此二者与PREM模型的建模差异均在东部地区和西北地区较大而在青藏高原地区较小.由图8可以看出，当0.01°<θ<1°时，AK135模型与PREM模型的负荷格林函数是完全分离的，而STW105模型与PREM模型的负荷格林函数存在交叉的情况，因此在相同的差异分布下，AK135模型与PREM模型的建模差异大于STW105模型与PREM模型的建模差异.然而在使用CRUST1.0模型进行负荷位移建模时，考虑了浅层地壳结构的局部差异，其负荷格林函数与PREM模型的负荷格林函数在东部地区的差异较小，在西北地区存在一定差异，在青藏高原地区的差异较为显著.因此尽管东部地区和西北地区的气压变化量大，建模差异同样较小.青藏高原地区的气压变化量小，一维地球模型之间的差异尚不足以发现明显的负荷位移变化，而较大的负荷格林函数差异使CRUST1.0模型可以得到与PREM相差超过0.4 mm的建模结果.

图8 垂向负荷格林函数红、绿、蓝三种虚线分别表示PREM、AK135和STW105模型的负荷格林函数值，黄色虚线表示CRUST1.0模型在中国陆区2.5°×2.5°各格网点的结果.1012(aθ) m per kg load表示每千克负荷引起的形变乘以1012(a θ)，单位是米，其中标准化量a为地球半径6.371×106 m，θ为负荷点与形变点间的角距，下同.Fig.8 Load Green′s functions for vertical displacement The red, green and blue dotted lines represent the load Green′s functions values of the AK135, PREM and STW105 models, respectively,the yellow dotted line represents the results of the CRUST1.0 model at 2.5°×2.5° grid points in Chinese continent. The normalizations a=earth′s radius, 6.371×106 m and θ=distance from load in radians, the same below.

对于一维地球模型，由于考虑的是全球或区域的综合影响，得到的是全球或区域平均的地壳结构参数，在进行负荷位移建模时无法考虑局部特殊的地壳结构.相关研究表明，负荷响应强烈依赖地壳和上地幔的弹性结构(毛伟建,1984)且局部负荷格林函数与岩层类型及厚度相关(Dill et al., 2015).为了进一步分析负荷格林函数差异与地壳结构之间的关系，本文通过图9展示了中国陆区上地壳层、中地壳层和下地壳层的厚度分布以及PREM模型与CRUST1.0模型的负荷格林函数在不同角距下的差值分布.由图9a—c可以看出中国陆区不同层级的地壳厚度呈阶梯状分布.其中，对于上地壳层，青藏高原地区的地壳厚度最大，西北地区和中部地区次之，东部地区最小，这与中国地形的三级阶梯分布比较相似；对于中地壳层，中部地区和西部地区的地壳厚度较大而东部地区的地壳厚度较小；对于下地壳层，青藏高原地区的地壳厚度最大，东部地区次之，西北地区和中部地区最小.图9d—f展示了当角距θ为0.16°、0.25°及0.50°时PREM模型与CRUST1.0模型的负荷格林函数差值分布.其中，青藏高原地区的差值最大，西北地区和中部地区次之，东部地区的差值最小.通过对比发现，负荷格林函数在近场的差值分布与上地壳厚度分布十分接近，因此负荷格林函数差异可能受上地壳厚度的影响较大.中国陆区地壳厚度分布不均，东部地区厚度小且变化平缓，西部地区厚度大但变化复杂，然而全球平均的一维地球模型无法反映局部特殊的地壳结构，并且由地球模型差异引起的负荷格林函数差异会产生较大的负荷位移建模误差，因此在对地壳结构复杂的区域进行负荷位移建模时建议使用考虑区域差异的三维地球模型如CRUST1.0.

图9 垂向负荷格林函数差值及不同地壳层厚度的分布(a、b、c) 分别为中国陆区上地壳层、中地壳层和下地壳层的地壳厚度分布； (d、e、f) 分别为当θ=0.16°、θ=0.25°和θ=0.50°时PREM模型与CRUST1.0模型在中国陆区的负荷格林函数差值分布.Fig.9 The distribution of vertical load Green′s function difference and different crustal thickness(a,b,c) The crustal thickness distribution of the upper, middle and lower crustal layers of the Chinese continent, respectively; (d,e,f) The distribution of load Green′s function difference between PREM and CRUST1.0 in Chinese continent when angular distance θ=0.16°, θ=0.25° and θ=0.50°.

4 结论

为了研究不同地球模型对大气负荷位移建模的影响，本文用全球地表气压再分析数据NCEP R-1和不同地球模型的负荷格林函数计算了中国区域2011—2017年大气负荷垂向位移；比较了中国陆区243个GPS测站的垂向残差与大气负荷垂向位移的一致性及在进行大气负荷效应改正后WRMS的变化；比较了不同地球模型的负荷格林函数差异及其对负荷位移计算的影响；分析了PREM与一维地球模型AK135、STW105和三维地球模型CRUST1.0的负荷位移建模差异及其对青藏高原地区测站GPS残差RMS的改正效果；讨论了地壳厚度与负荷格林函数的关系.主要研究结论如下：

(1) 大气负荷引起的垂向地表形变在我国呈现区域性分布，其中影响最大的是我国的华中地区、华东北部地区以及西北北部地区，周年振幅达到5 mm，其次是华北地区、西北西部地区及华南中部地区，周年振幅在4 mm左右，影响最小的是青藏高原地区，周年振幅约为2 mm.通过向GPS垂向残差中加入大气负荷效应改正，平均可使残差WRMS降低10.59%.

(2) 建模差异除了负荷格林函数外还与负荷量的大小有关.PREM与一维地球模型的大气负荷位移建模差异较小，但存在区域分化，即东部地区和西北地区差异大，青藏高原地区差异小.另外，不同同地球模型的位移格林函数在远场(θ>1°)基本一致，相比于其他一维模型，PREM模型在近场区域的位移格林函数更接近CRUST1.0模型的结果.因此，对于一维模型，建议采用PREM.

(3) 在青藏高原地区，两类模型的负荷位移建模结果明显不同.一维地球模型的大气负荷位移建模结果对GPS残差的改正效果相当.CRUST1.0模型对该区域GPS残差的改正效果较好，但与其他误差和影响相比并不占优.不同地球模型的负荷格林函数在近场的差异最大，而近场的负荷格林函数与上地壳厚度的关系较大.中国陆区的地壳厚度分布不均，仅使用一维地球模型并不足以体现其差异.因此在进行负荷位移建模时，地壳结构的局部差异的影响不容忽视.

致谢感谢国家地震科学数据中心提供的GPS坐标时间序列数据，感谢地震联合研究机构提供的地球模型参数，感谢Martens博士提供的LoadDef软件源码，感谢van Dam教授提供的全球海陆掩膜数据，感谢美国国家环境预报中心和美国国家大气研究中心提供的全球地表气压再分析数据，感谢全球地球物理流体中心提供的大气负荷位移全球格网数据以及本文实验中画图所用的GMT软件及其相关开发者们在此一并感谢.