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

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

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

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

1 累計降水量

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

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

2.1 缺測率統(tǒng)計

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

（1）

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

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

設(shè)計累計降水量的質(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)的累計降水量，與原文件進(jìn)行對比，這是通過時間一致性的原理對降水進(jìn)行質(zhì)控。另外，對比站點(diǎn)及其附近站點(diǎn)的累計降水量，通過空間一致性原理對降水量進(jìn)行質(zhì)控。質(zhì)控總體分以下步驟。第一步：進(jìn)行時間一致性檢查，如圖1所示。檢查時間步長內(nèi)的逐小時降水量是否完整。當(dāng)遍歷各時次文件時，對當(dāng)前時次的前3h/6h/12h/24h內(nèi)的逐小時降水量數(shù)據(jù)進(jìn)行檢查，如果有缺失，判斷當(dāng)前時次是否為定時觀測時次——基準(zhǔn)站、基本站人工定時觀測次數(shù)為每日5次（08時、11時、14時、17時、20時），一般站人工定時觀測次數(shù)為每日3次（08時、14時、20時）。如果為定時觀測時次，就通過質(zhì)控系統(tǒng)向臺站自動下發(fā)疑誤，質(zhì)控碼置為1。如果為非定時觀測時次，就直接將當(dāng)前時次文件內(nèi)相應(yīng)累計降水量數(shù)據(jù)作缺，質(zhì)控碼置為8[3]，常見質(zhì)控碼及含義見表1。

第二步：計算累計降水量。如果當(dāng)前時次的前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)臺站確認(rèn)反饋后可替換原文件中累計降水量。

第三步：對比計算值TPn與原文件累計值是否一致。原文件累計值通常為臺站端本地計算后上傳，有時僅修改錯誤的小時降水量，而忽略了累計降水量。如果計算值TPn與原文件累計值一致，那么判定原文件中累計降水量質(zhì)控結(jié)果為正確，質(zhì)控碼置為0。如果不一致，那么判定原文件累計降水量可疑，質(zhì)控碼置為1，向臺站下發(fā)相應(yīng)疑誤并提供計算值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)，再計算站點(diǎn)間距離，保留附近5km內(nèi)的站點(diǎn)[6]。這里計算的鄰近站點(diǎn)可以保存為常用表，不僅可用于降水量的空間一致性檢驗(yàn)，而且還可適用于其他多種氣象要素，例如地溫、氣溫、氣壓等的檢驗(yàn)或插值等。

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

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

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

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

地面站資料儲存在數(shù)據(jù)庫中，需要連接數(shù)據(jù)庫選取所需時段內(nèi)的降水?dāng)?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ù)。將查詢結(jié)果保存為DataFrame，有利于后續(xù)質(zhì)控的處理。

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

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

3.3 時間一致性檢查

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

經(jīng)過完整性檢查，符合條件的時次可以進(jìn)行累計降水量的計算。累計降水量p，得到計算值TP。判斷計算值TP與原文件中相應(yīng)時間步長的累計降水量是否一致。具體實(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ù)中的累計降水量；rain_hour為相應(yīng)時間步長的累計降水量列名（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，向臺站下發(fā)疑誤。如果相等，那么質(zhì)控碼為0，相應(yīng)累積降水量置為正確。

3.4 空間一致性檢查

經(jīng)過時間一致性檢查后，質(zhì)控碼為1的數(shù)據(jù)下發(fā)臺站質(zhì)疑，質(zhì)控碼為8的數(shù)據(jù)直接作缺測處理。數(shù)據(jù)質(zhì)控碼為0和4的累計降水量數(shù)據(jù)通過了時間一致性檢驗(yàn)，可進(jìn)行進(jìn)一步的空間一致性檢查。首先，挑選被檢點(diǎn)附近5km內(nèi)的站點(diǎn)，為減少計算量，篩選被檢點(diǎn)附近0.05°×0.05°的站點(diǎn)，經(jīng)緯度上1°約為110km，0.05°約為5.5km。其次，計算在此范圍內(nèi)的站點(diǎn)與被檢點(diǎn)的實(shí)際距離，使用geopy庫中的great_circle包，其計算站點(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)可能會導(dǎo)致站點(diǎn)間的降水相關(guān)性較小，而缺失其參考價值。針對一些地處甘孜州、阿壩州及涼山州等偏遠(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)時間步長的累計降水量值，其中最大值作為上限，最小值作下限，判斷被檢點(diǎn)的值是否在此范圍內(nèi)。如果不在范圍內(nèi)，那么質(zhì)控碼為1，疑誤下發(fā)臺站。如果在范圍內(nèi)，那么質(zhì)控碼為0。

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

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

再執(zhí)行時間一致性檢查，提出的疑誤如下。1）時間步長內(nèi)的逐小時降水正常，但3h/6h/12h/24h累計降水為缺測，針對這類疑誤，可應(yīng)用質(zhì)控程序的計算值補(bǔ)充缺失的累計降水量。2）質(zhì)控程序計算值與原文件累計降水量不一致，針對這一類疑誤，需要向臺站再次核實(shí)時間步長內(nèi)的逐小時降水量是否正確。3）時間步長內(nèi)的逐小時降水有缺測，但原文件有累計值，針對這類疑誤，需要向臺站核實(shí)累計值是否為人工觀測，如果不是，就按照質(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ì)控結(jié)果為0和4的累計降水量進(jìn)行空間一致性檢查，提出的疑誤描述為“超過鄰近站點(diǎn)降水范圍”，質(zhì)控后出現(xiàn)疑誤109條。經(jīng)時空一致性檢查后，2023年1月1日01時至2024年1月1日00時的四川省內(nèi)國家站地面逐小時文件總共提出疑誤1156條。這表明除每日8h、20時的12h累計降水量外，還有大量累計降水量存在疑誤但未得到有效質(zhì)控。

5 結(jié)語

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

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

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

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