基于時序柵格的暴雨事件提取與追蹤方法

楊光輝， 薛存金， 劉敬一， 鄔群勇， 伍程斌

1.福州大學(xué)衛(wèi)星空間信息技術(shù)國家地方聯(lián)合工程中心，福州350116

2.中國科學(xué)院數(shù)字地球重點(diǎn)實(shí)驗(yàn)室，北京100094

3.中國科學(xué)院空天信息研究院，北京100094

綜合對地觀測和眾源技術(shù)為獲取長時序、高時空分辨率的降水時序柵格數(shù)據(jù)集提供了技術(shù)支撐[1].高時空分辨率的暴雨事件集可以為獲取暴雨影響區(qū)域的時空分布、挖掘暴雨事件的時空模式、開展面向個體的精準(zhǔn)暴雨災(zāi)害預(yù)警預(yù)報服務(wù)提供重要數(shù)據(jù)基礎(chǔ)[2].因此，基于時序柵格數(shù)據(jù)集并結(jié)合暴雨事件的時空認(rèn)知方法提出暴雨事件提取方法，進(jìn)而構(gòu)建高時空分辨率的暴雨事件集，這對暴雨災(zāi)害及城市內(nèi)澇的監(jiān)測和預(yù)測具有重要意義[3].

目前，暴雨事件的提取方法主要分為3 類：1）基于站點(diǎn)日降水?dāng)?shù)據(jù)，即根據(jù)氣象上暴雨定義直接提取閾值，如文獻(xiàn)[4]以站點(diǎn)日降水?dāng)?shù)據(jù)提取暴雨天數(shù)，分析了中國近40年暴雨發(fā)生頻率的年代際時空變化特征；文獻(xiàn)[5]根據(jù)站點(diǎn)日降水?dāng)?shù)據(jù)分析了暴雨過程和暴雨日的年、季、月氣候分布及變化特征；文獻(xiàn)[6]基于日降水?dāng)?shù)據(jù)統(tǒng)計(jì)年代際暴雨雨量、雨日、雨強(qiáng)來分析中國1951—2010年的年代際暴雨時空變化格局.2）對柵格數(shù)據(jù)進(jìn)行日累加處理，設(shè)定特定閾值進(jìn)行暴雨事件提取，如文獻(xiàn)[7]對全球降水測量(global precitation measurement,GPM)與熱帶降雨測量任務(wù)(tropical rainfall measuring mission, TRMM)降雨數(shù)據(jù)進(jìn)行日累加后提取暴雨柵格；文獻(xiàn)[8]以日累加方式處理GPM 降雨數(shù)據(jù)后進(jìn)行適用性評價.3）基于雷達(dá)反射率閾值開展時刻狀態(tài)的暴雨事件提取，并利用暴雨的移動特征實(shí)現(xiàn)暴雨事件的追蹤，如文獻(xiàn)[9]提出用暴雨識別、跟蹤、分析和預(yù)報(thumderstorm identification, tracking,analysis, and nowcasting, TITAN)方法將空間上反射率超過一定閾值的連續(xù)區(qū)域作為暴雨區(qū)域進(jìn)行識別和匹配追蹤；文獻(xiàn)[10]基于雙閾值遲滯法識別暴雨，避免了只有幾個像素的暴雨區(qū)域；文獻(xiàn)[11]基于TITAN 在識別暴雨過程中使用多閾值分割技術(shù)（hierarchical threshold segmentation, HTS），減少了暴雨過程識別中錯誤的暴雨對象合并情況.

上述方法在暴雨提取分析方面取得了一系列成果.基于站點(diǎn)的降水?dāng)?shù)據(jù)在局部尺度上（站點(diǎn)位置）具有較高的觀測精度，但在偏遠(yuǎn)地區(qū)或地形復(fù)雜地區(qū)，因?yàn)橛^測站點(diǎn)分布稀疏，所以很難提供精度較高的降水空間分布特征[12].將時間分辨率高的柵格降水?dāng)?shù)據(jù)先進(jìn)行日累加再進(jìn)行提取分析則會損失原有數(shù)據(jù)精度，以致無法提取整個暴雨事件的發(fā)展細(xì)節(jié).基于TITAN 的暴雨事件識別與追蹤算法利用了暴雨事件的移動特征，為暴雨事件的提取與重構(gòu)提供了思路，但在追蹤過程中只考慮了暴雨路徑起止位置的分裂合并，因而無法識別暴雨路徑中間分裂合并行為特征.針對以上問題，本文設(shè)計(jì)了一種基于時序柵格的暴雨事件提取與追蹤方法，實(shí)現(xiàn)了對暴雨事件的完整提取.

1 暴雨事件提取與追蹤方法

暴雨事件是指降水量超過一定閾值且在時空中連續(xù)變化的降雨事件.暴雨事件提取與追蹤方法（rainstorm event extracting and tracking method, REETM）基于長時間序列柵格數(shù)據(jù)集，依據(jù)暴雨在時空上連續(xù)漸變的過程特性[13-15]開展暴雨事件提取與追蹤.REETM 主要包括以下3 個步驟：

步驟1在時間維度上提取暴雨柵格，依據(jù)降水?dāng)?shù)據(jù)劃分暴雨等級；

步驟2在空間維度上提取暴雨對象，連接空間上孤立的暴雨柵格并統(tǒng)計(jì)暴雨信息；

步驟3暴雨事件時空追蹤，在時空維度上追蹤暴雨并重構(gòu)暴雨事件.

1.1 時間維度暴雨提取

基于長時間序列降水量柵格數(shù)據(jù)集進(jìn)行時間維度的暴雨提取處理，提取出暴雨級別信息，獲取柵格暴雨分級數(shù)據(jù)集，如圖1 所示.

圖1 時間維度上的暴雨提取Figure 1 Rainstorm extraction in a temporal domain

具體步驟如下：

步驟1取出每個柵格像元所在位置的所有圖像柵格值，形成一個按照時間順序排列的一維數(shù)組.

步驟2遍歷這個一維數(shù)組，尋找極大值點(diǎn)，如果數(shù)組任意兩個不為0 的位置時間差值小于閾值t，則認(rèn)為降水連續(xù).以極大值點(diǎn)為中心，向前、向后搜尋降水連續(xù)的序列.

步驟3統(tǒng)計(jì)降水連續(xù)的序列范圍內(nèi)的累計(jì)降水量，如果序列時間長度超過24 h，則以極大值點(diǎn)前12 h 這一時刻為統(tǒng)計(jì)起始點(diǎn)，以極大值點(diǎn)后12 h 這一時刻為統(tǒng)計(jì)結(jié)束點(diǎn).依據(jù)中國氣象局制定的《降水量等級》[16]對該暴雨序列進(jìn)行分級，分級依據(jù)如表1 所示.

表1 暴雨等級劃分Table 1 Rainstorm classification

1.2 空間維度暴雨提取

在時間維度暴雨提取結(jié)果的基礎(chǔ)上進(jìn)行空間維度暴雨提取.在時間維度上的暴雨提取保證了暴雨事件在時間維度的連續(xù)性；在空間維度上的暴雨提取不但保證了暴雨事件在空間維度的相鄰性，而且能提取出暴雨對象及其信息，如圖2 所示.

圖2 空間維度上的暴雨提取Figure 2 Rainstorm connection in a spatial domain

具體步驟如下：

步驟1采用結(jié)點(diǎn)搜索法[17]連接孤立的暴雨柵格，將柵格暴雨分級數(shù)據(jù)集轉(zhuǎn)換為矢量暴雨數(shù)據(jù)集.

步驟2遍歷矢量暴雨數(shù)據(jù)集中每一個內(nèi)部連通矢量多邊形即暴雨對象，從原始柵格數(shù)據(jù)集和柵格暴雨分級數(shù)據(jù)集中尋找出該暴雨對象覆蓋范圍內(nèi)的所有柵格像元.

步驟3統(tǒng)計(jì)暴雨對象覆蓋范圍內(nèi)所有柵格像元信息以獲取暴雨對象信息，包括時間信息、空間信息、屬性信息.時間信息即該暴雨對象的發(fā)生時刻，空間信息包括暴雨對象空間范圍、面積和周長，屬性信息包括暴雨強(qiáng)度、降水體積、平均降水量、最大降水量、最小降水量.

步驟4保存所有暴雨對象，得到完整的矢量暴雨數(shù)據(jù)集.

1.3 暴雨事件時空追蹤

暴雨事件時空追蹤的核心思想如下：一個暴雨事件在時空中是連續(xù)變化的，它由多個時刻狀態(tài)組成，且每個時刻狀態(tài)包含一個或多個暴雨對象，于是可以通過追蹤相鄰時刻空間范圍拓?fù)湎嘟坏谋┯陮ο髞韺?shí)現(xiàn)暴雨事件的追蹤.

暴雨事件時空追蹤的流程與步驟如下：同一暴雨事件在短時間內(nèi)移動距離較短，空間范圍存在重疊.因此，對于相鄰時刻暴雨對象，若其空間范圍拓?fù)湎嘟?，則認(rèn)為暴雨對象屬于同一個暴雨事件，如圖3 所示.

圖3 暴雨事件時空追蹤Figure 3 Rainstorm tracking in spatiotemporal domain

在圖3 中，綠色多邊形表示暴雨對象，T1時刻暴雨對象1 與T2時刻暴雨對象3 在空間上拓?fù)湎嘟?，則認(rèn)為暴雨對象1 與暴雨對象3 存在關(guān)系1→3，即暴雨對象1 發(fā)展為暴雨對象3，且暴雨對象1 與暴雨對象3 屬于同一個暴雨事件；同理找出所有時刻屬于同一個暴雨事件的暴雨對象，并統(tǒng)計(jì)暴雨事件信息.具體算法步驟如下：

步驟1讀取矢量暴雨數(shù)據(jù)集，取出每一個暴雨對象并將其存儲到暴雨對象集合P 中.

步驟2從第1 個時刻開始追蹤暴雨對象，令t=1.

步驟3遍歷暴雨對象集合P 中t 時刻每一個暴雨對象pi，如果集合P 中t+1 時刻存在暴雨對象pj，使得pi的空間范圍與pj的空間范圍拓?fù)湎嘟?，則記錄暴雨對象關(guān)系pi→pj，并將該暴雨對象關(guān)系存儲到暴雨對象關(guān)系集合R 中；若pi與t+1 時刻多個暴雨對象空間范圍拓?fù)湎嘟?，則記錄pi與每個相交暴雨對象的關(guān)系并存儲到集合R 中.

步驟4令t 自增1，如果t 為最后時刻，則執(zhí)行步驟5；否則執(zhí)行步驟3.

步驟5若P 中不存在暴雨對象，則執(zhí)行步驟8；否則創(chuàng)建暴雨對象集合S 以存儲屬于同一個暴雨事件的所有暴雨對象，創(chuàng)建暴雨對象關(guān)系集合R'以存儲該暴雨事件中所有暴雨對象關(guān)系.創(chuàng)建暴雨對象集合P'和P''，從暴雨對象集合P 中取出一個暴雨對象存儲到P'中.

步驟6遍歷P'中每一個暴雨對象p，若R 中存在關(guān)系p →P'，則從R 中取出該暴雨關(guān)系儲存到暴雨對象關(guān)系集合R'中，并從P 中取出暴雨對象P'存儲到P''中.

步驟7將P'中所有暴雨對象取出后存儲到S 中.若P''中存在暴雨對象，則將P''中所有暴雨對象取出儲存到P'中，執(zhí)行步驟6；否則保存暴雨對象S 中的所有暴雨對象并統(tǒng)計(jì)該暴雨事件信息（暴雨級別、平均降水量、降水體積等）后予以保存.保存屬于該暴雨事件的對象關(guān)系集合R'中的所有暴雨對象關(guān)系，執(zhí)行步驟5.

步驟8所有暴雨事件提取完畢，則算法結(jié)束.

在以上算法中需要特別指出的是：暴雨事件是一個在時空上變化的實(shí)體，其平均降水量尚未有定義，這里暫且采用的計(jì)算公式為

式中，Ravg為暴雨事件平均降水量，n 為暴雨事件中的暴雨對象數(shù)目，Ai為第i 個暴雨對象的面積，Ri為第i 個暴雨對象的平均降水量.

2 實(shí)例驗(yàn)證

本文以GPM IMERG Final Products（簡稱為GPM 產(chǎn)品數(shù)據(jù)集）作為實(shí)例數(shù)據(jù)開展暴雨事件提取工作.GPM 產(chǎn)品數(shù)據(jù)集來源于NASA（disc.gsfc.nasa.gov），其空間分辨率為0.1?，時間分辨率為0.5 h；驗(yàn)證數(shù)據(jù)集包括暴雨預(yù)警信息（www.weather.com.cn）、49 個氣象臺站的日降水?dāng)?shù)據(jù)集和2016-07-01 的天氣雷達(dá)數(shù)據(jù)集（data.cma.cn）；研究區(qū)域?yàn)?0?E～140?E,0?～60?N，時間范圍為2016-06-01～2016-09-30.

將時間維度暴雨提取過程中最大降水連續(xù)間隔時間t 設(shè)置為3 h，提取出49 003 個暴雨事件，其中暴雨級別事件有48 937 個，大暴雨級別事件有66 個，特大暴雨級別事件有0 個.

2.1 與日累加法的對比分析

采用探測率（probability of detection, POD）nPOD統(tǒng)計(jì)指標(biāo)評價衛(wèi)星降水?dāng)?shù)據(jù)對暴雨事件的探測能力，該指標(biāo)表示探測正確的數(shù)目占實(shí)際數(shù)目的比例，其值越高表示暴雨事件提取方法的漏報率越低.表2 對比了日累加法與本文提出的REETM 方法的暴雨探測率，其中D0為暴雨真實(shí)天數(shù)，D1為采用日累加方法探測得到的暴雨天數(shù)，D2為采用本文方法探測得到的暴雨天數(shù).通過計(jì)算得出本文方法的探測率平均比日累加法高12.97%，原因是日累加法只統(tǒng)計(jì)每日固定起止時間累計(jì)降水量來識別暴雨，而本文方法統(tǒng)計(jì)累計(jì)降水量的起止時間是動態(tài)變化的，可以探測出固定時間段內(nèi)探測不到的暴雨，所以探測率較高.

表2 日累加法與REETM 的暴雨探測率對比Table 2 Comparison of rainstorm POD between daily accumulation method and REETM

2.2 與暴雨預(yù)警信息的對比分析

選取關(guān)鍵時刻展示一個暴雨事件的演變過程如圖4 所示，紫色多邊形為本方法提取出的暴雨事件覆蓋區(qū)域.該暴雨事件的起止時間是2016-07-18～2016-07-21，起源于湖南和陜西，經(jīng)過發(fā)展合并分裂等過程影響全國多個地區(qū)，由東北地區(qū)離開中國國境并逐漸消亡.2016-07-18 8:00 陜西出現(xiàn)暴雨，由西向東移動并逐漸擴(kuò)大，途徑陜西北部、山西大部分區(qū)域、河北大部分區(qū)域與河南北部.同時，湖南區(qū)域暴雨，由西南向東北移動并逐漸擴(kuò)大，途徑湖南北部、湖北中部和河南南部，于2016-07-20 0:00 在河南中部與起源于陜西的暴雨合并，然后繼續(xù)向東北移動，于2016-07-20 2:00 分裂為兩部分，一部分逐漸縮小，于2016-07-20 22:00 在河南境內(nèi)消亡；另一部分繼續(xù)向東北移動，途徑河南、河北、山東、北京、天津和遼寧等地，于2016-07-21 13:00 逐漸離開中國國境并消亡.將本暴雨事件與多地氣象臺預(yù)警信息進(jìn)行對比，可以發(fā)現(xiàn)本文提出的REETM 方法所提取的暴雨事件在發(fā)展過程中與氣象臺發(fā)布暴雨預(yù)警在時間和空間上基本吻合.

2.3 與雷達(dá)數(shù)據(jù)集的對比分析

圖5 將本文方法提取的暴雨事件圖像與雷達(dá)反射率圖像進(jìn)行對比，其中上側(cè)4 幅圖為本方法提取出的暴雨事件，下側(cè)4 幅圖為對應(yīng)時刻雷達(dá)反射率圖像.可以看出本文方法提取出的暴雨事件與雷達(dá)圖像中基本反射率較高的區(qū)域基本吻合，故本文方法提取的暴雨事件位置與實(shí)際暴雨發(fā)生位置基本吻合.

3 結(jié) 論

針對長時間序列、高時空分辨率的降水量柵格數(shù)據(jù)集，本文從時間維度、空間維度和時空追蹤這3 個層次設(shè)計(jì)了一種暴雨事件提取與追蹤方法.時間維度的暴雨提取保證了暴雨在時間上的連續(xù)性，空間維度的暴雨提取保證了暴雨在空間維度的相鄰性，暴雨時空追蹤保證了暴雨在時空上的連續(xù)性和一致性.通過實(shí)驗(yàn)例證與傳統(tǒng)提取方法的對比分析，驗(yàn)證了本文提取方法REETM 的正確性和可行性，并得出了以下3 個主要結(jié)論：

圖4 2016-07-18～2016-07-21 某暴雨事件Figure 4 A rainstorm event from July 18 to 21, 2016

圖5 暴雨事件圖像與雷達(dá)反射率圖像Figure 5 Rainstorm events images and radar reflectivity images

1）本文方法從時間、空間、時空3 個維度上進(jìn)行設(shè)計(jì)，實(shí)現(xiàn)了暴雨從產(chǎn)生發(fā)展到消亡事件的提取與追蹤；

2）相較于日累加暴雨提取方法，本文方法的探測率平均提高12.97%；

3）依據(jù)暴雨預(yù)警信息和雷達(dá)圖像驗(yàn)證了本文方法提取暴雨事件的有效性和可行性.

本文提取的暴雨事件清楚地體現(xiàn)了暴雨在整個發(fā)展過程中的空間分布以及動態(tài)變化過程，對暴雨的分析及預(yù)測有著重要的意義.比如，暴雨在某個時刻發(fā)生了分裂，可能由大氣氣流突然變化造成.然而，對于暴雨對象運(yùn)動過快以及前后時刻空間覆蓋范圍沒有重疊的情況，采用本文方法提取的準(zhǔn)確度會受到影響.因此，今后的研究重點(diǎn)在于結(jié)合暴雨運(yùn)動路徑對暴雨位置進(jìn)行預(yù)測，以完善優(yōu)化算法.