中國地基GNSS/MET 水汽產品質量控制及與再分析產品的對比評估

遠芳 廖捷 周自江

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

1 引言

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

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

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

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

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

2 數據

2.1 地基GNSS/MET 水汽資料

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

圖1 地基全球衛星導航系統氣象GNSS/MET 臺站分布，紅點代表氣象局（CMA）觀測站點，藍點代表中國大陸構造環境監測網絡CMONOC 觀測站點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 觀測資料與中國第一代全球大氣再分析產品CRA（China’s first-generation global atmosphere reanalysis，http://idata.cma/idata/web/fact/toTechReport2 [2021-10-08]）以及另外四套再分析資料提供的PWV 要素或整層積分水汽含量TCW（Total Column Water）進行了對比評估：ERA?nterim（Dee et al., 2011）、ERA5（Hersbach et al.,2020）、JRA55（Kobayashi et al., 2015）和NCEPDOE AM?P-??（Kanamitsu et al., 2002）。上述幾類再分析產品均未同化GNSS/MET 水汽資料。評估時取液態水的密度為1×103kg m-3，將TCW 換算成PWV，采用雙線性插值的方法將不同分辨率的再分析格點資料插值到最近的GNSS/MET 站點上，再用公式（1）和公式（2）計算觀測資料與再分析資料的偏差Bias 和均方根誤差RMSE。

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

CQC 方法包括質量檢查環節和綜合決策環節兩部分（圖2），其中質量檢查部分包括界限值檢查、臨近點檢查、濾波檢查等7 個檢查模塊。綜合決策環節包括權重評分系統和最終判斷算法。

圖2 GNSS/MET 水汽產品質量控制流程Fig. 2 Quality control flow of the GNSS/MET data

3.1 質量檢查

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

3.1.1 界限值檢查

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

3.1.2 臨近點檢查

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

表1 不同標記比例下的參數Table 1 Parameters of different flag rate

3.1.3 低通濾波檢查

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

式中，Fi是i點的加權滑動平均值（i=1, 2, 3,···,N，N為觀測數據量），hj是權重系數，d是數據與滑動窗口內均值的絕對偏差，n是滑動窗口內非缺測數據量。

3.1.4 鄰近站檢查

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

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

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

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

（3）空間分布方面，以目標站為中心，參考站盡可能位于四個不同象限內。對于滿足距離條件和相關性條件的待選參考站，分別在目標站的四個象限內按相關系數大小排序，目標站的第1 至第12 個參考站依次在四個象限內挑選，即第1、5 和9 個待選參考站是目標站為原點的第一象限中相關系數最高的三個臺站，第2、6 和10 個待選參考站是目標站為原點的第二象限中相關系數最高的三個臺站，以此類推。若某個象限中待選參考站不足則跳過該象限，在下一個象限中進行選擇。總共不超過12 個參考站。

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

圖3 山西汾陽站（紅色十字，站號53769）及其周圍站點，其中紅點為所選的汾陽站的鄰近站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 距平值檢查

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

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

3.1.6 峰/谷值檢查

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

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

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

3.2 綜合決策算法

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

表2 各檢查模塊權重Table 2 Weights of the checks

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

為了用戶使用方便同時也考慮當前數據業務質量控制碼的相關規范，經過綜合判斷之后我們用DMA 給數據標注最終質量控制碼，如表3 所示。

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

3.3 質量控制結果分析

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

圖4 山西永和站（站號53852）2019 年9 月9 日12:00（圖中橫坐標0 時刻）及前后24 h 大氣可降水量PWV 變化（黑色實線）及相關QC 閾值（PS5 參數。粉線：臨近點檢查，藍線：濾波檢查，橙線：鄰近站檢查，紫線：距平值檢查，綠點：背景場檢查，深紅線：峰/谷值檢查）。灰色區域是能夠PS5 參數條件下通過CQC 的數據取值范圍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 不同錯誤類型所對應的數據比例。從圖中可以看到，未通過質控的ZTD 數據中最多的三類組合依次是T+F（1.205）、A+F（1.24）和T+A（1.145），而PWV數據中最多的是T+F（1.205）、A+N（1.19）和F+N（1.25）（字母含義見表3）。對于ZTD，接近一半的疑誤數據是由時間一致性的兩項檢查（T+F）標記出來，而對于PWV，兩項空間一致性檢查的權重相對于ZTD 有所增加，這可能是由于從ZTD 解算PWV 的過程中用到了地面氣壓、氣溫和高度等參數，可能將局地因素引入PWV，比起ZTD，某個點的PWV 值偏離周圍臺站的可能性增加。

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

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

圖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

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

圖7 2018 年全國臺站ZTD 質量控制前（藍線）后（紅線）不同百分位值的分布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為某個時刻觀測樣本數。

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

圖8 （a）2019 年4 月1～7 日18:00 質量控制前（紅點加藍點）、后（藍點）PWV 觀測值（Obs）與背景場（CRA）對比的散點圖；（b）2018 年未通過綜合質量控制的觀測場（Obs）與背景場（CRA）的散點圖。圖中黑色虛線是觀測與背景場的線性擬合結果，灰色虛線是若觀測值與背景場相等時線性擬合結果。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 觀測與背景場對比的數據量和均方根誤差（RMSE）的時間序列（為保證資料穩定性，這里只挑選完整性超過40%的臺站與再分析資料進行對比，另外為了更獨立評估算法效果，圖9 中的質量控制關閉了背景場檢查）。結合表4可以看到，質量控制在一定程度上降低了觀測與再分析的均方根誤差，在2018 年1～5 月以及2019年4、5 月份原始RMSE 波動比較明顯時，質量控制效果也格外顯著，表明本文提出的方案能夠有效穩定資料質量。

圖9 2018～2019 年每日00:00 質量控制前（黑線）后（紅線）用于與背景場進行比較的（a）PWV 的數據量和（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 質量控制前后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 年經過質量控制的PWV觀測與多套再分析資料進行了對比。為了更客觀比較不同再分析資料，這里的質量控制關閉了背景場檢查。圖10 是觀測與背景場的平均偏差（Bias）和均方根誤差（RMSE）的時間序列（進行了5 點平滑）。從整體來看，春季和秋冬季二者偏差（O-B）基本在零線上下波動（圖10a），其中冬季略為負偏差，春秋季略為正偏差。在水汽更充沛的夏季出現較為明顯的正偏差，表明再分析資料中夏半年整層大氣含水量低于觀測。CRA、ERA?nterim 和ERA5 三套再分析資料與觀測結果更為一致，優于JRA55 和NCEP2 再分析結果，特別在夏半年。表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（藍線）、ERA5（黑線）、JRA55（黃線）和NCEP2（綠線））數據的偏差（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 再分析資料對比的空間分布，可以看到冬季全國大部分地區再分析資料相對于觀測以偏濕為主，四川和華南幾個省份略有偏干（圖11a 和c）。整體而言CRA 和ERA5 的Bias 分布較為一致（圖11e），CRA 在西北地區有少量正偏差，表明冬季CRA 在西北地區的含水量比ERA5 略低一些。有個別臺站觀測與再分析的差異較為顯著（主要位于江西地區），提示該地區可能存在數據觀測、處理或傳輸過程的錯誤，原因有待進一步分析。另外從RMSE的分布（圖11b 和d）可以看到各省之間存在較為明顯的差異，表明不同省份之間觀測資料質量也有差別。圖12 是夏季（7～9 月）觀測與CRA 和ERA5 再分析資料的對比（O-B），可以看到夏季CRA 和ERA5 再分析在太行山脈一線、華南、青藏高原以東地區均存在偏干（圖12a 和c），CRA在西北地區偏干比ERA5 更明顯（圖12e）。兩套再分析RMSE 的分布基本一致（圖12b、d 和f），CRA 在西北地區RMSE 更大一些；此外不同省份之間RMSE 的差異在夏季表現的更加明顯，可能與各省觀測所用儀器以及數據處理方式有關。綜上所述，CRA 和ERA5 兩套再分析與觀測的偏差空間分布基本一致，都在水汽較多的時段（夏半年）和空間范圍內（南方地區）整層大氣含水量略低于觀測，對干旱少雨中國西部地區模擬的水汽含量也略低，其中CRA 在西北地區的PWV 比ERA5 更低。

圖11 2018 年1～3 月每日00:00 GNSS/MET PWV 觀測數據分別與CRA、ERA5 再分析數據的對比分析：（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 數據Fig. 12 As in Fig.11, but for July, August and September at 0000 UTC everyday of 2018

5 結論與討論

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

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

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

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

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