大規(guī)模降雨監(jiān)測數(shù)據(jù)異常識別方法

田濟揚，劉含影，劉榮華，丁留謙，劉 宇

（1．中國水利水電科學(xué)研究院，北京 100038；2．水利部防洪抗旱減災(zāi)工程技術(shù)研究中心，北京 100038）

1 研究背景

降雨監(jiān)測是水文監(jiān)測的重要組成部分，是暴雨洪水災(zāi)害防御工作的耳目和參謀[1-2]。21世紀(jì)以來，水利部門加大了對自動監(jiān)測站建設(shè)的支持力度，特別是經(jīng)過山洪災(zāi)害防治項目建設(shè)，全國山洪災(zāi)害自動監(jiān)測站點達13．2萬個，自動雨量站網(wǎng)平均密度為38 km2/站，是2006年（6000站）的22倍，最小報汛時段縮短到10 min，數(shù)據(jù)信息量增加100余倍，極大地縮小了降雨監(jiān)測盲區(qū)，有力支撐了水旱災(zāi)害防御工作[3-5]。但由于部分測站建設(shè)標(biāo)準(zhǔn)偏低，位于山丘區(qū)的測站運維困難，數(shù)據(jù)質(zhì)量難以得到保障，常常出現(xiàn)冒大數(shù)、缺測、漏測等情況，且測站出現(xiàn)問題具有很強的隨機性，完全拋棄某一測站不切實際。

為使所有測站的監(jiān)測數(shù)據(jù)得到有效利用，需從眾多雨量監(jiān)測站中找出不同時段內(nèi)監(jiān)測數(shù)據(jù)準(zhǔn)確的站點，剔除數(shù)據(jù)質(zhì)量存在問題的站點。目前常用的異常站點識別方法有Hampel法、肖維勒準(zhǔn)則、格拉布斯準(zhǔn)則、K-medoids聚類法、離群點監(jiān)測等。Person等[6]提出了Hampel過濾方法，用中位數(shù)和中位絕對偏差代替對異常點敏感的平均數(shù)和標(biāo)準(zhǔn)差，使異常數(shù)據(jù)處理更有效。肖維勒準(zhǔn)則的置信概率偏低，且概率分布不均勻，對異常數(shù)據(jù)的剔除存在誤判的可能性[7]。格拉布斯準(zhǔn)則適用于樣本數(shù)量在[3，100]區(qū)間的異常數(shù)據(jù)判別[8]。在逐時雨量異常數(shù)據(jù)判別時，聚類分析契合了逐時降雨數(shù)據(jù)量大且變異性強的特點，但數(shù)據(jù)離散程度較大且區(qū)段劃分具有極強的主觀性[9]。離群點檢測的優(yōu)點是無需知道數(shù)據(jù)集的分布模型和分布參數(shù)，適用于存在某種距離度量手段的任何維度特征空間，其缺點是需要進行大量的距離計算與比較，對計算資源要求高[10]。

業(yè)務(wù)應(yīng)用對異常站點識別方法的高效性、穩(wěn)定性和可靠性均提出很高的要求。本文以福建全省雨量監(jiān)測站為例，基于Hampel法、格拉布斯準(zhǔn)則及周邊測站法等方法建立起遞進式異常站點篩查體系，并引入K-d tree（K-dimension tree）高級數(shù)據(jù)結(jié)構(gòu)和并行計算提高計算效率，以期為異常站點的快速識別、控制雨量監(jiān)測數(shù)據(jù)質(zhì)量提供參考。

2 研究區(qū)與數(shù)據(jù)

福建省位于我國東南沿海，山丘區(qū)面積占80%以上，且氣象條件復(fù)雜多變，受臺風(fēng)影響較大，暴雨洪澇災(zāi)害頻發(fā)[5]。為有效應(yīng)對強降水及其可能帶來的洪澇災(zāi)害，福建在全省范圍內(nèi)建設(shè)了大量的雨量監(jiān)測站，站網(wǎng)密度達到25 km2/站，居全國前列。研究選用雨量監(jiān)測站的數(shù)據(jù)來源于福建省水利廳。按照2020年的統(tǒng)計數(shù)據(jù)，福建全省發(fā)揮雨量監(jiān)測作用的各類雨量監(jiān)測站5234個，其中雨量站、河道水文站、水庫水文站、河道水位站個數(shù)分別為3268、119、1176、671，測站分布如圖1所示。氣象部門雨量監(jiān)測站報訊時間間隔為10 min，其余各類測站報訊時間間隔多為1 h。為便于異常站點識別，研究采用的雨量監(jiān)測數(shù)據(jù)時間間隔統(tǒng)一為1 h。自2010年山洪災(zāi)害防治項目實施以來，福建省雨量監(jiān)測站數(shù)量迅速增加，至2015年數(shù)量基本穩(wěn)定，后續(xù)每年仍有新建測站，但數(shù)量較少，大部分為更新改造。因此，研究選用福建全省雨量監(jiān)測站降雨數(shù)據(jù)的時間序列為2015—2020年。

圖1 福建省地形及站點分布圖

3 異常站點識別方法

遞進式異常站點篩查體系由初步判定基準(zhǔn)站、周邊測站分析和雷達輔助校驗三部分構(gòu)成（見圖2）。首先，基于長序列降雨觀測資料初步判定基準(zhǔn)站，采用Hampel法和改進的格拉布斯準(zhǔn)則初步篩選出雨量監(jiān)測較為穩(wěn)定、數(shù)據(jù)質(zhì)量相對較高的測站，稱為基準(zhǔn)站，并剔除存在明顯問題的測站；再以初步篩選出的基準(zhǔn)站為基準(zhǔn)，采用周邊測站分析法對逐小時的降雨監(jiān)測數(shù)據(jù)進行異常識別；最后通過雷達輔助校驗，對篩選出的異常站點做進一步驗證。

圖2 異常識別方法流程圖

3．1 初步判定基準(zhǔn)站

3．1．1 Hampel法 Hampel法可用于異常極端值的判別，其基本原理是對給定的數(shù)據(jù)集假設(shè)一個分布和概率模型，然后根據(jù)假設(shè)采用不一致檢驗對數(shù)據(jù)系列進行處理[11-13]。基于長序列測站的年雨量值，利用Hampel法對單一測站監(jiān)測數(shù)據(jù)的異常年份進行識別，方法如下：

式中：Xi為數(shù)據(jù)序列X中的某一值；Median為X的中位數(shù)；MAD（median absolute deviation）為數(shù)據(jù)集Y的中位數(shù)；X＝｛x1，x2，…，xn｝，是測站年降雨量數(shù)據(jù)序列；Y＝｛y1，y2，…，yn｝＝｛x1-Median，x2-Median，…，xn-Median｝。當(dāng) Zi（i＝1，2，…，n）值大于 2．24時，則判定 Xi為異常點，i為該測站的異常年份[14]。但考慮到大量測站建設(shè)年代較近，數(shù)據(jù)序列較短，僅通過Hampel法從時間維度上判定測站監(jiān)測數(shù)據(jù)的異常年份，可靠性還不夠。

3．1．2 改進的格拉布斯準(zhǔn)則 格拉布斯準(zhǔn)則適用于測量次數(shù)較少的情況（3≤n＜100），可一次性求出多個異常值。其基本原理是判斷可疑值取舍的過程中，將正態(tài)分布中平均值和方差這兩個最重要的參數(shù)引進來，提高判斷的準(zhǔn)確性[15-17]?？紤]到本文單站年累積雨量的時間序列較短，為了提高異常值判定的準(zhǔn)確性，需借助周邊測站從空間維度再做判斷，格拉布斯準(zhǔn)則判別法較為適用。改進的格拉布斯準(zhǔn)則是將原準(zhǔn)則公式中的平均值用中位數(shù)代替，可有效消除同側(cè)異常值的屏蔽效應(yīng)，更為穩(wěn)健[18]。具體方法如下。

先將通過Hampel法初步判定為某年異常的測站選出，在其周圍以20 km為半徑劃定區(qū)域，區(qū)域內(nèi)測站約50個左右，區(qū)域內(nèi)所有站點的年雨量值構(gòu)成樣本，通過從小至大排序為樣本序列X＝（x1，x2，…，xn），統(tǒng)計臨界系數(shù) G（a，n）的值G0（查臨界值表獲得），a為顯著性水平，本文取 a為0．05，并計算 G1、Gn：

式中：n為測站數(shù)量；X中為樣本中位數(shù)；σ為標(biāo)準(zhǔn)差。

若 G1≥Gn且G1＞G0，則判定x1為異常值并予以剔除；若Gn≥G1且Gn＞G0，則xn為異常值并予以剔除；若G1＜G0且Gn＜G0，則不存在異常值。若存在異常值，則剔除之后用剩余站點的年雨量值重新計算，重復(fù)上述步驟，直到無異常值為止。

為提高計算效率，縮短判斷待評估測站與周邊站點距離關(guān)系時程序的運行時長，在使用格拉布斯準(zhǔn)則對疑似異常站點進行判斷時，根據(jù)面積、測站分布、降雨空間分布情況等因素，將福建省分為7個區(qū)域，程序在查找20 km范圍內(nèi)測站時，僅在其中一個區(qū)域內(nèi)查找[19-20]。已經(jīng)判斷為異常的測站在后續(xù)判別待測站時不作為周邊測站參與比較，減少部分異常值對判別效果的不良影響。通過Hampel法和格拉布斯準(zhǔn)則共同判斷后，則完成基準(zhǔn)站初步判定。

3．2 周邊測站分析采用周邊測站分析時，優(yōu)先選用基準(zhǔn)站與待評估測站進行同時段雨量比較，基準(zhǔn)站距離較遠的（超過某一閾值），則選用距離待評估測站較近且經(jīng)基準(zhǔn)站評估后已確定為正常站點的非基準(zhǔn)站，與待評估測站進行同時段雨量比較。為了避免初步確定的基準(zhǔn)站在某一時刻降雨監(jiān)測出現(xiàn)問題，在采用基準(zhǔn)站評估待評估測站時，待評估測站也包括基準(zhǔn)站。通過比較待評估測站雨量與同時段周邊基準(zhǔn)站（或已評估合格的非基準(zhǔn)站）的平均雨量，判斷待評估測站是否異常。周邊測站分析僅對1 h（或3 h、6 h）雨量超過10 mm、12 h雨量超過15 mm、24 h雨量超過25 mm的測站進行評估。評估時分別比較待測站不同時間段（1 h、3 h、6 h、12 h、24 h）雨量值與周邊站點相應(yīng)時段的平均雨量值，當(dāng)雨量值相差超過一個等級時，則判定為異常站點。參考氣象部門對降雨等級的劃分規(guī)范，結(jié)合福建省的降雨特點，將雨量等級劃分如表1所示。

表1 雨量等級表

為選取合適的評估范圍，計算距離待評估測站5 km、10 km、15 km、20 km、25 km、30 km時，待評估測站的周邊站點平均數(shù)量分別為4個、13個、30個、50個、78個、115個，分析了距離選取對異常站點識別準(zhǔn)確率及計算時長的影響，具體見圖3?？紤]到福建省雨量監(jiān)測站網(wǎng)密度約為25 km2/站，且準(zhǔn)確率在15 km時達到最大，計算時長相對適中，故采用周邊測站分析時，分析范圍的半徑設(shè)為15 km，分析范圍內(nèi)站點平均數(shù)量30個。

圖3 距離與計算時長/準(zhǔn)確率的對應(yīng)關(guān)系

3．3 雷達輔助校驗經(jīng)過基準(zhǔn)站初步判定和周邊測站分析后，已完成異常站點的篩查，但人工校驗時發(fā)現(xiàn)，雨區(qū)與非雨區(qū)邊界、雨強差異較大的雨區(qū)邊界處，報訊正常的測站也容易被誤認為異常站點，因此仍需進一步驗證篩查結(jié)果。盡管雷達降雨反演的精度受反演算法等因素的影響大，但利用雷達反射率仍能有效判定雷達探測覆蓋范圍是否產(chǎn)生降雨或判斷降雨量級，能夠充分反映某一時段降雨的空間分布特征[21-23]。雷達輔助校驗方法包括：一是利用雷達低層仰角反射率超過20 dBZ閾值時即發(fā)生降雨的條件，驗證雨區(qū)與非雨區(qū)邊界處測站判定結(jié)果；二是通過雷達基數(shù)據(jù)反演降雨強度，并與測站降雨量級進行對比，驗證測站是否異常；三是通過反射率的空間變化梯度，驗證雨強差異較大的雨區(qū)邊界處測站的判定結(jié)果。

3．4 大規(guī)模降雨監(jiān)測數(shù)據(jù)處理算法福建省雨量監(jiān)測站數(shù)量多、密度大。為了能夠達到實時計算并判斷異常站點的目的，研究采用K-d tree高級數(shù)據(jù)結(jié)構(gòu)和并行計算方法，大幅提高計算效率。經(jīng)測試，全省站點進行一次異常識別，計算時長約5～8 min。

3．4．1 K-d tree高級數(shù)據(jù)結(jié)構(gòu) K-d tree高級數(shù)據(jù)結(jié)構(gòu)主要用于優(yōu)化一定范圍內(nèi)周邊站點的篩選過程。K-d tree是對數(shù)據(jù)點在 k維空間（如二維（x，y），三維（x，y，z），k維（x，y，z，…））中劃分的一種數(shù)據(jù)結(jié)構(gòu)。為了能有效的找到最近鄰，K-d tree將整個空間劃分為幾個小部分，且K-d tree索引的空間劃分不會出現(xiàn)區(qū)域重疊現(xiàn)象，更適合作為計算環(huán)境中的上層全局索引，從而在多維查詢過程中快速發(fā)現(xiàn)包含查詢結(jié)果的局部數(shù)據(jù)節(jié)點[24-25]。

3．4．2 并行計算 針對多個測站同時進行異常值識別這一問題，研究采用 CUDA（Compute Unified Device Architecture）平臺，使用支持CUDA的GPU（Graphics Processing Unit）進行并行編程，與傳統(tǒng)的GPU相比，CUDA的GPU能獲得相比同期CPU（Central Processing Unit）高幾倍乃至十幾倍的提速，大幅縮短計算時長[26-29]。該方法是對多個獨立站點同時進行異常識別，把來自不同站點的指令在不同的處理器上同時執(zhí)行，提高計算效率[30]。Hampel法和格拉布斯準(zhǔn)則所用數(shù)據(jù)為年尺度，不進行實時計算，計算量相對較小，對計算的時效性要求不高，僅在格拉布斯準(zhǔn)則法分區(qū)計算時采用并行算法。周邊測站法所用數(shù)據(jù)為逐小時降雨數(shù)據(jù)，數(shù)據(jù)量大且需要實時計算，對計算的時效性要求很高，因此周邊測站法全程采用并行計算將不同站點分配到不同的處理器進行同時計算。雷達輔助校驗的計算量取決于篩查出的異常站點數(shù)量，盡管時效性要求高，但計算量較小，未采用并行計算。

目前福建全省發(fā)揮雨量監(jiān)測作用的各類雨量監(jiān)測站5234個，某一時刻站點降雨數(shù)據(jù)量最多有5234條，以計算量最大的周邊測站分析法為例，按分析范圍15 km計，某一待檢測站周邊站點降雨數(shù)據(jù)約30條，采用支持CUDA的GPU并行計算可將計算時長控制在5～8 min。

4 結(jié)果分析

4．1 基準(zhǔn)站初步判定結(jié)果經(jīng)過Hampel法和格拉布斯準(zhǔn)則判別法對2015—2020年異常站點進行判定發(fā)現(xiàn)，2015年的異常站點數(shù)量最多，占比11．5%，之后異常站點逐年減少，至2020年，異常站點數(shù)量占比僅5．18%。主要原因是汛后逐年加強了測站的運行維護，測站監(jiān)測數(shù)據(jù)的質(zhì)量有所提升。異常站點中雨量站、水庫站的比例偏高，2015—2020年雨量站異常站點在所有異常站點中的占比分別為42．83%、44．13%、48．25%、47．83%、50．84%、47．97%，水庫水文站異常站點在所有異常站點中的占比分別為 35．33%、35．20%、31．88%、29．83%、25．18%、34．69%。分析計算 2015—2020年每年各類異常站點在全省相應(yīng)類型站點中的占比，發(fā)現(xiàn)水庫水文站、河道水位站的占比更高，各類異常站點的空間分布情況見圖5，異常站點數(shù)量及占比見圖6。為檢驗基準(zhǔn)站初步判定結(jié)果，將2015—2020年每年判定的異常站點與實際異常站點進行對比，基準(zhǔn)站判定的準(zhǔn)確率為95．4%，詳見表2，其中異常站點的類型主要有三種：一是年降雨總量為0的站點（站點不報數(shù)）；二是年降雨總量遠小于周邊測站的站點（站點年內(nèi)長時期不報數(shù)）；三是年降雨總量極大且遠高于周邊測站的站點。典型實例如圖7所示。

表2 判定異常站和實際異常站點數(shù)量 （單位：個）

圖5 2015—2020年異常站點分布圖

圖6 2015—2020年異常站點數(shù)量及占比情況

圖7 異常站點類型

4．2 基于周邊測站分析的異常站點識別結(jié)果由于每年的基準(zhǔn)站會發(fā)生一定的變化，因此在周邊測站分析時取前一年的基準(zhǔn)站為基準(zhǔn)開展周邊測站分析，如對2017年逐小時降雨數(shù)據(jù)進行異常識別時，選取基于2016年降雨資料判定的基準(zhǔn)站。分別選取2016—2020年每年的6月、7月和8月的1日8點、10日14點、30日20點的異常值識別結(jié)果進行人工驗證，異常識別正確率絕大多數(shù)超過了90%（見表3）。為進一步驗證方法的可用性，以2021年6月27—28日福建全省的雨量監(jiān)測站進行了異常站點實時判別，結(jié)果如圖8和圖9所示。

圖8 6月27日、28日4個時刻異常站點分布情況

圖9 6月27日、28日4個時刻異常站點數(shù)量及占比情況

表3 雨量異常值識別正確率 （單位：%）

通過周邊測站分析法對2021年6月27日和28日的逐小時降雨數(shù)據(jù)進行判定并選取共計8個時刻的結(jié)果進行統(tǒng)計分析。結(jié)果表明6月28日15∶00的異常站點數(shù)量最多，占全省站點的11．6%，6月27日3∶00的異常站點數(shù)量最少，占全省站點的3．8%，6月27日和28日21點的異常站點數(shù)量均較多，主要原因是福建全省在6月 27日19∶00—21∶00和6月 28日9∶00—15∶00降雨范圍較廣，報訊站較多，更容易發(fā)現(xiàn)異常站點。對于不同類型測站，雨量站異常站點占全部異常站點比例依然最大，其次是水庫水文站、河道水位站，占比最小的是河道水文站，各類異常雨量監(jiān)測站的數(shù)量情況見圖9（a）。分析計算8個時刻各類異常站點在全省相應(yīng)類型站點中的占比，發(fā)現(xiàn)雨量站、河道水文站的占比相對較高，具體情況見圖9（b）。

經(jīng)周邊測站分析法檢測出的異常站點主要分為三類：一是實際發(fā)生降雨但未報數(shù)（降雨值為0）的站點；二是雨量值明顯小于周邊測站的站點；三是降雨值極大且明顯高于周邊測站的站點（冒大數(shù)）。在三類異常站點中，處于不同等級降雨區(qū)域分界處的站點較難識別，典型實例如圖10所示。

圖10 異常站點類型

4．3 異常站點識別結(jié)果校驗與準(zhǔn)確率經(jīng)過基準(zhǔn)站初步判定和周邊測站分析后，已初步完成異常站點的篩查，但人工校驗時發(fā)現(xiàn)，處于雨區(qū)與非雨區(qū)邊界、雨強差異較大的雨區(qū)邊界的正常站點也容易被判定為異常值，因此仍需借助雷達反射率進一步驗證。在2021年6月28日11點的異常結(jié)果中選取屬于這類情況的4個站點，將站點雨量圖與雷達回波圖疊加進行人工判斷，具體情況見圖11。盡管通過上述方法將該4個站點判定為異常，但經(jīng)過雷達回波圖校驗發(fā)現(xiàn)，10∶00—11∶00，川里村站1 h雨量為0．4 mm，光澤西關(guān)站1 h雨量為1．5 mm，楓林村站1 h雨量為0．5 mm，王村站1h雨量為22．2 mm。10∶00—11∶00每間隔6 m in完成一次提掃并獲得一次雷達反射率，計算1 h內(nèi)4個站點位置處對應(yīng)的雷達反射率均值分別為14 dBz、25 dBz、27 dBz、20 dBz。對比雷達反射率和站點小時雨量可知，川里村站、光澤西關(guān)站無異常，楓林村站的雨量值偏小，而王村站的雨量值偏大。

圖11 異常站點識別結(jié)果校驗

考慮到福建全省在6月27日19∶00—21∶00和6月28日09∶00—15∶00降雨范圍較廣，報訊站較多，所以選擇如下6個時刻的異常結(jié)果進行雷達回波的校驗，具體情況如下表4所示。雷達校驗前異常識別結(jié)果的平均準(zhǔn)確率為89%，經(jīng)雷達校驗后平均準(zhǔn)確率提升至95%，這表明雷達輔助校驗方法非常適用于應(yīng)對處于雨區(qū)與非雨區(qū)邊界、雨強差異較大的雨區(qū)邊界的正常站點被錯誤判斷為異常站點的情況。

表4 6個時刻雷達校驗前后福建省雨量監(jiān)測站異常識別結(jié)果準(zhǔn)確率 （單位：%）

5 結(jié)論

本研究利用福建省5234個雨量站2015—2021年實測降雨資料，基于Hampel法、格拉布斯準(zhǔn)則、周邊測站分析法和雷達輔助校驗等方法構(gòu)建了遞進式異常站點篩查體系，對降雨監(jiān)測數(shù)據(jù)進行了異常識別。主要結(jié)論如下：

（1）基于測站年累積雨量，采用Hampel法和格拉布斯準(zhǔn)則對異常站點進行識別，結(jié)果表明異常站點類型主要有測站全年不報數(shù)、測站年內(nèi)長時期不報數(shù)、年降雨總量遠高于周邊測站3類。在2015—2020年，福建省雨量監(jiān)測站數(shù)據(jù)質(zhì)量不斷提升，異常站點數(shù)量逐年減少。其中，2015年異常站點數(shù)量最多，占比為11．5%；2020年異常站點數(shù)量最少，占比為5．18%。此外，福建省內(nèi)雨量站和水庫數(shù)量較其他類型站點明顯偏多，異常站點數(shù)量也較其他類型異常站點更多，2015—2020年雨量站異常站點占全部異常站點比例為42．83%～50．84%，水庫水文站異常站點占全部異常站點比例為25．18%～35．33%。

（2）采用周邊測站分析法對2016—2020年間9個時刻全省雨量監(jiān)測站進行異常識別，異常識別正確率超過90%，并以2021年6月27日和28日逐小時降雨監(jiān)測數(shù)據(jù)進行了進一步驗證，結(jié)果表明異常站點類型主要有實際發(fā)生降雨但未報數(shù)、明顯小于實際雨量、冒大數(shù)3類，處于不同等級降雨區(qū)域分界處更難準(zhǔn)確識別。當(dāng)降雨范圍較大時，更容易發(fā)現(xiàn)異常站點，雨量站異常站點占全部異常站點的比例最高，各類異常站點在全省相應(yīng)類型站點中，雨量站異常站點的占比也最高。

（3）通過雷達輔助校驗，6個典型時刻的異常識別平均準(zhǔn)確率由雷達校驗前的89%提升至95%，表明雷達輔助校驗方法非常適用于應(yīng)對處于雨區(qū)與非雨區(qū)邊界、雨強差異較大的雨區(qū)邊界的正常站點被錯誤判斷為異常站點的情況。

（4）本研究采用K-d tree高級數(shù)據(jù)結(jié)構(gòu)和并行計算方法，大幅提高了計算效率。經(jīng)測試，全省站點進行一次異常識別，計算時長約5～8 min。