山西省濁漳河北源流域洪水峰量演變規(guī)律及聯(lián)合分布研究

宋佳佳 張芳齊

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

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

中圖法分類號(hào)：TV122.5 文獻(xiàn)標(biāo)志碼：A DOI：10.15974/j.cnki.slsdkb.2021.10.001

文章編號(hào)：1006 - 0081（2021）10 - 0006 - 06

0 引 言

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

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

1 研究區(qū)概況

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

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

2.1 數(shù)據(jù)來(lái)源

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

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

2.2.1 可靠性

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

2.2.2 一致性

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

2.2.3 代表性

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

2.3 研究方法

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

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

2.3.2 邊緣分布函數(shù)

P-Ⅲ型曲線是一條一端有限一端無(wú)限的不對(duì)稱單峰、正偏曲線，被引入水文統(tǒng)計(jì)中后，中國(guó)基本采用P-Ⅲ型曲線，其概率密度函數(shù)見式（1）。本文選用矩法估計(jì)，得出P-Ⅲ型曲線的參數(shù)[5]，假設(shè)一隨機(jī)變量[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ù)估計(jì)

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

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

假設(shè)二維隨機(jī)變量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)指標(biāo)法進(jìn)行Copula函數(shù)參數(shù)估計(jì)，結(jié)果如表1所示。

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

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

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

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

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

式中：[Pc（i）]為理論聯(lián)合頻率值，[Po（i）]為經(jīng)驗(yàn)聯(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)期來(lái)描述洪水事件的嚴(yán)重性，洪水峰量聯(lián)合分布的重現(xiàn)期包括[D>d]或[S>s]和[D>d]且[S>s]兩種情況。

3 洪水峰量突變及顯著性水平檢驗(yàn)

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

1日洪量在1960～1964年處于增長(zhǎng)波動(dòng)變化，其他年份均呈現(xiàn)波動(dòng)減少變化。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]萬(wàn)m3，[Cv=0.94]，[Cs/Cv=2.16]。同時(shí)，采用Kolmogorov-Smirnov進(jìn)行檢驗(yàn)，洪峰流量和1日洪量K-S統(tǒng)計(jì)檢驗(yàn)值分別為0.063 3和0.133 5，顯著水平0.05對(duì)應(yīng)的臨界值是0.172 3，檢驗(yàn)結(jié)果均為0，通過(guò)檢驗(yàn)，可認(rèn)為洪峰流量和1日洪量均服從P-Ⅲ型分布。繪制洪峰流量和1日洪量的頻率累積曲線，見圖5。

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

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

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

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

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

5 結(jié) 論

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

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

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

（4）采用均方根誤差評(píng)定Gumbel-Hougaard，Clayton和Frank Copula三種Copula函數(shù)擬合結(jié)果。結(jié)果表明，濁漳河北源石棧道水文站適合運(yùn)用Gumbel-Hougaard Copula進(jìn)行洪水峰量特征分析，擬合誤差僅為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)期作為實(shí)際重現(xiàn)期的區(qū)間估計(jì)。

參考文獻(xiàn)：

[1] 唐建. 基于MK檢驗(yàn)的天山山區(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] 熊立華， 郭生練. 兩變量極值分布在洪水頻率分析中的應(yīng)用研究[J]. 長(zhǎng)江科學(xué)院院報(bào)，2004，21（2）：35-37.

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

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

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

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

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

[9] 郭生練，閆寶偉，肖義，等. Copula函數(shù)在多變量水文分析計(jì)算中的應(yīng)用及研究進(jìn)展[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