基于Copula函數的頭道拐水沙關系頻率分析研究

李弘瑞 李新杰 張紅濤

摘 要：為了掌握黃河中游水沙關系的變化特點和變化趨勢，基于黃河頭道拐站1958—1990年33 a和1998—2017年20 a的水沙資料，通過P-Ⅲ型曲線分別建立年徑流量與年輸沙量的頻率曲線，對年徑流量和年輸沙量進行豐平枯劃分;基于Copula函數對水沙關系進行變異診斷，構建年徑流量-年輸沙量聯合分布模型并通過擬合優度檢驗選取最合適的Copula函數，計算水沙豐平枯組合遭遇頻率。結果表明：頭道拐水文站在1958—1990年和1998—2017年的年徑流量和年輸沙量呈下降趨勢，1990年水沙關系發生變異;1958—1990年水沙關系存在一個變異點（1965年），1998—2017年水沙關系不存在明顯的變異點，整體水沙關系存在一定的平穩性;1966—1990年和1998—2017年兩個時期的水沙豐平枯組合遭遇頻率中豐枯同頻高于豐枯異頻，1966—1990年的年徑流量和年輸沙量在較大值處相關性較強，1998—2017年年徑流量和年輸沙量整體相關性較強，水量和沙量的豐平枯變化趨勢基本一致。

關鍵詞：水沙關系;聯合分布;變異診斷;遭遇頻率;頭道拐水文站;黃河

Abstract： In order to grasp the characteristics and trends of the relationship between runoff and sediment at the upper and middle reaches of the Yellow River， based on the water and sediment volume data of the Yellow River Toudaoguai Station in the 33 years from 1958 to 1990 and the past 20 years from 1998 to 2017， the P-Ⅲ curve was adopted to establish the runoff and sediment discharge frequency and divide annual runoff and annual sediment transport into hight， medium and low flow. It carried out mutation diagnosis of runoff and sediment relationship， and joint distribution model of annual runoff and annual sediment transport was constructed based on copula function by the goodness of fit test. Finally， it calculated the encounter frequency of wet-normal-dry runoff and sediment combination. The results show that the annual runoff and sediment discharge of Toudaoguai Hydrological Station in 1958-1990 and 1998-2017 show a downward trend and the relationship between runoff and sediment is remarkable especially after 1990; There is a variation point （1965） in the relationship between water and sediment during 1958-1990. There is no obvious variability point in the water and sediment relationship between 1998 and 2017 and there is a certain stability in the overall runoff and sediment relationship; the synchronous frequency in the frequency of the combination of wet-normal-dry runoff and sediment is greater than that the asynchronous frequency during 1966-1990 and 1998-2017， the annual runoff and annual sediment transport in the previous period have a strong correlation at the annual sediment load is larger. The correlation between annual runoff and annual sediment transport in the latter period is stronger and the variation trend of annual runoff volume has similar variation tendency to the annual sediment.

Key words： runoff and sediment relationship; joint distribution; mutation diagnosis; encounter frequency; Toudaoguai Hydrological Station; Yellow River

1 引 言

隨著黃河上游水利工程的運行和中游寧蒙河段灌區取用水等人類活動的影響，進入黃河中游的水沙過程發生了很大變化，對黃河中下游河道及水庫沖淤、流域生態環境等產生了重要影響[1]。因此，科學揭示水沙關系變化規律對黃河中下游水資源管理和高效利用具有重要意義。

針對黃河流域水沙關系，有關學者開展了相關研究[2-4]。金鑫等基于水沙序列的理論頻率曲線提出了黃河中游不同水文站的水沙頻率組合劃分標準[5];李新杰等從保證率的角度建立了潼關站1961—2018年58 a的年徑流量和輸沙量的頻率曲線，并確定了潼關站來水來沙豐平枯的劃分標準[6];馬雁等基于黃河上游蘭州站2014—2018年的水沙數據，總結了年徑流量和年輸沙量的關系和演變規律[7];考慮到環境變化和人類活動對流域來水來沙量的影響，李艷玲等結合滑動窗口和Copula函數，提出基于滑動Copula函數的降水和徑流關系變異診斷法[8];郭愛軍等提出滑動相關系數法用于涇河流域水沙關系的變異診斷，并基于Copula函數分析了徑河流域的水沙關系演變特征[9];姚曼飛等基于涇河干支流8個水文站實測水沙數據，采用Copula函數構建了水沙聯合分布模型，并計算了各個水文站的水沙豐平枯遭遇頻率[10]。

鑒于頭道拐站水沙關系對黃河中下游水沙調控的重要作用，以及已有研究多集中于徑流量和輸沙量特點或水沙關系單方面性質變化，對水沙組合遭遇情況分析稍有不足，筆者統計頭道拐站1958—1990年和1998—2017年共計53 a的年徑流量和年輸沙量原始數據，通過建立水沙頻率曲線，診斷水沙關系的變異特性，計算水沙豐平枯組合的遭遇頻率，基于二維Copula函數建立二維水沙聯合分布函數，分析黃河頭道拐站水沙關系整體的變化情況及豐平枯遭遇規律，以期為流域水沙調控提供技術支撐。

2 研究數據與方法

2.1 研究數據

頭道拐水文站位于內蒙古托克托縣雙河鎮，是黃河上游的出口站，也是中游萬家寨水庫的入口站，其徑流主要來自蘭州站以上，泥沙主要來自蘭州至頭道拐區間支流，其水沙變化直接影響萬家寨水庫的水沙調度情況[11-12]。本研究選取頭道拐水文站1958—1990年和1998—2017年兩個時期的年徑流量和年輸沙量數據，資料來源于黃河流域水文年鑒。

2.2 研究方法

2.2.1 頻率曲線

通過對年徑流量和年輸沙量序列的計算得到基本參數、CV、CS，進而確定參數α、β和b的值，即確定P-Ⅲ型分布函數。

本文運用適線法對頻率曲線進行尋優，通過P-Ⅲ型頻率分布建立年徑流量和年輸沙量的理論頻率曲線，應用數學期望公式計算年徑流量和年輸沙量的經驗頻率，計算理論頻率和經驗頻率擬合時的確定性系數R2，其計算公式如下[15]：

2.2.2 Copula函數

Copula函數是把隨機變量X1、X2、…、XN的聯合分布函數與各自的邊緣分布函數U1、U2、…、UN相連接的連接函數，即函數C（u1，u2，…，uN）[16]。

2.2.3 Copula函數類型及參數估計

Copula函數在將邊緣分布進行聯合分布時包含了變量所有的相依信息，在建立聯合分布函數的過程中保持信息的真實度[17]，更好地反映出水沙序列邊緣分布之間的關系變化。Copula函數在總體上大致分為3種：正態型、t型和Archimedean型。在水文統計中，主要使用Archimedean型Copula函數對水沙序列進行聯合分布，具體的表達式及參數范圍見表1[17]。

2.2.4 Copula函數擬合優度檢驗

為了檢驗Copula函數對于水沙序列聯合分布的擬合優度，引入經驗Copula函數，比較經驗累計頻率與理論累計頻率的擬合度。將樣本值代入Copula頻率公式計算出經驗累計頻率，同樣將樣本值代入3類Archimedean型Copula函數計算理論累計頻率。計算確定性系數R2并進行比較，再根據離差平方和OLS最小準則與AIC信息準則做進一步的擬合優度檢驗[18]。

通過擬合優度檢驗選取最合適的Copula函數建立水沙序列的聯合分布，從而得到接近實際值的理論分布。

2.2.5 變異診斷

本文基于Copula函數對水沙序列中水沙關系在時間分布上的變異情況進行診斷[21]。構造基于Copula函數的極大似然對數統計量公式如下：

原序列經過一次診斷后，若發現變異點則依據變異點劃分為兩個子序列，對子序列進行二次診斷，重復以上操作直至子序列的長度小于25。若原序列長度小于25，一次診斷即可。

對原水沙序列進行變異診斷是為了發現水沙序列在時間分布上是否存在平穩性，若存在，對原序列的水沙關系進行整體分析;若不存在，則根據變異點將原序列劃分為若干個子序列，對每個子序列的水沙關系獨立分析。

2.2.6 水沙豐平枯組合遭遇頻率劃分

在《水文基本術語和符號標準》[23]中，將河流的豐平枯年頻率劃分為5個級別，通過不同的保證率確定豐平枯年的頻率劃分標準，見表2[24]。

由于徑流量和輸沙量之間相依度高、相關性強，一般認為輸沙量頻率分布與徑流量頻率分布類似，因此輸沙量的豐平枯年的頻率劃分也可以采用表2的標準[25]。通過建立滿足P-Ⅲ型分布的年徑流量和年輸沙量頻率分布，依據頻率劃分標準確定豐平枯年年徑流量和年輸沙量的劃分標準，把年徑流量和年輸沙量均劃分為5個級別，水沙豐平枯組合具體劃分標準見表3。

3 數據分析

3.1 確定水沙頻率曲線

本文選取頭道拐站1958—1990年和1998—2017年兩個時段的年徑流量和年輸沙量序列，通過P-Ⅲ型分布函數建立年徑流量和年輸沙量序列的理論頻率分布，同經驗頻率分布進行擬合，在滿足確定性系數R2≥0.95的前提下，得到擬合度較高的年徑流量和年輸沙量的理論頻率曲線，見圖1和圖2。兩個時段年徑流量和年輸沙量頻率變化情況見圖3。

3.2 水沙聯合分布及擬合優度檢驗

在確定合適的邊緣分布后，作為變量輸入Clayton、Gumbel和Frank Copula函數以及經驗Copula函數，分別建立年徑流量和年輸沙量的聯合分布與經驗聯合分布并進行擬合，根據確定性系數R2、OLS和AIC信息準則3個擬合優度檢驗指標選取最合適的Copula函數，3類Archimedean型Copula函數的參數及擬合優度檢驗值的計算結果見表4。

根據表4的檢驗值可知：1958—1990年的聯合分布Clayton Copula的R2值最大、OLS和AIC的值均最小，即Clayton Copula的擬合優度最佳，表明Clayton Copula函數是用于頭道拐站1958—1990年水沙聯合分布最合適的Copula函數;1998—2017年的聯合分布中Frank Copula的R2值最大、OLS和AIC的值均最小，說明Frank Copula的擬合優度最佳，明顯優于其他兩類，Frank Copula函數是用于頭道拐站1998—2017年水沙聯合分布最合適的Copula函數。

3.3 水沙關系的變異診斷

基于Gumbel Copula函數對1958—1990年和1998—2017年水沙關系進行變異診斷（一次診斷），見圖4。

從圖4中可知，兩時段水沙關系的一次診斷中統計量Zn最大值（41.174 4）出現在1990年，在顯著水平α為0.05時Zn遠大于閾值（3.2），故存在變異點（1990年）。說明1958—1990年和1998—2017年兩時段水沙關系存在差異。

分別基于Clayton Copula函數和Frank Copula函數對1958—1990年和1998—2017年年徑流量和年輸沙量關系進行變異診斷，見圖5。

從圖5中可知，在1958—1990年的一次診斷中統計量Zn最大值（5.933 7）出現在1965年，大于顯著水平α=0.05時閾值（3.2），故1965年為變異點。在二次診斷中統計量Zn最大值（1.695 2）小于閾值，故1966—1990年水沙關系不存在變異點，具有一定的平穩性;在1998—2017年的一次診斷中統計量Zn最大值（0.203 9）小于閾值，因此認為在1998—2017年水沙關系不存在變異點，具有一定的平穩性。

頭道拐1998—2017年水沙邊緣分布的基于Frank Copula的聯合分布如圖8所示。從圖8（a）中可見，1998—2017年水沙規律性較強，服從固定分布，年徑流量和輸沙量概率密度函數呈現尾部相關的特點，上尾部和下尾部相關性最強且具有對稱性，表明年徑流量的極大或極小值處對年輸沙量有較大影響;從圖8（b）中可見，Frank Copula的分布函數具有對稱性，則水沙聯合分布的累計頻率具有對稱性。

3.4 計算遭遇頻率

根據建立的頭道拐站年徑流量和年輸沙量的P-Ⅲ型頻率分布，確定頭道拐站水沙豐平枯組合中年徑流量（X）和年輸沙量（Y）的劃分標準。選取1966—1990年和1998—2017年兩段時期計算，這兩段時期已被驗證水沙關系存在一定的平穩性，見表5～表8。

通過Clayton Copula建立水沙聯合分布計算頭道拐站1966—1990年水沙豐平枯組合遭遇頻率，見表9。

由表9可知，頭道拐站1966—1990年水沙豐平枯組合遭遇頻率中，水沙同頻組合（同豐、同偏豐、同平、同偏枯和同枯）頻率總和為0.645 6，水沙異頻組合頻率為0.354 4，同頻組合頻率明顯高于異頻組合;水沙異頻頻率的高頻部分主要集中于偏豐水平沙組合、平水偏豐沙組合、平水偏枯沙組合和偏枯水平沙組合，其頻率均大于0.04;總體而言，同偏豐組合（0.186 7）最大，同平組合（0.137 4）次之，水沙豐枯相反組合的頻率最小，說明頭道拐站1966—1990年年徑流量和年輸沙量頻率的關聯性強，且在年徑流量和年輸沙量偏豐時關聯性最強。

通過Frank Copula建立水沙聯合分布計算頭道拐站1998—2017年水沙豐平枯組合遭遇頻率，見表10。由表10可知，頭道拐站1998—2017年水沙豐平枯組合遭遇頻率中，同頻組合頻率總和為0.646 1，異頻組合頻率總和為0.353 9，則同頻組合頻率明顯高于異頻頻率;異頻頻率的高頻部分主要集中于豐水偏豐沙組合、偏豐水豐沙組合、偏豐水平沙組合、平水偏豐沙組合、平水偏枯沙組合、偏枯水平沙組合、偏枯水枯沙組合和枯水偏枯沙組合，其頻率均大于0.04;總體而言，同偏豐（同偏枯）組合（0.162 2）最大，同平組合（0.156 8）次之，水沙同豐（水沙同枯組）組合（0.082 5）第三，而枯水豐沙組合（豐水枯沙組合）最小，說明頭道拐站1998—2017年年徑流量和年輸沙量頻率的關聯性強，水沙狀態完全相反的組合的遭遇頻率極小。

4 結 論

以頭道拐站1958—1990年和1998—2017年兩個時期的年徑流量和年輸沙量序列為研究對象，利用P-Ⅲ頻率曲線探究水沙豐平枯組合遭遇頻率，基于Copula函數構建二維水沙聯合分布，分析水沙聯合分布的累計頻率，得到以下結論。

（1）頭道拐站年徑流量和年輸沙量呈下降趨勢，在1990年水沙關系發生變異。1958—1990年水沙序列中存在一個變異點（1965年），1966—1990年和1998—2017年序列水沙關系未發生明顯變異，具有一定平穩性。

（2）由基于Copula建立的水沙聯合分布的概率密度函數和分布函數可知，1958—1965年和1966—1990年水沙子序列概率密度函數上尾部相關性較強，即水沙序列中年徑流量的較大值處變化趨勢同年輸沙量變化趨勢關聯性較強;1998—2017年水沙子序列的概率密度函數上尾部和下尾部相關性較強，即水沙序列中年徑流量的高低變化趨勢同年輸沙量高低變化趨勢關聯性較強。

（3）由水沙豐平枯組合遭遇頻率分析可得，在1966—1990年和1998—2017兩個時段的水沙組合中水沙同頻明顯高于水沙異頻，其中水沙異頻組合發生的概率較小，頭道拐站1966—1990年年徑流量和年輸沙量在平豐時相關性強，1998—2017年年徑流量和年輸沙量整體相關性較強，水沙豐平枯變化趨勢基本一致。

參考文獻：

[1] 胡春宏.黃河水沙變化與治理方略研究[J].水力發電學報，2016，35（10）：1-11.

[2] 趙陽，胡春宏，張曉明，等.近70年黃河流域水沙情勢及其成因分析[J].農業工程學報，2018，34（21）：112-119.

[3] 高宗軍，馮國平.黃河水沙變化趨勢及成因分析[J].地下水，2020，42（1）：147-151.

[4] 高航，姚文藝，張曉華.黃河上中游近期水沙變化分析[J].華北水利水電學院學報，2009，30（5）：8-12.

[5] 金鑫，郝振純，張金良.黃河中游水沙頻率關系研究[J].泥沙研究，2006，31（3）：6-13.

[6] 李新杰，郜國明，朱亮，等.黃河潼關站水沙關系頻率及協調度分析研究[J].人民黃河，2020，42（5）：47-51.

[7] 馬雁，賈生海，張彥洪.黃河上游蘭州水文站2004—2018年水沙特性分析[J].水利規劃與設計，2020（4）：52-54，89.

[8] 李艷玲，暢建霞，黃強，等.基于滑動Copula函數的降水和徑流關系變異診斷[J].水力發電學報，2014，33（6）：20-24，60.

[9] 郭愛軍，黃強，暢建霞，等.基于Copula函數的涇河流域水沙關系演變特征分析[J].自然資源學報，2015，30（4）：673-683.

[10] 姚曼飛，黨素珍，孟美麗，等.基于Copula函數的涇河流域水沙豐枯遭遇頻率分析[J].水土保持研究，2019，26（1）：192-196，202.

[11] 冉大川，姚文藝，張攀，等.黃河頭道拐站水沙來源空間分布及其影響因素[J].泥沙研究，2015，40（1）：42-48.

[12] 任智慧，王婷，曲少軍.萬家寨水庫庫區沖淤特點分析[C]//中國大壩工程學會.水庫大壩高質量建設與綠色發展中國大壩工程學會2018學術年會論文集.北京：中國大壩工程學會，2018：172-177.

[13] 雷冠軍，王文川，殷峻暹，等.P-Ⅲ型曲線參數估計方法研究綜述[J].人民黃河，2017，39（10）：1-7.

[14] 黃繼文.P-Ⅲ型分布頻率分析在Excel中的實現及應用[J].水資源研究，2006，27（4）：7-9.

[15] 郭生練，葉守澤.論水文計算中的經驗頻率公式[J].武漢水利電力學院學報，1992，25（2）：38-45.

[16] PATTON A J. A Review of Copula Models for Economic Time Series[J].Journal of Multivariate Analysis，2012，110：4-18.

[17] 杜懿，麻榮永.不同Copula函數在洪水峰量聯合分布中的應用比較[J].水力發電，2018，44（12）：24-26，58.

[18] 李天元，郭生練，羅啟華，等.雙參數Copula函數在洪水聯合分布中的應用研究[J].水文，2011，31（5）：24-28，46.

[19] 賈玉紅，宋松柏.陜北地區年降水量頻率分布參數估算研究[J].水資源與水工程學報，2012，23（5）：48-50.

[20] 宋喜芳，李建平，胡希遠.模型選擇信息量準則AIC及其在方差分析中的應用[J].西北農林科技大學學報（自然科學版），2009，37（2）：88-92.

[21] HUANG S Z， LI P， HUANG Q， et al. Copula-Based Identification of the Non-Stationarity of the Relation Between Runoff and Sediment Load[J].International Journal of Sediment Research， 2017，32（2）：221-230.

[22] COSTA D A D. Copula Inference for Finance and Insurance[D].Zurich，Switzer-Land：ETH，2004：123-141.

[23] 中華人民共和國水利部.水文基本術語和符號標準：GB/T 50095—2014[S].北京：中國計劃出版社，2014：2-51.

[24] 徐宇程，朱首賢，張文靜，等.長江大通站徑流量的豐平枯水年劃分探討[J].長江科學院院報，2018，35（6）：19-23.

[25] 林沫，劉穎，叢遠飛，等.松遼流域主要河流水沙規律分析[J].東北水利水電，2009，27（12）：40-42，72.

【責任編輯 張 帥】