基于Copula似然比檢驗凌汛多變量序列變點診斷

李蓉蓉 于坤霞 熊斌

摘 要：在氣候變化和人類活動影響下，凌汛期的水文情勢分析對防凌減災研究有重要意義。以黃河寧蒙河段三湖河口水文站的開河期最大流量為研究對象，由相關性分析得到其影響變量——開河期最高水位和冰厚，利用Mann-Kendall和Pettitt檢驗進行單變量序列變點檢驗，在單變量變點檢驗結果的基礎上采用Copula似然比檢驗（CLR）方法對多變量凌汛序列進行變點診斷分析。結果表明：對于單變量序列，開河期最大流量無顯著突變，開河期最高水位在1968年發生突變，冰厚在1987年發生突變;對于三變量Pair-Copula結構，tree 1結構沒有發生顯著變化，tree 2結構在1996年發生突變。經過對比發現，變點前后變量間的相關性有顯著差異。造成突變的可能原因是河槽持續淤積導致封凍期和開河期水位逐年抬升。

關鍵詞：開河期最大流量;Copula函數;相關結構;變點診斷;黃河寧蒙河段

中圖分類號：P333;TV882.1 文獻標志碼：A

doi：10.3969/j.issn.1000-1379.2021.08.006

引用格式：李蓉蓉，于坤霞，熊斌.基于Copula似然比檢驗凌汛多變量序列變點診斷[J].人民黃河，2021，43（8）：33-38.

Abstract： Under the influence of climate change and human activities， the analysis of the Hydrological situation during ice flood period is of great significance to the study of prevention and mitigation of ice flood disaster. In this paper， the peak break-up discharge of Sanhuhekou Hydrological Station in Ning-Meng reach of the Yellow River was selected as the research object， its effect variables （the peak break-up water level and average ice thickness） were obtained through correlation analysis. Mann-Kendall and Pettitt tests were used for change-point tests of univariate variables， on the basis of these results， Copula-based likelihood ratio test （CLR） method was used for change-point detection of the multivariate series. The results show that for univariate series， there is no significant mutation for the peak break-up discharge， the peak break-up water level mutates in 1968 and average ice thickness mutates in 1987. For Pair-Copula structure， the tree 1 does not change significantly， while the tree 2 mutates in 1996. After comparison， it finds out that the correlation between variables before and after the change-point is significantly different. The possible reason for the mutation is that the continuous silting of the channel which leads to the keep rising of water level during the freeze-up period and the break-up period.

Key words： peak break-up discharge; Copula function; dependence structure; change-point detection; Ning-Meng reach of the Yellow River

1 引 言

凌汛是因冰凌對水流的阻水作用而導致江河水位抬升的一種自然現象，一般發生于冬季的封河期和春季的開河期。通常所說的凌汛水位包括封河期水位和開河期水位，凌汛流量指封凍期流量和開河期流量。凌汛水位、流量作為凌汛事件最直接觀測的水文變量，對于凌汛研究及凌洪減災有關鍵影響，其中凌汛極值水位、流量的影響最顯著。有學者對凌汛極值水位、流量做了相關研究。黃文元等[1]應用灰色自記憶模型對黃河內蒙古段巴彥高勒、三湖河口、頭道拐3個水文站的封河期最高水位建立模型并進行預測，得到較高的預報及擬合合格率。王富強等[2]統計分析了黃河寧蒙河段石嘴山、巴彥高勒、三湖河口和頭道拐4個水文站1951—2000年和2001—2010年的凌汛時間特征，并對開河期最大流量和槽蓄水增量進行系統分析，發現寧蒙河段開河期最大流量在2000年之后呈普遍減小趨勢。

受氣候變化和人類活動影響，世界很多地區的水文序列呈現出明顯的突變點或變化趨勢[3-6]。變化環境下寧蒙河段的凌汛水文變量發生了變化。李超群等[7]對寧蒙河段石嘴山、巴彥高勒、三湖河口、頭道拐4個水文站的凌汛水位與凌汛流量進行分時段統計，表明內蒙古河段的封、開河期最高水位上升，特別是巴彥高勒和三湖河口站的凌汛期最高水位上升明顯。水位上升將導致凌洪安全隱患的增加，因此水文序列的變異診斷顯得越來越重要。目前對于單變量水文序列的變點檢驗方法包括Mann-Kendall方法[8-10]、有序聚類分析法[10]、貝葉斯推斷方法[11-12]等。多變量聯合分布變點分析方法近來才引起人們的關注。Chebana等[13]采用非參數Mann-Kendall檢驗和Spearman秩相關檢驗檢測多變量水文序列的變化趨勢;Ben等[14]采用多元線性回歸中多變點檢驗的貝葉斯方法檢測雙變量序列相關結構中的變化點。這些方法可以檢驗多變量時間序列的變化趨勢或變化點，但沒有明確區分邊緣分布和多維聯合分布這兩種不同類型的變化。Copula似然比檢驗（CLR）方法避免了這一限制。該方法由Dias等[15]結合似然比檢驗理論與Copulas在2004年提出，并用于檢驗多變量金融序列相關結構的變化點。水文學中，Xiong等[16]采用Cramér-von Mises（CvM）方法和CLR方法分別檢驗多變量水文序列邊緣分布和相關結構中的變點，并將建立的模型結構用于中國漢江流域上游由年最大日流量、年最大3 d洪量和年最大15 d洪量組成的三變量洪水序列的變點分析。Huang等[17]采用CLR方法對渭河流域的降雨徑流雙變量序列進行了變點檢驗。

由上述研究可以看出，針對多維聯合分布的CLR變點檢驗方法的應用集中于洪水要素與降雨徑流研究，在凌汛方面的應用較少。凌汛的形成機制和影響因素與洪水及降雨徑流不相同，構建凌汛因子間的多維聯合分布模型，并對變化環境下的凌汛水文變量序列進行變異診斷顯得十分必要。因此，本文以黃河寧蒙河段三湖河口水文站的開河期最大流量為研究對象，通過相關性分析得到研究對象的影響因子，檢驗單變量序列的變化趨勢及突變點，在此基礎上利用藤Copula函數構建凌汛多變量序列間的多維聯合分布模型，并采用CLR方法對聯合分布模型相關結構進行變點檢驗，以此為多維聯合分布模型構建及變點分析在凌汛方面的應用提供參考。

2 方 法

2.1 Copula函數

Copula函數的核心是Sklar定理[18]。它是將多個隨機變量X1，X2，…，Xd的聯合分布F（X1，X2，…，Xd）與它們各自的邊緣分布F1（x1），F2（x2），…，Fd（xd）連接起來的函數，定義域為[0，1]。它的優點是可以將聯合分布分為邊緣分布和變量間的相關結構兩個獨立部分分別處理，且對邊緣分布的分布類型沒有要求，即任意邊緣分布經過Copula 函數連接都可以構建聯合分布，并且在轉換過程中不會產生信息失真[19]。Copula有很多函數族[20]，本文選用水文領域應用較多的Clayton Copula、Gumbel Copula和Frank Copula[19]作為備選函數構建聯合分布，并采用AIC信息準則進行Copula聯合分布優選[21]。

高維分布的構建通常選用Pair-Copula函數，將多維聯合密度分解為d（d-1）/2（d為聯合分布的維度）個二維Copula密度函數。Pair-Copula主要分為C藤和D藤兩種類型。C藤適用于多變量中存在一個關鍵變量的情況，D藤適用于多變量間關系相對獨立的情況[22]。三維情況下C藤有3種不同的分解形式，包含3個變量、2個tree結構和3條邊（edges）[22-23]，對應的聯合概率密度函數為式（1）（其中x1為關鍵變量）。

2.2 邊緣分布推求

用于單變量水文序列趨勢和變點檢驗的方法分別是Mann-Kendall檢驗[24]和Pettitt檢驗[25]。在d維隨機向量情況下，根據單變量序列Xm （m=1，2，…，d）變點分析的結果，相應的邊緣概率由經驗概率公式進行估計[16]。

如果單變量水文序列不存在變點，則邊緣分布為

2.3 CLR方法

CLR（Copula-based Likelihood-Ratio test）方法[16]即Copula似然比檢驗方法。假設Copula函數類型不變，則CLR方法通過測量給定的Copula參數變化來檢驗多變量相關結構的變化點。

似然比統計量Z1/2n的漸近分布由式（7）給出。

根據數據序列長度，選定顯著性水平為10%，由式（7）計算統計量Zn的邊界值為7.2。

3 研究區域及數據

黃河寧蒙河段是寧夏、內蒙古河段的統稱，位于黃河流域的最北端、黃河上游的末端。寧蒙河段干流在北緯37°17′—40°51′之間，河段高程在1 000 m以上，全長1 237 km，其中寧夏河段397 km、內蒙古河段830 km。寧蒙河段的氣候類型為大陸性氣候，年平均降水量155～366 mm，降水的年內和年際變化較大。河流由低緯度流向高緯度，冬季干燥、寒冷、漫長，極端氣溫為1988年1月1日頭道拐站的-39 ℃。河道形態復雜，淺灘彎道較多，主流擺動不定。寧蒙河段上游（寧夏河段）氣溫比下游（內蒙古河段）氣溫高，冬季下游先于上游封河，而在春季則是上游先開河，這將導致冰塞、冰壩等凌情的出現。在每年的11月到次年3月通常會發生凌汛，有時冰情嚴重導致凌災發生，因此黃河寧蒙河段是黃河防凌的重點河段。

寧蒙河段防凌的重要控制站包括石嘴山、巴彥高勒、三湖河口、頭道拐水文站等，其分布見圖1。

三湖河口水文站建于1950年8月，位于內蒙古自治區烏拉特前旗公廟鎮三湖河口村（東經108°46′，北緯40°37′），集水面積約為35萬km2。受地形影響，黃河在三湖河口站測驗河段上游形成一個大的轉折。該站設在這一大轉折的下游相對穩定的河道上，自設站以來，測驗河段始終處于變動之中。該站處于干旱半干旱地區，多年平均降水量237.8 mm，主要集中在6—9月;多年平均徑流量223.5億m3，約占黃河徑流量的39%;多年平均輸沙量1.08億t，約占黃河輸沙量的7%。自上游龍羊峽、劉家峽水庫運用以來，受兩水庫徑流調節的影響，汛期基本無較大洪水發生。

本文數據序列包括三湖河口開河期最大流量、開河期最高水位、封凍天數、封凍期來水量、累計負氣溫、最大槽蓄水增量、封凍期最高水位、冰厚，各序列對應年份均為1958—2005年。

4 結果與分析

三湖河口開河期最大流量的可能影響因素包括開河期最高水位、封凍天數、封凍期來水量、累計負氣溫、最大槽蓄水增量、封凍期最高水位、冰厚等。它們之間的Kendall秩相關檢驗結果見表1。

表1表明：開河期最大流量（D）與開河期最高水位（L）和冰厚（I）相關。進一步對L和I做相關性分析（τ=0.142，pτ=0.160）可知二者之間不相關。

4.1 單變量趨勢及變點檢驗

利用Mann-Kendall和Pettitt方法分別對D、L、I進行趨勢和變點檢驗，結果見表2。

由表2可知，在10%的顯著性水平下，開河期最高水位和冰厚呈顯著減小趨勢，而開河期最大流量無顯著變化。Pettitt檢驗的結果表明，開河期最高水位和冰厚存在明顯的突變點，它們對應的突變年份分別為1968年和1987年。考慮黃河上游大型水庫的修建及運行時間（劉家峽水庫1968年投入運行，龍羊峽水庫1986年開始蓄水運行，青銅峽水庫1978年竣工，《黃河劉家峽水庫凌期水量調度暫行辦法》（防總國汛[1989]22號）頒布時間為1989年），突變發生的原因可能與水庫調蓄有關。

基于變點檢驗結果的經驗頻率曲線見圖2～圖4。可以看出，開河期最高水位和冰厚變點前后序列的經驗頻率存在明顯差異，說明變點檢驗結果合理，開河期最高水位和冰厚變點后的值均小于變點前的，這可能是水庫調度及河道淤積所導致的[7]。

4.2 三變量C藤模型構建及變點檢驗

4.2.1 模型構建

開河期最大流量（D）的影響因素為開河期最高水位（L）和冰厚（I）。由于L和I不相關，因此在此以D為關鍵變量，構建D、L和I三者間的C藤Copula模型，模型結構見圖5。

基于4.1節單變量變點檢驗結果，利用2.2節介紹的方法推求D、L和I三變量的邊緣分布。采用Clayton、Gumbel和Frank Copula函數分別對C藤Copula模型tree 1結構中（D，L）和（D，I）的聯合分布進行極大似然參數估計，并用AIC信息準則比較3種Copula函數的擬合程度，從而選出最優Copula聯合分布類型，見表3。

由表3可知，對于（D，L），Clayton Copula對應的AIC值最小，因此最優Copula聯合分布為Clayton Copula;對于（D，I），Frank Copula對應的AIC值最小，因此最優Copula聯合分布為Frank Copula。

利用優選出的Copula函數類型構建C藤Copula模型tree 2結構，分析相關性發現兩條件分布（D與L的條件分布、D與I的條件分布）存在負相關關系。為了避免Copula類型差異對變點檢驗結果的影響，本文選用Frank Copula構建多變量序列的聯合分布。

4.2.2 相關結構變點檢驗

三變量Pair-Copula的結構分為兩個層次，每一層次的Copula參數描述了對應兩個變量之間的關系，因此多變量序列的變點檢驗可以分解為雙變量問題，包括兩個步驟：①檢驗tree 1結構（θ12、θ13）的變化點;②檢驗tree 2結構（θ23|1）的變化點。利用CLR方法對聯合分布相關結構進行變點檢驗，tree 1結構變點檢驗結果見圖6和圖7，tree 2結構變點檢驗結果見圖8。

圖6和圖7表明，開河期最大流量與開河期最高水位、開河期最大流量與冰厚間的聯合分布相關結構均不存在顯著變點。

由圖8可知，tree 2結構存在一個突變點，對應年份為1996年。龍羊峽、劉家峽水庫的聯合調蓄使洪峰流量和汛期水量削減，降低了洪水的輸沙能力，加上上游來水減少、支流來沙量增多使得內蒙古河道自1987年以來主河槽持續淤積，且1994年、1998年毛不浪孔兌、西柳溝、罕臺川相繼發生洪水，挾帶了大量泥沙進入黃河干流，以致內蒙古河段的淤積量增大[26]。1998年后，內蒙古河段從上至下巴彥高勒、三湖河口、頭道拐3個斷面的主槽均淤積抬升[27]。河道淤積使河槽變得寬淺散亂，輸冰輸沙能力下降，導致封、開河期水位偏高，從而在寬淺河段和河道轉彎處形成冰凌堵塞，1998—1999年封河期三湖河口段的封河水位為1 020.57 m，達到歷史較高水位。因此，1996年tree 2結構發生突變的原因可能是河道淤積。

4.2.3 相關結構變點前后對比分析

計算1996年前后（D，L）、（D，I）、（D，L，I）相關結構的Kendall相關系數（τ）與Copula聯合分布參數（θ），結果見表4。

由表4可知，變點后（D，L）和（D，I）變量間相關性增加，（D，L，I）變量間的相關性由正相關轉變為負相關。顯然，（D，L）和（D，I）相關性的變化改變了（D，L，I）整體的相關性結構。這說明受河道淤積與水利工程修建等的影響，開河期最大流量及其影響因素間的相關結構發生了較大變化。這一結論可為黃河寧蒙河段防凌工作及水利工程規劃實施和流域管理等提供參考。

5 結 語

以黃河寧蒙河段三湖河口水文站的開河期最大流量為研究對象，通過相關性分析得到研究對象的影響因素，在單變量變點檢驗結果的基礎上，采用C藤Copula函數構建多維聯合分布，并利用CLR方法對聯合分布相關結構進行變點檢驗。研究表明：河道淤積使得三湖河口水文站的凌情特征發生變化，導致開河期最高水位在1968年發生突變，冰厚在1987年發生突變;以開河期最大流量為關鍵變量的C藤Copula相關結構在1996年發生突變，這可能與內蒙古河段的主槽淤積抬升有關。

凌汛的影響因素復雜多變，在以后的研究中可以加入更多的可能影響因子（考慮更多氣象因子及除水電站工程外的人類活動因素）分析其對凌汛的影響，選擇并比較更多Copula函數構建模型，以期選出較適合凌汛研究的一整套模型結構。同時，探索新的研究方法對多維（三維及以上）聯合分布相關結構進行變點檢驗，并對比變點前后的序列特征。

參考文獻：

[1] 黃文元，樊忠成.灰色自記憶模型在黃河凌汛水位預報中的應用[J].人民黃河，2013，35（2）：3-5.

[2] 王富強，王雷.近10年黃河寧蒙河段凌情特征分析[J].南水北調與水利科技，2014，12（4）：21-24.

[3] DOUGLAS E M， VOGEL R M， KROLL C N. Trends in Floods and Low Flows in the United States： Impact of Spatial Correlation[J]. Journal of Hydrology， 2000， 240： 90-105.

[4] XIONG L H， GUO S L. Trend Test and Change-Point Detection for the Annual Discharge Series of the Yangtze River at the Yichang Hydrological Station[J]. Hydrological Sciences Journal， 2004， 49（1）： 99-112.

[5] VILLARINI G， SMITH J A， NAPOLITANO F. Nonstationary Modeling of a Long Record of Rainfall and Temperature over Rome[J]. Advances in Water Resources， 2010， 33（10）： 1256-1267.

[6] JEON J J， SUNG J H， CHUNG E S. Abrupt Change Point Detectionof Annual Maximum Precipitation Using Fused Lasso[J]. Journal of Hydrology， 2016， 538： 831-841.

[7] 李超群，劉紅珍.黃河內蒙古河段凌情特征及變化研究[J].人民黃河，2015，37（3）：36-39.

[8] KENDALL M G. Rank Correlation Methods [M]. 4th ed. London： Charles Griffin， 1970： 19-33.

[9] YUE S， PILON P， CAVADIAS G. Power of the Mann-Kendall and Spermans Rho Tests for Detecting Monotonic Trendsin Hydrological Series[J]. Journal of Hydrology， 2002， 259： 254-271.

[10] 張利茹，王興澤，王國慶，等.變化環境下水文資料序列的可靠性與一致性分析[J].水文，2015，35（2）：39-43.

[11] PERREAULT L， BERNIER J， BOBE B， et al. Bayesian Change-Point Analysisin Hydrometeorological Time Series：Part 1. The Normal Model Revisited[J]. Journal of Hydrology， 2000， 235（3）： 221-241.

[12] 熊立華，周芬，肖義，等.水文時間序列變點分析的貝葉斯方法[J].水電能源科學，2003，21（4）：39-41.

[13] CHEBANA F， OUARDA T B M J， DUONG T C. Testing for Multivariate Trendsin Hydrologic Frequency Analysis[J]. Journal of Hydrology， 2013， 486： 519-530.

[14] BEN A M A， CHEBANA F， OUARDA T B M J， et al. Dependence Evolutionof Hydrological Characteristics， Applied to Floods in a Climate Change Context in Quebec[J]. Journal of Hydrology， 2014， 519： 148-163.

[15] DIAS A， EMBRECHTS P. Change-Point Analysisfor Dependence Structures in Finance and Insurance[G]// Risk Measures for the 21st Century. Chichester： Wiley Finance Series， 2004： 321-335.

[16] XIONG L H， JIANG C， Xu C Y， et al. A Framework of Change-Point Detection for Multivariate Hydrological Series[J]. Water Resources Research， 2015， 51（10）： 8198-8217.

[17] HUANG S Z， HUANG Q， CHANG J X， et al. Variations in Precipitation and Runoff from a Multivariate Perspective in the Wei River Basin， China[J]. Quaternary International， 2017， 440： 30-39.

[18] NELSEN R B. An Introduction to Copulas[M]. New York： Springer，2006： 14-20.

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

[20] GENEST C， FAVRE A C. Everything You Always Wantedyo Know About Copula Modeling but Were Afraid to Ask[J]. Journal of Hydrologic Engineering， 2007， 12（4）： 347-368.

[21] 侯蕓蕓，宋松柏.基于Copula函數的洪峰洪量聯合分布研究[J].人民黃河，2010，32（11）：39-41.

[22] 宋松柏，蔡煥杰，金菊良，等.Copulas函數及其在水文中的應用[M].北京：科學出版社，2012：208-222.

[23] AAS K， CZADO C， FRIGESSI A， et al.Pair-Copula Constructions of Multiple Dependence[J]. Insurance： Mathematics and Economics， 2009， 44（2）： 182-198.

[24] LIBISELLER C AND GRIMVALL A. Performance of Partial Mann-Kendall Testsfor Trend Detection in the Presence of Covariates[J]. Environmetrics， 2002， 13： 71-84.

[25] PETTITT A N.A Non-Parametric Approach to the Change Point Problem[J]. Applied Statistics， 1979， 28： 126-135.

[26] 葉春江.黃河寧蒙河段河道整治的實踐與研究[D].西安：西安理工大學，2003：12-14.

[27] 金雙彥，翟家瑞，蔣昕暉，等.凌汛期頭道拐小流量過程變化及影響因素分析[J].人民黃河，2012，34（12）：27-29.

【責任編輯 許立新】