申龍杰 章傳銀 徐鵬飛 楊 康

1 蘭州交通大學(xué)測繪與地理信息學(xué)院，蘭州市安寧西路88號，730070 2 中國測繪科學(xué)研究院，北京市蓮花池西路28號，100830 3 甘肅省地理國情監(jiān)測工程實驗室，蘭州市安寧西路88號，730070 4 地理國情監(jiān)測技術(shù)應(yīng)用國家地方聯(lián)合工程研究中心，蘭州市安寧西路88號，730070

地球表層系統(tǒng)中的環(huán)境變化會引起固體地球發(fā)生形變，表現(xiàn)在地面垂直形變、地面重力、地傾斜與大地水準(zhǔn)面變化[1-2]等方面。大氣負(fù)荷是地表環(huán)境負(fù)荷的主要影響因素之一，是研究GNSS時間序列周期信號特征及地表負(fù)荷形變場不可忽略的因素。近年來，許多學(xué)者在研究固體地球形變和GNSS站時間序列分析中均已考慮大氣壓變化的負(fù)荷影響。Martens等[3]分析GPS基站位置的大氣壓力負(fù)荷表明，季節(jié)性變形存在差異；Yue等[4]結(jié)合GNSS基準(zhǔn)站分析不同地區(qū)大氣變化與地殼形變的相關(guān)性；安旭偉等[5]對祁連山南麓區(qū)域的大氣負(fù)荷動力學(xué)效應(yīng)進(jìn)行研究；王海濤等[6]結(jié)合地表溫度變化與大氣壓變化，綜合分析溫度變化與大氣壓變化對基準(zhǔn)站垂向形變的影響規(guī)律；劉宇等[7]計算浙江地區(qū)的負(fù)荷形變場和重力場變化發(fā)現(xiàn)，該區(qū)域的大氣負(fù)荷影響達(dá)cm級，該研究雖然結(jié)合了區(qū)域內(nèi)部氣象站數(shù)據(jù)，但氣象站分布離散且稀疏，內(nèi)插格網(wǎng)化時會引入噪聲。上述研究表明，在高精度高分辨率區(qū)域進(jìn)行大氣負(fù)荷形變場解算時，結(jié)合區(qū)域高分辨率大氣壓格網(wǎng)數(shù)據(jù)的相關(guān)研究較少，因此本文利用中國氣象網(wǎng)滇西地區(qū)高分辨率的大氣壓數(shù)據(jù)進(jìn)一步精化該區(qū)域的大氣負(fù)荷影響，提高大地測量精度。

本文以滇西地區(qū)為例，引入局部重力場逼近的移去-恢復(fù)思想，基于負(fù)荷球諧系數(shù)和區(qū)域負(fù)荷格林函數(shù)，結(jié)合滇西地區(qū)的全球大氣壓變化數(shù)據(jù)和區(qū)域高分辨率大氣壓變化數(shù)據(jù)，分析大氣壓變化對滇西地區(qū)的負(fù)荷影響。

1 理論與方法

大氣壓變化引起的負(fù)荷形變場和時變重力場變化由負(fù)荷位和負(fù)荷附加位兩個部分組成[8]，常用的計算方法有球諧系數(shù)法和格林函數(shù)積分法等。

1.1 全球負(fù)荷形變場的球諧分析

地球表層大氣壓變化為非潮汐運(yùn)動，可用地面等效水高表示，地面上任意一點(R,θ,λ)處的等效水高h(yuǎn)w可表示為規(guī)格化負(fù)荷球諧級數(shù)[8]：

hw(R,θ,λ)=

(1)

由負(fù)荷形變理論[1]可知，地面及地球外部的高程異常負(fù)荷影響可表示為：

(2)

地面重力負(fù)荷影響計算公式為：

(3)

大地高負(fù)荷影響計算公式為：

(4)

1.2 區(qū)域負(fù)荷影響的格林函數(shù)積分公式

地面測站位置的變化也可用負(fù)荷數(shù)(負(fù)荷格林函數(shù))進(jìn)行表示[9]：

(5)

式中，ΔΘ(φ,λ)為地面測站位置的變化量，ψ為地面測站點和流動點(φ′,λ′)之間的球面角距，G(ψ)為負(fù)荷格林函數(shù)，ρw為水的密度，S為地球表面積。

格林函數(shù)積分法計算的負(fù)荷位[10]為：

(6)

格林函數(shù)積分法計算的負(fù)荷附加位為：

(7)

式(6)、(7)中，G為地球平均半徑，M為地球總質(zhì)量，k′n為n階位負(fù)荷數(shù)。

地面重力的負(fù)荷格林函數(shù)總影響[10]可表示為：

(8)

大地高的負(fù)荷格林函數(shù)總影響計算公式為：

(9)

式(8)、(9)中，當(dāng)ΔΘ(φ,λ)為地面重力的變化量時，G(ψ)為重力格林函數(shù)；當(dāng)ΔΘ(φ,λ)為地面大地高的變化量時，G(ψ)為徑向格林函數(shù)；Pn(cosψ)為n階Legendre函數(shù)。

2 研究區(qū)與實驗數(shù)據(jù)處理

2.1 研究區(qū)概況

本文將滇西地區(qū)作為研究區(qū)，包括保山、大理兩個地級市及周邊區(qū)域(97.5°～101.5°E，23.5°～27°N)。該區(qū)地貌以高山為主，怒江、瀾滄江、金沙江等河流流經(jīng)該區(qū)。受地形地貌多樣性影響，雨季局部地區(qū)暴雨頻發(fā)。

2.2 實驗數(shù)據(jù)

2.2.1 全球大氣壓變化數(shù)據(jù)

全球大氣壓數(shù)據(jù)采用歐洲中期天氣預(yù)報中心(European Centre for Medium-Range Weather Forecasts)的全球氣候和天氣第五代再分析產(chǎn)品(http:∥cds.climate.copernicus.eu)，該產(chǎn)品在其前身基礎(chǔ)上采用先進(jìn)的數(shù)據(jù)同化和模式系統(tǒng)， 可用來估計更為準(zhǔn)確的大氣狀況。本次實驗下載全球范圍數(shù)據(jù)，空間分辨率為0.25°× 0.25°，時間分辨率為h，時間跨度為2018-01-01 00:00～2020-12-31 23:00。利用MATLAB程序提取原始NC文件中的surface press(表面壓力)，將每天24 h數(shù)據(jù)取平均值作為日值數(shù)據(jù)。圖1為2018-01-01全球大氣壓數(shù)據(jù)分布圖。

圖1 全球大氣壓數(shù)據(jù)分布Fig.1 Distribution of global atmospheric pressure data

2.2.2 區(qū)域高分辨率大氣壓變化數(shù)據(jù)

區(qū)域高分辨率大氣壓數(shù)據(jù)采用國家氣象科學(xué)數(shù)據(jù)中心中國氣象局陸面數(shù)據(jù)同化系統(tǒng)(CLDAS-V2.0)近實時產(chǎn)品數(shù)據(jù)集(http:∥data.cma.cn)，范圍為96°～103°E、22°～29°N，空間分辨率為0.062 5°×0.062 5°，時間分辨率為d，時間跨度為2018-01-01～2020-12-31。CLDAS-V2.0產(chǎn)品利用多種來源于地面、衛(wèi)星的觀測資料，范圍覆蓋亞洲區(qū)域，在中國區(qū)域質(zhì)量優(yōu)于國際同類產(chǎn)品，且時空分辨率更高?？紤]到格林積分函數(shù)計算負(fù)荷形變時是對區(qū)域進(jìn)行卷積積分，在選擇區(qū)域高分辨率大氣壓變化數(shù)據(jù)范圍時適當(dāng)擴(kuò)展區(qū)域范圍。圖2為2018-01-01區(qū)域高分辨率大氣壓變數(shù)據(jù)分布。

圖2 區(qū)域高分辨率大氣壓數(shù)據(jù)分布Fig.2 Distribution of regional high resolution atmospheric pressure data

2.2.3 CORS站大地高變化數(shù)據(jù)

通過解算滇西地區(qū)內(nèi)72個CORS站每天的觀測數(shù)據(jù)，并扣除固體潮、海潮等潮汐影響，得到CORS站大地高變化時間序列，時間跨度為2018-01-01～2020-12-31，時間分辨率為d。

2.3 數(shù)據(jù)處理流程

本文將全球大氣壓數(shù)據(jù)和區(qū)域高分辨率大氣壓數(shù)據(jù)轉(zhuǎn)換為數(shù)字格網(wǎng)模型，然后將其轉(zhuǎn)換為等效水高格網(wǎng)模型(地表大氣負(fù)荷質(zhì)量變化)，再采用距離反比插值方法將全球大氣等效水高格網(wǎng)的空間分辨率由0.25°轉(zhuǎn)為0.5°，區(qū)域高分辨率等效水高格網(wǎng)的空間分辨率由0.062 5°轉(zhuǎn)為30″。由區(qū)域高分辨率大氣壓變化數(shù)據(jù)，采用基于全球負(fù)荷形變場球諧系數(shù)模型的移去-恢復(fù)法，計算區(qū)域高精度大氣負(fù)荷形變場和時變重力場。計算流程為：

1)對全球等效水高變化格網(wǎng)進(jìn)行球諧分析展開，得到 360 階次負(fù)荷變化球諧系數(shù)模型，并通過球諧綜合計算得到滇西地區(qū)參考等效水高和負(fù)荷影響模型值格網(wǎng)；

2)將區(qū)域等效水高變化格網(wǎng)數(shù)據(jù)移去參考等效水高，獲得區(qū)域等效水高變化殘差格網(wǎng)數(shù)據(jù)，并利用負(fù)荷格林函數(shù)積分理論計算區(qū)域殘余負(fù)荷影響模型值格網(wǎng)；

3)將全球負(fù)荷影響模型值格網(wǎng)與區(qū)域殘余負(fù)荷影響模型值格網(wǎng)結(jié)合，恢復(fù)為滇西地區(qū)區(qū)域高精度負(fù)荷影響格網(wǎng)模型。

3 實驗結(jié)果與分析

為了統(tǒng)一基準(zhǔn)，本文將全球大氣壓變化格網(wǎng)和區(qū)域高分辨率大氣壓變化格網(wǎng)的日值數(shù)據(jù)按照GPS周進(jìn)行周平均，得到周值文件，并將2018年52個周值文件的平均值作為基準(zhǔn)扣除，得到大氣壓周變化格網(wǎng)。然后根據(jù)移去-恢復(fù)法計算大氣壓變化對滇西地區(qū)高精度負(fù)荷形變場及時變重力場的影響，計算得到空間分辨率為30″、時間分辨率為周。高精度負(fù)荷形變場和時變重力場可以更好地反映大氣壓變化對滇西地區(qū)的影響。

3.1 地殼形變

地殼垂直形變是地殼在高程方向的變化量，即大地高變化。表1為大氣壓變化對滇西地區(qū)高程異常、大地高等大地測量參數(shù)的統(tǒng)計信息?？梢钥闯?，大氣壓變化對地殼高程方向的影響較大，年變化幅度達(dá)cm級；而對地殼水平方向的影響較小，年變化幅度最大僅為0.944 mm。其中，大氣負(fù)荷對高程異常的影響等于其對正常高影響與大地高影響之差。

表1 大氣負(fù)荷地殼形變信息統(tǒng)計

由于不同時間內(nèi)大氣壓變化存在差異，因此本文選取年變化幅度最大的2019年中每個月第1周的大氣壓變化對大地高的影響進(jìn)行分析。圖3為2019年大氣壓變化對滇西地區(qū)大地高的影響，可以看出，大氣壓變化對大地高的影響具有十分明顯的季節(jié)性變化規(guī)律。由春入夏過程中，滇西地區(qū)大地高變化緩慢抬升，6月份大地高變化量達(dá)到最大值；而從秋到冬，大地高變化又逐漸下降，12月份大地高變化量達(dá)到最小值。

為更好地分析大氣壓負(fù)荷對滇西地區(qū)的影響，選取滇西區(qū)域內(nèi)YNLC(臨滄)、YNLJ(麗江)、YNRL(瑞麗)、YNSD(施甸)、YNYA(姚安)、YAYS(永善)6個CORS站進(jìn)行分析。從原始大地高時間序列中扣除均值并進(jìn)行粗差探測，移去線性項和常數(shù)項，得到非線性大地高時間序列，并按照GPS周進(jìn)行周平均，得到周值數(shù)據(jù)，結(jié)果如圖4所示?？梢钥闯觯嵛鞯貐^(qū)CORS站大地高變化具有十分明顯的年、半年周期規(guī)律，變化幅度達(dá)到-30～30 mm。

圖4 CORS站大地高變化時間序列Fig.4 Time series of geodetic height variation at CORS stations

將利用移去-恢復(fù)法計算的大地高負(fù)荷影響采用距離反比插值方法，得到6個CORS站點處大氣壓變化對大地高的負(fù)荷影響，結(jié)果如圖5所示。對比圖4與圖5可知，大氣負(fù)荷對CORS站大地高的影響與CORS站大地高變化的趨勢一致，且變化幅度為-8～6 mm。圖5中黑色虛線為區(qū)域高分辨率大氣壓數(shù)據(jù)與全球大氣壓模型數(shù)據(jù)計算的大地高負(fù)荷影響的差值，可以看出，兩者差異較小，這可能與研究區(qū)地理位置有關(guān)。由于大氣壓變化與海拔相關(guān)，我國大氣壓變化引起的負(fù)荷影響，南方比北方小，西部比東部小，本文研究區(qū)位于西南地區(qū)，因此計算的負(fù)荷影響較小。

圖5 CORS站大氣負(fù)荷時間序列Fig.5 Time series of atmospheric load at CORS stations

為進(jìn)一步說明在分析GNSS時間序列周期信號特征及地表環(huán)境負(fù)荷形變場時大氣負(fù)荷改正的重要性，本文通過分析扣除大氣負(fù)荷前后CORS站大地高時間序列的均方根(RMS)來反映大氣負(fù)荷影響對CORS站大地高的貢獻(xiàn)。表2為計算結(jié)果，在扣除大氣負(fù)荷后，CORS站非線性大地高時間序列的RMS值減小，大氣負(fù)荷的貢獻(xiàn)可達(dá)8.62%。這一方面明確了大氣壓變化引起的大地高變化在CORS站大地高時間序列中的貢獻(xiàn)率，不同CORS站的大地高變化對大氣負(fù)荷的敏感性有所不同；另一方面也說明，在研究CORS站非線性運(yùn)動中，大氣負(fù)荷改正具有必要性。

表2 大氣負(fù)荷對CORS站大地高時間序列的貢獻(xiàn)

3.2 時變重力場

本文主要從地面重力、擾動重力等方面來說明大氣壓變化對時變重力場的影響。大氣壓變化對地面重力與擾動重力的影響包括負(fù)荷位引起的直接影響和負(fù)荷附加位引起的間接影響。計算過程中將地面重力和擾動重力的直接影響與間接影響相加，得到大氣壓變化對地面重力和擾動重力的總負(fù)荷影響；垂線偏差為垂線偏差南向與垂線偏差西向的矢量和。為了解大氣負(fù)荷對滇西地區(qū)時變重力場的影響，對時變重力場參數(shù)進(jìn)行統(tǒng)計分析，通過計算其年變化振幅來了解其影響的量級規(guī)律(表3)。由表3可知，大氣壓變化對地面重力的影響較大，年變化振幅約為8 μGal；對垂線偏差的影響較小，年變化振幅最大僅為2.4 ms。

表3 大氣負(fù)荷時變重力場信息統(tǒng)計

圖6為2019-01～12第1周大氣壓變化對滇西地區(qū)地面重力的影響分布。可以看出，地面重力與大氣壓成反比，大氣壓越大，所引起的地面重力影響越小，且大氣負(fù)荷對地面重力的影響同樣具有年周期規(guī)律。規(guī)律顯示，大氣負(fù)荷對地面重力的影響在春、夏季具有增大趨勢，在秋、冬季則表現(xiàn)為減小趨勢。滇西地區(qū)雨季暴雨頻發(fā)，主要集中在6～8月，該段時間大氣壓較小，對地面重力的影響較大，最大可達(dá)3.8 μGal。

圖6 大氣壓變化對滇西地區(qū)地面重力的影響Fig.6 Influence of atmospheric pressure change on surface gravity in western Yunnan

4 結(jié) 語

本文以滇西地區(qū)為研究區(qū)，根據(jù)地球重力場理論，引入移去-恢復(fù)思想，結(jié)合全球大氣壓變化數(shù)據(jù)和區(qū)域高分辨率大氣壓變化數(shù)據(jù)，基于負(fù)荷球諧系數(shù)和區(qū)域負(fù)荷格林函數(shù)方法，得到大氣壓變化對滇西地區(qū)高精度負(fù)荷形變場及時變重力場的影響。結(jié)果表明，大氣壓變化引起的負(fù)荷影響具有十分明顯的年周期規(guī)律。大氣壓變化對地殼形變中垂直形變的影響在1 a內(nèi)大于10 mm，且大氣壓變化引起的大地高變化在CORS站大地高時間序列中的貢獻(xiàn)率占8%左右；對水平方向的負(fù)荷影響較小。大氣壓變化對時變重力場中地面重力的影響約為8 μGal，且與大氣變化呈現(xiàn)明顯的負(fù)相關(guān)性。本文研究結(jié)果對負(fù)荷形變場的精化具有一定的參考意義。