謝劭峰 潘清瑩 黃良珂 朱 葛 王義杰 魏朋志

1 桂林理工大學(xué)測繪地理信息學(xué)院，桂林市雁山街319號，541006 2 廣西空間信息與測繪重點(diǎn)實(shí)驗(yàn)室，桂林市雁山街319號，541006

當(dāng)無線電信號通過對流層時(shí)，傳播路徑會發(fā)生彎曲而產(chǎn)生對流層延遲，這種延遲無法通過雙頻改正的方法消除[1]，需要采用模型改正法消除。天頂對流層延遲(zenith total delay, ZTD)由天頂靜力學(xué)延遲(zenith hydrostatic delay, ZHD)和天頂濕延遲(zenith wet delay, ZWD)組成。隨著大氣再分析資料積分計(jì)算的對流層延遲在GNSS精密定位和GNSS大氣反演中的廣泛應(yīng)用，提高對流層延遲信息的精度具有重要意義。

大氣再分析資料的格網(wǎng)點(diǎn)高程與探空站、GNSS 站等用戶位置高程不一致，且ZTD/ZWD在垂直方向的變化遠(yuǎn)大于水平方向，因此 ZTD/ZWD 垂直剖面模型是實(shí)現(xiàn)高精度對流層格網(wǎng)產(chǎn)品空間插值的關(guān)鍵。部分學(xué)者基于二次多項(xiàng)式、負(fù)指數(shù)函數(shù)或高斯函數(shù)構(gòu)建對流層垂直剖面模型[2-7]，所得的區(qū)域或全球?qū)α鲗友舆t垂直剖面模型具有各自的優(yōu)勢，但仍然存在時(shí)空分辨率低、所需參數(shù)較多、部分區(qū)域適用性較差等不足，中國區(qū)域高精度、高分辨率、實(shí)時(shí)ZTD和ZWD垂直剖面模型仍然缺乏。

本文利用高時(shí)空分辨率的MERRA-2再分析資料積分計(jì)算的分層ZTD、ZWD剖面信息，基于負(fù)指數(shù)函數(shù)表達(dá)ZTD、ZWD垂直剖面，詳細(xì)分析ZTD、ZWD高程縮放因子的時(shí)空分布特性，為構(gòu)建中國區(qū)域高精度ZTD、ZWD垂直剖面模型提供參考。

1 數(shù)據(jù)來源

MERRA-2是由美國宇航局(NASA)推出的高時(shí)空分辨率再分析資料(https:∥goldsmr4.gesdisc.eosdis.nasa.gov/data/MERRA2)，可提供地表氣象參數(shù)和分層氣象參數(shù)，水平分辨率均為0.5°×0.625°(緯度×經(jīng)度)，分層氣象參數(shù)的垂直分辨率為72層和42層標(biāo)準(zhǔn)氣壓層。本文選取42層分層數(shù)據(jù)，地表氣象參數(shù)的時(shí)間分辨率為1 h，分層氣象參數(shù)的時(shí)間分辨率為6 h。因需同時(shí)用到地表氣象參數(shù)和分層氣象參數(shù)，為統(tǒng)一時(shí)間分辨率，地表氣象參數(shù)采樣率取6 h。根據(jù)MERRA-2大氣再分析資料提供的分層數(shù)據(jù)，通過積分法可得到ZTD、ZWD在每個(gè)格網(wǎng)點(diǎn)的垂直剖面，表達(dá)式如下[8]：

(1)

(2)

(3)

(4)

(5)

ZTD=ZHD+ZWD

(6)

式中，e為水汽壓(單位hPa)，Sh為比濕，P為大氣壓(單位hPa)，N為總折射率，Nw為濕折射率，T為溫度(單位K)，H為高程(單位km)，hL為大氣資料積分計(jì)算的最底層高度，htop為大氣積分資料積分計(jì)算的最頂層高度，k1=77.604 K/Pa、k2=64.79 K/Pa、k3=375 463 K2/hPa、k′2=22.97 K/Pa為常數(shù)系數(shù)。

由于大氣再分析資料MERRA-2頂層還有殘余大氣，需要利用Saastamoinen模型計(jì)算殘余大氣延遲值，并將其附加到格網(wǎng)點(diǎn)每一層的積分結(jié)果上。對流層層頂?shù)膶α鲗訚裱舆t較小，可忽略不計(jì)，因此僅計(jì)算ZHD即可：

ZHDSaas=

(7)

式中，Ptop為對流層頂層大氣壓(單位hPa)，φ為緯度，ZHDSaas為利用Saastamoinen模型計(jì)算的MERRA-2層頂?shù)臍堄郱HD值。為便于積分計(jì)算，對式(4)和式(5)進(jìn)行離散化：

(8)

(9)

式中，Ni、Nwi和ΔHi分別為第i層的大氣總折射率、濕折射率和大氣厚度，n為積分大氣層數(shù)。

Huang等[9]利用全球無線電探空數(shù)據(jù)和國際全球?qū)Ш叫l(wèi)星系統(tǒng)服務(wù)的精確ZTD產(chǎn)品對大氣再分析資料MERRA-2得到的ZTD和ZWD的性能進(jìn)行評估。結(jié)果表明，全球范圍的MERRA-2再分析資料得到的ZTD和ZWD值在中國區(qū)域具有較高的精度和較好的穩(wěn)定性，因此可將MERRA-2再分析資料作為中國區(qū)域ZTD/ZWD垂直剖面格網(wǎng)模型構(gòu)建的數(shù)據(jù)源。本文選取2014～2016年70°～135°E、15°～55°N范圍內(nèi)MERRA-2大氣再分析分層數(shù)據(jù)及對應(yīng)的地表數(shù)據(jù)進(jìn)行積分計(jì)算，得到該范圍內(nèi)各格網(wǎng)點(diǎn)ZTD、ZWD的垂直剖面資料。

2 ZWD、ZTD高程縮放因子獲取

用二次項(xiàng)、高斯函數(shù)表示ZTD、ZWD在垂直方向的變化[2-7]，用指數(shù)函數(shù)表示ZTD、ZWD與高程之間的關(guān)系[4,6-7]：

ZTD(h)=a1·eb1·h

(10)

ZWD(h)=a2·eb2·h

(11)

式中，e為自然對數(shù)的底數(shù)，h為高程，ZTD(h)、ZWD(h)分別為高程h處的ZTD、ZWD值，a1、a2、b1、b2為模型參數(shù)。

為進(jìn)一步驗(yàn)證指數(shù)函數(shù)的適用性，選取2015-01-01 00:00(UTC)中國區(qū)域3個(gè)具有代表性的MERRA-2格網(wǎng)點(diǎn)，根據(jù)式(1)～(9)分別積分計(jì)算不同等壓層面的ZTD、ZWD信息，并通過負(fù)指數(shù)函數(shù)分別對ZTD、ZWD進(jìn)行擬合，結(jié)果如圖1和圖2所示。

圖1 ZTD值在高程上的變化與指數(shù)函數(shù)擬合結(jié)果及其擬合殘差Fig.1 The changes of ZTD values with elevation andfitting results of exponential function and fitting residual

圖2 ZWD值在高程上的變化與指數(shù)函數(shù)擬合結(jié)果及其擬合殘差Fig.2 The change of ZWD value with elevation and fitting results of exponential function and fitting residual

從圖1和圖2可以看出，利用負(fù)指數(shù)函數(shù)可較好地表達(dá)ZTD、ZWD在高程上的變化，ZTD擬合偏差在5.35 mm以內(nèi)，擬合RMS誤差在18.72 mm以內(nèi)，ZWD擬合偏差在1.99 mm以內(nèi)，擬合RMS誤差在8.25 mm以內(nèi)，表現(xiàn)出較高的擬合精度，因此可采用負(fù)指數(shù)函數(shù)來表達(dá)ZTD、ZWD的垂直剖面。位置A、B處的ZTD、ZWD可分別用式(12)～(15)表示，通過推導(dǎo)可得式(16)和式(17)：

ZTDA=a1·eb1·HA

(12)

ZTDB=a1·eb1·HB

(13)

ZWDA=a2·eb2·HA

(14)

ZWDB=a2·eb2·HB

(15)

ZTDA=ZTDB·eb1(HA-HB)

(16)

ZWDA=ZWDB·eb2(HA-HB)

(17)

式中，ZTDA、ZTDB、ZWDA、ZWDB分別為位置A、B處ZTD、ZWD值，HA、HB分別為位置A、B處高程。由于b1、b2較小，為方便表示，對式(16)和式(17)進(jìn)行改寫：

(18)

(19)

式中，ZTDt和ZWDt分別為目標(biāo)高程Ht處ZTD和ZWD值，ZTDr和ZWDr分別為參考高程Hr處ZTD和ZWD值，Hs和Hw分別為ZTD和ZWD的高程縮放因子(單位km)。本文中參考高程處ZTD、ZWD值為利用MERRA-2地表氣象參數(shù)計(jì)算得到的ZTD、ZWD值，目標(biāo)高程處ZTD、ZWD值為利用MERRA-2分層氣象參數(shù)計(jì)算得到的ZTD、ZWD值，最終可得到Hs和Hw。

3 ZWD、ZTD高程縮放因子的時(shí)空特性分析

圖3和圖4分別為2015年不同季節(jié)ZTD、ZWD高程縮放因子在中國區(qū)域的分布情況。從圖3可以看出，ZTD高程縮放因子最小值7.55 km出現(xiàn)在低緯度地區(qū)，最大值8.23 km出現(xiàn)在西部地區(qū)。春季和秋季ZTD高程縮放因子分布狀況相似，且秋季低緯度地區(qū)季均ZTD高程縮放因子較??；夏季ZTD高程縮放因子最大值出現(xiàn)在40°～50°N附近，呈現(xiàn)向低緯度區(qū)域減小的趨勢；冬季ZTD高程縮放因子最大值出現(xiàn)在20°～30°N附近，呈現(xiàn)向中緯度區(qū)域減小的趨勢。從4季尺度來看，中緯度地區(qū)ZTD高程縮放因子先增大后減小，而低緯度地區(qū)先減小后增大。由圖4可知，ZWD高程縮放因子最小值、最大值均出現(xiàn)在低緯度地區(qū)，最小值為1.35 km，最大值為3.21 km。春季和秋季ZWD高程縮放因子分布狀況相似，且秋季低緯度地區(qū)季均ZWD高程縮放因子較?。幌募竞投綵WD高程縮放因子呈現(xiàn)相反的分布趨勢，夏季最大值出現(xiàn)在20°～30°N和40°N附近，呈現(xiàn)向中緯度區(qū)域減小的趨勢，冬季最大值出現(xiàn)在30°～40°N附近，呈現(xiàn)向低緯度區(qū)域減小的趨勢。由圖3(b)可知，青藏高原地區(qū)(20°00′～39°47′N，73°19′～104°47′E)ZTD高程縮放因子與其附近地區(qū)相比存在明顯差異；圖4(b)中青藏高原地區(qū)ZWD高程縮放因子與其附近地區(qū)相比也相對較小，這可能是由于青藏高原平均海拔在4 000 m以上、夏季溫度低于其他地區(qū)的緣故。由圖4可知，從春季到夏季，中國區(qū)域ZWD高程縮放因子呈現(xiàn)增大的趨勢，而從夏季至冬季呈現(xiàn)減小的趨勢。值得注意的是，40°N附近、75°～105°E區(qū)域(即塔里木盆地)ZTD、ZWD高程縮放因子全年均處于較大值，可能是由于該區(qū)域位于內(nèi)陸，且為非季風(fēng)性氣候，受地形因素影響，全年較干旱。

圖3 2015年中國區(qū)域ZTD高程縮放因子季均值分布Fig.3 Distribution of seasonal mean value of ZTD height scale factor in China in 2015

圖4 2015年中國區(qū)域ZWD高程縮放因子季均值分布Fig.4 Distribution of seasonal mean value of ZWD height scale factor in China in 2015

綜上分析可知，ZTD、ZWD高程縮放因子均呈現(xiàn)一定的季節(jié)性分布差異，我國南北跨緯度較大，氣溫呈現(xiàn)一定的差異，因此可推斷ZTD、ZWD高程縮放因子的分布與地理位置有關(guān)，且具有季節(jié)性。ZTD、ZWD高程縮放因子的時(shí)空分布存在差異，因此需分別構(gòu)建模型。

對流層延遲具有顯著的年周期和半年周期變化，為探究ZTD、ZWD高程縮放因子是否具有周期性，本文隨機(jī)選取4個(gè)格網(wǎng)點(diǎn)進(jìn)行分析，圖5和圖6為其Hs、Hw時(shí)間序列及傅里葉頻譜分析。

由圖5可知，Hs時(shí)間序列在中國區(qū)域呈現(xiàn)出明顯的年周期和半年周期特性，波動范圍為7.4～8.5 km，變化范圍相對穩(wěn)定；相對于另外3個(gè)格網(wǎng)點(diǎn)，格網(wǎng)點(diǎn)(20°N, 80°E)表現(xiàn)出更為明顯的年周期和半年周期特性。

圖5 ZTD高程縮放因子時(shí)間序列及傅里葉頻譜分析Fig.5 Time series and Fourier spectrum analysis of ZTD height scale factor

由圖6可知，Hw在中國區(qū)域也表現(xiàn)出顯著的年周期和半年周期特性，尤其在格網(wǎng)點(diǎn)(20°N, 80°E)和格網(wǎng)點(diǎn)(20°N, 105°E)處；在格網(wǎng)點(diǎn)(40°N, 80°E)處Hw變化范圍相對較大，為1.5～6 km，另外3個(gè)格網(wǎng)點(diǎn)Hw變化范圍相對較小，為1～4 km。

圖6 ZWD的高程縮放因子的時(shí)間序列及傅里葉頻譜分析Fig.6 Time series and Fourier spectrum analysis of ZWD height scale factor

分別計(jì)算2014～2016年中國區(qū)域MERRA-2再分析資料計(jì)算的每個(gè)格網(wǎng)點(diǎn)Hs、Hw的年均值、年周期振幅及半年周期振幅值，以進(jìn)一步分析Hs、Hw在中國區(qū)域的分布特性，結(jié)果見圖7和圖8。

圖7 ZTD高程縮放因子年均值、年周期振幅、半年周期振幅Fig.7 The annual mean value, annual periodic amplitude, and semi-annual periodic amplitudeof the ZTD height scale factor

圖8 ZWD高程縮放因子年均值、年周期振幅、半年周期振幅Fig.8 The annual mean value, annual periodic amplitude, and semi-annual periodic amplitude of the ZWD height scale factor

由圖7可知，Hs在中國西部到中部的三角形區(qū)域存在較大的年均值，可能是由于該地區(qū)以高原為主，平均海拔在4 000 m以上，地形起伏較大。年周期振幅在中緯度地區(qū)較小，從中緯度地區(qū)向西南部、東北部呈現(xiàn)增大趨勢，且低緯度地區(qū)西南部數(shù)值較大?？赡苁怯捎谥袊途暥鹊貐^(qū)屬于亞熱帶季風(fēng)氣候，水汽變化較其他區(qū)域大，進(jìn)而導(dǎo)致Hs的年周期振幅較大。中國西南地區(qū)和東北地區(qū)存在相對較大的半年周期振幅，西南地區(qū)數(shù)值較大，而西北地區(qū)至東南地區(qū)半年周期振幅較小。

由圖8可知，Hw在中國西北部和東南部存在較大的年均值，其中，中國西北部為塔里木盆地，處于非季風(fēng)區(qū)，位于內(nèi)陸，終年受大陸氣團(tuán)控制，冬季嚴(yán)寒，夏季高溫，常年少雨；而中國西南部處于非季風(fēng)區(qū)，位于沿海地區(qū)，該區(qū)域包含溫帶季風(fēng)氣候和亞熱帶季風(fēng)氣候，受冬夏季風(fēng)交替控制，水汽充足，冬季寒冷干燥，夏季高溫多雨。由此可見，ZWD高程縮放因子受多種因素影響。年周期振幅在緯度上表現(xiàn)為由低緯度向高緯度逐漸遞減，低緯度地區(qū)存在較大值，中緯度地區(qū)存在較小值。可能是由于中國低緯度地區(qū)屬于亞熱帶季風(fēng)氣候，水汽變化較其他區(qū)域大，進(jìn)而導(dǎo)致Hw的年周期振幅較大。中國西南部和東北部存在相對較大的半年周期振幅。

綜上，為保證ZTD、ZWD垂直剖面模型的精度，在構(gòu)建中國區(qū)域ZTD、ZWD垂直模型時(shí)需顧及Hs、Hw的年周期和半年周期變化。

利用皮爾遜相關(guān)性分析法分析Hs、Hw的年均值與經(jīng)緯度的相關(guān)性，結(jié)果見表1。由表1可知，Hs與經(jīng)度呈負(fù)相關(guān)，相關(guān)性顯著，與緯度呈正相關(guān)，相關(guān)性顯著；Hw與經(jīng)度相關(guān)性弱，與緯度呈負(fù)相關(guān)，相關(guān)性顯著；Hs與Hw呈正相關(guān)。

表1 相關(guān)性分析

4 結(jié) 語

針對中國區(qū)域缺乏實(shí)時(shí)、高時(shí)空分辨率的ZTD、ZWD垂直剖面模型，本文利用2014～2016年MERRA-2再分析數(shù)據(jù)，基于負(fù)指數(shù)函數(shù)，詳細(xì)分析ZTD、ZWD高程縮放因子在中國區(qū)域的時(shí)空分布特性。結(jié)果表明：

1)夏季Hs最大值出現(xiàn)在40°～50°N附近，呈現(xiàn)向低緯度區(qū)域減小的趨勢；冬季Hs最大值出現(xiàn)在20°～30°N附近，呈現(xiàn)向中緯度地區(qū)減小的趨勢。

2)Hs在中國西部到中部的三角形區(qū)域存在較大的年均值，中緯度地區(qū)先增大后減小，而低緯度地區(qū)先減小后增大，Hs年周期振幅、半年周期振幅具有顯著的空間分布特性。

3)Hw的最值出現(xiàn)在低緯度地區(qū)，夏季Hw最大值出現(xiàn)在20°～30°N附近及40°N附近，呈現(xiàn)向中緯度區(qū)域減小的趨勢；冬季Hw最大值出現(xiàn)在30°～40°N附近，呈現(xiàn)向低緯度地區(qū)減小的趨勢。

4)Hw在中國區(qū)域西北部和東南部存在較大的年均值，Hw的年周期振幅、半年周期振幅也具有顯著的空間分布特性。

綜上可知，在構(gòu)建中國區(qū)域ZTD、ZWD垂直剖面模型時(shí)需顧及Hs、Hw的年周期、半年周期特性，且Hs、Hw都與經(jīng)緯度具有一定的相關(guān)性。研究結(jié)果可為ZTD、ZWD垂直剖面模型的構(gòu)建提供參考。在后續(xù)工作中可進(jìn)一步研究構(gòu)建中國區(qū)域ZTD、ZWD垂直剖面模型，以實(shí)現(xiàn)高精度的對流層格網(wǎng)產(chǎn)品空間插值。