本文采用的大气数据是由欧洲中期天气预报中心生成的ECMWF Interim全球三维大气再分析资料,包括2002—2017年共192个月的观测数据。ECMWF Interim提供的37个大气等压层的高程、位势高度、比湿、气压、气温等要素时间分辨率为1个月,水平空间分辨率为0.125°×0.125°。其中,地表及各规定等压层的Tesseroid单元体大气密度可利用各要素根据流体静力平衡方程计算得到。

根据地表负荷函数理论(Farrell,1972),青藏高原南缘拉萨(LHAS)、日喀则(XZRK)、仲巴(XZZB)和噶尔(XZGE)4个绝对重力测站(图1)大气负荷重力效应可用其质量和负荷格林函数的褶积积分计算得到:

L(θ,λ)=_sρ(θ,λ)h(θ,λ)G(ψ)ds(1)

式中:θ与λ分别表示绝对重力测站的纬度和经度; ρ表示积分面元的大气密度; h表示积分面元的高程; ψ表示观测点到积分面元的角距; G(ψ)表示重力负荷格林函数,可表示为:

G(ψ)=G^N(ψ)+G^D(ψ)=(g)/(M)∑^∞_n=0nP_n(cosψ)+(g)/(M)∑^∞_n=0[-(n+1)k_n+2h_n]P_n(cosψ)(2)

式中:g表示地球平均重力; M表示地球质量; P_n(cosψ)表示勒让德函数; k_n和h_n表示负荷勒夫数; G^N(ψ)表示牛顿引力项直接效应; G^D(ψ)表示弹性地球产生的形变以及由此使地球内部质量重新分布而引起的弹性项间接效应。

在式(1)的实际计算过程中,本文将大气数据划分为近区(ψ<1°)、中区(1°<ψ<20°)和远区(ψ>20°)3个部分进行计算。对于近区的Tesseroid单元体大气质量,本文利用泰勒级数展开二阶近似计算其牛顿引力项直接效应(Neumeyer et al,2004; Heck,Seitz,2007)。

图1 青藏高原南缘绝对重力测站分布
Fig.1 Distribution of absolute gravity stations in southern margin of Qinghai-Tibet Plateau

本文采用的GRACE时变重力场模型是由德克萨斯大学空间中心CSR提供的Level 2 RL06版本GSM和GAC数据,包括2002年8月至2017年6月共160个月的观测数据。其中,GRACE GSM数据在解算过程中已扣除了非潮汐大气和海洋质量的影响,因此其反映的主要是陆地水质量重新分布引起的变化。对于GRACE时变重力场模型,本文对其一阶项(Swenson et al,2008)和C₂₀(Cheng,Tapley,2004)进行了替换,采用P5M11去相关滤波(Chen et al,2007)和300 km高斯滤波(Wahr et al,1998)的组合滤波方法以压制高阶项球谐系数误差,绝对重力测站重力异常和垂直位移的反演公式可分别表示为:

Δg(θ,λ)=(GM)/(R2)∑nmaxn=0(n-1)∑nm=0P^-nm(cosθ)

[ΔC_nmcos(mλ)+ΔS_nmsin(mλ)](3)

Δr(θ,λ)=R∑nmaxn=1(hn)/(1+kn)∑nm=0P^-nm(cosθ)

[ΔC_nmcos(mλ)+ΔS_nmsin(mλ)](4)

式中:θ和λ分别表示余纬和经度; R表示地球平均半径; ΔC和ΔS分别表示GRACE月时变重力场模型球谐系数; P^-_nm(cosθ)表示正则化缔合勒让德函数。

利用GAMIT/GLOBK和OQCA软件处理2002—2017年青藏高原南缘地区“陆态网络”GPS连续站观测数据。除了LHAS测站以外,另外3个绝对重力测站与GPS连续站均为并址观测。对GPS观测数据的具体处理步骤为:①联合青藏高原南缘及全球均匀分布的IGS站进行GAMIT基线解算,得到测站坐标单日松弛解。解算过程主要的参数设置为:天线相位中心改正采用绝对天线相位中心模型; 对流层延迟映射函数采用GMF模型; 全球固体潮、海潮模型改正遵循IERS 2004协议,但未进行非潮汐大气和海洋改正。②利用QOCA软件综合多天的单日松弛解平差计算得到ITRF2008框架下测站坐标的时间序列。经过以上处理,对GPS连续站垂直位移时间序列中由于天线更换、强震同震位移等影响造成的阶跃进行了修复,并剔除了大于3倍误差的异常值。

对于绝对重力测站上的重力异常和地壳垂直位移时间序列,本文重点分析其趋势及周年变化特征。假定其时间序列主要包含了线性趋势项β₁、周年项β₂,然后基于最小二乘原理对时间序列进

行拟合分析:

Δg(θ,λ,t)=β₀(θ,λ)+β₁(θ,λ)t+β₂cos[ω₁t+α₁(θ,λ)](5)

式中:β₀,β₁,α₂,β₁为求解的参数; ω₁=2π。

为了分析时间序列的趋势异常特征,采用递归最小二乘方法对线性趋势参数的稳定性做判断,并利用递归残差累积和检验识别出时间序列中结构断点的分布位置。如果时间序列存在结构断点,则对时间序列进行分段线性拟合。

基于2002—2017年的ECMWF Interim数据,采用负荷格林函数方法计算得到了LHAS,XZRK,XZZB和XZGE这4个测站的大气负荷重力效应时间序列(图2),并分别提取其周年振幅及趋势变化特征(表1)。基于ECMWF Interim内插计算得到的4个测站的气压数据,本文还采用大气导纳值(-0.303 7 μGal/hPa)方法(罗少聪,2003)计算得到4个测站的大气负荷重力效应时间序列。从图2可以看出,采用负荷格林函数与大气导纳值方法得到的4个测站的大气负荷重力效应较为一致,

图2 4个测站的大气负荷重力效应时间序列
Fig.2 Time series of atmospheric gravity effect at four stations

差值标准差小于0.4 μGal。从图2和表1可以看出,4个测站的时间序列均反映出青藏高原南缘的大气负荷重力效应存在明显的季节性变化特征,一般在夏季和冬季分别达到最小值和最大值,周年振幅分别为1.0,0.9,1.4和1.2 μGal。4个测站的时间序列均没有表现出明显的长期趋势及趋势异常特征,说明该地区的大气重力效应主要表现为周年变化特征。

基于2002—2017年的GRACE GSM数据,采用卫星重力反演方法计算得到了LHAS,XZRK,XZZB和XZGE这4个测站的陆地水负荷重力效应时间序列(图3),并提取其周年振幅及趋势变化特征(表1)。从图3可以看出,4个测站的时间序列均反映出青藏高原南缘的陆地水负荷重力

表1 青藏高原南缘4个测站非构造重力效应
Tab.1 The non-tectonic gravity effect of four stations in Southern Margin of Qinghai-Tibet Plateau

效应存在较为明显的季节性变化特征,LHAS,XZRK,XZZB测站一般在春季和秋季分别达到最小值和最大值,而XZGE测站一般在夏季和冬季分别达到最小值和最大值,4个测站周年振幅分别为2.0,2.6,2.0和1.6 μGal。另外,LHAS和XZRK测站的时间序列表现出了较为明显的长期趋势特征,速率分别为-0.6和-0.4 μGal/a,而XZZB和XZGE测站时间序列的长期趋势相对较小,速率分别为-0.1和-0.2 μGal/a。4个测站的时间序列均存在明显的结构断点,不同时间跨度的GRACE陆地水负荷重力效应趋势变化呈现出明显的差异性特征,说明该地区陆地水负荷重力效应存在较为显著的年际变化特征,这可能与土壤水、积雪、冰川和地下水等多种陆地水文过程有关。

基于2002—2017年的GRACE GSM&GAC和GPS数据,计算得到了4个测站的地壳垂直位移时间序列(图4)。从图4中可以看出,4个测站的垂直位移时间序列均表现出明显的季节性变化特征,二者给出的周年相位较为一致,但是由GRACE数据得出的周年振幅量级明显偏小,分别为6.0,6.5,5.2和1.9 mm,而由GPS数据得出的周年振幅分别为6.5,9.2,7.0和3.6 mm,采用布格模型(-0.2 μGal/mm)换算成重力效应

图3 4个测站的陆地水负荷重力效应时间序列
Fig.3 Time series of hydrological gravity effect at four stations

图4 4个测站的地壳垂直位移时间序列
Fig.4 Time series of crustal vertical displacement at four stations

相当于1.3,1.8,1.4和0.7 μGal(表1)。由GRACE数据得到4个测站的垂直位移时间序列长期趋势分别为0.4,0.3,0.2和0.2 mm/a,与之相比GPS得到的时间序列具有更加明显的长期变化趋势,其速率分别为1.3,0.9,1.4和0.5 mm/a,引起的重力效应相当于-0.3,-0.2,-0.3和-0.1 μGal/a(表1)。

对于青藏高原南缘的绝对重力测站而言,GRACE GSM&GAC反演得到的地壳垂直位移时间序列主要反映了大气和陆地水等大尺度地表流体的负荷效应,与之相比GPS观测得到的结果则更容易受到测站周边局部区域流体质量负荷效应的影响,这可能是GRACE垂直位移时间序列量级明显偏小的主要原因。以上比较分析结果也从侧面反映出青藏高原南缘的绝对重力测站局部地区存在着显著的地表流体质量变化,而目前的大气再分析资料ECMWF Interim以及GRACE GSM数据还不能有效反映其非构造重力效应。