累計降水量質(zhì)控方法設計及應用

摘 要：氣象臺站每小時上傳的觀測資料中包括3h/6h/12h/24h累計降水量，但目前質(zhì)控系統(tǒng)僅對每日8時、20時的12h累計降水量進行質(zhì)控，無法滿足用戶對各類累計降水數(shù)據(jù)的需求。本文基于時空一致性原理設計了一種針對累計降水數(shù)據(jù)的質(zhì)控方法，應用于2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件，完成了數(shù)據(jù)完整性檢查并且根據(jù)場景提出相應疑誤，同時實現(xiàn)了錯誤數(shù)據(jù)的自動更正功能。結果表明，該質(zhì)控方法能有效質(zhì)控地面站累計降水量數(shù)據(jù)，可以作為當前質(zhì)控系統(tǒng)的有力補充。

關鍵詞：降水；累計降水量；質(zhì)量控制方法；自動化

中圖分類號：TP 399" " 文獻標志碼：A

累計降水量作為重要的氣象要素，使用十分廣泛，它是指一段時間內(nèi)某地接收到的雨水總量。不同的應用場景下，用戶對累計降水量的時間步長需求不同。例如氣候研究中，用戶常使用的ERA5、CFSR、CRA等再分析資料主要包括6h累計降水。研究降水的日變化、月變化、年變化、年際變化時，用戶需要更長時間步長的累計降水量。在天氣研究中，用戶常使用的衛(wèi)星、雷達定量降水估計產(chǎn)品主要包括1h/3h/6h/12h/24h累計降水[1]。當對不同的氣象資料進行對比時，首先需要統(tǒng)一時間步長。而地面資料的質(zhì)控主要集中在1h降水和每日8、20時的12h累計降水，對其他時次的累計降水量質(zhì)控不足[2]。用戶在使用時需要自行計算所需時間步長的累計降水，由于算法的不同或用戶對降水數(shù)據(jù)的了解不足，可能會出現(xiàn)錯誤，影響用戶的后續(xù)結果。并且用戶自行計算的數(shù)據(jù)無法得到實時有效的質(zhì)控，也無法與臺站核實數(shù)據(jù)的真實性。因此，提供經(jīng)過質(zhì)控的多種不同時間步長的累計降水量數(shù)據(jù)，對用戶更友好，用戶可以直接選擇需要的累計降水量數(shù)據(jù)，提高研究效率。

1 累計降水量

地面站累計降水量是指一段時間內(nèi)的降水總量，通常由小時降水量累計，僅在設備異常時采用人工觀測代替[3]。目前累計降水的時間步長包括1h/3h/6h/12h/24h，臺站上傳的逐小時文件均包括這些累計降水數(shù)據(jù)。然而，目前省級質(zhì)控程序僅對1h降水量和每日8、20時的12h累計降水進行質(zhì)控，其他時次的累計降水量未得到有效質(zhì)控[4]。本文使用2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件中的降水量作為試驗數(shù)據(jù)，應用于設計的質(zhì)控程序中。

2 累計降水量質(zhì)控邏輯構架

2.1 缺測率統(tǒng)計

統(tǒng)計數(shù)據(jù)的缺測率有利于用戶掌握數(shù)據(jù)集的整體概況，對缺測時次較多的臺站來說，可以及時提醒其檢查設備儀器是否出現(xiàn)故障等，保障地面觀測平穩(wěn)運行。先統(tǒng)計時間步長分別為1h/3h/6h/12h/24h的降水量文件出現(xiàn)時次，再計算被統(tǒng)計的時間段內(nèi)應收時次。缺測率如公式（1）所示。

（1）

式中：Rt為已收時次；Rall為總時次；R為缺測率，用百分數(shù)表示。

2.2 累計降水量質(zhì)控

設計累計降水量的質(zhì)控方法，首先考慮累計降水量的計算方法。累計降水量是指單站在一段時間內(nèi)的降水量。臺站每小時均會上傳該小時內(nèi)的降水量，例如臺站在13：00時上傳的文件中的小時降水量是指12：01-13：00內(nèi)該站的降水量。在13：00時上傳的文件中的3h累計降水量是指11-13時上傳的3個文件中的逐小時累計降水量累計值，即10：01-13：00內(nèi)該站的累計降水量。6h/12h/24h累計降水量類似，分別指7：01-13：00、1：01-13：00、前一日13：01-13：00。因此，要計算累計降水量，需要保證時間步長內(nèi)的數(shù)據(jù)完整。然后計算時間步長內(nèi)的累計降水量，與原文件進行對比，這是通過時間一致性的原理對降水進行質(zhì)控。另外，對比站點及其附近站點的累計降水量，通過空間一致性原理對降水量進行質(zhì)控。質(zhì)控總體分以下步驟。第一步：進行時間一致性檢查，如圖1所示。檢查時間步長內(nèi)的逐小時降水量是否完整。當遍歷各時次文件時，對當前時次的前3h/6h/12h/24h內(nèi)的逐小時降水量數(shù)據(jù)進行檢查，如果有缺失，判斷當前時次是否為定時觀測時次——基準站、基本站人工定時觀測次數(shù)為每日5次（08時、11時、14時、17時、20時），一般站人工定時觀測次數(shù)為每日3次（08時、14時、20時）。如果為定時觀測時次，就通過質(zhì)控系統(tǒng)向臺站自動下發(fā)疑誤，質(zhì)控碼置為1。如果為非定時觀測時次，就直接將當前時次文件內(nèi)相應累計降水量數(shù)據(jù)作缺，質(zhì)控碼置為8[3]，常見質(zhì)控碼及含義見表1。

第二步：計算累計降水量。如果當前時次的前3h/6h/12h/24h內(nèi)的逐小時降水量數(shù)據(jù)完整，計算累計降水量，如公式（2）所示。

（2）

式中：TPn為累計降水量；Pi為時間步長內(nèi)的逐小時降水量，i=1，2，...，n；n為時間步長，n=3，6，12，24。

計算值TPn作為參考值保存并提供給數(shù)據(jù)處理人員，經(jīng)臺站確認反饋后可替換原文件中累計降水量。

第三步：對比計算值TPn與原文件累計值是否一致。原文件累計值通常為臺站端本地計算后上傳，有時僅修改錯誤的小時降水量，而忽略了累計降水量。如果計算值TPn與原文件累計值一致，那么判定原文件中累計降水量質(zhì)控結果為正確，質(zhì)控碼置為0。如果不一致，那么判定原文件累計降水量可疑，質(zhì)控碼置為1，向臺站下發(fā)相應疑誤并提供計算值TPn。

第四步：空間一致性檢查（如圖2所示）。查找被檢站點附近5km的鄰近站，站點距離如公式（3）所示[5]。

?R=Rc?arccos[sinφ1sinφ2+cosφ1cosφ2（λ2-λ1）] （3）

式中：?R為站點間距離；Rc為地球半徑，Rc=6371km；φ1、λ1為被檢站點的緯度、經(jīng)度；φ2、λ2為鄰近站點的緯度、經(jīng)度。

為節(jié)省運算量，首先根據(jù)經(jīng)緯度篩選被檢點附近0.05°×

0.05°范圍內(nèi)的站點，再計算站點間距離，保留附近5km內(nèi)的站點[6]。這里計算的鄰近站點可以保存為常用表，不僅可用于降水量的空間一致性檢驗，而且還可適用于其他多種氣象要素，例如地溫、氣溫、氣壓等的檢驗或插值等。

第五步：查詢5km內(nèi)鄰近站點相應時間步長累計降水量的質(zhì)控碼，質(zhì)控碼為0或4的臨近站中的最大值和最小值可以作為被檢站點的上限與下限值范圍參考，如果被檢站點值在范圍內(nèi)，那么質(zhì)控結果為正確，質(zhì)控碼置為0。如果不在范圍內(nèi)，那么質(zhì)控碼置為1，提出疑誤。如果鄰近站點累計降水量質(zhì)控碼為1、8，就不作為參考站點。

3 基于python的質(zhì)控模塊設計與實現(xiàn)

按照2的邏輯架構，設計以下3個模塊以實現(xiàn)質(zhì)控功能。

3.1 連接數(shù)據(jù)庫

地面站資料儲存在數(shù)據(jù)庫中，需要連接數(shù)據(jù)庫選取所需時段內(nèi)的降水數(shù)據(jù)?；赑ython的cx_Oracle，使用用戶名及密碼連接數(shù)據(jù)庫。執(zhí)行以下SQL（Structured Query Language）語句選取2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件中的1h降水量和3h/6h/12h/24h累計降水量數(shù)據(jù)：

sql=\"select station，time，rain_1h，rain_3h，rain_6h，rain_12h，rain_24h from OBS where timegt;=TO_DATE（'2023/01/01 01：L00：00'，'yyyy/mm/dd hh：00：00'）and timelt;=TO_DATE（'2024/01/01 00：00：00'，'yyyy/mm/dd hh：00：00'）

其中，rain_1h/rain_3h/rain_6h/rain_12h/rain_24h是指原始地面文件的1h/3h/6h/12h/24h累計降水量數(shù)據(jù)。將查詢結果保存為DataFrame，有利于后續(xù)質(zhì)控的處理。

3.2 整體數(shù)據(jù)完整性檢查

檢查所有地面站的所有時次文件中1h/3h/6h/12h/24h累計降水量是否完整。統(tǒng)計缺測時次，得到缺測率。本功能使用datetime庫計算總時次，使用Pandas的數(shù)據(jù)結構DataFrame統(tǒng)計已收時次，其差值為缺測時次，可計算各步長累計降水量缺測率。具體實現(xiàn)如下。

startime=datetime.datetime.strptime（'2022/12/31 13：00：00'，'%Y/%m/%d %H：00：00）

endtime=datetime.datetime.strptime（'2023/12/31 12：00：00'，'%Y/%m/%d %H：00：00）

duration=endtime-startime + datetime.timedelta（hours=1）

duration_hours=duration.days*24 + duration.seconds/3600

R=1 - （len（df）/duration_hours）

其中，df為已收數(shù)據(jù)；len（df）為已收時次；duration_hours為2023年1月1日01時至2024年1月1日00時的總時數(shù)，即應收總時次；R為統(tǒng)計的缺測率，可以用百分數(shù)表示。

3.3 時間一致性檢查

分臺站依時次讀取地面站降水量數(shù)據(jù)，以讀取時次作為被檢驗時次。檢查被檢臺站的被檢驗時次前3h/6h/12h/24h內(nèi)降水數(shù)據(jù)是否完整，本功能使用datetime庫生成時間戳，檢查原文件是否包含該時次過去3h/6h/12h/24h內(nèi)的所有時間戳。具體實現(xiàn)如下。

For i in range（1，hours）：

time_end=time

time_fore=time_end - datetime.timedelta（hours=i）

if （df['time'] == time_fore）==：

break

else：

p = df['time']['rain_1h']

其中hours的值分別為3/6/12/24，df['rain']為原文件中被檢時次的前3h/6h/12h/24h內(nèi)的逐小時降水量。針對缺測的逐小時累計降水量，其累計值也置為缺測保存。檢查該被檢時次是否為定時觀測時次，如果是，那么質(zhì)控碼置為1，向臺站下發(fā)疑誤。如果為非定時觀測時次，那么直接修改原文件該時次累計值為缺測，質(zhì)控碼置為8。

經(jīng)過完整性檢查，符合條件的時次可以進行累計降水量的計算。累計降水量p，得到計算值TP。判斷計算值TP與原文件中相應時間步長的累計降水量是否一致。具體實現(xiàn)如下。

rain_sum=df[‘time’][‘rain_hour’]

if np.allclose（TP，rain_sum，atol=0.1）==：

code=0

else：

code=1

其中，rain_sum為原數(shù)據(jù)中的累計降水量；rain_hour為相應時間步長的累計降水量列名（3.1中提到）；np.allclose用于判斷數(shù)組是否在誤差范圍內(nèi)相等；atol為絕對容忍度，在這里表示TP和rain_sum在0.1的誤差范圍內(nèi)相等；誤差設置為0.1是由于地面降水量觀測值的精度為0.1mm；code為質(zhì)控碼，如果TP與rain_sum不相等，那么質(zhì)控碼為1，向臺站下發(fā)疑誤。如果相等，那么質(zhì)控碼為0，相應累積降水量置為正確。

3.4 空間一致性檢查

經(jīng)過時間一致性檢查后，質(zhì)控碼為1的數(shù)據(jù)下發(fā)臺站質(zhì)疑，質(zhì)控碼為8的數(shù)據(jù)直接作缺測處理。數(shù)據(jù)質(zhì)控碼為0和4的累計降水量數(shù)據(jù)通過了時間一致性檢驗，可進行進一步的空間一致性檢查。首先，挑選被檢點附近5km內(nèi)的站點，為減少計算量，篩選被檢點附近0.05°×0.05°的站點，經(jīng)緯度上1°約為110km，0.05°約為5.5km。其次，計算在此范圍內(nèi)的站點與被檢點的實際距離，使用geopy庫中的great_circle包，其計算站點間的原理如公式（3）所示，具體實現(xiàn)如下。

Apoint=（Alat，Alon）

Bpoint=（Blat，Blon）

distance=great_circle（Apoint，Bpoint）.km

其中，Apoint為被檢點；Bpoint為可能鄰近站。

將鄰近站的距離范圍控制在5km內(nèi)是由于降水具有局地性和陣性的特點，距離過遠可能會導致站點間的降水相關性較小，而缺失其參考價值。針對一些地處甘孜州、阿壩州及涼山州等偏遠地區(qū)的站點，平均站點距離超過10km，其5km內(nèi)可能沒有鄰近站點，這些站點無法進行空間一致性檢查?；蛘呤琴|(zhì)控碼為1、8，本身數(shù)據(jù)可疑的站點也無法作為參考站點。針對符合條件的鄰近站點，選取其相應時間步長的累計降水量值，其中最大值作為上限，最小值作下限，判斷被檢點的值是否在此范圍內(nèi)。如果不在范圍內(nèi)，那么質(zhì)控碼為1，疑誤下發(fā)臺站。如果在范圍內(nèi)，那么質(zhì)控碼為0。

4 應用結果

利用2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件，應用該質(zhì)控程序，缺測率統(tǒng)計結果如圖3所示。國家站整體數(shù)據(jù)質(zhì)量較高，僅部分站點由設備故障等原因出現(xiàn)累計降水量缺測，因此圖3省略了沒有出現(xiàn)缺測的站點。出現(xiàn)缺測的站點缺測率最高不超過0.6%，主要集中在0.25%以下。由于1h降水量和3h/6h/12h/24h累計降水量的內(nèi)在聯(lián)系，因此各時間步長的缺測率分布較為一致。

再執(zhí)行時間一致性檢查，提出的疑誤如下。1）時間步長內(nèi)的逐小時降水正常，但3h/6h/12h/24h累計降水為缺測，針對這類疑誤，可應用質(zhì)控程序的計算值補充缺失的累計降水量。2）質(zhì)控程序計算值與原文件累計降水量不一致，針對這一類疑誤，需要向臺站再次核實時間步長內(nèi)的逐小時降水量是否正確。3）時間步長內(nèi)的逐小時降水有缺測，但原文件有累計值，針對這類疑誤，需要向臺站核實累計值是否為人工觀測，如果不是，就按照質(zhì)控程序缺測處理。對2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件質(zhì)控后，第一類疑誤出現(xiàn)226條，第二類疑誤出現(xiàn)534條，第三類疑誤出現(xiàn)287條，見表2。

對經(jīng)過時間一致性檢查且質(zhì)控結果為0和4的累計降水量進行空間一致性檢查，提出的疑誤描述為“超過鄰近站點降水范圍”，質(zhì)控后出現(xiàn)疑誤109條。經(jīng)時空一致性檢查后，2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件總共提出疑誤1156條。這表明除每日8h、20時的12h累計降水量外，還有大量累計降水量存在疑誤但未得到有效質(zhì)控。

5 結語

本文利用Python設計了一種基于時空一致性的累計降水量質(zhì)控方法，有效彌補了當前業(yè)務中質(zhì)控程序的不足，有利于提升觀測數(shù)據(jù)質(zhì)量，滿足用戶使用數(shù)據(jù)時的多種需求。應用結果表明，當前還有大量累計降水量數(shù)據(jù)未得到有效質(zhì)控。該質(zhì)控程序能夠高效完成數(shù)據(jù)完整性檢查及質(zhì)控流程，可作為補充質(zhì)控手段應用于實際業(yè)務中。目前，該質(zhì)控方法主要應用于國家站，除國家站外，四川省內(nèi)還有數(shù)以倍計的區(qū)域站，涉及的數(shù)據(jù)量更多。同樣，目前僅對區(qū)域站每日8時、20時的12h累計降水量進行質(zhì)控，有的區(qū)域站僅上傳逐小時降水數(shù)據(jù)，這一類站沒有累計降水量數(shù)據(jù)，也不對其累計降水量進行質(zhì)控。因此，下一步將該質(zhì)控方法推廣到區(qū)域站，針對僅上傳逐小時降水數(shù)據(jù)的站還可補齊其不同時間步長的累計降水量，并且由于數(shù)據(jù)量的增加，還可以進一步優(yōu)化算法，節(jié)省質(zhì)控時間。

參考文獻

[1]陳姝宜，高慶九，譚桂容，等.CRA40再分析資料在“21·7”極端強降水研究中適用性評估[J].氣象科學，2024，44（1）：82-93.

[2]蓋曉波，陳明，宮麗艷，等.氣象資料業(yè)務系統(tǒng)MDOS降水類疑誤數(shù)據(jù)處理[J].農(nóng)業(yè)與技術，2016，36（16）：209-209.

[3]中國氣象局.地面氣象觀測規(guī)范[M].北京：氣象出版社，2003.

[4]宋中玲，張勇.雙套自動氣象站數(shù)據(jù)質(zhì)控及參數(shù)設置應用探討[J].自動化儀表，2020，41（8）：106-110.

[5]文圣豐，于翔，張亞平.國產(chǎn)平臺下機載矢量地圖設計[J].航空計算技術，2023，53（2）：88-91，95.

[6]徐慧琳，趙鑫，于波，等.一種多分辨率特征提取紅外圖像語義分割算法[J].紅外技術，2024，46（5）：556-564.