基于CRA40產品的對流層延遲計算及對比分析

周要宗，樓益棟，張衛(wèi)星，梁 宏，施 闖，吳 迪，曹云昌

1. 武漢大學衛(wèi)星導航定位技術研究中心，湖北 武漢 430079； 2. 中國氣象局氣象探測中心，北京 100081； 3. 北京航空航天大學衛(wèi)星導航與移動通信融合技術工業(yè)和信息化部重點實驗室，北京 100191； 4. 南寧師范大學地理與海洋研究院，廣西 南寧 530001

對流層延遲是全球衛(wèi)星導航系統(tǒng)(GNSS)、甚長基線干涉測量(very long baseline interferometry，VLBI)、衛(wèi)星激光測距(satellite laser ranging，SLR)、衛(wèi)星測高等空間大地測量技術實現(xiàn)高精度參數(shù)估計的主要誤差源之一，在數(shù)據(jù)處理中可利用外部數(shù)據(jù)產品估算對流層延遲并予以改正[1-3]。全球大氣再分析產品，如歐洲中尺度天氣預報中心(European Centre for Medium-range Weather Forecasts，ECMWF)再分析中間產品(ECMWF reanalysis interim，ERA-Interim，簡稱“ERAI”)和第5代產品(the fifth generation ECMWF reanalysis，ERA5)、美國國家環(huán)境預報中心再分析產品(The United States National Centers for environmental prediction reanalysis，NCEP)，通過同化多源歷史氣象觀測資料，具有精度高、時空覆蓋完整的優(yōu)點[4]，被廣泛用于對流層天頂延遲(zenith path delay，ZPD)和斜路徑延遲(slant path delay，SPD)的計算及建模[5-8]，發(fā)布了一系列ZPD經驗模型(如SHAO系列模型[9]、IGGtrop系列模型[10]、TropGrid系列模型[11]、ITG模型[12]、GTrop模型[13]和HGPT模型[14])和映射函數(shù)模型(如UNB-VMF1[15]、VMF3[16]和GFZ-VMF1[17])。

2013年11月，中國氣象局啟動了我國的全球大氣再分析計劃，總體目標是建成我國第1代全球大氣再分析業(yè)務系統(tǒng)，并建成跨度40 a(1979—2018)的中國第1代全球大氣/陸面再分析產品(China 's first generation 40 a global atmosphere and land reanalysis，CRA40)[18-19]，質量超過國際第2代，在中國區(qū)域接近或達到國際第3代大氣再分析水平(http:∥data.cma.cn/analysis/cra40)。2020年12月CRA40產品正式上線并業(yè)務化發(fā)布[20]。探討利用CRA40進行對流層延遲計算的方法，并評估CRA40用于對流層延遲計算的精度，尤其是在中國區(qū)域的精度，具有重要的意義。文獻[21]針對CRA40在中國區(qū)域GNSS水汽反演中的應用，對2016年全年CRA40和ERA5的氣溫、氣壓、大氣水汽加權平均溫度和天頂總延遲(zenith total delay，ZTD)進行了較為系統(tǒng)的比較評估，結果表明以GNSS ZTD為參考，CRA40 ZTD精度為1.35 cm，總體較ERA5差。但文獻[21]主要討論的是中國區(qū)域的天頂延遲結果，并未開展全球范圍內包括斜路徑延遲在內的對流層延遲的系統(tǒng)評估。

本文主要對基于CRA40計算的對流層延遲開展系統(tǒng)的評估工作，首先討論基于CRA40采用射線追蹤技術計算對流層延遲的方法，然后分別在231個國際GNSS服務(international GNSS service，IGS)站和213個中國大陸構造環(huán)境監(jiān)測網絡(crustal movement observation network of China，CMONOC)站處對ERAI、ERA5和CRA40對流層延遲計算精度(ZPD和5°高度角SPD)進行系統(tǒng)比較評估，并對中國區(qū)域CRA40對流層濕延遲計算精度的變化特征進行分析。

1 數(shù)據(jù)與方法

1.1 CRA40產品及預處理

CRA40(1979—2018)是由中國氣象局國家氣象信息中心牽頭，聯(lián)合多家單位開發(fā)的，整體經歷了兩個階段，即準備階段(2014—2016年)和實施階段(2016—2019年)[20]。準備階段主要完成歷史氣象觀測資料的準備和模式同化系統(tǒng)的建設[22-23]，實施階段主要完成CRA-Interim(2006年7月—2016年12月)產品的制作和評估及CRA40產品的研制[24-26]。CRA40產品是分10段進行研制的，每段4.5 a，段與段之間有半年重疊，以保證產品的連續(xù)性。2019年8月全部40 a CRA40再分析產品研制完成，并于2020年12月正式上線業(yè)務化[20]。本文采用CRA40氣壓層產品進行對流層延遲計算，產品模式同化系統(tǒng)為GFS/GSI-3Dvar，時間分辨率為6 h。產品水平格網采用高斯經緯度投影，經向為均勻網格劃分，分辨率為0.312 5°，緯向為非均勻網格劃分，分辨率在0.312 5°左右。選用參數(shù)包括位勢高、溫度和比濕，位勢高和溫度參數(shù)氣壓層層數(shù)為47，層頂氣壓為1 hPa(42 km左右)，比濕參數(shù)氣壓層層數(shù)為37，層頂氣壓為100 hPa(16 km左右)。

相關研究表明：同一再分析產品隨著水平分辨率的提高，射線追蹤的計算效率會顯著降低，而對射線追蹤ZTD計算精度的影響并不顯著[7]，因此綜合考慮射線追蹤計算效率和精度，本文采用NCEP官網提供的wgrib2工具利用雙線性插值方法將再分析資料的水平分辨率統(tǒng)一重采樣為1°，來開展后續(xù)的再分析資料評估試驗。此外，CRA40比濕參數(shù)僅有37層，與位勢高和溫度參數(shù)氣壓層層數(shù)不一致，因此需要對比濕參數(shù)層頂100 hPa以上氣壓層(1、2、3、5、7、10、20、30、50、70 hPa)的比濕參數(shù)進行處理，得到氣壓層層數(shù)為47的比濕參數(shù)?？蛇x的處理方法主要有3種：方法1認為100 hPa(16 km)以上水汽已經相當稀疏，因此，可以簡單地將100 hPa以上氣壓層的比濕參數(shù)值設為0 g·kg-1；方法2由現(xiàn)有層頂100 hPa和相鄰氣壓層(如125、150、175 hPa等)的比濕參數(shù)來線性外推層頂以上的比濕參數(shù)；方法3直接采用100 hPa的比濕參數(shù)值作為其他層頂氣壓層的數(shù)值。

為評估以上方法的差異，圖1(a)給出了2018年7月19日UTC 15時BZRG測站基于ERA5計算的天頂濕延遲(zenith wet delay，ZWD)隨高度的變化曲線，圖1(b)給出了2018年7月ERA5月平均產品全球比濕參數(shù)均值隨氣壓(37層)的變化曲線。由圖1可知，16 km(100 hPa)以上ZWD數(shù)值在亞毫米量級，100 hPa以上比濕參數(shù)并不為0，幾乎處于一個不變值，因此方法1可行但不嚴密，而方法3相對更為合理。相比之下，方法2利用層頂100 hPa和相鄰氣壓層比濕來線性外推頂層比濕，可能會得到小于0 g·kg-1的比濕。綜合考慮計算的可行性和合理性，本文采用方法3處理頂層比濕參數(shù)。經預處理后，CRA40產品時間和水平分辨率分別為6 h和1.0°，氣壓層層數(shù)為47，層頂氣壓為1 hPa(42 km)。

圖1 ZWD隨高度和比濕隨氣壓的變化曲線Fig.1 Changing curves of ZWD with height and specific humidity with pressure

1.2 ERA5和ERAI產品

ERA5(1950年至今)是ECMWF發(fā)布的最新一代全球大氣再分析產品，產品模式同化方法為Ensemble of 4D-Var，時間延遲在5 d左右，時間分辨率為1 h，水平分辨率為31 km，氣壓層產品垂直層數(shù)為37[27-28]。ERAI(1979—2019年)是ERA5的上一代全球大氣再分析產品，產品模式同化方法為4D-Var，時間分辨率為6 h，水平分辨率為80 km，氣壓層產品垂直層數(shù)為37[29]。本文選用ERA5和ERAI氣壓層產品進行對流層延遲計算，參數(shù)包括位勢、溫度和比濕，產品時間和水平分辨率與CRA40重采樣產品保持一致，分別為6 h和1.0°，層頂氣壓為1 hPa。

1.3 基于射線追蹤技術的對流層延遲計算方法

組合使用再分析產品(0～42 km)和美國標準大氣模型(COESA 1976)(42～84 km)來實現(xiàn)參數(shù)的中性大氣層全覆蓋[30]。依據(jù)文獻[31]給出的方法，采用EGM2008模型[32]將再分析位勢高/位勢轉換成大地高。利用文獻[33]給出的方法將再分析比濕轉換為水汽壓。采用文獻[34]給出的大氣垂直分層方法(0～2 km：10 m間隔；2～6 km：20 m間隔；6～16 km：50 m間隔；16～36 km：100 m間隔；36～84 km：500 m間隔)來進行垂直大氣細化分層。分別使用線性、指數(shù)和指數(shù)插值方法將溫度、氣壓和水汽壓由等壓層內插到細化高度層[30]。以溫度、氣壓和水汽壓為輸入，通過式(1)來計算靜力學折射率(nh)、濕折射率(nw)和總折射率(n)[35]

(1)

(2)

分別以6 h和1.0°時空分辨率的CRA40(47層)、ERA5(37層)和ERAI(37層)再分析溫度、比濕和位勢高/位勢為輸入，采用前述的對流層延遲射線追蹤方法，計算2018年全年全球231個IGS站處(圖2(a))和國內213個CMONOC站處(圖2(b))的ZPD和5°高度角、16個方位角(0°～337.5°，間隔22.5°)下的SPD。取16個方位角的SPD均值作為5°高度角的SPD。最終得到3種再分析產品在所有選定測站處2018年時間分辨率為6 h的ZPD和5°高度角SPD產品。

圖2 IGS和CMONOC站點分布Fig.2 Distribution of selected IGS and CMONOC stations

2 對流層天頂總延遲精度評估

事后GNSS ZTD未被CRA40、ERAI和ERA5同化[26-29]，因此可以將其作為獨立參考，來評估基于再分析資料計算的ZTD精度。從IGS官方FTP(gdc.cddis.eosdis.nasa.gov)下載2018年全年選定231個IGS站處時間分辨率為5 min的ZTD產品，產品統(tǒng)計精度可達4 mm[39](https:∥www.igs.org/products#about)。采用PANDA軟件和文獻[40]所示的數(shù)據(jù)處理策略對所有選定213個CMONOC測站的GNSS觀測數(shù)據(jù)進行事后精密單點定位(precise point positioning，PPP)處理[41]，得到選定CMONOC測站2018年全年時間分辨率為1 h的ZTD產品，產品統(tǒng)計精度優(yōu)于5 mm[40]。

以GNSS ZTD為參考，統(tǒng)計2018年全年每個選定測站基于再分資料計算的ZTD誤差偏差(Bias)和均方根(root mean square，RMS)。圖3給出了全球范圍(IGS站網)和中國區(qū)域(CMONOC站網)基于3種再分析產品的ZTD Bias地理分布和直方圖，圖3(a)、(b)和(c)給出全球范圍CRA40、ERA5和ERAI的ZTD Bias分布，圖3(d)、(e)和(f)給出相應的統(tǒng)計直方圖，圖3(g)、(h)和(i)給出中國區(qū)域CRA40、ERA5和ERAI的ZTD Bias分布，而圖3(j)、(k)和(l)則給出相應的統(tǒng)計直方圖。在全球范圍，ERA5 ZTD Bias平均最小，均值為0.20 cm(圖3(b)、(e))。ERAI ZTD Bias所有測站均值為0.34 cm，ZTD偏大的測站主要分布在歐洲、中亞、北美等地區(qū)(圖3(c)、(f))。相比之下，CRA40 ZTD Bias最大，均值為0.48 cm，主要出現(xiàn)在歐洲、北美及赤道區(qū)域(圖3(a)、(d))。

圖3 ZTD差值Bias分布和直方圖Fig.3 Distribution and histogram of ZTD Bias

在中國區(qū)域，ERA5 ZTD Bias仍然最小，均值為0.02 cm(圖3(h)、(k))。ERAI ZTD Bias次之，均值為0.11 cm(圖3(i)、(l))。CRA40 ZTD Bias同樣最大，平均為0.28 cm，川滇地區(qū)、新疆北部的ZTD正偏差較為明顯，負偏差則主要出現(xiàn)在新疆南部區(qū)域(圖3(g)、(j))。再分析ZTD整體偏差均為正，表明再分析ZTD與GNSS ZTD間可能存在系統(tǒng)性偏差(再分析ZTD數(shù)值偏大)，文獻[42—43]也給出了類似的結果。

除了平均偏差外，圖4也相應地給出了基于3種再分析資料計算的ZTD RMS地理分布和直方圖，圖4(a)、(b)和(c)分別給出全球范圍CRA40、ERA5和ERAI的ZTD RMS分布，圖4(d)、(e)和(f)分別給出相應的統(tǒng)計直方圖，圖4(g)、(h)和(i)分別給出中國區(qū)域CRA40、ERA5和ERAI的ZTD RMS分布，而圖4(j)、(k)和(l)則給出相應的統(tǒng)計直方圖。在全球范圍，再分析ZTD RMS與測站緯度相關，赤道及其附近區(qū)域RMS明顯偏大，這主要是受赤道及其附近區(qū)域含量豐富且變化劇烈水汽的影響[44]。ERA5 ZTD RMS整體最小，平均為1.16 cm(圖4(b)、(e))。CRA40 ZTD RMS次之，平均為1.39 cm，其中赤道及其附近區(qū)域相比ERA5 RMS明顯更大(圖4(a)、(d))。ERAI ZTD RMS整體最大，平均為1.47 cm，RMS較大的測站主要分布在赤道及其附近區(qū)域以及中東、中亞等地區(qū)(圖4(c)、(f))。

在中國區(qū)域，ERA5 ZTD RMS仍然最好，均值為1.16 cm(圖4(h)、(k))。ERAI ZTD RMS次之，川滇地區(qū)以及新疆、河北、山東、遼寧等地均相對ERA5明顯變差，均值為1.39 cm(圖4(i)、(l))。CRA40 ZTD RMS相對ERAI略差，主要集中在川滇地區(qū)和新疆，均值為1.41 cm(圖4(g)、(j))。整體統(tǒng)計上，ERA5 ZTD RMS明顯小于CRA40和ERAI，而CRA40和ERAI ZTD RMS比較接近，表明ERA5 ZTD精度顯著優(yōu)于ERAI和CRA40，而ERAI和CRA40 ZTD精度相近。

圖4 ZTD差值RMS分布和直方圖Fig.4 Distribution and histogram of ZTD RMS

3 對流層斜路徑延遲精度評估

由于斜路徑延遲缺乏可靠的參考，無法直接對基于再分析資料計算的SPD絕對精度進行評估，而對流層天頂總延遲精度評估結果顯示ERA5 ZTD精度最高，理論上ERA5計算的ZPD和SPD精度也最高，因此就以ERA5計算的ZPD和SPD為參考，來評估CRA40和ERAI計算的ZPD和SPD精度。首先計算2018年全年所有選定測站CRA40和ERAI ZPD和SPD同ERA5 ZPD和SPD的差異Bias和RMS，然后統(tǒng)計出差異Bias和RMS的最小值(min)、最大值(max)和均值(mean)，結果見表1。

由表1可知，5°高度角SPD差異Bias和RMS約是相應ZPD差異RMS的10倍，這是因為5°高度角SPD和ZPD比值是在10.2左右[45]。CRA40和ERA5對流層延遲差異主要體現(xiàn)在濕延遲部分。全球范圍(IGS站網)CRA40和ERA5對流層延遲(ZPD和SPD)差異RMS分別為1.04(ZTD)、0.13(ZHD)、1.04(ZWD)、10.83(STD)、1.37(SHD)和10.80 cm(SWD)，中國區(qū)域(CMONOC站網)為1.18、0.14、1.16、12.30、1.49和12.13 cm。以ERA5對流層延遲為參考，在全球范圍CRA40同ERAI的對流層延遲計算精度基本相當，而在中國區(qū)域ERAI對流層延遲計算精度略優(yōu)于CRA40，可能原因是CRA40在中國區(qū)域同化了一些國內特有的氣象觀測數(shù)據(jù)，導致在中國區(qū)域CRA40與ERA5在同化數(shù)據(jù)源上有更大差異。

表1 再分析ZPD和SPD差異Bias和RMS統(tǒng)計Tab.1 Statistical bias and RMS between reanalysis-based ZPD and SPD cm

4 中國區(qū)域CRA40對流層濕延遲精度變化特征

在中國區(qū)域CRA40對流層總延遲(ZTD和STD)精度主要取決于濕延遲部分(ZWD和SWD)計算精度，因此重點對中國區(qū)域CRA40濕延遲精度的變化特征進行分析。從中國區(qū)域選取4個代表性IGS測站，即BJFS(北京房山，溫帶季風氣候)、JFNG(武漢九峰，亞熱帶季風氣候)、LHAZ(拉薩，高原山地氣候)和URUM(烏魯木齊，溫帶大陸性氣候)，來分析2018年CRA40 ZWD和SWD精度變化特征，結果如圖5所示。圖5中，橫軸表示年積日，縱軸(左，藍色)表示ZWD/SWD差值，藍色曲線表示CRA40 ZWD/SWD同ERA5 ZWD/SWD的差值時序，縱軸(右，紅色)表示ZWD/SWD，紅色曲線表示CRA40 ZWD/SWD時序。

由圖5可知，同一站點的ZWD精度變化相位同SWD一致、變化振幅約是SWD的1/10，而不同站點間ZWD和SWD精度變化相位和振幅存在明顯差異，主要與站點所處氣候類型相關。同屬季風氣候類型的BJFS和JFNG站，ZWD和SWD精度均有明顯的季節(jié)變化特征：夏季精度最差，冬季精度最好，且溫帶BJFS站的ZWD和SWD精度變化振幅明顯小于亞熱帶JFNG站，尤其是在冬季(圖5(a)、(b)、(c)和(d))。高原山地氣候的LHAZ站，ZWD和SWD精度全年較為穩(wěn)定，基本沒有季節(jié)變化(圖5(e)、(f))。溫帶大陸性氣候的URUM站，ZWD和SWD精度也存在一定的季節(jié)變化特征，但變化振幅相對BJFS和JFNG站要小得多(圖5(g)、(h))。

圖5 2018年BJFS、JFNG、LHAZ和URUM站CRA40 ZWD和SWD差值時序Fig.5 Time series of CRA40 ZWD and SWD differences at BJFS, JFNG, LHAZ and URUM stations in 2018

實際上，CRA40 ZWD和SWD精度直接與對流層水汽的含量及變化幅度相關，對流層水汽含量越豐富、變化幅度越劇烈，ZWD和SWD精度就越差，反之亦然。中國區(qū)域季風氣候區(qū)對流層水汽含量及變化幅度夏季大于冬季，因此季風氣候區(qū)ZWD和SWD精度夏季差于冬季，高原山地氣候區(qū)對流層水汽含量較低、變化幅度較小，因此高原山地氣候區(qū)ZWD和SWD精度全年較為穩(wěn)定。

5 結論與展望

本文探討了基于CRA40采用射線追蹤技術計算對流層延遲的方法，在全球范圍和中國區(qū)域對CRA40、ERAI和ERA5對流層延遲(ZPD和5°高度角SPD)的精度進行了系統(tǒng)比較評估，并對中國區(qū)域CRA40 ZWD和SWD精度變化規(guī)律進行了分析。結果表明：

(1) 以GNSS ZTD為參考，CRA40 ZTD統(tǒng)計Bias在全球范圍和中國區(qū)域分別為0.48和0.28 cm；統(tǒng)計RMS約為1.40 cm，在全球范圍略優(yōu)于ERAI，在中國區(qū)域基本與ERAI相當，整體上較ERA5差。

(2) 以ERA5 STD為參考，CRA40 STD統(tǒng)計Bias在全球范圍和中國區(qū)域分別為2.84和2.58 cm；統(tǒng)計RMS在全球范圍和中國區(qū)域分別為10.83和12.30 cm，同ERAI無明顯差異。

(3) 中國區(qū)域CRA40 ZWD和SWD精度變化規(guī)律與氣候類型相關，非季風氣候區(qū)精度全年較為穩(wěn)定，季風氣候區(qū)冬季精度明顯優(yōu)于夏季。

結果證實國產自主可控CRA40產品對流層延遲射線追蹤精度基本與第3代全球大氣再分析ERAI相當，可以作為可靠的數(shù)據(jù)源來進行對流層延遲射線追蹤、建模等相關研究。本文僅對基于CRA40計算的對流層延遲的精度進行了評估，下一步將利用CRA40對流層延遲進行映射函數(shù)、水平梯度建模及北斗/GNSS定位應用研究。

致謝：感謝ECMWF提供的ERAI和ERA5數(shù)據(jù)及中國氣象局提供的CRA40數(shù)據(jù)；感謝IGS提供的IGS ZTD產品和中國地震局提供的GNSS觀測數(shù)據(jù)。