基于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.

【责任编辑 许立新】