中國地基GNSS/MET 水汽產(chǎn)品質(zhì)量控制及與再分析產(chǎn)品的對比評估

遠(yuǎn)芳 廖捷 周自江

國家氣象信息中心, 北京 100081

1 引言

水汽是大氣中重要的變量，它既在地球氣候系統(tǒng)的能量和水循環(huán)中扮演關(guān)鍵的角色，也是災(zāi)害性天氣形成和演變中的重要因子。地基全球衛(wèi)星導(dǎo)航系統(tǒng)氣象（GNSS/MET）水汽資料利用放置在地面上的接收機(jī)測量GNSS 衛(wèi)星信號穿過大氣層到達(dá)地面時所引起的時間延遲量并進(jìn)一步反演出天頂方向整層大氣或信號斜路徑上的水汽累積量。地基GNSS/MET 水汽探測可以獲取高時效、高時空分辨率的大氣水汽場，有助于更準(zhǔn)確地分析天氣系統(tǒng)的演變特征。隨著地基GNSS/MET 水汽觀測站網(wǎng)不斷加密和資料傳輸業(yè)務(wù)化程度不斷提高，該資料已廣泛用于天氣、氣候特征分析（梁宏等, 2006;Fujita et al., 2012）以及與衛(wèi)星反演水汽、探空、再分析等資料的對比評估（Liu et al., 2006; Wang et al., 2007; Zhang et al., 2018; 梁宏等, 2012）。有研究表明，大氣可降水量PWV（Perceptible Water Vapor）資料在短時間內(nèi)的快速增加與降水有密切關(guān) 系（Seco et al., 2012; Shoji, 2013; Yao et al.,2017），另外隨著地基水汽資料的積累，已有研究人員開始利用該資料開展氣候變化研究（van Malderen et al., 2014; Wang et al., 2016）。資 料 同化與數(shù)值預(yù)報是地基GNSS/MET 水汽資料重要業(yè)務(wù)應(yīng)用方向，同化地基GNSS/MET 資料能夠提高模式對水汽相關(guān)要素的預(yù)報效果。在美國，NOAA 自1998 年開始評估該資料在天氣預(yù)報同化/模式系統(tǒng)的效果，測試表明同化PWV 資料對濕度的預(yù)報有一定的正效果（Gutman et al., 2004）。英國氣象局自2007 年開始在其北大西洋和歐洲數(shù)值預(yù)報模式中同化了天頂總延遲（Zenith Total Delay，簡稱ZTD）要素，改進(jìn)了相對濕度和云的預(yù)報（Bennitt and Jupp, 2012）。法國氣象局2006 年6月起在業(yè)務(wù)系統(tǒng)Arpege 中同化ZTD 資料，試驗表明同化ZTD 能夠改進(jìn)天氣尺度環(huán)流和降水評分（Poli et al. 2007)。仲躋芹等（2017）的研究也表明，同化ZTD 可以有效提升預(yù)報系統(tǒng)的降水預(yù)報效果，特別是在無探空資料參加同化的預(yù)報時次，同時也發(fā)現(xiàn)同化ZTD 的效果優(yōu)于同化PWV 的效果。英國氣象局統(tǒng)計了單位觀測數(shù)據(jù)量在同化中的影響力，結(jié)果發(fā)現(xiàn)GNSS/MET 資料排名第二（Jones et al., 2020）。

有諸多原因會影響地基GNSS 水汽產(chǎn)品的質(zhì)量，衛(wèi)星相關(guān)誤差（如軌道誤差、衛(wèi)星鐘差、衛(wèi)星仰角過低等），信號傳播相關(guān)誤差（如電離層誤差、多路徑效應(yīng)等），接收機(jī)相關(guān)誤差（如接收機(jī)鐘差、設(shè)備故障等）以及觀測過程相關(guān)誤差（如觀測環(huán)境噪聲、障礙物對天線的遮擋等）都可能造成資料出現(xiàn)各種錯誤（Wang et al., 2007; Bock et al., 2016）。另外在原始資料解算過程中需要用到投影函數(shù)和地面氣壓、氣溫等變量，若數(shù)據(jù)傳輸過程中出現(xiàn)要素不全或數(shù)據(jù)文件打包壓縮過程出現(xiàn)失誤或數(shù)據(jù)上傳不及時等問題也會影響數(shù)據(jù)質(zhì)量。因此在進(jìn)行應(yīng)用之前對資料進(jìn)行質(zhì)量控制是必不可少的步驟。Wang et al.（2007）利用PWV 與探空和微波輻射計資料對比之前采用了離群值檢查，選取4 倍標(biāo)準(zhǔn)差作為閾值，剔除了不到0.1%的數(shù)據(jù)。英國氣象局的資料進(jìn)入同化前剔除了接收機(jī)高度與背景場模式地形相差300 m 以上或與模式背景場相差超過55 mm 的ZTD 資料（Bennitt and Jupp, 2012）。法國氣象局同化ZTD 資料時采用了初猜場質(zhì)量控制檢查，剔除了大約不到3%的數(shù)據(jù)（Poli et al.,2007）。李昊睿等（2014）在同化前對PWV 進(jìn)行了界限值檢查并剔除了與背景場相差較大的數(shù)據(jù)（絕對差值超過6 mm）。仲躋芹等（2017）也基于模式背景差設(shè)計了由多個檢查步驟組成的質(zhì)控算法，各月份未通過質(zhì)控的數(shù)據(jù)比例約為2%～8%。

總體而言，現(xiàn)有的研究在開展天氣與氣候分析之前以基本的離群值檢查（即剔除偏離均值3 或4倍標(biāo)準(zhǔn)差的粗大誤差）為主，較少進(jìn)行嚴(yán)格的質(zhì)量控制。同化前的質(zhì)量控制算法更細(xì)致，但是被剔除的數(shù)據(jù)與模式密切相關(guān)，例如剔除與模式地形相差較大的臺站數(shù)據(jù)（Bennitt and Jupp, 2012）等。此外相對于數(shù)據(jù)本身的質(zhì)量而言，同化前的質(zhì)量控制主要關(guān)注與背景場的偏差。Bock et al.（2016）建立了一套獨(dú)立于模式的質(zhì)量控制方案，但是其中多個檢查步驟用到了G?SP 軟件（Zumberge et al.,1997）相關(guān)的參數(shù)，對其他軟件（例如GAM?T；Herring et al., 2010）并不完全適用。

氣象資料的質(zhì)量控制算法基本都帶有統(tǒng)計特征，比較常用的方法是將觀測值與某個參考值進(jìn)行比較，超出特定的閾值范圍則認(rèn)為該觀測值是疑誤數(shù)據(jù)（WMO, 1993），參考值可以是鄰近站或相鄰時間點(diǎn)上的值或另一個觀測平臺或數(shù)值模式結(jié)果（Bennitt and Jupp, 2012）。閾值的選取方法一般包括統(tǒng)計信度、偏離平均值的程度或預(yù)設(shè)標(biāo)記比例等（Durre et al., 2008）。作為統(tǒng)計結(jié)果，幾乎所有算法都面臨第一類錯誤（棄真）和第二類錯誤（存?zhèn)危┑娘L(fēng)險，所以需要充分評估兩類錯誤后綜合確定閾值。本文參考了Graybeal et al.（2004）提出的預(yù)設(shè)標(biāo)記比例的思路，即每一個檢查都對預(yù)設(shè)比例的數(shù)據(jù)進(jìn)行標(biāo)記（例如5%或1%）?？紤]到這種方法必然會導(dǎo)致棄真的問題（Durre et al., 2008），本文引入綜合決策算法，定位被多項檢查反復(fù)標(biāo)記的數(shù)據(jù)，并最終判斷這些數(shù)據(jù)是否正確、可疑或錯誤，以達(dá)到降低誤判、提高判斷準(zhǔn)確率的效果。

本文第2 部分介紹本研究用的觀測和再分析數(shù)據(jù)，第3 部分介紹質(zhì)量控制方法及評估結(jié)果，第4部分利用質(zhì)量控制后的數(shù)據(jù)開展幾套再分析資料的對比評估，最后是結(jié)論與討論部分。

2 數(shù)據(jù)

2.1 地基GNSS/MET 水汽資料

本文采用全國1254 站2016～2019 年的小時觀測數(shù)據(jù)作為樣本來統(tǒng)計質(zhì)量控制所用參數(shù)以及評估質(zhì)量控制方案效果。Liang et al.（2015）介紹了原始觀測資料的收集與處理流程。圖1 是臺站分布，這些臺站中245 個來自中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)（CMONOC，藍(lán)點(diǎn)），1019 個來自中國氣象局觀測站網(wǎng)。從圖中可以看到氣象局臺站主要密集分布在我國東部地區(qū)，空間分布不均勻，不同省份之間有較大差異；CMONOC 臺站空間分布相對均勻，我國中西部地區(qū)主要以CMONOC 臺站為主。

圖1 地基全球衛(wèi)星導(dǎo)航系統(tǒng)氣象GNSS/MET 臺站分布，紅點(diǎn)代表氣象局（CMA）觀測站點(diǎn)，藍(lán)點(diǎn)代表中國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)CMONOC 觀測站點(diǎn)Fig. 1 Distribution of GNSS/MET (Ground-based Navigation Satellite System/METeorology) sites. The red dots denote the CMA (China Meteorological Administration) sites, whereas the blue dots denote the CMONOC (Crustal Movement Observation Network Of China) sites

2.2 再分析資料

本文用PWV 觀測資料與中國第一代全球大氣再分析產(chǎn)品CRA（China’s first-generation global atmosphere reanalysis，http://idata.cma/idata/web/fact/toTechReport2 [2021-10-08]）以及另外四套再分析資料提供的PWV 要素或整層積分水汽含量TCW（Total Column Water）進(jìn)行了對比評估：ERA?nterim（Dee et al., 2011）、ERA5（Hersbach et al.,2020）、JRA55（Kobayashi et al., 2015）和NCEPDOE AM?P-??（Kanamitsu et al., 2002）。上述幾類再分析產(chǎn)品均未同化GNSS/MET 水汽資料。評估時取液態(tài)水的密度為1×103kg m-3，將TCW 換算成PWV，采用雙線性插值的方法將不同分辨率的再分析格點(diǎn)資料插值到最近的GNSS/MET 站點(diǎn)上，再用公式（1）和公式（2）計算觀測資料與再分析資料的偏差Bias 和均方根誤差RMSE。

3 綜合質(zhì)量控制（CQC）方法與效果評估

CQC 方法包括質(zhì)量檢查環(huán)節(jié)和綜合決策環(huán)節(jié)兩部分（圖2），其中質(zhì)量檢查部分包括界限值檢查、臨近點(diǎn)檢查、濾波檢查等7 個檢查模塊。綜合決策環(huán)節(jié)包括權(quán)重評分系統(tǒng)和最終判斷算法。

圖2 GNSS/MET 水汽產(chǎn)品質(zhì)量控制流程Fig. 2 Quality control flow of the GNSS/MET data

3.1 質(zhì)量檢查

質(zhì)量檢查包括界限值檢查、考察時間一致性的臨近點(diǎn)檢查和低通濾波檢查，考察空間一致性的鄰近站檢查、距平值檢查和峰谷值檢查，以及針對同化應(yīng)用的基于背景場的粗大誤差檢查。本文中“閾值”指的是某個觀測值合理的取值范圍的上下限，而“參數(shù)”則是與閾值的選取相關(guān)的統(tǒng)計量。CQC 中除界限值和峰谷值檢查以外的幾項檢查采用的是預(yù)設(shè)標(biāo)記比例的方法（Graybeal et al.,2004），這里我們對10%、5%、1%和0.1%四個不同標(biāo)記比例（以下稱為PS10、PS5、PS1 和PS0.1，PS 為Parameter Scheme 的簡稱）進(jìn)行評估后，最終選擇PS1 作為質(zhì)量控制參數(shù)。

3.1.1 界限值檢查

界限值檢查通常是質(zhì)量控制算法的第一步，目的是檢查數(shù)據(jù)是否在該要素觀測值允許的物理范圍內(nèi)（儀器、邏輯等）。ZTD 的界限值范圍是[1000.0 mm, 3000.0 mm]，PWV 的界限值范圍是[0 mm, 100 mm]（Jones et al., 2020）。

3.1.2 臨近點(diǎn)檢查

作為小時值數(shù)據(jù)，相鄰兩個時次數(shù)據(jù)之間的變化量應(yīng)該在合理的范圍內(nèi)。為了確定這個范圍，我們計算了所有臺站的時間序列上每個非缺測點(diǎn)與其前一個非缺測點(diǎn)的差值絕對值 γti=|Vi-Vi-1|的概率分布函數(shù)（i=1, 2, 3,···,N，N為觀測數(shù)據(jù)量），以ZTD 為例，90%、95%、99%和99.9%的函數(shù)值對應(yīng)的γti值分別是14 mm、18 mm、27.5 mm和37 mm（表1），即γti值超過上述參數(shù)的比例分別為總數(shù)據(jù)的10%、5%、1%和0.1%。當(dāng)待檢數(shù)據(jù)與前后兩個時次數(shù)據(jù)的絕對偏差γti均超過表1中PS1 對應(yīng)的參數(shù)時，視為不通過臨近點(diǎn)檢查。

表1 不同標(biāo)記比例下的參數(shù)Table 1 Parameters of different flag rate

3.1.3 低通濾波檢查

低通濾波檢查與臨近點(diǎn)檢查都是針對單個臺站，目的是找出時間序列上的離群點(diǎn)。針對所有臺站的時間序列，選取7 小時時間窗，對觀測序列進(jìn)行加權(quán)滑動平均[見公式（3）和公式（4）]得到其低通濾波值Fi，計算每個非缺測點(diǎn)Vi與Fi的差值絕對值γfi=|Vi-Fi|的概率分布函數(shù)，確定不同概率分布函數(shù)值所對應(yīng)的參數(shù)PF（表1）。當(dāng)待檢查數(shù)據(jù)與其濾波值的差值絕對值γfi超過表1 中PS1 對應(yīng)的參數(shù)時，視為未通過該步檢查。

式中，F(xiàn)i是i點(diǎn)的加權(quán)滑動平均值（i=1, 2, 3,···,N，N為觀測數(shù)據(jù)量），hj是權(quán)重系數(shù)，d是數(shù)據(jù)與滑動窗口內(nèi)均值的絕對偏差，n是滑動窗口內(nèi)非缺測數(shù)據(jù)量。

3.1.4 鄰近站檢查

地基水汽資料有較好的空間一致性，在地形相差不大的情況下，某個站的某次觀測值與其周圍鄰近參考站的差值應(yīng)該保持在合理范圍內(nèi)。鄰近站檢查的目標(biāo)是找出那些數(shù)值上偏離周圍臺站的觀測。對于某時刻的空間場，計算每個點(diǎn)Vi與其鄰近站平均值Mi的差值絕對值γni=|Vi-Mi|的概率分布函數(shù)，確定不同概率分布函數(shù)值所對應(yīng)的參數(shù)PN（表1）之后，選擇PS1 作為質(zhì)量控制參數(shù),對于某個目標(biāo)點(diǎn)Vi，取為[Mi-PN，Mi+PN]作為Vi的閾值，超出閾值范圍視為未通過該步檢查。

本研究中鄰近參考站選取方法如下：

（1）空間距離方面，目標(biāo)站和待選參考站空間距離不超過200 km，海拔2500 m 以下兩站高度差不超過200 m，海拔2500 m 以上兩站高度差不超過500 m。

（2）相關(guān)性方面，以2016～2019 年共四年數(shù)據(jù)為統(tǒng)計樣本，要求目標(biāo)站和待選參考站相關(guān)系數(shù)超過0.8，同時要求待選參考站樣本量超過6500，并且與目標(biāo)站在同一觀測時間能匹配的樣本量超過2000。有些臺站自身或周圍臺站序列較短或缺測較多則可能無法取得參考站。

（3）空間分布方面，以目標(biāo)站為中心，參考站盡可能位于四個不同象限內(nèi)。對于滿足距離條件和相關(guān)性條件的待選參考站，分別在目標(biāo)站的四個象限內(nèi)按相關(guān)系數(shù)大小排序，目標(biāo)站的第1 至第12 個參考站依次在四個象限內(nèi)挑選，即第1、5 和9 個待選參考站是目標(biāo)站為原點(diǎn)的第一象限中相關(guān)系數(shù)最高的三個臺站，第2、6 和10 個待選參考站是目標(biāo)站為原點(diǎn)的第二象限中相關(guān)系數(shù)最高的三個臺站，以此類推。若某個象限中待選參考站不足則跳過該象限，在下一個象限中進(jìn)行選擇?？偣膊怀^12 個參考站。

按照上述條件，在臺站相對密集的東部地區(qū)，多數(shù)臺站都能在四個不同象限找到鄰近參考站。圖3 是山西汾陽站（站號53679，紅色十字）的12 個參考站（紅點(diǎn)）分布，可以看到鄰近參考站相對均勻的分布在目標(biāo)站周圍。圖中有些空間距離較近的臺站由于時間序列過短或缺測數(shù)據(jù)太多而無法成為目標(biāo)站的參考站。

圖3 山西汾陽站（紅色十字，站號53769）及其周圍站點(diǎn)，其中紅點(diǎn)為所選的汾陽站的鄰近站Fig. 3 Fengyang station of Shanxi Province (red cross, station ?D:53769) and nearby stations. The red dots denote the selected reference stations

3.1.5 距平值檢查

距平值檢查也是針對空間場的檢查，旨在剔除空間距平場中的離群點(diǎn)。首先計算各臺站各月的多年平均值Mv，隨后計算各臺站數(shù)據(jù)Vi相對于Mv的距平ANOi并進(jìn)行升序排列，得到上下四分位值（Q75和Q25，即75%和25%），然后計算四分位間距外的距平值與Q75和Q25的差值 γai（公式6）的概率密度函數(shù)，

其中，i=1, 2, 3,···,N，N為某時刻的觀測數(shù)據(jù)量。這里比較 ANOi與上下四分位而不是中位數(shù)或平均值的差異，是因為不同時刻的 A NOi可能不是正態(tài)分布，例如某時刻有大范圍強(qiáng)降水發(fā)生時，空間場上正距平會明顯多于負(fù)距平。計算選取PS1 時 γa對應(yīng)的參數(shù)PA（表1）之后，對于某個目標(biāo)點(diǎn)，取[Mv-Q25-PA，Q75+PA+Mv]作為閾值。如超出閾值，視為未通過該步檢查。

3.1.6 峰/谷值檢查

對于大量樣本的評估發(fā)現(xiàn)地基GNSS/MET 水汽數(shù)據(jù)的異常值經(jīng)常表現(xiàn)為一個空間場中的孤立極值，正常情況下這些極值可以通過對比該點(diǎn)與時間和空間上的臨/鄰近點(diǎn)進(jìn)行檢查，但是對歷史資料的分析發(fā)現(xiàn)有部分時次、部分臺站數(shù)據(jù)完整性較差，難以找到參考點(diǎn)對目標(biāo)數(shù)據(jù)進(jìn)行判斷，因此本文在CQC 算法中增加了峰/谷值檢查，考察某個時刻的整個空間場中最大和最小的三個點(diǎn)是否超出特定范圍。將某時刻空間場中所有N個非缺測數(shù)據(jù)進(jìn)行降序排列Vi（i=1, 2, 3,···,N，N為觀測數(shù)據(jù)量），評估最大（峰值，i=1, 2, 3）和最?。ü戎?，i=N-2,N-1,N）的三個點(diǎn)偏離其他點(diǎn)的程度。這里我們參考了Houchi et al.（2015）的方法，比較相鄰兩組數(shù)據(jù)的差異，并利用經(jīng)驗試錯的方法對參數(shù)進(jìn)行調(diào)整，最終確定Di，計算公式為對峰值，從V3開始，若D3超過8.0，則將V3作為該時刻的閾值，空間場中超過或等于V3的數(shù)據(jù)會被標(biāo)記出來，若D3未超過8.0，則考察D2。對谷值從VN-2開始，也進(jìn)行類似處理。

3.1.7 基于背景場的粗大誤差檢查

地基水汽資料在同化和預(yù)報中有重要意義，本文也提供了針對背景場的基本質(zhì)量控制，剔除與背景場偏差較大的觀測點(diǎn)。對于某時刻的空間場，計算每個點(diǎn)（O）與背景場（B）的差值絕對值γbi=|Oi-Bi|的概率分布函數(shù)，確定不同概率分布函數(shù)值所對應(yīng)的參數(shù)PB（表1）之后，對于某個目標(biāo)點(diǎn)Vi，取 [Bi-PB，Bi+PB]作為Vi的閾值，對超出閾值的數(shù)據(jù)進(jìn)行標(biāo)記。

3.2 綜合決策算法

綜合決策算法包括兩部分：權(quán)重評分系統(tǒng)和最終判斷算法（Decision Making Algorithm，簡稱DMA）。首先，根據(jù)CQC 方案的設(shè)計思路給出相應(yīng)權(quán)重（表2）。若某數(shù)據(jù)通過了某項質(zhì)量控制檢查，則該檢查權(quán)重為0，若未通過，則權(quán)重為表2所示。權(quán)重的選取遵循以下原則：（1）若某個目標(biāo)觀測數(shù)據(jù)未通過“界限值檢查”和“峰/谷值檢查”，則可以認(rèn)定該數(shù)據(jù)為“錯誤”，這兩個檢查的權(quán)重設(shè)定為1.5～2.0 之間。（2）對于采用預(yù)設(shè)標(biāo)記比例的檢查（即除界限值檢查和峰谷值檢查的其他五項檢查），任意一項檢查的權(quán)重低于1.0，任意兩項檢查權(quán)重相加低于1.5，任意三項檢查權(quán)重相加高于1.5?；谏鲜鰲l件，五項檢查的權(quán)重設(shè)置為0.5～1.0 之間。（3）在給定范圍內(nèi)對各檢查模塊的權(quán)重進(jìn)行微調(diào)，使得任意項檢查權(quán)重值相加的結(jié)果互不相同，以保留質(zhì)量控制過程信息。

表2 各檢查模塊權(quán)重Table 2 Weights of the checks

每個觀測數(shù)據(jù)經(jīng)過質(zhì)量控制后會得到一個質(zhì)控評分（SQ，公式13），為其經(jīng)過的檢查權(quán)重之和。數(shù)據(jù)經(jīng)過質(zhì)量檢查后可能出現(xiàn)的結(jié)果共64 種情況，權(quán)重值的設(shè)計使得不同組合的SQ不會重復(fù)，通過分析SQ能夠反推某個數(shù)據(jù)未通過哪些檢查。例如，某數(shù)據(jù)的SQ是2.87，則表示它未通過濾波檢查、峰/谷值檢查和背景場檢查（0.65+1.53+0.69=2.87）。

為了用戶使用方便同時也考慮當(dāng)前數(shù)據(jù)業(yè)務(wù)質(zhì)量控制碼的相關(guān)規(guī)范，經(jīng)過綜合判斷之后我們用DMA 給數(shù)據(jù)標(biāo)注最終質(zhì)量控制碼，如表3 所示。

表3 最終判斷算法Table 3 Decision making algorithm

3.3 質(zhì)量控制結(jié)果分析

以2019 年數(shù)據(jù)為樣本，選取PS1 為各檢查模塊的預(yù)設(shè)標(biāo)記參數(shù)，經(jīng)過綜合判斷后未通過質(zhì)控的ZTD 數(shù)據(jù)比例為0.125%，PWV 數(shù)據(jù)比例為0.129%。我們同時分析了標(biāo)記比例分別為PS10、PS5 和PS1 時未通過質(zhì)量控制的數(shù)據(jù)比例。當(dāng)標(biāo)記比例為PS0.1 時，未通過質(zhì)控的ZTD 數(shù)據(jù)比例為0.043%，PWV 數(shù)據(jù)比例為0.046%，存在較多漏判數(shù)據(jù)。當(dāng)標(biāo)記比例為PS5 時，未通過質(zhì)控的ZTD數(shù)據(jù)比例為0.836%，PWV 數(shù)據(jù)比例為0.841%。通過PS1 方案但是未通過PS5 方案的數(shù)據(jù)多在閾值邊界，對其正確性的判斷較為困難，可能存在誤判的風(fēng)險。以山西永和站為例，圖4 中橫坐標(biāo)0 時刻是2019 年9 月9 日12:00（協(xié)調(diào)世界時，下同）的觀測，該點(diǎn)未通過濾波檢查和鄰近站檢查?？梢钥吹皆擖c(diǎn)確實是一個極大值點(diǎn)，也是前后24 h 內(nèi)的最大值，但是該點(diǎn)未表現(xiàn)出明顯的“離群”特征。在實時數(shù)據(jù)處理業(yè)務(wù)或者在研制基礎(chǔ)數(shù)據(jù)集時，由于數(shù)據(jù)用途廣泛，應(yīng)該采取較為“保守”的質(zhì)量控制算法以盡可能保留更多數(shù)據(jù)，因此本文選擇了PS1 作為標(biāo)記的參數(shù)。

圖4 山西永和站（站號53852）2019 年9 月9 日12:00（圖中橫坐標(biāo)0 時刻）及前后24 h 大氣可降水量PWV 變化（黑色實線）及相關(guān)QC 閾值（PS5 參數(shù)。粉線：臨近點(diǎn)檢查，藍(lán)線：濾波檢查，橙線：鄰近站檢查，紫線：距平值檢查，綠點(diǎn)：背景場檢查，深紅線：峰/谷值檢查）。灰色區(qū)域是能夠PS5 參數(shù)條件下通過CQC 的數(shù)據(jù)取值范圍Fig. 4 GNSS/MET precipitable water vapor PWV data (black dotted line) 24 h before and after 1200 UTC on September 9, 2019, from Yonghe station in Shanxi Province (station ?D 53852) and the associated QC thresholds of buddy check (pink line), low-pass filter check (blue line),neighboring station check (orange line), anomaly check (purple line), background check (green dots), peak–valley value check (carmine line) (PS5 parameters; the gray shaded area denotes the range of the data that can pass CQC under the PS5 parameter conditions)

圖5 給出PS1 條件下，ZTD 和PWV 不同錯誤類型所對應(yīng)的數(shù)據(jù)比例。從圖中可以看到，未通過質(zhì)控的ZTD 數(shù)據(jù)中最多的三類組合依次是T+F（1.205）、A+F（1.24）和T+A（1.145），而PWV數(shù)據(jù)中最多的是T+F（1.205）、A+N（1.19）和F+N（1.25）（字母含義見表3）。對于ZTD，接近一半的疑誤數(shù)據(jù)是由時間一致性的兩項檢查（T+F）標(biāo)記出來，而對于PWV，兩項空間一致性檢查的權(quán)重相對于ZTD 有所增加，這可能是由于從ZTD 解算PWV 的過程中用到了地面氣壓、氣溫和高度等參數(shù)，可能將局地因素引入PWV，比起ZTD，某個點(diǎn)的PWV 值偏離周圍臺站的可能性增加。

圖5 不同錯誤類型所占比例Fig. 5 Proportions of different error types

圖6 是多個錯誤類型的個例，可以看到質(zhì)量控制算法能夠通過不同檢查的組合剔除離群點(diǎn)（圖中0 時刻的點(diǎn)）。對質(zhì)控結(jié)果的分析發(fā)現(xiàn)，未通過兩項質(zhì)控而被檢出的疑誤觀測數(shù)據(jù)通常超出閾值（圖中陰影）但與閾值較為接近（圖6b、c、d），我們通過對大量個例的人工分析認(rèn)為在PS1 條件下將檢出的數(shù)據(jù)判為疑誤是較為合理的。

圖6 同圖4，但為（a）湖北宜昌（站號57461）2019 年8 月27 日12:00，（b）河北灤平（站號54420）2019 年7 月13 日12:00，（c）重慶巫山（站號57349）2019 年7 月3 日05:00，（d）江蘇徐州（站號58027）2019 年7 月15 日06:00 PWV 觀測和PS1 條件下的閾值范圍Fig. 6 Same as Fig.4, but for (a) Yichang station in Hubei Province (station ?D 57461), (b) Luanping station in Hebei Province (station ?D 54420), (c)Wushan station in Chongqing Province (station ?D 57349), and (d) Xuzhou station in Jiangsu Province (station ?D 58027). The gray shaded area denotes the range of the data that can pass CQC under the PS1 parameter conditions

分析發(fā)現(xiàn)，觀測數(shù)據(jù)中的離群點(diǎn)經(jīng)常以異常大值的形式出現(xiàn)。為了驗證質(zhì)控算法的有效性，我們參考Houchi et al.（2015）采用的百分位方法對質(zhì)量控制前后的觀測數(shù)據(jù)進(jìn)行了分析。將每個時刻所有臺站的觀測值按照從小到大的順序排列，利用公式（14）計算百分位值。圖7 是ZTD 2018 年質(zhì)量控制前后各時刻的空間場中多個百分位的變化曲線，這幾條最外層的百分位曲線基本能夠展現(xiàn)大值的分布，其中圖7a 中“100%”代表了各時刻觀測場中的最大值?？梢钥吹劫|(zhì)量控制后絕大部分異常大值被剔除，但是圖7a 中仍有極個別異常大值由于觀測數(shù)據(jù)量太少、CQC 算法難以發(fā)揮作用而無法剔除。

圖7 2018 年全國臺站ZTD 質(zhì)量控制前（藍(lán)線）后（紅線）不同百分位值的分布Fig. 7 Time series of the ZTD percentiles of all sites before (blue line) and after (red line) quality control in 2018

式中，i=1, 2, 3,···,n，n為某個時刻觀測樣本數(shù)。

在缺乏“真值”的情況下，利用數(shù)值模式結(jié)果對觀測資料進(jìn)行評估是驗證資料質(zhì)量的重要手段（ 王丹等, 2020）。目前全球多套再分析資料均未同化地基GNSS/MET 水汽資料，因此觀測與再分析資料作為兩個相對獨(dú)立的數(shù)據(jù)源可以充分利用各自優(yōu)勢進(jìn)行對比分析。本文采用CRA 作為背景場來評估質(zhì)量控制的效果。為了更客觀地進(jìn)行對比，這里的質(zhì)量控制結(jié)果關(guān)閉了背景場檢查。圖8 給出的是PS1 條件下觀測與背景場對比的效果，從圖8a可以看到被剔除的點(diǎn)觀測與背景場的值均存在較大差異，質(zhì)控前后觀測與背景場的相關(guān)系數(shù)從0.940提高到0.964。圖8b 則是2018 年未通過CQC 的觀測PWV 與背景場的對比，圖中黑色虛線是二者的擬合曲線（公式16），可以看到其較為明顯的偏離了對角線。觀測值與背景場的相關(guān)系數(shù)為0.742，明顯低于圖8a 所示的正常情況。

圖8 （a）2019 年4 月1～7 日18:00 質(zhì)量控制前（紅點(diǎn)加藍(lán)點(diǎn)）、后（藍(lán)點(diǎn)）PWV 觀測值（Obs）與背景場（CRA）對比的散點(diǎn)圖；（b）2018 年未通過綜合質(zhì)量控制的觀測場（Obs）與背景場（CRA）的散點(diǎn)圖。圖中黑色虛線是觀測與背景場的線性擬合結(jié)果，灰色虛線是若觀測值與背景場相等時線性擬合結(jié)果。Fig. 8 Comparison of PWV between observation and background field (CRA): (a) The red dots denote the error data detected by the QC algorithm at 1800 UTC from April 1 to 7, 2019; (b) blue dots denote the error data detected by the QC algorithm in 2018 (the black dashed line denotes the linear fitting, the gray dashed line denotes the linear fitting when the observations equal to the background data)

圖9 是PWV 觀測與背景場對比的數(shù)據(jù)量和均方根誤差（RMSE）的時間序列（為保證資料穩(wěn)定性，這里只挑選完整性超過40%的臺站與再分析資料進(jìn)行對比，另外為了更獨(dú)立評估算法效果，圖9 中的質(zhì)量控制關(guān)閉了背景場檢查）。結(jié)合表4可以看到，質(zhì)量控制在一定程度上降低了觀測與再分析的均方根誤差，在2018 年1～5 月以及2019年4、5 月份原始RMSE 波動比較明顯時，質(zhì)量控制效果也格外顯著，表明本文提出的方案能夠有效穩(wěn)定資料質(zhì)量。

圖9 2018～2019 年每日00:00 質(zhì)量控制前（黑線）后（紅線）用于與背景場進(jìn)行比較的（a）PWV 的數(shù)據(jù)量和（b）觀測與背景場的均方根誤差（RMSE）時間序列Fig. 9 Data quality of (a) CRA for PWV and (b) RMSE of observation and background fields at 0000 UTC everyday of 2018–2019：Before (black)and after (red) quality control

表4 2018～2019 年每日00:00 質(zhì)量控制前后PWV 觀測與CRA 的均方根誤差（RMSE）Table 4 RMSE of observation and CRA for PWV before and after quality control at 0000 UTC everyday of 2018–2019

4 與再分析資料的對比評估

本文用2018～2019 年經(jīng)過質(zhì)量控制的PWV觀測與多套再分析資料進(jìn)行了對比。為了更客觀比較不同再分析資料，這里的質(zhì)量控制關(guān)閉了背景場檢查。圖10 是觀測與背景場的平均偏差（Bias）和均方根誤差（RMSE）的時間序列（進(jìn)行了5 點(diǎn)平滑）。從整體來看，春季和秋冬季二者偏差（O-B）基本在零線上下波動（圖10a），其中冬季略為負(fù)偏差，春秋季略為正偏差。在水汽更充沛的夏季出現(xiàn)較為明顯的正偏差，表明再分析資料中夏半年整層大氣含水量低于觀測。CRA、ERA?nterim 和ERA5 三套再分析資料與觀測結(jié)果更為一致，優(yōu)于JRA55 和NCEP2 再分析結(jié)果，特別在夏半年。表5 給出觀測與CRA、ERA-?nterim、和ERA5 三個再分析資料的平均偏差和均方根誤差，可以看到CRA 的全國平均偏差為0.633 mm，均方根誤差為3.650 mm，介于ERA-?nterim 和ERA5 之間。

圖10 2018～2019 年00:00 GNSS/MET PWV 觀測與幾套再分析（CRA（紅線）、ERA-?nterim（藍(lán)線）、ERA5（黑線）、JRA55（黃線）和NCEP2（綠線））數(shù)據(jù)的偏差（Bias, a）和均方根誤差（RMSE, b）時間序列Fig. 10 Time series of (a) Bias and (b) RMSE of PWV observation and reanalysis: CRA (red), ERA-?nterim (blue), ERA5 (black), JRA55 (yellow)and NCEP2 (green) at 0000 UTC of 2018–2019

表5 2018 年GNSS/MET PWV 00:00 觀測資料與CRA、ERA-Interim 和ERA5 再分析資料的平均偏差（Bias）和均方根誤差（RMSE）Table 5 Averaged Bias and RMSE of PWV observation and reanalysis of CRA, ERA-Interim and ERA5 for the year 2018 at 0000 UTC

圖11 是冬季（1～3 月）PWV 觀測與CRA 和ERA5 再分析資料對比的空間分布，可以看到冬季全國大部分地區(qū)再分析資料相對于觀測以偏濕為主，四川和華南幾個省份略有偏干（圖11a 和c）。整體而言CRA 和ERA5 的Bias 分布較為一致（圖11e），CRA 在西北地區(qū)有少量正偏差，表明冬季CRA 在西北地區(qū)的含水量比ERA5 略低一些。有個別臺站觀測與再分析的差異較為顯著（主要位于江西地區(qū)），提示該地區(qū)可能存在數(shù)據(jù)觀測、處理或傳輸過程的錯誤，原因有待進(jìn)一步分析。另外從RMSE的分布（圖11b 和d）可以看到各省之間存在較為明顯的差異，表明不同省份之間觀測資料質(zhì)量也有差別。圖12 是夏季（7～9 月）觀測與CRA 和ERA5 再分析資料的對比（O-B），可以看到夏季CRA 和ERA5 再分析在太行山脈一線、華南、青藏高原以東地區(qū)均存在偏干（圖12a 和c），CRA在西北地區(qū)偏干比ERA5 更明顯（圖12e）。兩套再分析RMSE 的分布基本一致（圖12b、d 和f），CRA 在西北地區(qū)RMSE 更大一些；此外不同省份之間RMSE 的差異在夏季表現(xiàn)的更加明顯，可能與各省觀測所用儀器以及數(shù)據(jù)處理方式有關(guān)。綜上所述，CRA 和ERA5 兩套再分析與觀測的偏差空間分布基本一致，都在水汽較多的時段（夏半年）和空間范圍內(nèi)（南方地區(qū)）整層大氣含水量略低于觀測，對干旱少雨中國西部地區(qū)模擬的水汽含量也略低，其中CRA 在西北地區(qū)的PWV 比ERA5 更低。

圖11 2018 年1～3 月每日00:00 GNSS/MET PWV 觀測數(shù)據(jù)分別與CRA、ERA5 再分析數(shù)據(jù)的對比分析：（a）PWV 與CRA 的偏差（BiasC）；（b）PWV 與ERA5 的 偏 差（BiasE）；（c）PWV 與CRA 的 均 方 根 誤 差（RMSEC）；（d）PWV 與ERA5 的 均 方 根 誤 差（RMSEE）；（e）BiasC 與BiasE 的差值；（f）RMSEC 與RMSEE 的差值Fig. 11 Comparison of observation and reanalysis data: (a) Bias between PWV observation and CRA (BiasC); (b) Bias of PWV observation and ERA5 (BiasE); (c) RMSE of PWV observation and CRA (RMSEC); (d) RMSE between PWV observation and ERA5 (RMSEE); (e) difference between BiasC and BiasE; (f) difference between RMSEC and RMSEE for January to March at 0000 UTC everyday of 2018

圖12 同圖11，但為2018 年7～9 月每日00:00 數(shù)據(jù)Fig. 12 As in Fig.11, but for July, August and September at 0000 UTC everyday of 2018

5 結(jié)論與討論

地基GNSS/MET 水汽資料有高時效、高時空分辨率的等優(yōu)點(diǎn)，在資料同化、降水特別是強(qiáng)降水預(yù)報等方面有重要作用，但是多種原因都可能造成資料出現(xiàn)各種錯誤，在資料投入科研業(yè)務(wù)應(yīng)用前將錯誤數(shù)據(jù)剔除是必不可少的步驟。本文介紹了針對GNSS/MET 地基水汽資料的綜合質(zhì)量控制（CQC）算法。CQC 算法由質(zhì)量檢查環(huán)節(jié)和綜合決策環(huán)節(jié)兩部分組成，質(zhì)量檢查環(huán)節(jié)首先用界限值檢查剔除超出物理允許值范圍的數(shù)據(jù)，隨后針對時間一致性進(jìn)行臨近點(diǎn)檢查和低通濾波檢查，針對空間一致性進(jìn)行鄰近站檢查、距平值檢查和峰谷值檢查，以及針對同化應(yīng)用的基于背景場的粗大誤差檢查。

本文采用預(yù)設(shè)標(biāo)記比例的方法獲得了質(zhì)量控制過程中所需要的大量參數(shù)，但是如前文所述，該方法必然會造成第一類錯誤（“棄真”），因此引入綜合決策算法來定位被各檢查模塊反復(fù)標(biāo)記的數(shù)據(jù)以減少誤判，并最終判斷是否正確、可疑或錯誤。本文選用PS1 方案統(tǒng)計確定了各個質(zhì)量控制檢查過程所需的參數(shù)和閾值。評估表明，在假設(shè)每項檢查未通過檢查的比例約為1%的前提下，對各檢查結(jié)果進(jìn)行綜合判識后，最終未通過質(zhì)量控制的ZTD數(shù)據(jù)比例為0.125%，PWV 數(shù)據(jù)比例為0.129%。本文通過個例分析、百分位分析以及與再分析資料的對比評估等三個方面證明質(zhì)量控制算法較為合理，能有效剔除粗大誤差。

隨近些年觀測手段的不斷發(fā)展，新型觀測資料越來越多，但是新型資料的質(zhì)量控制方法、特別是質(zhì)量控制所需參數(shù)往往比較匱乏，本文提出的“多個檢查模塊+預(yù)設(shè)標(biāo)記比例+綜合決策”的質(zhì)量控制思路可以供讀者參考借鑒，開展對這些新型資料的清洗。

在實時應(yīng)用時，按照Liang et al.（2015）所介紹的資料處理流程，目前我國GNSS/MET 資料業(yè)務(wù)上采用組網(wǎng)處理的方法，各臺站觀測文件基本上傳完成后再開展解算ZTD 和PWV 的過程，因此在質(zhì)量控制時空間一致性相關(guān)的檢查（鄰近站檢查、距平值檢查、峰/谷值檢查）都可以正常進(jìn)行。時間一致性相關(guān)的兩個檢查中，濾波檢查需要用到當(dāng)前時刻之后三小時的觀測值，在實時業(yè)務(wù)中無法開展；臨近點(diǎn)檢查可在本文3.1 節(jié)的基礎(chǔ)上略做調(diào)整，改為考察當(dāng)前時刻與前一時刻的差值是否超出閾值。此外，在實際業(yè)務(wù)應(yīng)用時，也可以根據(jù)需要選取表1 中所示的更嚴(yán)格的閾值。

利用質(zhì)量控制后的地基水汽觀測數(shù)據(jù)與CRA、ERA-?nterim 和ERA5 等五套再分析資料進(jìn)行了對比評估。評估表明，幾套再分析資料整層大氣含水量在冬季整體略高于觀測，夏季則明顯低于觀測；空間上在中國南方地區(qū)和西部地區(qū)模擬的水汽含量低于觀測，這種情況在夏半年更加明顯。相對于觀測，CRA 的平均偏差（O-B）為0.633 mm，均方根誤差為3.650 mm，介于ERA-?nterim 和ERA5 之間，明顯優(yōu)于JRA55 和NCEP2 再分析結(jié)果。本文的評估采用的是各套再分析提供的PWV 或TCW（整層水汽含量）要素，不同模式地形差異、各高度層水汽含量差異等因素造成的影響都包含在其中，造成圖10 所示差別的原因還有待于進(jìn)一步分析。