累計(jì)降水量質(zhì)控方法設(shè)計(jì)及應(yīng)用

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

關(guān)鍵詞：降水；累計(jì)降水量；質(zhì)量控制方法；自動(dòng)化

中圖分類號：TP 399" " 文獻(xiàn)標(biāo)志碼：A

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

1 累計(jì)降水量

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

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

2.1 缺測率統(tǒng)計(jì)

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

（1）

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

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

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

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

（2）

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

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

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

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

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

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

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

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

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

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

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

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

地面站資料儲(chǔ)存在數(shù)據(jù)庫中，需要連接數(shù)據(jù)庫選取所需時(shí)段內(nèi)的降水?dāng)?shù)據(jù)?；赑ython的cx_Oracle，使用用戶名及密碼連接數(shù)據(jù)庫。執(zhí)行以下SQL（Structured Query Language）語句選取2023年1月1日01時(shí)至2024年1月1日00時(shí)的四川省內(nèi)國家站地面逐小時(shí)文件中的1h降水量和3h/6h/12h/24h累計(jì)降水量數(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累計(jì)降水量數(shù)據(jù)。將查詢結(jié)果保存為DataFrame，有利于后續(xù)質(zhì)控的處理。

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

檢查所有地面站的所有時(shí)次文件中1h/3h/6h/12h/24h累計(jì)降水量是否完整。統(tǒng)計(jì)缺測時(shí)次，得到缺測率。本功能使用datetime庫計(jì)算總時(shí)次，使用Pandas的數(shù)據(jù)結(jié)構(gòu)DataFrame統(tǒng)計(jì)已收時(shí)次，其差值為缺測時(shí)次，可計(jì)算各步長累計(jì)降水量缺測率。具體實(shí)現(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）為已收時(shí)次；duration_hours為2023年1月1日01時(shí)至2024年1月1日00時(shí)的總時(shí)數(shù)，即應(yīng)收總時(shí)次；R為統(tǒng)計(jì)的缺測率，可以用百分?jǐn)?shù)表示。

3.3 時(shí)間一致性檢查

分臺(tái)站依時(shí)次讀取地面站降水量數(shù)據(jù)，以讀取時(shí)次作為被檢驗(yàn)時(shí)次。檢查被檢臺(tái)站的被檢驗(yàn)時(shí)次前3h/6h/12h/24h內(nèi)降水?dāng)?shù)據(jù)是否完整，本功能使用datetime庫生成時(shí)間戳，檢查原文件是否包含該時(shí)次過去3h/6h/12h/24h內(nèi)的所有時(shí)間戳。具體實(shí)現(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']為原文件中被檢時(shí)次的前3h/6h/12h/24h內(nèi)的逐小時(shí)降水量。針對缺測的逐小時(shí)累計(jì)降水量，其累計(jì)值也置為缺測保存。檢查該被檢時(shí)次是否為定時(shí)觀測時(shí)次，如果是，那么質(zhì)控碼置為1，向臺(tái)站下發(fā)疑誤。如果為非定時(shí)觀測時(shí)次，那么直接修改原文件該時(shí)次累計(jì)值為缺測，質(zhì)控碼置為8。

經(jīng)過完整性檢查，符合條件的時(shí)次可以進(jìn)行累計(jì)降水量的計(jì)算。累計(jì)降水量p，得到計(jì)算值TP。判斷計(jì)算值TP與原文件中相應(yīng)時(shí)間步長的累計(jì)降水量是否一致。具體實(shí)現(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ù)中的累計(jì)降水量；rain_hour為相應(yīng)時(shí)間步長的累計(jì)降水量列名（3.1中提到）；np.allclose用于判斷數(shù)組是否在誤差范圍內(nèi)相等；atol為絕對容忍度，在這里表示TP和rain_sum在0.1的誤差范圍內(nèi)相等；誤差設(shè)置為0.1是由于地面降水量觀測值的精度為0.1mm；code為質(zhì)控碼，如果TP與rain_sum不相等，那么質(zhì)控碼為1，向臺(tái)站下發(fā)疑誤。如果相等，那么質(zhì)控碼為0，相應(yīng)累積降水量置為正確。

3.4 空間一致性檢查

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

Apoint=（Alat，Alon）

Bpoint=（Blat，Blon）

distance=great_circle（Apoint，Bpoint）.km

其中，Apoint為被檢點(diǎn)；Bpoint為可能鄰近站。

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

4 應(yīng)用結(jié)果

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

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

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

5 結(jié)語

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

參考文獻(xiàn)

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

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

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

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

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

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