淮河流域暴雨危險度時空特征分析

吉中會 ，張蝶，于小兵

(1.南京信息工程大學管理工程學院，南京210044；
2.南京信息工程大學氣象災害預報預警與評估協同創新中心，南京210044；
3.南京信息工程大學風險治理與應急決策研究院，南京210044)

在氣候變暖的大背景下，全球暴雨的時空格局正在發生重大變化。中國飽受洪澇災害侵擾，各流域暴雨強度和頻次變化趨勢并不一致(胡畔等，2021；
謝五三和王勝，2010)?；春恿饔蛭挥谖覈鴸|部(111°55′—121°20′E，30°55′—36°20′N)，地跨河南、安徽、江蘇及山東四省，屬于東亞季風氣候區，是中國南北氣候的過渡地帶。根據氣候歷史資料記載，淮河流域在

1470—1991 年期間發生一般等級洪澇災害的年次數占總年數的30.8%，平均約3 a一次；
嚴重(包括極嚴重)的洪澇災害發生頻次為14.9%，約7 a 發生一次(溫克剛和丁一匯，2008)。由于上游落差大，中游平緩、支流水系復雜眾多，以及歷史上黃河長期奪淮入海等原因，淮河流域極易洪泛(劉倪和葉金印，2014)。因此，研究淮河流域的暴雨危險度，對于有效防治暴雨洪澇、降低災害風險，切實保障人民群眾生命財產安全和社會經濟持續健康發展具有重大意義。

暴雨危險度是對暴雨致災因子危險性程度的客觀量化，國內外學者在區域暴雨危險度研究方法方面已經取得了一定的成果。在單因素量化分析方面，降水指標通常是對其趨勢性和周期性等的分析(董旭光等，2014；
孫燕和朱偉軍，2013；
閔愛榮等，2016)。趨勢分析主要采用時間序列分析方法，如趨勢擬合、平滑法等(黃埡飛等，2021；
Shao et al.，2017)；
周期分析則主要采用各類譜分析法(魏鳳英，2007)，比較經典的有傅里葉變換(劉瑞鑫，2019)、小波分析(潘雅婧等，2012；
劉毅等，2020)等。傅里葉變換是一種重要的時間序列數據處理方法，經過處理的數據中包含各種頻率成分，但是丟失了數據的時間信息。小波變換可以將數據的頻域信息和時域信息結合起來，克服了傅里葉變換的不足，但是由于小波基底選擇不同，會對結果產生顯著影響(郭建平和趙立龍，2021)。經驗模態分解(Empirical Mode Decomposition，EMD)作為一種復雜數據分析方法，最早由Huang等(1998)提出，該方法既可以處理非線性、非平穩數據，也可以處理線性的、平穩的數據，具有良好的自適應屬性，能夠比其他方法更好地反映時間序列數據的物理意義(楊永鋒和吳亞峰，2013)。因此，EMD 方法在暴雨危險度的單因素分析中具有較好的應用潛力。在多因素綜合量化分析方面，主要通過多指標加權綜合方法實現對暴雨危險度的綜合表達。例如，王晶等(2019)采用層次分析法考慮了降水強度、最多持續降水量、最長持續降水日、高程及河網密度，評估了淮河流域中上游地區洪澇災害的致災因子危險性；
嚴志濤等(2021)采用加權求和方法考慮降水強度、最大持續降水量和最長持續降水天數指標，評估了西南四省的暴雨洪澇危險性。這些方法主要是將不同的影響因子進行相加或相乘，來代表總的危險性程度，存在明顯不足。Copulas函數在刻畫多變量之間相依結構方面具有明顯優勢，給暴雨危險度的多指標綜合提供了新的思路。Copulas 函數最早于1959年被提出，受計算能力限制，直到20世紀90年代末才被廣泛地應用在金融領域，它以著名的Sklar定理為基礎(Nelsen，2007；
Rodriguez，2007)，目前在氣象災害領域的應用逐步興起，例如在沙塵暴、水文災害等領域中的應用(李寧等，2013；
Tosunoglu and Can，2016)。

本文利用淮河流域及其周邊站點降水數據，嘗試采用EMD 方法來分析反映暴雨危險度的單指標時間序列數據，再選擇Copulas 函數方法分析淮河流域暴雨的綜合危險度，探討淮河流域暴雨危險度的時空變化特征。

1.1 資料來源

考慮到暴雨危險度分析的空間完整性，本研究區涵蓋了整個淮河流域及其周邊的部分區域。暴雨災害風險的主要致災因子為降水，參考相關文獻(王晶等，2019；
嚴志濤等，2021；
鄭艷等，2018；
張正濤等，2014)，選擇淮河流域及其周邊站點的歷年年暴雨日數(日降水量≥50 mm的日數)作為頻數指標，各站點歷年的年最大日降水量作為強度指標來分析淮河流域暴雨危險度的時空規律。數據來源于中國氣象局國家氣象信息中心提供的1951—2017 年淮河流域的降水日值數據集。通過對數據缺測值、異常值、均一性等方面的檢驗，最終選取淮河流域及其周邊共68個站點(圖1)的逐日降水數據進行分析。其中部分站觀測數據起始年份晚于1951年(所分析的時間序列最少不短于57 a，即最晚不遲于1960 年)，本文基于各站的觀測數據對淮河流域暴雨危險度進行時空特征分析。

圖1 研究區概況及氣象站點分布Fig.1 Study area and meteorological stations

1.2 研究方法

1.2.1 暴雨危險度的周期及趨勢分析

EMD 方法的實質是將數據中存在的不同尺度或周期下的波動或變化趨勢由高頻向低頻逐級提取出來，產生一系列不同尺度的分量，每個分量被稱為一個特征模態函數(Intrinsic Mode Function，IMF)(Huang et al.，1998)。其中，IMF需要滿足以下兩個條件：

(1)極值的個數和過零點的個數必須相等或相差一個；

(2)在任何一點上，由局部極大值定義的上包絡線和由局部極小值定義的下包絡線的平均值為零。

EMD方法具體分解步驟為：

(1)對于任意給定的時間序列數據x(t)，求解其所有的極值點(極大值和極小值)。用三次樣條函數連接所有極大值點形成上包絡線u(t)，連接所有極小值點形成下包絡線b(t)；

(2)計算上下包絡線的均值m(t)，然后用原始數據x(t)減去m(t)得到新的序列h(t)，計算方法為

(3) 將h(t)作為新的數據序列，重復上述兩個步驟，直到h(t)滿足IMF 的兩個條件，就得到了第一個IMF，記為IMF1，其包含了原始數據信息的最高頻率成分；

(4)重復上述步驟，直到殘余序列rn(t)是一個單調函數，即不能再分解出IMF 分量為止，這時原始序列x(t)可以用所有IMF(IMF)與殘余序列res.(r)之和來表達，其公式為式(3)

其中，j為IMF的序號；
n為IMF的個數；
rn(t)為殘余序列，代表原始數據中的平均趨勢(楊永鋒和吳亞峰，2013)。

綜上可以發現，EMD方法可以將由不同時間尺度或周期模式的原始時間序列數據進行分解，用簡單、有規律的IMF分量來表達，以便分析判斷原始序列數據的時間演變規律。

1.2.2 暴雨危險度的聯合概率及重現期分析

Copulas 函數在多因子聯合概率分析中具有顯著優勢，它不需要假設原始變量服從一致的邊緣分布，函數形式靈活多樣，建立的聯合函數能更接近實際情況；
能夠捕捉變量之間的相依結構，所構建的函數不受變量偏態分布和異常值影響；
模型描述的相依結構是一個完整分布，而不僅僅是平均狀態，可以揭示變量在不同分位數上相依性強度的變化(Ji and Wan，2021)。由于作為暴雨危險度致災因子的暴雨日數和年最大日降雨量兩指標之間可能存在線性相關關系，也有可能存在非線性或非對稱關系，且各變量有可能服從不同的邊緣分布，因此，嘗試采用Copulas函數方法對研究區所選站點兩指標數據進行聯合分析，通過描述兩變量之間的相依關系，構建相應的聯合分布函數，綜合客觀地探討兩變量共同影響下的暴雨危險度時空規律。

Copulas函數構建聯合分布的主要步驟為：

(1)確定所選的時間序列數據指標的邊緣分布及相應的參數；

(2) 度量所選指標變量之間的相依性(一般采用Pearson 相關系數(r)、Kendall 秩相關系數(τ)或Spearman 相關系數(ρ)進行計算)；

(3)選擇不同類型的Copulas函數并進行參數估計；

(4)通過擬合優度檢驗，選擇最優的Copula 函數，建立聯合概率分布。

Copula函數類型豐富，最常見的有橢圓簇(tCopula、Gaussian Copula等)和阿基米德簇(Frank、Clayton 及Gumbel Copula 等)(Nelsen，2007)。本文考慮這五種常見的Copulas函數類型，邊緣分布函數的參數估計采用最大似然法(Karmakar and Simonovic，2009)，通過Kendall 秩相關系數τ描述變量之間的相依性(Klein et al.，2011)，結合離差平方RMSE 和一般最小二乘法(Ordinary Least Squares，OLS)最小準則進行擬合優度檢驗，來選擇最優的Copula函數類型(劉雪琴等，2012)，最小OLS計算經驗與理論概率之間的距離(SOLS)的公式為

其中，n為樣本量，PE和PT分別代表計算的經驗概率和理論概率。

聯合概率反映的是兩變量共同作用下暴雨危險度的概率分布，五種Copulas 函數模型的計算公式可參考文獻(謝中華，2015)。由于最優Copula 函數計算所得的聯合概率是一個完整的概率分布，為進一步分析各站點聯合概率大小即危險性程度，可以通過統計概率分布的特征值，即均值、分位數、極值(最大值/最小值)等指標來反映區域暴雨的危險度。

暴雨聯合重現期是指每出現一次暴雨的平均時間間隔，是災害長期預警預測的重要參考指標(姚莉等，2010)。多變量的聯合重現期包括同現重現期(and)和非同現重現期(or)兩種類型，其公式分別為式(5)和式(6)。

其中，T(x,y)為聯合重現期，L為暴雨發生的平均間隔時間，P代表聯合概率，FX(x)和FY(y)為兩個單變量的邊緣分布函數，C為Copula聯合概率分布。

在所構建的最優Copula函數模型的基礎之上，計算淮河流域整體站點平均的聯合概率及重現期，以及各站點的暴雨綜合危險度，并分析其時空變化特征。

2.1 淮河流域暴雨危險度的時間變化特征

2.1.1 單指標趨勢和周期的時間序列分析

將所選站點的暴雨頻率和強度兩指標數據進行站點平均，獲取反映淮河流域暴雨危險度單指標的時間序列，通過EMD 方法進行分析，結果如圖2 所示。從圖2 可見，暴雨頻率指標分成了3 個分量和殘余序列，而暴雨強度指標則被分解成4 個分量和殘余序列。兩指標的IMF1 與原始序列的相符程度最高，包含了最高的頻率成分，但是周期和趨勢情況均不易直接判讀；
IMF2振幅隨著時間推移均呈現逐漸變小的趨勢，且周期波動規律也越來越弱。殘余序列(res.)反映的是原始序列的平均趨勢，可以看出暴雨頻率在1980年左右出現明顯轉折，由原來顯著下降轉變為顯著上升趨勢；
暴雨強度則一直呈現出顯著的上升趨勢。通過進一步對1980—2017 年的殘余序列進行線性擬合，得到兩指標的平均趨勢函數分別為y=0.009 6x+2.221和y=0.078 2x+86.266，判定系數(R)高達0.977 和0.996，說明其線性上升趨勢非常顯著。因此，通過對兩個指標EMD結果的分析發現：淮河流域的暴雨頻率的周期波動規律逐步減弱，說明流域降水系統的穩定性在下降，是可能出現極端暴雨的重要特征；
暴雨頻率和強度在最近40 a(1980年之后)均呈顯著的線性上升趨勢，說明淮河流域的暴雨頻次和強度均在不斷增強，特大暴雨增多。

圖2 年暴雨日數(a1—a5,單位:d)和日最大降雨量(b1—b6,單位:mm)的EMD結果Fig.2 EMD results of the indexes of(a1-a5)the annual rainstorm days(unit:d)and(b1-b6)the maximum daily rainfall(unit:mm)

根據各模態極值點的位置和數量求算平均周期(即求算相鄰極值點之間的平均時間長度)，將由極大值點求得的平均周期與由極小值點求得的平均周期再進行平均，作為各分量的平均周期。求得暴雨頻率序列IMF1—IMF3的平均周期分別為2.6 a，6.7 a，14.3 a，暴雨強度序列IMF1—IMF4 的平均周期分別為2.4 a，5.5 a，9.1 a，27 a，兩者對應關系較好，尤其是前兩個分量序列的周期，分別約在2 a和6 a，這與已有研究中ENSO 在年際變化上的信號具有一定相關性(于文金等，2016)。各分量的方差貢獻率可反映其周期波動情況(表1)，兩指標IMF1 的方差貢獻率分別為73.51%和63.48%，說明淮河流域的暴雨頻率和強度的主周期均為準2 a；
IMF1與IMF2的方差累計貢獻率分別達到97.34%和82.62%，說明準6 a也是淮河流域暴雨的主要周期。

表1 年暴雨日數和日最大降雨量各分量對應的方差貢獻率(單位:%)Table 1 Variance contribution rate(unit:%)corresponding to each component(IMF)of the annual rainstorm days and the maximum daily rainfall.

2.1.2 聯合概率及重現期的時間序列分析

將站點平均的頻率和強度指標進行Copula 聯合概率建模，構建反映淮河流域暴雨危險度的Copula模型。采用Normal、Lognormal、Gamma、Weibull、Exponential、Extreme Value、Generalized Extreme Value、Logistic 等常見的概率分布類型對兩變量進行邊緣分布擬合，通過對頻率分布直方圖的直觀判斷，以及Anderson-Darling(AD)檢驗，兩變量的最優邊緣分布均為一般極值分布(Generalized Extreme Value)。在此基礎上，計算Kendall 秩相關系數τ為0.602 7，通過了0.01的顯著性水平檢驗，說明兩變量之間存在較高的正向關系，并通過Kendallτ與Copula函數的參數之間的關系(Nelsen，2007)，構建五種Copulas 函數模型，結合RMSE 和OLS 準則進行擬合優度檢驗，從中選擇出最優的函數類型為Gaussian Copula(表2)，相應的聯合累積概率分布如圖3所示。

表2 五種Copula函數模型的擬合優度檢驗結果Table 2 Test results of fit goodness for five Copula function models.

圖3 暴雨危險性聯合累積概率分布圖Fig.3 Joint cumulative probability distribution of the torrential rain hazard

從圖3可以判斷出不同頻率和強度組合下的暴雨危險度的聯合累積概率。從年暴雨日數指標來看，全流域平均3 d及以內的暴雨頻次占總暴雨頻次的90%以上，超過3 d的以上的暴雨頻次很少，并隨著天數的增加暴雨頻次增加緩慢；
從年最大日降水量指標來看，最大日降水量小于100 mm 的情況占總的暴雨強度的90%以上，且超過100 mm后，隨著最大日降水量的增加，暴雨強度增加顯著變緩。

根據聯合概率分布函數，進一步求算淮河流域暴雨的聯合重現期(式(5)和式(6))，結果如圖4 所示?？梢钥闯?，圖4a 中暴雨的同現重現期(and)，隨著兩指標值的逐步增加而增大，且增大速率變快(等值線密度增大)，說明隨著暴雨頻次和強度的同時增加，暴雨的危險度快速增大，重現期高值出現的頻率增加，即極端暴雨事件顯著增多。圖4b 暴雨的非同現重現期(or)也表現出與同現重現期相同的趨勢，即隨著兩指標的增加，危險度增加速度加快，但是整體重現期值較小，說明兩因素非同時作用下暴雨危險度相對較低。

圖4 暴雨危險性的聯合同現重現期“and”(a)與聯合非同現重現期“or”(b)Fig.4 The joint return periods of the heavy rain hazard for(a)“and”and (b)“or”scenario

2.2 淮河流域暴雨危險度的空間變化特征

2.2.1 單指標平均趨勢的空間變化特征分析

分別對流域各站點的兩指標數據進行經驗模態分解，并對反映平均趨勢的殘差項進行線性擬合，通過F檢驗判定線性趨勢的顯著性，結果顯示：68 個站點中頻率指標和強度指標分別有52和53個站點的線性平均趨勢通過了α=0.05 顯著性檢驗。將通過顯著性水平檢驗的站點的線性趨勢的回歸系數通過Arc-Map 中的含障礙的核插值(Kernel Interpolation with Barriers，KIB)方法展布到空間上(圖5)，已有研究表明KIB方法在水文和氣象要素的空間插值應用中具有良好效果(Yan and Sun，2009)。從圖5 中可見，淮河流域暴雨兩個指標平均變化趨勢的高值空間格局表現基本一致，上升趨勢最大的地區主要分布在淮河流域的西南部分，其中暴雨頻次上升趨勢最大的區域為河南信陽和安徽六安，暴雨強度上升趨勢最大的為安徽六安。兩指標對應的上升和下降趨勢的區域面積比例相當，河南信陽、駐馬店，安徽六安、阜陽、淮南、天長、滁州、蚌埠，江蘇揚州、興化、淮陰、宿遷、徐州，以及山東棗莊和臨沂南部等地區的暴雨頻次均呈增加趨勢，暴雨頻次下降趨勢主要分布在淮河流域北部的河南鄭州、開封，山東菏澤、濟寧北部、日照，以及江蘇鹽城的東北部。暴雨強度下降趨勢最大的區域分布在河南平頂山、鄭州、開封、許昌，及江蘇鹽城等。流域內上升與下降趨勢明顯區域之間過渡地帶的變化趨勢平緩，在0值線附近波動。

圖5 暴雨頻率(a)和強度(b)平均變化趨勢的空間分布Fig.5 Spatial distribution of the average change trends of the heavy rain(a)frequency and(b)intensity

2.2.2 綜合危險度的空間變化特征分析

對流域各站點分別求算兩指標聯合概率，根據RMSE 和OLS 最小準則進行擬合優度檢驗，選取最優Copula 函數模型(圖6a)并進行危險性程度的概率估計。由于聯合概率計算的是各站點的整個概率分布，其特征值中的均值能夠反映暴雨危險性程度的平均狀態，而極值可以反映暴雨的極端狀態，因此，在計算的聯合概率的基礎上，統計了各站點相應的特征值：均值、最小值與最大值(圖6b—d)作為暴雨綜合危險度的表征，并采用KIB 方法展布到整個流域。圖6a 可見，流域大多數站點的暴雨危險度聯合概率模型的最優Copula類型為Clayton，其次為Frank，第三為tcopula，其他兩種類型均沒有選入。圖6b—d為聯合概率的空間分布情況，可以看出三個特征值與單指標平均變化趨勢的高低值(圖5)的空間分布存在顯著差異。聯合概率均值和最小值的空間格局基本一致，對應高值出現在流域的北部和西北地區，最低值出現在流域的東南部，并呈現由東南向西北方向逐步升高的趨勢；
聯合概率最大值分布則比較分散，整個流域的南部及東北部均出現了最大值的高值。結合單指標平均變化趨勢的空間分析，可以發現：淮河流域西北區域的部分城市，如鄭州、開封、平頂山北部等城市雖然暴雨頻次和強度呈現下降趨勢，但是聯合概率估計的暴雨綜合危險度均值和最小值卻很高。此外，淮河流域中南部的河南開封、商丘，安徽淮北、亳州、淮南、六安一帶是聯合概率最大值的高值所在，暴雨綜合危險度極高。

圖6 各站點最優Copula類型(a)及危險性聯合概率均值(b)、最大值(c)和最小值(d)空間分布Fig.6(a)Optimal Copula type for each station,and the spatial distributions of the joint probability for(b)mean, (c)maximum and(d)minimum of rainstorm hazard

本文選擇年暴雨日數和年最大日降水量作為暴雨危險度的頻次和強度指標，采用EMD方法分析淮河流域單指標時間序列數據的平均變化趨勢及周期，采用Copulas 函數方法估算暴雨的綜合危險度及重現期。通過站點平均數據，分析流域暴雨危險度在單指標和兩指標聯合情況下的時間序列特征；
再通過對各站點的數據分別計算，分析流域暴雨危險度單指標的平均變化趨勢和綜合危險度的空間分布規律。主要結論如下：

(1)流域整體暴雨頻次呈現先下降后上升的趨勢，轉折點出現在1980年左右；
暴雨強度則一直呈現顯著的線性上升趨勢。暴雨頻率和強度的主周期均為準2 a周期，6 a也是淮河流域暴雨的主要周期。

(2)流域整體暴雨危險度的最優聯合概率模型為Gaussian Copula。全流域90%以上的暴雨頻次在3 d及以內且暴雨強度小于100 mm。超過該閾值后，暴雨頻率和強度的增加顯著變緩。暴雨災害的同現重現期值和非同現重現期，均隨兩單指標值的增加而增大，且增大速率逐漸變快，但是同現重現期高值較高，出現極端暴雨的概率更大；
而非同現重現期值相對較小，說明兩因素非同時作用下的暴雨危險度相對較小。

(3)暴雨頻率和強度指標的平均變化趨勢的高值空間格局表現基本一致，上升趨勢最大的地區主要分布淮河流域的西南部分。兩個指標上升和下降趨勢的區域面積比相當，下降趨勢主要分布在淮河流域北部。暴雨強度下降趨勢最大的區域分布在流域的西北部，以及江蘇鹽城的沿海地區。

(4)流域大多數站點的暴雨危險度聯合概率模型的最優類型為Clayton Copula，綜合危險度的三個特征值與單指標平均變化趨勢高低值的空間分布存在顯著差異。危險度的均值和最小值的空間格局基本一致，對應高值出現在流域的西北地區，最低值出現在流域的東南部，并呈現由東南向西北方向逐步升高的趨勢；
聯合概率的最大值分布則比較分散，高值出現在流域的中南部地區。

流域西北的部分城市地區雖然暴雨頻次和強度呈現下降趨勢，但是暴雨綜合危險度均值和最小值卻很高，政府部門在作防災減災決策時需要充分考慮暴雨的綜合危險度。流域中南部地區為暴雨危險度最大值的高值所在，氣象部門在暴雨災害的預警工作中尤其需要注意這些地區可能發生的極端暴雨事件。

本文嘗試采用EMD和Copulas函數方法對暴雨危險度的單指標和綜合指標開展時空特征分析。降水作為暴雨災害風險評估中的致災因子危險度的主要指標，除了暴雨日數和最大日降水量之外，在以后的研究中還可以增加考慮暴雨持續時間、暴雨事件頻數等相關危險性指標。此外，Copula函數還可以拓展到三維甚至是高維，在未來研究中可以進一步挖掘暴雨危險性指標之間的關系，考慮更多類型的Copula 函數，以期獲得更為客觀準確的區域暴雨危險度的時空分布規律。