李長(zhǎng)海， 鐘萍， 姜中山， 湯苗， 楊興海， 由曉文

西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院， 成都 611756

0 引言

環(huán)境負(fù)載如大氣、海洋、陸地水等的變化會(huì)使固體地球表面產(chǎn)生明顯的位移變化，并同步引起重力場(chǎng)變化.環(huán)境負(fù)載引起的地表形變的準(zhǔn)確計(jì)算有利于相關(guān)地球物理信號(hào)的研究(Van Dam et al., 1994; 王敏等，2005; 袁林果等，2008).相比于其他環(huán)境負(fù)載，地表氣壓的觀測(cè)數(shù)據(jù)、數(shù)值模式和同化方案都較為完善(趙天保等, 2010)，并且當(dāng)前大氣負(fù)荷效應(yīng)的研究也日趨成熟.相關(guān)研究表明，大氣負(fù)荷引起的垂向地表形變具有明顯的季節(jié)性特征且在中國(guó)地區(qū)的變化幅度最大超過(guò)20 mm(張?jiān)娪竦? 2006; 周江存和孫和平, 2007; Yue et al., 2020)，對(duì)GPS站點(diǎn)位移進(jìn)行大氣負(fù)荷修正可以明顯減小數(shù)據(jù)殘差(Van Dam et al., 1994; Tregoning and Watson, 2009; 姜衛(wèi)平等, 2016).因此，對(duì)于地殼構(gòu)造運(yùn)動(dòng)的研究及大地參考框架的建立，精確改正大氣負(fù)荷效應(yīng)是十分必要的.

負(fù)荷位移建模方法主要有格林函數(shù)法和球諧函數(shù)法，兩者在數(shù)學(xué)上是等價(jià)的，其精度在誤差水平上也是一致的(沈迎春等, 2017).格林函數(shù)法是應(yīng)用卷積計(jì)算全球或區(qū)域格點(diǎn)負(fù)荷對(duì)測(cè)站的影響，在負(fù)荷位移建模中被廣泛采用(Farrell, 1972; van Dam et al., 1994; Dill et al., 2015; Li et al., 2020).早期研究采用一維模型來(lái)表征地球內(nèi)部結(jié)構(gòu)及密度和速度等特征(Dziewonski and Anderson, 1981)，當(dāng)前進(jìn)行負(fù)荷位移建模也普遍采用一維地球模型.但是隨著理論研究的深入和觀測(cè)精度的不斷提高，建立更加完整和精確的負(fù)荷模型愈發(fā)重要.馮銳(1985)的研究表明中國(guó)陸區(qū)地殼厚度分布不均，東部變化平緩，西部變化復(fù)雜.基于對(duì)各種地球物理數(shù)據(jù)的統(tǒng)計(jì)分析，Laske等(2013)提出了三維地球模型CRUST1.0，該模型最深可達(dá)上地幔，考慮了地殼的區(qū)域差異.Wang等(2012)計(jì)算了不同地球模型(PREM、iasp91、ak135+crust2.0)的負(fù)荷勒夫數(shù)，這些勒夫數(shù)在超過(guò)200階后差異較大.他們認(rèn)為這可能是由于模型參數(shù)在220 km深度處不連續(xù)和不同模型在地殼中彈性參數(shù)差異的影響所致，并建議當(dāng)使用格林函數(shù)角距Θ<1°時(shí)，用CRUST2.0的地殼結(jié)構(gòu)代替PREM或ak135的上層.賈路路等(2014)的研究表明，全球平均的地球模型不能很好地反映中國(guó)陸區(qū)特別是青藏高原地區(qū)的地殼結(jié)構(gòu).Dill等(2015)采用局部負(fù)荷格林函數(shù)來(lái)計(jì)算負(fù)荷形變，結(jié)果表明，區(qū)域地殼結(jié)構(gòu)對(duì)負(fù)荷位移建模的影響不容忽視.

中國(guó)陸區(qū)的地殼結(jié)構(gòu)復(fù)雜，是研究地球模型差異對(duì)負(fù)荷建模影響的天然試驗(yàn)場(chǎng).為了研究地球模型差異對(duì)利用地表氣壓數(shù)據(jù)估算中國(guó)陸區(qū)地表垂向形變的影響，本文基于不同地球模型計(jì)算了對(duì)應(yīng)的負(fù)荷格林函數(shù)，并與全球地表氣壓再分析數(shù)據(jù)估算中國(guó)陸區(qū)大氣負(fù)荷垂向形變，對(duì)比分析不同地球模型的負(fù)荷位移建模差異.結(jié)合研究區(qū)域內(nèi)GPS參考站的去線性時(shí)間序列數(shù)據(jù)，根據(jù)大氣負(fù)荷修正后GPS殘差WRMS的改正情況來(lái)評(píng)估不同地球模型的建模差異.

1 數(shù)據(jù)

1.1 地球模型

本文討論了四種地球模型，包括一維地球模型PREM、STW105、AK135和三維地球模型CRUST1.0.以上地球模型的模型參數(shù)(不同地表深度下的縱波波速VP、橫波波速VS和巖石密度ρ)由地震聯(lián)合研究機(jī)構(gòu)(Incorporated Research Institutions for Seismology，IRIS)提供.其中CRUST1.0模型由冰層、水層、沉積層和地殼層組成.在計(jì)算時(shí)，本文將CRUST1.0模型地幔層及以下的結(jié)構(gòu)用PREM模型中的對(duì)應(yīng)部分補(bǔ)充.不同地球模型的模型參數(shù)如圖1所示，當(dāng)深度超過(guò)1000 km后，不同地球模型的模型參數(shù)基本一致，為了便于展示，我們僅繪制了1000 km以上的部分.

圖1 不同地球模型的模型參數(shù)縱坐標(biāo)表示地球內(nèi)部一點(diǎn)至地表的深度.橫坐標(biāo)表示該深度下的模型參數(shù)值.紅、綠、藍(lán)三種實(shí)線分別表示AK135、PREM和STW105模型的模型參數(shù)(包括縱波波速VP，橫波波速VS和密度ρ)，黃色虛線表示CRUST1.0模型在中國(guó)陸區(qū)3°×3°各格網(wǎng)點(diǎn)的模型參數(shù).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殘差時(shí)間序列

為了正確分析大氣負(fù)荷效應(yīng)對(duì)中國(guó)陸區(qū)GPS測(cè)站的影響，本文采用國(guó)家地震科學(xué)數(shù)據(jù)中心提供的單天解去線性項(xiàng)的精密坐標(biāo)時(shí)間序列(http:∥www.eqdsc.com/data/pgv-sjxl.htm)，該數(shù)據(jù)使用麻省理工學(xué)院提供的GAMIT/GLOBK軟件進(jìn)行解算，采用IGS提供的精密星歷信息，采用GMF模型改正先驗(yàn)對(duì)流層延遲，使用IERS03模型改正固體潮，使用otl_FES2004全球海潮模型改正海潮負(fù)荷效應(yīng).數(shù)據(jù)時(shí)間跨度為2011—2017年.由于本文獲取的GPS數(shù)據(jù)已經(jīng)扣除了部分地球物理信號(hào)(固體潮、海潮負(fù)荷)，因此文中將其稱為GPS殘差.為了獲得更可靠的統(tǒng)計(jì)結(jié)果，本文剔除了GPS殘差時(shí)間序列中超過(guò)三倍中誤差的歷元數(shù)據(jù)，篩選出了數(shù)據(jù)完整率大于70%的243個(gè)測(cè)站.

1.3 地表氣壓數(shù)據(jù)

本文采用美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心(National Centers for Environmental Prediction, NCEP)和美國(guó)國(guó)家大氣研究中心(National Center for Atmospheric Research, NCAR)提供的全球地表氣壓再分析數(shù)據(jù)NCEP R-1(https:∥psl.noaa.gov/data/gridded/data.ncep.reanalysis.surface.html)計(jì)算大氣負(fù)荷引起的地表形變，數(shù)據(jù)時(shí)間跨度為2011—2017年.該數(shù)據(jù)的空間分辨率為2.5°×2.5°，時(shí)間分辨率為6 h，有3～5天的數(shù)據(jù)延遲.由于原GPS數(shù)據(jù)未扣除大氣潮汐效應(yīng)，為了保持一致，本文在計(jì)算氣壓變化量時(shí)未扣除地表氣壓數(shù)據(jù)中的潮汐項(xiàng)(如S1、S2).

2 計(jì)算方法

2.1 負(fù)荷格林函數(shù)計(jì)算

本文通過(guò)負(fù)荷格林函數(shù)與全球格網(wǎng)點(diǎn)負(fù)荷卷積計(jì)算得到大氣負(fù)荷引起的地表形變.負(fù)荷格林函數(shù)表示單位點(diǎn)負(fù)荷變化引起的地表形變，其中垂向位移格林函數(shù)可以通過(guò)下式計(jì)算：

(1)

其中，Gr表示垂向位移格林函數(shù)，Θ表示負(fù)荷源與形變點(diǎn)之間的角距離，RE為地球半徑，me為地球質(zhì)量，n為計(jì)算階數(shù)，hn為負(fù)荷勒夫數(shù)，Pn為勒讓德多項(xiàng)式.以地表至6371 km深度的地球模型參數(shù)作為輸入數(shù)據(jù)，本文使用Martens等(2019)提供的LoadDef程序中的run_ln計(jì)算了PREM、STW105、AK135和CRUST1.0模型10000階次的負(fù)荷勒夫數(shù)，并將該結(jié)果作為輸入數(shù)據(jù)，使用run_gf計(jì)算了相應(yīng)的負(fù)荷格林函數(shù).對(duì)于CRUST1.0地球模型，本文用CRUST1.0模型每個(gè)2.5°×2.5°網(wǎng)格點(diǎn)的地殼層參數(shù)代替PREM模型的地殼層參數(shù)，計(jì)算了全球陸地2.5°×2.5°的結(jié)果，即每個(gè)格網(wǎng)點(diǎn)的地球模型所對(duì)應(yīng)的負(fù)荷勒夫數(shù)和負(fù)荷格林函數(shù).

2.2 大氣負(fù)荷位移計(jì)算

本文采用格林函數(shù)方法(Farrell, 1972; van Dam et al., 1994; Dill et al., 2015; Li et al., 2020)計(jì)算地表氣壓變化引起的垂向地表形變，公式如下：

(2)

其中，u表示形變點(diǎn)r′處的垂向位移，ΔP表示負(fù)荷源r處的氣壓變化量，g表示重力加速度，Ω表示積分區(qū)域.考慮到海洋的逆氣壓效應(yīng)，需要對(duì)海洋區(qū)域的氣壓進(jìn)行反變氣壓改正，本文采用Spiridonov和Vinogradova(2019)、Li等(2020)和van Dam and Ray(2010)的改正策略，即僅對(duì)陸地區(qū)域進(jìn)行積分計(jì)算，并使用van Dam提供的海陸掩膜數(shù)據(jù).在計(jì)算氣壓變化時(shí)，本文取20年(1991—2010年)的氣壓均值作為參考?xì)鈮?van Dam and Ray, 2010).需要說(shuō)明的是，一維地球模型的負(fù)荷格林函數(shù)僅與負(fù)荷點(diǎn)和形變點(diǎn)之間的角距有關(guān)，而用CRUST1.0模型替換PREM模型的地殼結(jié)構(gòu)后，該地球模型的負(fù)荷格林函數(shù)還與負(fù)荷源所處的位置有關(guān)，本文用CRUST1.0地球模型的負(fù)荷格林函數(shù)進(jìn)行負(fù)荷位移建模時(shí)，參照了Dill等(2015)的計(jì)算方法，公式如下：

(3)

與(2)式不同的是，(3)式中的Gr為負(fù)荷源r所處格網(wǎng)的局部負(fù)荷格林函數(shù).在實(shí)際計(jì)算時(shí)，兩種方法的原理是相同的.

2.3 負(fù)荷修正效果評(píng)估

本文采用WRMS分析法(Jiang et al., 2013)，通過(guò)扣除GPS時(shí)間序列中的大氣負(fù)荷效應(yīng)，統(tǒng)計(jì)扣除前后GPS殘差的加權(quán)均方根值(weighted root mean square, WRMS)的相對(duì)變化量來(lái)評(píng)估大氣負(fù)荷效應(yīng)修正對(duì)測(cè)站坐標(biāo)時(shí)間序列的修正效果.令ω(i)=1/sigu(i)2，WRMS相對(duì)變化量可通過(guò)下式計(jì)算：

WRMS(diff percentage)%=

(4)

WRMS(var)=

(5)

式中var(i)和sigu(i)分別表示i時(shí)刻的變量及其不確定度；Dgps表示GPS垂向殘差；Dload表示大氣負(fù)荷引起的垂向位移；ndat表示測(cè)站的觀測(cè)值數(shù)量.計(jì)算不同地球模型的大氣負(fù)荷效應(yīng)改正所對(duì)應(yīng)的WRMS值時(shí)只需替換對(duì)應(yīng)的Dload分量即可.由式(4)—(5)可知，WRMS相對(duì)變化量大于0表示大氣負(fù)荷校正可以降低GPS時(shí)間序列的RMS.

2.4 周年振幅估計(jì)

本文根據(jù)

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

(6)

采用最小二乘方法對(duì)周年振幅進(jìn)行估計(jì).式中，y表示垂向位移，ti為以年為單位的時(shí)間歷元，a為常數(shù)項(xiàng)，b為線性速度項(xiàng)，c和d為周年運(yùn)動(dòng)的系數(shù)項(xiàng)，vi代表殘差序列，A為擬合得到的周年振幅.

2.5 相關(guān)系數(shù)計(jì)算

本文采用Pearson相關(guān)系數(shù)(Correlation Coefficient,CC)來(lái)衡量不同數(shù)據(jù)之間的一致性，計(jì)算方法如下：

(7)

3 結(jié)果和討論

3.1 計(jì)算結(jié)果驗(yàn)證

為了驗(yàn)證本文大氣負(fù)荷位移計(jì)算結(jié)果的正確性，根據(jù)全球地球物理流體中心(Global geophysical fluids center, GGFC)提供的大氣負(fù)荷位移格網(wǎng)數(shù)據(jù)(van Dam and Ray, 2010)，本文計(jì)算了與之對(duì)應(yīng)的2.5°×2.5°中國(guó)陸區(qū)2011年大氣負(fù)荷垂向位移的時(shí)間序列并將二者進(jìn)行比較.由于GGFC計(jì)算大氣負(fù)荷位移時(shí)用低通濾波扣除了地表氣壓數(shù)據(jù)的高頻項(xiàng)，因此本文比較兩種結(jié)果時(shí)使用的是日均值.由圖2可以看出中國(guó)陸區(qū)大部分地區(qū)本文的計(jì)算結(jié)果和GGFC提供的產(chǎn)品精度相當(dāng).由于計(jì)算時(shí)采取的積分間隔不同，因此東部沿海地區(qū)的負(fù)荷位移建模結(jié)果存在略微差異，但并不影響不同地球模型間的比較.

圖2 本文計(jì)算結(jié)果與GGFC提供的產(chǎn)品對(duì)比(a) 本文用PREM模型計(jì)算的與GGFC提供的2011年中國(guó)陸區(qū)2.5°×2.5°各格網(wǎng)點(diǎn)大氣負(fù)荷位移時(shí)間序列差值的均方根; (b) 三個(gè)具有代表性的格網(wǎng)點(diǎn)的時(shí)間序列結(jié)果(包括本文計(jì)算結(jié)果、GGFC提供的產(chǎn)品及兩者的差值)及所有格網(wǎng)點(diǎn)的RMS統(tǒng)計(jì)結(jié)果.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 大氣負(fù)荷位移特征分析

為了對(duì)中國(guó)陸區(qū)大氣負(fù)荷垂向位移的空間分布有一個(gè)直觀的了解，本文首先用PREM模型計(jì)算了1°×1°中國(guó)陸區(qū)2011—2017年大氣負(fù)荷垂向位移的時(shí)間序列，并根據(jù)公式(7)采用最小二乘估計(jì)方法擬合得到了周年振幅.計(jì)算結(jié)果如圖3所示.從圖中可以發(fā)現(xiàn)大氣負(fù)荷效應(yīng)最大的是我國(guó)的華中北部地區(qū)、華東西北部地區(qū)及華北東南部地區(qū)，周年振幅可達(dá)5.0 mm以上.其次是華北中部地區(qū)、華中南部地區(qū)、西北北部地區(qū)及華東中部地區(qū)，周年振幅超過(guò)4.0 mm.影響最小的是青藏高原地區(qū)，周年振幅小于2.0 mm.該結(jié)果與相關(guān)研究(王敏等, 2005; Yue et al., 2020)的結(jié)論基本一致.大氣負(fù)荷效應(yīng)呈現(xiàn)復(fù)雜多樣的分布情況主要與中國(guó)陸區(qū)地形變化、跨越多個(gè)緯度帶以及季風(fēng)氣候有關(guān).一方面，氣壓的季節(jié)性變化受氣壓帶和季風(fēng)環(huán)流的影響.另一方面，受緯度、海拔等因素的影響，氣壓變化會(huì)隨著緯度的升高而增大(姜衛(wèi)平等，2016)，而隨著海拔升高、大氣密度減小、氣壓變化減小，負(fù)荷形變也隨之減小.

圖3 中國(guó)陸區(qū)大氣負(fù)荷垂向位移周年振幅圖中標(biāo)注了4個(gè)典型測(cè)站(SDJX、YNLJ、FJXP、XJJJ).Fig.3 Annual amplitude of PREM-derived ATML vertical displacement in China 4 typical stations (SDJX，YNLJ，F(xiàn)JXP，XJJJ) are marked in the figure.

為了分析大氣負(fù)荷效應(yīng)對(duì)中國(guó)陸區(qū)不同區(qū)域GPS測(cè)站的影響，本文用PREM地球模型計(jì)算了中國(guó)陸區(qū)243個(gè)GPS測(cè)站的大氣負(fù)荷垂向位移，通過(guò)WRMS分析法來(lái)評(píng)估大氣負(fù)荷效應(yīng)修正效果，并用Pearson相關(guān)系數(shù)衡量GPS垂向殘差與大氣負(fù)荷垂向形變之間的一致性.兩個(gè)數(shù)據(jù)集的相關(guān)性及WRMS相對(duì)修正量如圖4所示.統(tǒng)計(jì)結(jié)果表明，對(duì)于中國(guó)陸區(qū)的243個(gè)GPS測(cè)站，大氣負(fù)荷垂向位移與GPS垂向殘差的平均相關(guān)系數(shù)為0.44，其中強(qiáng)相關(guān)(0.6<|CC|<0.8)的測(cè)站有58個(gè)，主要分布在華北地區(qū)及西北北部.中等程度相關(guān)(0.4<|CC|<0.6)的測(cè)站有74個(gè)，主要分布在華中、華北和西北地區(qū).弱相關(guān)(0.2<|CC|<0.4)的測(cè)站有91個(gè)，主要分布在西南地區(qū)、西北南部以及華東地區(qū).極弱相關(guān)(|CC|<0.2)的測(cè)站有20個(gè)，主要分布在新疆西南部、西藏西部以及華南地區(qū).在加入大氣負(fù)荷效應(yīng)改正后，WRMS的平均相對(duì)改正量為10.59%.其中，222個(gè)測(cè)站在改正后GPS殘差的RMS得到不同程度的降低.效果最好的是NMEL站(二連浩特)，其WRMS相對(duì)改正量為39.59%.另外，也存在部分測(cè)站在大氣負(fù)荷修正后WRMS反而變大.當(dāng)數(shù)據(jù)一致性越高時(shí)，大氣負(fù)荷效應(yīng)在GPS殘差中的比例越高，WRMS的修正效果也越好.

圖4 大氣負(fù)荷垂向位移與GPS垂向殘差的關(guān)系(a) 大氣負(fù)荷垂向形變與GPS垂向殘差的相關(guān)性系數(shù)； (b) 大氣負(fù)荷效應(yīng)改正后各測(cè)站W(wǎng)RMS的修正百分比.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中標(biāo)注了4個(gè)典型測(cè)站SDJX、YNLJ、XJJJ和FJXP，下面通過(guò)位移時(shí)間序列進(jìn)一步分析GPS殘差中大氣負(fù)荷效應(yīng)的影響，結(jié)果如圖5所示.其中，嘉祥站(SDJX)位于華北地區(qū)，該測(cè)站的大氣負(fù)荷垂向位移與GPS垂向殘差具有很好的一致性，說(shuō)明在扣除固體潮和海潮等影響后該區(qū)域的非構(gòu)造形變以大氣負(fù)荷為主.另外，從圖中還可以看出，嘉祥站的GPS殘差時(shí)間序列在2014年1月和2015年1月前后存在明顯的跳動(dòng)，而大氣負(fù)荷位移時(shí)間序列的變化基本一致.這說(shuō)明在其他信號(hào)影響小或改正較好的情況下，GPS觀測(cè)對(duì)氣壓變化響應(yīng)十分靈敏，同時(shí)大氣負(fù)荷效應(yīng)也可以用來(lái)解釋GPS監(jiān)測(cè)到的地表垂向形變.麗江站(YNLJ)位于西南地區(qū)，受熱帶海洋氣團(tuán)和極地大陸氣團(tuán)的交替控制，屬于亞熱帶季風(fēng)氣候區(qū)，兼具海洋性氣候和大陸性氣候特征，干濕季節(jié)分明，其GPS殘差時(shí)間序列具有明顯的季節(jié)性且變化幅度遠(yuǎn)大于大氣負(fù)荷形變，該地區(qū)非構(gòu)造形變的成因除了氣壓變化以外，受地表水負(fù)荷、土壤水遷移等其他環(huán)境因素的影響較大，因此GPS殘差與大氣負(fù)荷形變的一致性較低.霞浦站(FJXP)位于東南沿海地區(qū)，受大氣負(fù)荷的影響較小，并且由于其他環(huán)境負(fù)荷及海潮模型不準(zhǔn)確的影響，其GPS殘差時(shí)間序列除了表現(xiàn)出明顯的季節(jié)性以外還存在較強(qiáng)的高頻噪聲，因此大氣負(fù)荷效應(yīng)修正對(duì)WRMS的改正較小.芨芨站(XJJJ)位于新疆東部，與前者相比，該區(qū)域年降水量少且遠(yuǎn)離海洋，受陸地水負(fù)荷及海潮負(fù)荷的影響小.從時(shí)間序列來(lái)看，兩者在峰谷值所處時(shí)刻基本一致，大氣負(fù)荷效應(yīng)可以解釋GPS觀測(cè)結(jié)果的主體部分.但是該站的GPS殘差數(shù)據(jù)中存在大量噪聲，因此與嘉祥站相比，芨芨站的數(shù)據(jù)一致性及WRMS改正量均不高，分別為0.57和7.83%.

圖5 部分典型測(cè)站的位移時(shí)間序列，其中紅線表示GPS垂向殘差，藍(lán)線表示由PREM模型計(jì)算的大氣負(fù)荷垂向形變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 不同地球模型的建模差異

為了研究地球模型對(duì)大氣負(fù)荷位移建模的影響，本文分別用三種地球模型(AK135、STW105和CRUST1.0)計(jì)算了2011—2017年中國(guó)陸區(qū)1°×1°的大氣負(fù)荷垂向位移，接著與PREM模型的計(jì)算結(jié)果進(jìn)行對(duì)比.圖6分別為2011年不同模型的負(fù)荷位移的最大差值.其中，PREM與AK135模型的負(fù)荷位移建模結(jié)果存在區(qū)域性差異，西北、東北、華中和華東地區(qū)的建模差異最大，超過(guò)0.2 mm，青藏高原地區(qū)和華南地區(qū)的建模差異較小，僅0.1 mm左右.PREM和STW105的負(fù)荷位移建模的差異分布與前者基本一致，但是后者的計(jì)算結(jié)果更接近，最大建模差異小于0.05 mm.PREM與CRUST1.0模型的建模結(jié)果也存在區(qū)域性差異，但與前兩者的分布不同，其中青藏高原地區(qū)的差異最大，超過(guò)0.4 mm，西北地區(qū)次之，在0.2 mm左右，東部地區(qū)最小，低于0.05 mm.

圖6 不同模型計(jì)算的大氣負(fù)荷位移在2011年的最大差值的絕對(duì)值(a、b、c) 分別為STW105、AK135和CRUST1.0與PREM地球模型的最大負(fù)荷位移差.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.

使用一維地球模型進(jìn)行負(fù)荷位移建模時(shí)無(wú)法考慮淺層地殼結(jié)構(gòu)的區(qū)域差異，單一的結(jié)構(gòu)使不同一維地球模型的建模差異具有相同的分布特征即主要集中在負(fù)荷較大的區(qū)域.青藏高原地區(qū)的大氣負(fù)荷較小，因此一維地球模型在該區(qū)域的負(fù)荷建模結(jié)果基本一致.CRUST1.0模型考慮了地殼結(jié)構(gòu)的區(qū)域差異，表現(xiàn)出與前者不同的差異分布，尤其是青藏高原地區(qū).為了比較不同地球模型在該區(qū)域的負(fù)荷位移建模精度，本文用AK135、STW105和CRUST1.0模型計(jì)算了青藏高原地區(qū)44個(gè)GPS測(cè)站2011—2017年的大氣負(fù)荷垂向位移及WRMS的相對(duì)改正量，并與PREM模型的改正效果進(jìn)行對(duì)比，如圖7所示.圖中的值為

圖7 WRMS相對(duì)改正量差值(a、b、c) 分別為STW105、AK135和CRUST1.0與PREM模型對(duì)青藏高原地區(qū)不同測(cè)站W(wǎng)RMS相對(duì)改正量之差.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)

其中，每個(gè)地球模型的WRMS(diff percentage)通過(guò)公式(4)計(jì)算得到.結(jié)果為正表示相比于PREM模型，使用其他地球模型計(jì)算的大氣負(fù)荷位移可以更好地降低GPS殘差的WRMS，結(jié)果為負(fù)則表示改正效果不如PREM模型.

由圖7可知，STW105和PREM模型對(duì)WRMS的改正效果相當(dāng).與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模型的改正效果明顯優(yōu)于PREM，改正最好的是拉薩站(LHAZ)，WRMS減小了0.93%，最差的是祁連站(QHQL)，WRMS增大了0.30%，平均降低了0.30%.由前面的分析可知，青藏高原地區(qū)的垂向地表形變受大氣負(fù)荷效應(yīng)的影響較小，受其他環(huán)境負(fù)荷的影響更大，因此不同模型對(duì)殘差WRMS的改正較小.但從該結(jié)果可以看出，考慮地殼淺層局部差異的CRUST1.0模型與一維地球模型的負(fù)荷位移建模結(jié)果存在明顯差異.但是對(duì)于不到1%的改正量來(lái)講，其效果雖然有，但相比于其他的誤差和影響來(lái)說(shuō)，三維模型并不占優(yōu).

3.4 負(fù)荷格林函數(shù)差異

負(fù)荷位移建模中不同地球模型的差異主要體現(xiàn)為負(fù)荷格林函數(shù)的差異，而由負(fù)荷格林函數(shù)差異引起的建模差異可以簡(jiǎn)化為下式：

(9)

其中,Δu為位移差，S為負(fù)荷作用區(qū)域的面積，ΔG為不同地球模型的負(fù)荷格林函數(shù)差值.

由該式可知，不同地球模型的建模差異主要與負(fù)荷作用面積、負(fù)荷格林函數(shù)差值以及氣壓變化量有關(guān).圖8為不同地球模型的標(biāo)準(zhǔn)化的垂向負(fù)荷格林函數(shù)，其中CRUST1.0模型表示的是中國(guó)陸區(qū)2.5°×2.5°的結(jié)果.可以發(fā)現(xiàn)，不同地球模型的負(fù)荷格林函數(shù)在遠(yuǎn)場(chǎng)基本一致，差異集中在近場(chǎng)(θ<1°).因此，在比較不同地球模型的大氣負(fù)荷位移建模結(jié)果時(shí)，主要考慮近場(chǎng)的氣壓變化量及負(fù)荷格林函數(shù)差異.當(dāng)近場(chǎng)氣壓變化量固定時(shí)，建模差異與負(fù)荷格林函數(shù)之差成正比.當(dāng)負(fù)荷格林函數(shù)之差固定時(shí)，建模差異與近場(chǎng)氣壓變化量成正比.對(duì)于STW105模型和AK135模型，由于東部地區(qū)和西北地區(qū)的氣壓變化量大，青藏高原地區(qū)的氣壓變化量小，因此二者與PREM模型的建模差異均在東部地區(qū)和西北地區(qū)較大而在青藏高原地區(qū)較小.由圖8可以看出，當(dāng)0.01°<θ<1°時(shí)，AK135模型與PREM模型的負(fù)荷格林函數(shù)是完全分離的，而STW105模型與PREM模型的負(fù)荷格林函數(shù)存在交叉的情況，因此在相同的差異分布下，AK135模型與PREM模型的建模差異大于STW105模型與PREM模型的建模差異.然而在使用CRUST1.0模型進(jìn)行負(fù)荷位移建模時(shí)，考慮了淺層地殼結(jié)構(gòu)的局部差異，其負(fù)荷格林函數(shù)與PREM模型的負(fù)荷格林函數(shù)在東部地區(qū)的差異較小，在西北地區(qū)存在一定差異，在青藏高原地區(qū)的差異較為顯著.因此盡管東部地區(qū)和西北地區(qū)的氣壓變化量大，建模差異同樣較小.青藏高原地區(qū)的氣壓變化量小，一維地球模型之間的差異尚不足以發(fā)現(xiàn)明顯的負(fù)荷位移變化，而較大的負(fù)荷格林函數(shù)差異使CRUST1.0模型可以得到與PREM相差超過(guò)0.4 mm的建模結(jié)果.

圖8 垂向負(fù)荷格林函數(shù)紅、綠、藍(lán)三種虛線分別表示PREM、AK135和STW105模型的負(fù)荷格林函數(shù)值，黃色虛線表示CRUST1.0模型在中國(guó)陸區(qū)2.5°×2.5°各格網(wǎng)點(diǎn)的結(jié)果.1012(aθ) m per kg load表示每千克負(fù)荷引起的形變乘以1012(a θ)，單位是米，其中標(biāo)準(zhǔn)化量a為地球半徑6.371×106 m，θ為負(fù)荷點(diǎn)與形變點(diǎn)間的角距，下同.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.

對(duì)于一維地球模型，由于考慮的是全球或區(qū)域的綜合影響，得到的是全球或區(qū)域平均的地殼結(jié)構(gòu)參數(shù)，在進(jìn)行負(fù)荷位移建模時(shí)無(wú)法考慮局部特殊的地殼結(jié)構(gòu).相關(guān)研究表明，負(fù)荷響應(yīng)強(qiáng)烈依賴地殼和上地幔的彈性結(jié)構(gòu)(毛偉建,1984)且局部負(fù)荷格林函數(shù)與巖層類(lèi)型及厚度相關(guān)(Dill et al., 2015).為了進(jìn)一步分析負(fù)荷格林函數(shù)差異與地殼結(jié)構(gòu)之間的關(guān)系，本文通過(guò)圖9展示了中國(guó)陸區(qū)上地殼層、中地殼層和下地殼層的厚度分布以及PREM模型與CRUST1.0模型的負(fù)荷格林函數(shù)在不同角距下的差值分布.由圖9a—c可以看出中國(guó)陸區(qū)不同層級(jí)的地殼厚度呈階梯狀分布.其中，對(duì)于上地殼層，青藏高原地區(qū)的地殼厚度最大，西北地區(qū)和中部地區(qū)次之，東部地區(qū)最小，這與中國(guó)地形的三級(jí)階梯分布比較相似；對(duì)于中地殼層，中部地區(qū)和西部地區(qū)的地殼厚度較大而東部地區(qū)的地殼厚度較??；對(duì)于下地殼層，青藏高原地區(qū)的地殼厚度最大，東部地區(qū)次之，西北地區(qū)和中部地區(qū)最小.圖9d—f展示了當(dāng)角距θ為0.16°、0.25°及0.50°時(shí)PREM模型與CRUST1.0模型的負(fù)荷格林函數(shù)差值分布.其中，青藏高原地區(qū)的差值最大，西北地區(qū)和中部地區(qū)次之，東部地區(qū)的差值最小.通過(guò)對(duì)比發(fā)現(xiàn)，負(fù)荷格林函數(shù)在近場(chǎng)的差值分布與上地殼厚度分布十分接近，因此負(fù)荷格林函數(shù)差異可能受上地殼厚度的影響較大.中國(guó)陸區(qū)地殼厚度分布不均，東部地區(qū)厚度小且變化平緩，西部地區(qū)厚度大但變化復(fù)雜，然而全球平均的一維地球模型無(wú)法反映局部特殊的地殼結(jié)構(gòu)，并且由地球模型差異引起的負(fù)荷格林函數(shù)差異會(huì)產(chǎn)生較大的負(fù)荷位移建模誤差，因此在對(duì)地殼結(jié)構(gòu)復(fù)雜的區(qū)域進(jìn)行負(fù)荷位移建模時(shí)建議使用考慮區(qū)域差異的三維地球模型如CRUST1.0.

圖9 垂向負(fù)荷格林函數(shù)差值及不同地殼層厚度的分布(a、b、c) 分別為中國(guó)陸區(qū)上地殼層、中地殼層和下地殼層的地殼厚度分布； (d、e、f) 分別為當(dāng)θ=0.16°、θ=0.25°和θ=0.50°時(shí)PREM模型與CRUST1.0模型在中國(guó)陸區(qū)的負(fù)荷格林函數(shù)差值分布.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 結(jié)論

為了研究不同地球模型對(duì)大氣負(fù)荷位移建模的影響，本文用全球地表氣壓再分析數(shù)據(jù)NCEP R-1和不同地球模型的負(fù)荷格林函數(shù)計(jì)算了中國(guó)區(qū)域2011—2017年大氣負(fù)荷垂向位移；比較了中國(guó)陸區(qū)243個(gè)GPS測(cè)站的垂向殘差與大氣負(fù)荷垂向位移的一致性及在進(jìn)行大氣負(fù)荷效應(yīng)改正后WRMS的變化；比較了不同地球模型的負(fù)荷格林函數(shù)差異及其對(duì)負(fù)荷位移計(jì)算的影響；分析了PREM與一維地球模型AK135、STW105和三維地球模型CRUST1.0的負(fù)荷位移建模差異及其對(duì)青藏高原地區(qū)測(cè)站GPS殘差RMS的改正效果；討論了地殼厚度與負(fù)荷格林函數(shù)的關(guān)系.主要研究結(jié)論如下：

(1) 大氣負(fù)荷引起的垂向地表形變?cè)谖覈?guó)呈現(xiàn)區(qū)域性分布，其中影響最大的是我國(guó)的華中地區(qū)、華東北部地區(qū)以及西北北部地區(qū)，周年振幅達(dá)到5 mm，其次是華北地區(qū)、西北西部地區(qū)及華南中部地區(qū)，周年振幅在4 mm左右，影響最小的是青藏高原地區(qū)，周年振幅約為2 mm.通過(guò)向GPS垂向殘差中加入大氣負(fù)荷效應(yīng)改正，平均可使殘差WRMS降低10.59%.

(2) 建模差異除了負(fù)荷格林函數(shù)外還與負(fù)荷量的大小有關(guān).PREM與一維地球模型的大氣負(fù)荷位移建模差異較小，但存在區(qū)域分化，即東部地區(qū)和西北地區(qū)差異大，青藏高原地區(qū)差異小.另外，不同同地球模型的位移格林函數(shù)在遠(yuǎn)場(chǎng)(θ>1°)基本一致，相比于其他一維模型，PREM模型在近場(chǎng)區(qū)域的位移格林函數(shù)更接近CRUST1.0模型的結(jié)果.因此，對(duì)于一維模型，建議采用PREM.

(3) 在青藏高原地區(qū)，兩類(lèi)模型的負(fù)荷位移建模結(jié)果明顯不同.一維地球模型的大氣負(fù)荷位移建模結(jié)果對(duì)GPS殘差的改正效果相當(dāng).CRUST1.0模型對(duì)該區(qū)域GPS殘差的改正效果較好，但與其他誤差和影響相比并不占優(yōu).不同地球模型的負(fù)荷格林函數(shù)在近場(chǎng)的差異最大，而近場(chǎng)的負(fù)荷格林函數(shù)與上地殼厚度的關(guān)系較大.中國(guó)陸區(qū)的地殼厚度分布不均，僅使用一維地球模型并不足以體現(xiàn)其差異.因此在進(jìn)行負(fù)荷位移建模時(shí)，地殼結(jié)構(gòu)的局部差異的影響不容忽視.

致謝感謝國(guó)家地震科學(xué)數(shù)據(jù)中心提供的GPS坐標(biāo)時(shí)間序列數(shù)據(jù)，感謝地震聯(lián)合研究機(jī)構(gòu)提供的地球模型參數(shù)，感謝Martens博士提供的LoadDef軟件源碼，感謝van Dam教授提供的全球海陸掩膜數(shù)據(jù)，感謝美國(guó)國(guó)家環(huán)境預(yù)報(bào)中心和美國(guó)國(guó)家大氣研究中心提供的全球地表氣壓再分析數(shù)據(jù)，感謝全球地球物理流體中心提供的大氣負(fù)荷位移全球格網(wǎng)數(shù)據(jù)以及本文實(shí)驗(yàn)中畫(huà)圖所用的GMT軟件及其相關(guān)開(kāi)發(fā)者們?cè)诖艘徊⒏兄x.