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

國(guó)家氣象信息中心, 北京 100081

1 引言

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

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

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

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

2 數(shù)據(jù)

2.1 地基GNSS/MET 水汽資料

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

圖1 地基全球衛(wèi)星導(dǎo)航系統(tǒng)氣象GNSS/MET 臺(tái)站分布，紅點(diǎn)代表氣象局（CMA）觀測(cè)站點(diǎn)，藍(lán)點(diǎn)代表中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)CMONOC 觀測(cè)站點(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 觀測(cè)資料與中國(guó)第一代全球大氣再分析產(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)行了對(duì)比評(píng)估：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 水汽資料。評(píng)估時(shí)取液態(tài)水的密度為1×103kg m-3，將TCW 換算成PWV，采用雙線性插值的方法將不同分辨率的再分析格點(diǎn)資料插值到最近的GNSS/MET 站點(diǎn)上，再用公式（1）和公式（2）計(jì)算觀測(cè)資料與再分析資料的偏差Bias 和均方根誤差RMSE。

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

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

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

3.1 質(zhì)量檢查

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

3.1.1 界限值檢查

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

3.1.2 臨近點(diǎn)檢查

作為小時(shí)值數(shù)據(jù)，相鄰兩個(gè)時(shí)次數(shù)據(jù)之間的變化量應(yīng)該在合理的范圍內(nèi)。為了確定這個(gè)范圍，我們計(jì)算了所有臺(tái)站的時(shí)間序列上每個(gè)非缺測(cè)點(diǎn)與其前一個(gè)非缺測(cè)點(diǎn)的差值絕對(duì)值 γti=|Vi-Vi-1|的概率分布函數(shù)（i=1, 2, 3,···,N，N為觀測(cè)數(shù)據(jù)量），以ZTD 為例，90%、95%、99%和99.9%的函數(shù)值對(duì)應(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ù)與前后兩個(gè)時(shí)次數(shù)據(jù)的絕對(duì)偏差γti均超過表1中PS1 對(duì)應(yīng)的參數(shù)時(shí)，視為不通過臨近點(diǎn)檢查。

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

3.1.3 低通濾波檢查

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

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

3.1.4 鄰近站檢查

地基水汽資料有較好的空間一致性，在地形相差不大的情況下，某個(gè)站的某次觀測(cè)值與其周圍鄰近參考站的差值應(yīng)該保持在合理范圍內(nèi)。鄰近站檢查的目標(biāo)是找出那些數(shù)值上偏離周圍臺(tái)站的觀測(cè)。對(duì)于某時(shí)刻的空間場(chǎng)，計(jì)算每個(gè)點(diǎn)Vi與其鄰近站平均值Mi的差值絕對(duì)值γni=|Vi-Mi|的概率分布函數(shù)，確定不同概率分布函數(shù)值所對(duì)應(yīng)的參數(shù)PN（表1）之后，選擇PS1 作為質(zhì)量控制參數(shù),對(duì)于某個(gè)目標(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)計(jì)樣本，要求目標(biāo)站和待選參考站相關(guān)系數(shù)超過0.8，同時(shí)要求待選參考站樣本量超過6500，并且與目標(biāo)站在同一觀測(cè)時(shí)間能匹配的樣本量超過2000。有些臺(tái)站自身或周圍臺(tái)站序列較短或缺測(cè)較多則可能無法取得參考站。

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

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

圖3 山西汾陽(yáng)站（紅色十字，站號(hào)53769）及其周圍站點(diǎn)，其中紅點(diǎn)為所選的汾陽(yáng)站的鄰近站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 距平值檢查

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

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

3.1.6 峰/谷值檢查

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

3.1.7 基于背景場(chǎng)的粗大誤差檢查

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

3.2 綜合決策算法

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

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

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

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

圖4 山西永和站（站號(hào)53852）2019 年9 月9 日12:00（圖中橫坐標(biāo)0 時(shí)刻）及前后24 h 大氣可降水量PWV 變化（黑色實(shí)線）及相關(guān)QC 閾值（PS5 參數(shù)。粉線：臨近點(diǎn)檢查，藍(lán)線：濾波檢查，橙線：鄰近站檢查，紫線：距平值檢查，綠點(diǎn)：背景場(chǎng)檢查，深紅線：峰/谷值檢查）。灰色區(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 不同錯(cuò)誤類型所對(duì)應(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）。對(duì)于ZTD，接近一半的疑誤數(shù)據(jù)是由時(shí)間一致性的兩項(xiàng)檢查（T+F）標(biāo)記出來，而對(duì)于PWV，兩項(xiàng)空間一致性檢查的權(quán)重相對(duì)于ZTD 有所增加，這可能是由于從ZTD 解算PWV 的過程中用到了地面氣壓、氣溫和高度等參數(shù)，可能將局地因素引入PWV，比起ZTD，某個(gè)點(diǎn)的PWV 值偏離周圍臺(tái)站的可能性增加。

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

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

圖6 同圖4，但為（a）湖北宜昌（站號(hào)57461）2019 年8 月27 日12:00，（b）河北灤平（站號(hào)54420）2019 年7 月13 日12:00，（c）重慶巫山（站號(hào)57349）2019 年7 月3 日05:00，（d）江蘇徐州（站號(hào)58027）2019 年7 月15 日06:00 PWV 觀測(cè)和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)，觀測(cè)數(shù)據(jù)中的離群點(diǎn)經(jīng)常以異常大值的形式出現(xiàn)。為了驗(yàn)證質(zhì)控算法的有效性，我們參考Houchi et al.（2015）采用的百分位方法對(duì)質(zhì)量控制前后的觀測(cè)數(shù)據(jù)進(jìn)行了分析。將每個(gè)時(shí)刻所有臺(tái)站的觀測(cè)值按照從小到大的順序排列，利用公式（14）計(jì)算百分位值。圖7 是ZTD 2018 年質(zhì)量控制前后各時(shí)刻的空間場(chǎng)中多個(gè)百分位的變化曲線，這幾條最外層的百分位曲線基本能夠展現(xiàn)大值的分布，其中圖7a 中“100%”代表了各時(shí)刻觀測(cè)場(chǎng)中的最大值?？梢钥吹劫|(zhì)量控制后絕大部分異常大值被剔除，但是圖7a 中仍有極個(gè)別異常大值由于觀測(cè)數(shù)據(jù)量太少、CQC 算法難以發(fā)揮作用而無法剔除。

圖7 2018 年全國(guó)臺(tái)站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為某個(gè)時(shí)刻觀測(cè)樣本數(shù)。

在缺乏“真值”的情況下，利用數(shù)值模式結(jié)果對(duì)觀測(cè)資料進(jìn)行評(píng)估是驗(yàn)證資料質(zhì)量的重要手段（ 王丹等, 2020）。目前全球多套再分析資料均未同化地基GNSS/MET 水汽資料，因此觀測(cè)與再分析資料作為兩個(gè)相對(duì)獨(dú)立的數(shù)據(jù)源可以充分利用各自優(yōu)勢(shì)進(jìn)行對(duì)比分析。本文采用CRA 作為背景場(chǎng)來評(píng)估質(zhì)量控制的效果。為了更客觀地進(jìn)行對(duì)比，這里的質(zhì)量控制結(jié)果關(guān)閉了背景場(chǎng)檢查。圖8 給出的是PS1 條件下觀測(cè)與背景場(chǎng)對(duì)比的效果，從圖8a可以看到被剔除的點(diǎn)觀測(cè)與背景場(chǎng)的值均存在較大差異，質(zhì)控前后觀測(cè)與背景場(chǎng)的相關(guān)系數(shù)從0.940提高到0.964。圖8b 則是2018 年未通過CQC 的觀測(cè)PWV 與背景場(chǎng)的對(duì)比，圖中黑色虛線是二者的擬合曲線（公式16），可以看到其較為明顯的偏離了對(duì)角線。觀測(cè)值與背景場(chǎng)的相關(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 觀測(cè)值（Obs）與背景場(chǎng)（CRA）對(duì)比的散點(diǎn)圖；（b）2018 年未通過綜合質(zhì)量控制的觀測(cè)場(chǎng)（Obs）與背景場(chǎng)（CRA）的散點(diǎn)圖。圖中黑色虛線是觀測(cè)與背景場(chǎng)的線性擬合結(jié)果，灰色虛線是若觀測(cè)值與背景場(chǎng)相等時(shí)線性擬合結(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 觀測(cè)與背景場(chǎng)對(duì)比的數(shù)據(jù)量和均方根誤差（RMSE）的時(shí)間序列（為保證資料穩(wěn)定性，這里只挑選完整性超過40%的臺(tái)站與再分析資料進(jìn)行對(duì)比，另外為了更獨(dú)立評(píng)估算法效果，圖9 中的質(zhì)量控制關(guān)閉了背景場(chǎng)檢查）。結(jié)合表4可以看到，質(zhì)量控制在一定程度上降低了觀測(cè)與再分析的均方根誤差，在2018 年1～5 月以及2019年4、5 月份原始RMSE 波動(dòng)比較明顯時(shí)，質(zhì)量控制效果也格外顯著，表明本文提出的方案能夠有效穩(wěn)定資料質(zhì)量。

圖9 2018～2019 年每日00:00 質(zhì)量控制前（黑線）后（紅線）用于與背景場(chǎng)進(jìn)行比較的（a）PWV 的數(shù)據(jù)量和（b）觀測(cè)與背景場(chǎng)的均方根誤差（RMSE）時(shí)間序列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 觀測(cè)與CRA 的均方根誤差（RMSE）Table 4 RMSE of observation and CRA for PWV before and after quality control at 0000 UTC everyday of 2018–2019

4 與再分析資料的對(duì)比評(píng)估

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

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

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

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

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

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

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