多變量洪水非平穩頻率分析方法研究
——以淮河流域蚌埠水文站為例

徐鵬程，劉楠楠，李 帆，仇建春，張昌盛，戴 笠

（1. 揚州大學水利科學與工程學院，江蘇 揚州 225009； 2. 江蘇省水利建設工程有限公司，江蘇 揚州 225009；3. 泰州市港航事業發展中心，江蘇 泰州 225300）

0 引 言

氣候變化和人類活動加劇將改變全球水文循環的現狀，對降水、蒸發、徑流、土壤濕度等關鍵水文氣象要素造成直接影響，導致流域性大洪水的發生風險增加，而這些正成為制約國民經濟的主要瓶頸［1-3］。現行的洪水頻率分析是建立在平穩性假設前提下的，這意味著影響洪水特征變量的環境控制因素諸如氣候要素、土地利用率在過去、現在和將來都遵循一成不變的物理機制［4］。而在氣候變化、城市化等因素的多重影響下，基于平穩性假設的洪水頻率分析方法的適用性受到了挑戰［5，6］。為了充分保證洪水設計值的可靠性，探索非平穩性條件下的洪水頻率分析方法以及工程設計值的推求顯得至關重要。

目前，非平穩水文頻率分析方法分為單變量非平穩和多變量非平穩分析方法。其中，單變量非平穩性已經逐步納入到洪水極值事件［7］、干旱極值事件［8］和極端降雨事件［9，10］頻率分析中。Rigby 和Stasinopoulos［11］提出了考慮位置、尺度、形狀參數的廣義可加模型（Generalised Additive Models for Location Scale and Shape， GAMLSS）， 其基于參數回歸模型定量解析了極值序列的統計參數與物理解釋變量間線性/非線性關系。GAMLSS模型涵蓋了多種隨機變量分布函數類型，克服了傳統的分布函數對于特定分布類型（超峰度、高度正偏/負偏和平頂峰度類型的分布）擬合過度問題［12］。Chen 等［13］基于臺灣地區降雨極值序列對比分析了單變量情景下3 種非平穩頻率（非參數離散小波變換（Discrete Wavelet Transform， DWT）、經驗模態分解法（Ensemble Empirical Mode Decomposition， EEMD）以及參數加權最小二乘法）分析方法的適用性。

傳統的洪水頻率分析主要聚焦于洪峰量級的分析，由于洪水事件（過程）并不能憑借單個特征變量（洪峰）得以完全概化，一次洪水過程往往包含洪峰、洪量、歷時等要素信息，是一個包含頻域、時域和空間域的復雜水文過程，僅從單變量角度無法有效刻畫峰高量大的流域型洪水風險。為此，近年來，一些水文研究也逐步關注到多變量水文序列的非平穩性問題［14］:Sarhadi 等［15］提出了一種貝葉斯動態Copula 模型有效模擬了多維水文氣象極端事件的混合連續與離散的多變量時變相依結構。Jiang 等［16］采用動態VineCopula 函數有效建立了西江流域的年最大日徑流-年最大3 日洪量-最大7 日-最大15 日洪量的四維聯合分布模型，為非平穩條件下多元水文極值設計提供了理論依據。Xu 等［17］基于動態Copula 函數量化分析了不同物理解釋變量對洪水設計值的影響，并提出了變化環境下多變量洪水風險的短期預估的一般框架。

淮河流域人口密集，人水爭地矛盾突出。因人類活動（城市化）、氣候變化引發的洪水極端事件［7］，極端降雨事件［18］都呈現一定的非平穩性， 但是目前對于淮河流域洪水頻率的非平穩研究主要聚焦于單變量角度開展的。過去幾十年以來，城市化對流域洪水特性的影響逐步成為氣候變化領域的研究熱點［19，20］。由于多種因素的相互作用，高度城市化地區的洪峰流量增幅較為明顯［21-23］。城市化地區不透水面積的增加很大程度上弱化了自然滲透量，導致暴雨徑流量顯著增加；同時，城市規劃設計層面考慮高效的人工排水溝替代天然河道，減少了徑流響應的滯后時間。為此，如何從非平穩性入手，定量解析城市化和氣候變化相關的驅動因子與多變量洪水風險的響應關系，是目前非平穩頻率分析領域亟待解決的關鍵科學問題。為此本文基于淮河流域蚌埠（吳家渡）水文站的洪峰（P）洪量（V）序列為研究對象，構建非平穩條件下P、V相依性結構的二維動態Copula 模型，定量分析氣候變化因子（極端降雨因子、長尺度氣候濤動因子）、城市化因子（不透水面積率（Percent of Impervious Surface Area， PISA））對單變量、多變量洪水風險的影響程度。

1 基于動態Copula 函數的非平穩風險模型構建

由于不受限于邊緣分布的形式，Copula 函數是一種備受推崇的水文多變量分析模型，主要用于擬合隨機（極值）變量之間的相依性結構。針對氣候變化和人類活動導致的單變量、多變量水文極值序列的非平穩性，本文擬構建如下的動態Copula 函數模型:假設［X，Y］為洪水-洪量的極值序列對， 基于時變Copula和邊緣分布的多變量聯合概率分布模型可以如下公式求得。

以上所有這些參數和時間變量或者物理解釋變量呈現線性或非線性函數關系；物理協變量因子集合主要包括氣候濤動因子（NAO、SOI、NINO3.4）、極值降雨量和城市化率；u和v是邊緣分布的累積概率值。在本研究中，非平穩多變量風險評估模型分為2 個個階段:①采用邊緣分布和Copula 為主的時變矩多變量模型，并結合對數似然比法（LR）和最小化AICc（修正化的赤池信息）準則，從分布參數角度定量識別洪水極值序列中潛在非平穩性；②基于不同非平穩分布模式下特定分位數的風險值的計算定量剝離城市化、氣候變化因素對多變量洪水風險的影響。

1.1 非平穩邊緣分布模型構建

本節討論5種常見的極值分布函數為原型如何構建以物理相關因子作為分布參數的解釋變量非平穩模型。表1 顯示了5種常用極值分布的概率密度函數（PDF）形式。由于篇幅限制，本處著重以GEV分布的位置參數（μ）為例展示非平穩條件下的時變分布模型的形式:

表1 5種潛在的邊緣分布的概率密度函數Tab.1 The probability density function （PDF） for the five potential marginal distribution

式中:Covk(t)為表示第k個潛在的時變物理變量序列｛大尺度氣候濤動因子（NAO， SOI， NINO3.4），氣象因子［極值面降雨量（Maxpre）］，城市化因子（PISA）｝。

式中采用了線性和非線性趨勢型非平穩特征。當然多種因子的復合影響，可以增加公式（2）中的協變量因子。考慮到不同分布對應的參數數量存在差異，為此對于三參數分布（GEV 分布和PIII 分布），潛在非平穩邊緣分布模型種類明顯多于兩參數分布的潛在非平穩模型種類。

1.2 非平穩Copula模型構建

為了模擬非平穩條件下洪峰和洪量之間的時變相依性結構，本節討論物理解釋變量作為參數協變量的動態Copula 的建立過程。在多變量水文頻率分析中，可以考慮如下多種潛在的二元Copula 函數作為動態Copula 模型的母函數:5 種單參數Copula（Joe、Frank、Gumbel、Clayton 和Gaussian）和5 種雙參數Copula［Clayton-Gumbel（BB1）、t、Joe-Gumbel（BB6）、Joe-Clayton（BB7）和Joe-Frank（BB8）］。以上Copula 模型的計算擬合可以通過R語言環境下CDVine包［24］中實現。

類似于非平穩邊緣分布模型的構建過程，非平穩Copula 模型同樣是通過構建參數與物理解釋變量間的線性函數關系，可以建立如下公式:

式中:Cov1(t)，Cov2(t)，…，Covk(t)為k個潛在的時變物理變量序列。

術前常規使用0.5%左氧氟沙星滴眼液每日4次，連用3 d，術后使用0.3%妥布霉素地塞米松滴眼液和0.1%玻璃酸鈉滴眼液每日4次，1周后停用妥布霉素地塞米松滴眼液，0.1%玻璃酸鈉滴眼液連續使用1個月。

1.3 基于LR-AICc準則的最優非平穩模型的優選

一旦構建了多種潛在的非平穩邊緣分布和Copula 函數模型，接下來需要綜合運用極大似然比準則（LR）和AICc準則獲得最優的非平穩/平穩模型。 假設ST0代表平穩模型，ST1代表特定的非平穩時變模型；同時LL0和LL1為平穩和非平穩模型的對數似然值，其對應的對數似然比值可以按如下公式求得:

考慮到可能存在多個非平穩性模型能夠通過以上的LR準則檢驗，需要計算各個非平穩模型的修正化赤池信息量［25］（AICc）:

式中:k代表非平穩模型對應的參數個數；n代表研究序列的樣本長度。AICc準則運用原則為:越小的AICc值代表了模型的擬合效果更好。

2 非平穩條件下單變量/多變量風險的計算

為了定量剝離城市化因子（PISA），氣候變化因子（Maxpre，NAO，SOI，NINO3.4）對洪水風險的影響，因此可以通過重現期得出單變量和多變量風險。

對于單變量情形，每個洪水極值（P或V）的單變量重現期為:

因而單變量風險:

式中:FX(x0)為邊緣分布的累積概率值（CDF）；RP為單變量重現期；UR為單變量風險（累積超越概率）。

AND情景［7］下的雙變量聯合重現期（JRP）可計算如下:

式中:FX(x0)，FY(y0)是洪水變量對應的邊緣分布累積概率值；JRP代表聯合重現期；BR是多變量風險（累積聯合超越概率）。

3 研究區域與數據

淮河流域位于東經111°55'～121°25'、北緯30°55'～36°36'，東西長約700 km，南北寬約400 km，面積27 萬km2。流域地跨湖北、河南、安徽、江蘇、山東五省，覆蓋40 個市、163 個縣（市、區）的全部或部分區域。淮河流域地貌具有類型復雜多樣、層次分明、平原地貌類型豐富的特點。流域西部、南部及東北部為山區和丘陵區，其余為平原、湖泊和洼地。淮河流域的植被具有明顯的過渡性特點。流域北部的植被屬暖溫帶落葉闊葉林與針葉松林混交林；中部低山丘陵區屬亞熱帶落葉闊葉林與常綠闊葉林混交林；南部山區為常綠闊葉林、落葉闊葉林與針葉松林混交林，并夾有竹林，山區腹部有部分原始森林。淮河流域地處我國南北氣候過渡地帶，氣候變化復雜，降雨時空分布不均。流域內眾多支流多為扇形網狀水系結構，洪水集流迅速。由于淮河流域位于多重過渡地帶的孕災區，特定的地理環境和下墊面條件、復雜多變的氣候，加之黃河奪淮和人類活動的不利影響，導致流域水、旱、風暴潮等自然災害的頻繁發生。淮河流域的洪水呈現高洪峰-高水位-長歷時的特征，易于發生全流域性大洪水。本研究選取淮河流域中游蚌埠（吳家渡）水文站（圖1）作為研究對象，基于站點1959/1/1-2016/12/31 的日徑流序列提取洪峰-洪量序列并構建動態Copula 函數風險分析模型。由于本文采用不透水面積率（PISA）作為物理協變量納入到非平穩模型的參數中，所使用的PISA（城市化率）序列可以從如下的網址鏈接中下載:http://data.ess.tsinghua.edu.cn/。PISA的逐年演變圖見圖2。

圖1 研究區域水文氣象站點分布圖Fig.1 Map of location of hydrometeorological stations in the study area

圖2 淮河中游子流域（蚌埠水文站所控制的子流域）的城市化范圍演變圖Fig.2 The evolution map of urbanization scope of Huaihe River middle reaches sub-river basin （controlled by Bengbu hydrological station）

多個長尺度氣候濤動因子序列［南方濤動因子（SOI），北大西洋濤動因子（NAO）和厄爾尼諾因子（NINO3.4）］由如下鏈接下載:https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries。同時基于蚌埠水文站控制的子流域邊界內的氣象站1959-2016年的日降雨序列提取子流域內的極值面雨量（Maxpre）序列作為非平穩模型潛在的氣候變化因子。

4 結果分析

4.1 非平穩優化結果

首先依據LR-AICc準則對潛在平穩/非平穩模型進行優化選擇，具體模型優化結果見表2。對于蚌埠水文站提取的洪峰序列，非平穩模型GEV3 其對應的LR檢驗的p值為0.003 2 顯著小于0.05 的閾值，表明洪峰序列的趨勢型非平穩性特征在5%顯著水平上能夠被檢測到；盡管GEV1-GEV2 模型的p值可以通過LR檢驗，但是考慮到GEV3模型的AICc值相對其他兩個非平穩模型更小，所以GEV3 是洪峰序列對應的最優邊緣分布模型，其中位置參數和面極值雨量（Maxpre）呈現二次型統計學關系（非線性），尺度參數與NAO 因子呈現線性關系，形狀參數和城市化率（PISA）呈現線性關系。按照同樣的方法，模型GA3為洪量序列的最優非平穩模型，而由于AIC值最小化準則選定平穩Frank Copula（FR0）模型為最優的Copula模型。

表2 最優的非平穩邊緣分布和Copula模型優選結果Tab.2 Results of best-fitted marginal distribution and Copula models

4.2 LR-AICc 準則對于趨勢型非平穩特征檢驗的有效性

分別對蚌埠水文站的洪峰、洪量進行非參數的單變量、多變量Mann-Kendall（M-K）趨勢型檢驗［26］，其檢驗結果為:ZQ=1.156，ZV= 0.689。由M-K 檢驗結果可知洪峰，洪量序列在5%顯著水平上趨勢不顯著。由于M-K 檢驗側重于變量均值的潛在的單調性趨勢［27］，對于某一分布（如GEV 分布），將存在這樣的特定情形:位置和尺度參數的反向趨勢會使得M-K檢驗結果為不顯著（|Z|＜1.96）， 這在一定程度上體現了本文采用的對數似然比（LR）檢測方法的必要性。為了探究引發上述非參數檢測法和LR檢測法結果的不一致性，圖3 繪制洪峰、洪量分布參數的隨時間的演變關系。由于GEV 分布的尺度、形狀參數（斜率分別為-0.196 7和-0.001 4）呈現下降趨勢性而其位置參數呈現上升趨勢性（斜率為6.474 1），由于分布參數間存在的反向趨勢性導致了洪峰序列在M-K 檢驗時趨勢不顯著而LR檢驗時趨勢性顯著；同樣的結果適用于洪量序列，其尺度和形狀參數也存在反向趨勢性（斜率分別為0.004 3，-0.003 9）。

圖3 分布參數趨勢的協同性對極值序列非平穩性檢測結果的影響Fig.3 The influence of the trend synergy of distribution parameters on the detection of nonstationarity of extreme value series

4.3 非平穩條件下多變量洪水設計值推求

由于洪峰、洪量序列存在一定的趨勢型非平穩特征，為此特定聯合重現期下的洪峰-洪量的設計分位數對存在一定的時變特性（圖4）。雖然為Q-V之間的相依性結構的最優擬合模型為FR0 模型，表明其沒有顯著非平穩性，但在50 年一遇的聯合重現期（50-JRP）水平（或圖5 中的BR=0.02）下，Q-V的設計分位數等高線是逐年變化的。由于樣本量為58（1959-2016），每年可以生成這兩個極值變量在50-JRP水平上的等值線，這將繪制58 條JRP等值線。為了直觀呈現設計分位數的時變性，1960 年、1970 年、1990 年、2000 年、2010 的JRP等值線如圖5 所示。圖5 中綜合考慮氣候變化、城市化因素的情形對應的50-JRP等值線表現出3 種時變特征:①1960-1970 年間， 洪峰、洪量的設計值呈增大趨勢；②1970-1990 年間呈現下降趨勢；③1990-2010年呈現逐漸增大趨勢。單一考慮氣候變化或者城市化因素都會出現一定程度上的高估或低估洪峰-洪量的設計分位數。

圖4 考慮不同因素影響下的時變聯合重現期等值線圖（JRP=50）Fig.4 Contour map of time-varying joint return period considering the influence of different factors（JRP=50）

圖5 城市化、氣候變化對洪水單變量、多變量風險的影響Fig.5 Impacts of urbanization and climate change on flood univariate and multivariate risks

4.4 不同協變量影響下的單變量、多變量風險值對比分析

為了定量剝離城市化（PISA）、氣候變化（Maxpre+NAO）對于單變量、多變量風險的影響，依據公式（7）和（9），考慮經驗頻率值分別為［70%，90%］（步長為0.01）對應洪峰、洪量分位數集合，獲得不同非平穩模型（GEV1-GEV3， GA1-GA3）下對應的單變量（UR）和多變量風險值（BR）。為了量化城市化對于單變量和多變量風險的影響，可以分別依據GEV3（綜合考慮了城市化+氣候變化因素）模型和GEV2模型（僅僅考慮氣候變化因素）中分別計算70%～90%分位數區間｛（Q70%～Q90%)， （V70%～V90%）｝對應的單變量和多變量風險值［22，28］。例如對于經驗分位數Q70%，對于GEV3和GEV1 可以分別計算兩個不同的單變量風險值（URGEV1-Q70%，URGEV3-Q70%）（如圖5中左上角回歸線圖）。接著可以對風險樣本序列（URGEV1，URGEV3）進行線性擬合，其中城市化對單變量風險的影響可以通過回歸曲線斜率和對角線函數斜率之間的相對差值百分率來定量表示（對角線斜率等于1）。由圖5可知，對于蚌埠水文站控制的子流域，城市化對于洪水單變量、多變量風險的影響明顯弱于氣候變化因素對于洪水風險值的影響。

5 分析與討論

由于頻率分布模型地建立是基于非平穩性假設基礎之上，為此采用有效的非平穩性診斷/識別方法是進行非平穩頻率分析的必要前提。非平穩性診斷是從數理統計角度判斷水文極值序列是否存在趨勢的特征。目前，常用的水文序列趨勢診斷方法是基于非參數法的Mann-Kendall（M-K）和Spearman。同時對于水文極值序列的趨勢性檢測在水文氣象極端事件領域開展了大量的研究，如洪水趨勢分析［29，30］、極端降雨趨勢分析［31，32］，干旱趨勢分析［3，33］。本文采用的LR-AICc趨勢型非平穩檢測方法相比于非參數的M-K檢測法，有效解決了序列的非單調性趨勢 （二次型的趨勢:例如先上升后下降型）和參數的反向趨勢性問題，為后續非平穩洪水頻率分析模型的構建奠定良好的基礎。

由于多變量水文序列的概率分布可以分解為單變量的邊緣概率分布和變量間的相依結構描述兩方面，非平穩性頻率模型構建同樣分為基于邊緣分布的非平穩模型和基于聯合分布的非平穩模型。盡管以物理因子作為洪水序列統計參數的解釋變量較時間變量能夠更好地描述與預測洪水序列中的非平穩性，但是洪水頻率分析最終還是以洪水風險的計算作為最終落腳點，同時定量量化各個驅動因子對于洪水設計值的影響程度有利于解析洪水過程的驅動演進機制，本研究采用不同物理協變量組合下不同非平穩模型對比分析的思路有利于定量剝離氣候變化因子、人類活動因子（城市化）對于洪水多變量風險的影響。

6 結 論

通過構建動態Copula 函數為核心的聯合分布模型，采用LR-AICc 結合的趨勢性檢測手段對淮河流域蚌埠水文站的洪峰-洪量序列進行非平穩條件下的多變量頻率分析，主要獲得如下結論。

（1）由于城市化、氣候變化因素的綜合影響下，蚌埠水文站的洪峰、洪量呈現了一定的趨勢型非平穩特征。其中以PISA因子、Maxpre以及NAO因子為分布參數的協變量非平穩廣義極值分布模型（GEV3）是洪峰序列的最優邊緣分布模型，而以PISA因子、Maxpre以及NAO因子為分布參數的協變量非平穩Gamma分布模型（GA3）是洪量序列的最優邊緣分布模型；平穩Frank Copula模型為Q-V相依性結構的最優聯合分布模型。

（2）由于邊緣分布參數的時變特性，50 年一遇聯合重現期水平下的洪峰-洪量分位數對呈現一定的時變特性。蚌埠水文站控制的子流域內，氣候變化因素是影響洪水單變量、多變量風險的主導因素。