宋佳佳 張芳齊

摘要：為研究山西省濁漳河北源流域洪水特征，采用濁漳河北源石棧道水文站1957～2016年共60 a逐日流量資料，對洪峰流量、1日洪量進行研究。通過Mann-Kendall非參數(shù)檢驗法，揭示了該水文站洪峰流量和1日洪量演變特征及規(guī)律，并通過Copula函數(shù)研究了濁漳河北源流域洪水峰量聯(lián)合分布情況。結(jié)果表明：洪峰流量和1日洪量均呈現(xiàn)明顯的下降趨勢;時間序列上洪峰流量在1971年后有顯著的下降趨勢，1日洪量在1974年后呈現(xiàn)顯著的下降趨勢;Copula函數(shù)可較好擬合濁漳河北源洪水峰量的相關(guān)關(guān)系。研究成果可為該流域水利工程規(guī)劃設(shè)計和風險評估提供科學依據(jù)。

關(guān)鍵詞：洪峰流量;1日洪量;Mann-Kendall;Copula函數(shù);濁漳河北源流域

中圖法分類號：TV122.5 文獻標志碼：A DOI：10.15974/j.cnki.slsdkb.2021.10.001

文章編號：1006 - 0081（2021）10 - 0006 - 06

0 引 言

在全球氣候變化和人類活動的影響下，流域內(nèi)洪水情勢發(fā)生變化，需對其變化和發(fā)展趨勢進行分析。洪水峰量是研究洪水的重要特征，它不僅反映洪水的強度，也預示防洪的級別，因此對洪水峰量的演變規(guī)律研究尤為重要[1]。洪水需要多個特征變量才能完整描述，如洪峰、峰現(xiàn)時間、時段洪量和洪水歷時等存在相關(guān)性的特征變量。多個特征變量的聯(lián)合分析相比于單變量分析能更好地描述洪水特性并將有利于建設(shè)安全的防洪工程系統(tǒng)。Copula函數(shù)是多變量聯(lián)合分布構(gòu)建理論與方法的重大突破，De Michele[2]采用Copula函數(shù)建立了降雨強度和降雨歷時的聯(lián)合概率分布。熊立華等[3]采用Copula函數(shù)構(gòu)建了長江流域某站點洪峰和洪量的聯(lián)合分布。近十幾年來，利用Copula函數(shù)研究洪水的成果眾多，主要集中在洪水峰量、洪水歷時等變量的聯(lián)合分布和重現(xiàn)期研究上[4]。多變量重現(xiàn)期以聯(lián)合重現(xiàn)期、同現(xiàn)重現(xiàn)期和Kendall重現(xiàn)期應用最為廣泛[5] ， 結(jié)構(gòu)荷載重現(xiàn)期[6]也引起了關(guān)注。

為較全面研究濁漳河北源流域洪水極端事件，需研究多變量的聯(lián)合分布，Copula函數(shù)是實現(xiàn)這種相關(guān)性分析的有效方法。本文以濁漳河北源流域為例，采用石棧道水文站實測水文資料，利用Mann-Kendall非參數(shù)檢驗法，揭示了濁漳河北源洪水峰量的演變特征及規(guī)律。洪水峰量均單獨服從P-Ⅲ型分布，通過Copula函數(shù)將濁漳北源流域洪水峰量結(jié)合起來建立聯(lián)合分布模型，為水利工程規(guī)劃設(shè)計和風險評估提供科學依據(jù)。

1 研究區(qū)概況

濁漳河屬漳衛(wèi)南運河水系，上游有北、南、西三大支流，稱濁漳北源、濁漳南源、濁漳西源。濁漳河北源是濁漳河三源之一，發(fā)源于山西省晉中市榆社縣社城鎮(zhèn)焦紅寺大牛村，在山西省長治市襄垣縣小峧村南匯入濁漳河干流（圖1）。濁漳北源流域總面積3 797 km2，主河道長135 km。

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)來源

本文數(shù)據(jù)來源為濁漳河北源河道榆社縣石棧道水文站，利用石棧道水文站1957～2016年共60 a逐日流量資料，采取年內(nèi)日最大選樣法選取洪峰流量、1 日洪量為研究對象，運用Mann-Kendall檢驗法對洪水峰量進行演變規(guī)律分析，并采用Copula函數(shù)對石棧道水文站站址處洪水峰量的聯(lián)合分布進行研究。

2.2 數(shù)據(jù)資料“三性分析”

2.2.1 可靠性

石棧道水文站位于晉中市榆社縣石棧道村濁漳河北源河道上，地理坐標為東經(jīng)112°58′、北緯37°04′，于1957年6月設(shè)立。該站屬于山西省水文水資源勘測局，為基本水文站。該站觀測編制水文資料均符合有關(guān)規(guī)程、規(guī)范要求，資料質(zhì)量可靠，精度滿足水文分析要求。

2.2.2 一致性

石棧道水文站上游約15 km處建有雙峰水庫，控制流域面積190.2 km2。雙峰水庫自1975年興建以來，共經(jīng)過兩次續(xù)建，于2008年完成全部工程，達到現(xiàn)狀規(guī)模。由于庫區(qū)移民問題，2017年才開始蓄水。本文選取石棧道水文站1957～2016年實測流量系列資料具有一致性。

2.2.3 代表性

石棧道水文站實測流量系列資料為1957～2016年，流量系列相對較長，包括了豐、平、枯年。繪制石棧道水文站模比差積曲線及洪峰流量和1 日洪量變化過程線（圖2～3）。石棧道水文站年最大洪峰流量和1日洪量變化規(guī)律基本一致，即大洪峰大洪量，小洪峰小洪量，但年最大洪峰流量和最大1日洪量不一定出現(xiàn)在同一日。石棧道水文站實測流量系列資料具有代表性。

2.3 研究方法

2.3.1 Mann-Kendall非參數(shù)檢驗法

本文采用Mann-Kendall非參數(shù)檢驗法研究水資源演變特征[7]，選擇95%的置信水平檢驗，取值為1.96。同時，利用Mann-Kendall非參數(shù)檢驗法進行突變識別。

2.3.2 邊緣分布函數(shù)

P-Ⅲ型曲線是一條一端有限一端無限的不對稱單峰、正偏曲線，被引入水文統(tǒng)計中后，中國基本采用P-Ⅲ型曲線，其概率密度函數(shù)見式（1）。本文選用矩法估計，得出P-Ⅲ型曲線的參數(shù)[5]，假設(shè)一隨機變量[X（x1，x2，...，xn）]：

[f（x）=βαΓ（α）（x-a0）α-1e-β（x-a0）] （1）

其中，? ? ? ? ? ? ? ? ? ? [α=4C2s]

[β=2xCvCs]

[a0=x1-2CvCs]

[x=1ni=1nxi]

[Cv=i=1nki-12n-1]

[Cs=i=1nki-13（n-3）C3v]

[ki=xix（i=1， 2， ...， n）]

式中：[α]，[β]，[a0]分別為P-Ⅲ型分布的形狀、尺度和位置參數(shù);[Γ（α）]為伽瑪函數(shù);[x]，[Cv]，[Cs]分別為系列均值、變差系數(shù)、偏態(tài)系數(shù)。

2.3.3 Copula函數(shù)與參數(shù)估計

通過洪水峰量來分析洪水頻率及重現(xiàn)期時，需計算洪水峰量兩者的聯(lián)合關(guān)系，水文領(lǐng)域中常利用Copula函數(shù)建立兩者的聯(lián)合分布[7-8]。常用的Copula函數(shù)分為Archimedean、橢圓型和二次型三大類[9] 。 其中，Archimedean Copula函數(shù)具有結(jié)構(gòu)簡便、計算相對簡單、可構(gòu)造形式多樣、適應性強等優(yōu)點，在實際應用中較為廣泛。分析洪水事件中最常用的Archimedean Copula 函數(shù)主要有Gumbel-Hougaard、Clayton和Frank Copula。

本文先計算出3種Copula函數(shù)的參數(shù)θ，然后采用Kolmogorov-Smirnov檢驗的統(tǒng)計值對Copula函數(shù)模型進行檢驗，選用均方根誤差（RMSE）確定擬合優(yōu)度最高的Copula函數(shù)。利用選取的最優(yōu)Copula函數(shù)分析洪水峰量的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期，可以更好地掌握洪水事件特征，并應用于設(shè)計洪水估算和風險分析。

假設(shè)二維隨機變量D和S，它們的邊緣分布函數(shù)是[u=FD（d）=P（D≤d）]，[v=FS（s）=P（S≤s）]，聯(lián)合分布函數(shù)為[FD，S（d， s）=P[D≤d， S≤s]]，則三者表示為

[u=FD（d）， v=FS（s）]，

[FD， S（d， s）=exp{-[（-lnu）θ+（-lnv）θ]1θ}]

[FD，S（d， s）=（u-θ+v-θ-1）-1θ]

[FD， S（d， s）=-1θln[1+（e-θu-1）（e-θv-1）（e-θ-1）]]

式中：θ為Copula函數(shù)的參數(shù)。

本研究采用相關(guān)指標法進行Copula函數(shù)參數(shù)估計，結(jié)果如表1所示。

二維Copula函數(shù)經(jīng)驗頻率計算公式如下：

[Po（i）=（mi-0.44）/（n+0.12）] （2）

式中：[mi]為聯(lián)合觀測樣本中滿足條件[D≤di]且[S≤si]的觀測個數(shù)，[n]為樣本容量

采用均方根誤差評定各種Copula函數(shù)擬合結(jié)果

[RRMSE=1ni=1nPc（i）-Po（i）2] （3）

式中：[Pc（i）]為理論聯(lián)合頻率值，[Po（i）]為經(jīng)驗聯(lián)合頻率值。

[TDS=E（L）PD>d?S>s = E（L）1-FD（d）-FS（s）+FD， S（d，s）] （4）

[TDS′=E（L）PD>d?S>s=E（L）1-FD， S（d， s）]（5）

式中：[TDS]為洪旱事件的同現(xiàn)重現(xiàn)期（[D>d]且[S>s]），[TDS′]為洪旱事件的聯(lián)合重現(xiàn)期（[D>d]或[S>s]），[E（L）]為洪旱間隔的期望值。

根據(jù)重現(xiàn)期來描述洪水事件的嚴重性，洪水峰量聯(lián)合分布的重現(xiàn)期包括[D>d]或[S>s]和[D>d]且[S>s]兩種情況。

3 洪水峰量突變及顯著性水平檢驗

根據(jù)M-K檢驗成果（圖4），洪峰流量在1960～1964年處于增長波動變化，其他年份均呈現(xiàn)波動減少變化。洪峰流量UF和UB曲線相交于1971年，且均落于0.05顯著性水平閾值內(nèi)，即洪峰流量在1971年出現(xiàn)突變，因正序列曲線在1976年突破0.05顯著性水平閾值，表明洪峰流量出現(xiàn)顯著性突變減少。

1日洪量在1960～1964年處于增長波動變化，其他年份均呈現(xiàn)波動減少變化。1日洪量UF和UB曲線相交于1974年，且均落于0.05顯著性水平閾值內(nèi)，即1日洪量在1974年出現(xiàn)突變，因正序列曲線在1985年突破0.05顯著性水平閾值，表明1日洪量出現(xiàn)顯著性突變減少。石棧道水文站1957～2016年洪水峰量平均值及UF均值見表2。

4 洪水峰量聯(lián)合分布研究

4.1 洪水特征變量邊緣分布

采用適線法得到石棧道水文站洪水峰量的頻率曲線。洪峰流量采用[Q均值=353.78] m3/s，[Cv=0.84]，[Cs/Cv=1.82];最大1日洪量采用[W均值=536.10]萬m3，[Cv=0.94]，[Cs/Cv=2.16]。同時，采用Kolmogorov-Smirnov進行檢驗，洪峰流量和1日洪量K-S統(tǒng)計檢驗值分別為0.063 3和0.133 5，顯著水平0.05對應的臨界值是0.172 3，檢驗結(jié)果均為0，通過檢驗，可認為洪峰流量和1日洪量均服從P-Ⅲ型分布。繪制洪峰流量和1日洪量的頻率累積曲線，見圖5。

4.2 Copula函數(shù)參數(shù)估計及擬合優(yōu)度檢驗

采用Gumbel-Hougaard、Clayton和Frank Copula三種Copula函數(shù)建立洪水峰量的聯(lián)合分布，分別運用相關(guān)指標法估計Copula函數(shù)參數(shù)，并采用均方根誤差評定各種Copula函數(shù)擬合結(jié)果。結(jié)果表明：濁漳河北源石棧道水文站適合運用Gumbel-Hougaard Copula進行洪水峰量特征分析，擬合誤差僅為0.300 5（表3）。

4.3 組合變量分布函數(shù)與單變量分布函數(shù)對比

（1）組合變量與單變量對應的重現(xiàn)期不同。1996年最大洪峰流量為951 m3/s，最大1日洪量為1 866萬m3，由P-Ⅲ曲線可知，洪峰流量為951 m3/s時，重現(xiàn)期為15 a，洪量為1 866萬m3時重現(xiàn)期為20 a。按G-H Copula函數(shù)洪水峰量聯(lián)合分布計算，同時發(fā)生最大洪峰流量為951 m3/s、最大1日洪量為1 866億m3大洪水的重現(xiàn)期為59 a。采用邊緣分布函數(shù)時，1996年場次洪水洪峰流量的重現(xiàn)期為15 a，洪量的重現(xiàn)期為20 a。采用G-H Copula聯(lián)合分布函數(shù)時，相同大小的洪峰或洪量任一出現(xiàn)的聯(lián)合重現(xiàn)期[TDS′]為13 a，相同大小的洪峰和洪量同時出現(xiàn)的同現(xiàn)重現(xiàn)期[TDS]為59 a。

（2）聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期可作為實際重現(xiàn)期的區(qū)間估計。選取不同的邊緣分布重現(xiàn)期得到聯(lián)合分布重現(xiàn)期，邊緣分布的重現(xiàn)期介于[TDS′]與[TDS]之間， 聯(lián)合分布的兩種重現(xiàn)期可以看作邊緣分布的兩種極端情況。可以根據(jù)聯(lián)合分布的重現(xiàn)期作為實際重現(xiàn)期的區(qū)間估計，當邊緣分布的重現(xiàn)期T為2 a時，石棧道水文站實際重現(xiàn)期為1.3～2.7 a，當邊緣分布的重現(xiàn)期為100 a時，石棧道站實際的重現(xiàn)期在74.6～142.3 a之間，見表4。

5 結(jié) 論

（1）根據(jù)洪水峰量的趨勢性分析，洪峰流量和1日洪量均呈現(xiàn)顯著的下降趨勢，洪峰流量在1976年以后開始表現(xiàn)出明顯的下降，1日洪量在1985年以后開始表現(xiàn)出明顯的下降。

（2）通過Mann-Kendall非參數(shù)檢驗法對石棧道水文站1957～2016年洪峰流量和1日洪量序列變化規(guī)律進行了分析。結(jié)果表明：濁漳河北源洪峰流量在1971年發(fā)生了突變，1日洪量在1974年發(fā)生了突變，均呈現(xiàn)顯著的下降趨勢。

（3）1957～2016年石棧道水文站洪峰流量和1日洪量變化過程線可看出，石棧道水文站洪峰流量和1日洪量變化規(guī)律基本一致，即大洪峰大洪量，小洪峰小洪量。通過矩法估計P-Ⅲ型分布函數(shù)的參數(shù)值，同時采用Kolmogorov-Smirnov進行檢驗，洪峰流量和1日洪量K-S統(tǒng)計檢驗值分別為0.063 3和0.133 5，顯著水平0.05對應的臨界值是0.172 3，檢驗結(jié)果均為0，通過檢驗，可認為洪峰流量和1日洪量均服從P-Ⅲ型分布。

（4）采用均方根誤差評定Gumbel-Hougaard，Clayton和Frank Copula三種Copula函數(shù)擬合結(jié)果。結(jié)果表明，濁漳河北源石棧道水文站適合運用Gumbel-Hougaard Copula進行洪水峰量特征分析，擬合誤差僅為0.300 5。

（5）基于G-H Copula函數(shù)研究分析洪峰、洪量聯(lián)合分布下的聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期。結(jié)果顯示，選取不同的邊緣分布重現(xiàn)期得到聯(lián)合分布重現(xiàn)期，邊緣分布的重現(xiàn)期介于聯(lián)合重現(xiàn)期與同現(xiàn)重現(xiàn)期之間。聯(lián)合分布的兩種重現(xiàn)期可以看作邊緣分布的兩種極端情況，可將聯(lián)合分布的重現(xiàn)期作為實際重現(xiàn)期的區(qū)間估計。

參考文獻：

[1] 唐建. 基于MK檢驗的天山山區(qū)近55年降水特征分析[J]. 甘肅水利水電技術(shù)，2019，7（1）：5-8.

[2] De MICHELE C，SALVADORI G.A Generalized Pareto intensity-duration model of storm rainfall exploiting 2-Copulas[J]. Journal of Geophysical Research： Atmospheres， 2003， 108（D2）： 4067.

[3] 熊立華， 郭生練. 兩變量極值分布在洪水頻率分析中的應用研究[J]. 長江科學院院報，2004，21（2）：35-37.

[4] 唐建. 基于Copula函數(shù)的洪峰流量與洪水歷時聯(lián)合分布研究[J]. 中國農(nóng)村水利水電，2017（2）：204-214.

[5] 何英，彭亮，鄭淑文，等. 基于Copula函數(shù)的葉爾羌河流域洪水要素聯(lián)合分布研究[J]. 中國農(nóng)村水利水電，2019（4）：74-79.

[6] 劉章君，郭生練，許新發(fā)，等. 兩變量洪水結(jié)構(gòu)荷載重現(xiàn)期與聯(lián)合設(shè)計值研究[J]. 水利學報，2018，49（8）： 956-965.

[7] 梁艷芹. 海河流域暴雨洪水演變趨勢分析[J]. 南水北調(diào)與水利科技，2014（3）：180-185.

[8] 楊興，趙玲玲，陳子燊，等. 中小流域洪峰流量與水位聯(lián)合分布的設(shè)計洪水分析[J]. 水電能源科學，2019，（8）：43-46.

[9] 郭生練，閆寶偉，肖義，等. Copula函數(shù)在多變量水文分析計算中的應用及研究進展[J]. 水文， 2008，28（3）：1-7.

（編輯：李 慧）

Analysis of evolution law and combined frequency of flood peak and volume in northern source of Zhuozhang River of Shanxi Province

SONG Jiajia ， ZHANG Fangqi

（Institute of Architectural Design and Research，Taiyuan University of Technology，Taiyuan 030024，China）

Abstract：To study the flood characteristics in the northern source of Zhuozhang River of Shanxi Province， the daily discharge data of 60 years from 1957 to 2016 at the Shizhandao Hydrological Station in the basin were used to analyze the flood peak and maximum 1 day flood volume. Through Mann-Kendall nonparametric test method， the characteristics and evolution laws of flood peak and maximum 1 day flood volume at the Shizhandao Hydrological Station were revealed， and the combined frequency of flood peak volume in the northern source of the Zhuozhang River was studied by Copula function. The results showed that the peak discharge and maximum 1 day flood volume were in an obvious decreasing trend; in the time series， the peak flood discharge has a significant decreasing trend after 1971， and the 1 day flood volume has a significant decreasing trend after 1974; Gumbel-Hougaard Copula function can well fit the correlation between the peak and volume of the flood in the northern source of the Zhuozhang River. The results can provide scientific basis for planning and risk assessment of water conservancy projects in the basin.

Key words： flood peak flow; maximum 1 day flood volume; Mann-Kendall test; Copula function;? northern source basin of Zhuozhang River