摘 要：
水文系列的变异性检验是非一致性水文频率分析过程中的重要环节，其中跳跃性变异是非一致性的重要表现形式，准确识别突变点对认识水文过程发生的变化及开展实际水文水资源工作具有重要意义。对于多维时间序列而言，可能会出现各单维变点不统一的情况，为此，将启发式分割算法(BG算法)与线性加权综合统计量结合，应用于多维水文时间序列的变异点检测，通过综合各维度的统计结果，确定最终的突变点。统计试验结果表明：该法有利于排除多维系统中虚假突变点或弱突变点的干扰，得到更为合理的检测结果。对汉江黄龙滩入库洪水极值系列变异性检测的应用结果表明，洪峰和最大7 d洪量的单变量检测突变点分别为1990年和1985年，通过综合统计量判定最终的突变点为1985年，为解决多维水文时间序列突变点的统一检测问题提供参考。
关键词：
非一致性;突变点检测;多维时间序列;启发式分割算法;
作者简介：
刘杨(1997—),女，硕士研究生，主要从事水文水资源研究。
*梁忠民(1962—),男，教授，博士，主要从事水文水资源研究。
基金：
国家自然科学基金项目(51709073);
中央高校基本科研费专项资金(2019B03214);
引用：
刘杨，梁忠民，罗序义，等． 多维时间序列突变点检测方法研究［J］. 水利水电技术( 中英文) ，2022，53( 5) : 65-72.
LIU Yang，LIANG Zhongmin，LUO Xuyi，et al． Study on change point detection method of multidimensional time series［J］. Water Ｒesources and Hydropower Engineering，2022，53( 5) : 65-72.
0 引 言
现行水文频率分析方法假设水文系列满足一致性，即水文极值系列的概率分布或统计规律在过去和未来保持不变。但随着全球气候变化及人类活动对水文过程影响的日益显著，众多研究表明该一致性假设遭到一定破坏，导致传统基于一致性假设的水文频率分析方法计算结果可靠性受到质疑,因此研究非一致性条件下水文频率分析方法具有重要意义。
水文系列的变异性检验是非一致性水文频率分析过程的重要环节，其中，跳跃性变异是水文样本系列非一致性的一种重要表现形式，一般定义为系列从某点(年份)开始突然或急剧地从某种状态变化到另外一种状态。因此，对跳跃性系列突变点(变异点)的识别和检测是关键所在。目前国内外水文学者提出了许多突变点识别和检测方法，如Lee-Heghinian法、Pettitt检验法、有序聚类法、滑动T检验法、滑动F检验法、滑动秩和法、R/S检验法、贝叶斯法和Mann-Kendall检验法等。一些研究者对各方法的特点、适用性、优缺点进行了比较分析，如雷红富等采用统计试验的方法比较分析了10种常用变异检验方法在3类变异情形下(均值变异类、CV变异类、CS变异类)的检验性能；管晓祥等比较了4种检验方法对于武隆站水文资料突变点检验的性能，结果表明，M-K突变点检验法性能最优；田小靖等选用9种水文序列突变点识别方法，比较了其对于长时间输沙数据序列的检验性能。但此类比较研究结果各异，尚未能获得在多数情形下均表现优异的方法，其原因是这些水文序列变异点检验方法的假设前提和适用条件各不相同，因此在实际应用中各种方法的结果以及性能往往会存在差异，需视数据的具体特性选用相对较适用的突变点检验方法。
上述方法一般针对单维序列，然而，水文事件实际上具有多方面的特征属性，例如洪水事件包含洪峰流量、时段洪量和洪水过程线三要素，干旱事件一般也需考虑干旱烈度、历时和范围等要素，因此传统的单变量分析无法全面反映水文事件的真实特性。为更全面描述水文事件的内在规律，近些年关于多变量极值系列的非一致性检验已受到越来越多的重视，常见的多变量跳跃性变异检验方法有自适应检验方法、基于Kolmogorov-Smirnov统计量的非参数检验法、基于Cramer-von Mises(CvM)统计的非参数检验法等。 AISSIA等基于多元线性回归中多变点检验的贝叶斯方法分析了具有相关性的两变量洪水序列的变点。XIONG等提出了一个基于Copula函数的多变量水文序列突变点的诊断框架，即首先通过对单个水文变量进行非一致性诊断，然后采用合适的Copula函数构建相关结构，在此基础上基于似然比检验对相关结构的突变点进行诊断。但这些方法尚未能很好地解释各维变点不统一的问题，即对于描述同一水文事件的多维时间序列，各单维度的变点检测结果有可能不一致，给后续的非一致性频率分析造成困难。
为此，本文将BERNAOLA-GALVAN等的启发式分割算法(以下简称BG算法)与高桢的线性加权综合统计量相结合，应用于多维水文时序的突变点检测；并设计了统计试验方案以说明该方法的合理性，通过黄龙滩水库入库洪水年极值系列的示例分析，进一步检验其可行性。
1 研究方法1.1 启发式分割算法
2001年，BERNAOLA-GALVAN等提出了启发式分割算法(BG算法),该方法基于滑动t检验的思想，且在分割时采用了二分的迭代算法，能将一个非平稳的时间序列分割成多个平稳的子序列，各序列的均值互不相同并代表不同的物理背景，大大减少了计算量，实用性较好。陈广才等将BG算法与10种常用突变点检验方法进行了比较，结果表明该方法能较好排除虚假突变点干扰，较准确识别突变点数目及其位置，检验性能优于其他方法，适用于水文变异分析。
设时间序列X(t)由N个点组成，从序列的左端向右端滑动选取分割点i,分别计算分割点左右两侧子序列的均值μ1(i),μ2(i)及标准差s1(i),s2(i),则分割点的合并偏差SD(i)可表示为
式中，N1、N2分别为分割点左右两侧子序列的长度。
采用t检验的统计值T(i)度量分割点两侧子序列均值的差异，计算公式为
对时间序列中的每个点重复上述计算，得到与之相对应的检验统计序列T(t)。其中，T值越大，表示该点两侧子序列均值差异越明显。计算T(t)中最大值Tmax对应的统计显著性P(Tmax),即事件T≤Tmax的概率为
一般情况下，P(Tmax)可近似表示为
式中，变量均为由蒙特卡洛模拟得到的经验公式 ,η=4.19lnN-11.54,v=N-2,δ=0.40,其中Ix(a,b)为不完全β函数，表达式为
实际应用时，首先设定一个临界值P0,如果P(Tmax)≥P0,则在该点处将X(t)分割成两个均值有一定差异的子序列，否则不进行分割。同理，对得到的新序列重复上述操作，直至子序列的长度小于等于l0(l0为最小分割尺度),则不再对其进行分割。通过上述步骤，可将原时间序列X(t)分割成若干均值不同的子序列，而分割点即为该序列的突变点。通常，P0的取值范围为0.5～0.95,l0≥25,通过改变给定的P0和l0,可以实现对序列不同尺度的变异检验。
1.2 多维时间序列突变点检测
多维时间序列的突变点检测与单维时间序列的分析方法不尽相同，考虑到不同维时间序列包含的信息对整体的重要程度有差别，为综合考虑各维度对整体的影响，本文首先基于启发式分割BG算法对多维时间序列中的每一维进行单独分析，即在每一维度上计算对应的统计序列T(t),然后综合各维度统计结果，最终确定多维系统的综合突变点。
对于各维统计结果的比较与综合，考虑两种思路：第一种是每一维度上分别确定突变点，通过显著性比较或专家经验等方式选择某一维度上的突变点作为整个多维时间序列最终的突变点，此方法考虑了不同维度对整体的影响，但是有可能导致突变点“过检验”的问题。第二种是将各个维度上相同时刻的统计量进行综合，得到一维的综合统计量，在此基础上确定综合突变点。为了避免单维时间序列上虚假突变点或弱突变点的干扰，且更好地综合各维检测结果，本文采用第二种思路确定最终的突变点，并在1.3节通过统计试验来说明该思路的合理性及可行性。
首先计算K维序列上每一点统计量的平均值Tavettave和最大值T maxt为
理论上，不同维度的统计量对整体的重要性可能是有所不同的，Tavettave可能掩盖了某一维度上的明显变化，而Tmaxttmax易受某一维度上噪声的影响，因此若对Tavettave和Tmaxttmax直接取平均可能会过于简化整体结构。为了解决上述问题，高桢提出对Ttave和Ttmax线性加权以表示综合统计量
通过上式得到的综合统计量序列，既能考虑不同维度上统计量的差异，又减少了噪声的干扰，可较好地反应多维时间序列的整体变异情况。加权因子α的取值可视实际情况进行调整，当时间序列混杂着大量的噪声时，可适当减小α,增强公式的平滑作用，此时α的取值范围为[0,0.5);若需强调不同维度上变化的差异，可考虑增大α,从(0.5,1]中进行取值。
1.3 统计试验方案
本文设计了一个统计试验，将利用上述检测方法得到的突变点与已知的真实突变点进行比较，并设计三种不同的情形来说明该方法在多维时间序列突变点检测中的合理性与可行性。由于P-Ⅲ分布对于我国大部分河流的水文资料拟合较好，因此，基于P-Ⅲ分布分别生成三种情形下的三维随机序列，令xab(t)=f(X0,CV,CS),其中，xab(t)为服从P-Ⅲ分布的随机数，a表示情形编号，b表示序列维数，t为生成的随机数个数，X0,CV,CS分别为均值、变差系数和偏态系数。
1.3.1 情形一
令该三维序列均发生显著变异，且突变点位于相同位置，即生成长度为2 000,理论突变点均为1 000的随机序列
采用上述方法进行突变点检测，结果如图1所示，其中图1(a)中折线为随机序列的过程线，竖线指示突变点所在位置，突变点前后水平线为前后子序列的均值；图1(b)为三个序列的统计量过程线。BG算法检测三个序列的突变点位置分别为1 000,1 001,1 002,基本符合理论突变点位置，产生的微小偏差是生成样本时的随机误差所导致的。分别取权重α=0.4,0.5,0.6,计算三维序列的变异性综合统计量如图2所示，结果显示，无论α取何值，突变点均位于1 001处，符合理论突变点位置。
图1 情形一样本序列及统计量过程线
图2 情形一综合统计量过程线
1.3.2 情形二
控制突变点位置不变，改变变异显著性，即令该三维序列的序列1发生显著性变异，序列2与序列3发生显著性较弱的跳跃性变异，三个序列均在1 000处发生变异，生成随机序列
检测结果如图3所示，图中说明参照情形1。BG算法检测三个序列的突变点位置分别为1 003,1 001,993,基本符合理论突变点位置。分别取权重α=0.4,0.5,0.6,计算三维序列的变异性综合统计量如图4所示，结果显示，无论α取何值，突变点均位于1 003处，基本符合理论突变点位置，且结果较为稳定，减小了单维序列的随机误差影响。
图3 情形二样本序列及统计量过程线
图4 情形二综合统计量过程线
1.3.3 情形三
改变变异显著性及突变点位置，即令该三维序列的序列1发生显著性变异，序列2与序列3发生显著性较弱的跳跃性变异，三个序列分别在1 000,600和1 400处发生变异，生成随机序列
检测结果如图5所示，图中说明参照情形1。BG算法检测三个序列的突变点位置分别为1 000,599,1 402,基本符合理论突变点位置。分别取权重α=0.4,0.5,0.6,计算三维序列的变异性综合统计量如图6所示，结果显示，无论α取何值，突变点均位于1 002处，基于变异性综合统计量得出的结论表明，该综合检测方法可以突出显著性较强的序列对多维时间序列整体的影响，而排除虚假突变点或显著性较弱的单序列突变点。
图5 情形三样本序列及统计量过程线
图6 情形三综合统计量过程线
1.3.4 结果分析
以上三种情形下的统计试验结果表明，无论哪种情形，无论α取何值，基于启发式分割BG算法的多维时间序列突变点检测方法找到的突变点位置与理想突变点位置基本相同，试验中产生的微小偏差是生成样本时的随机误差造成的，因此该方法对多维时间序列的变异检测是行之有效的。对于复杂的多维系统，该法通过计算综合变异性统计量，可以提供更加合理的统一结果，有利于排除虚假突变点或弱突变点的干扰，结论更为可靠。
2 实例分析2.1 研究区域及数据
本文研究对象为黄龙滩水库的入库年极值洪水序列。黄龙滩水库于1978年建成，位于堵河下游湖北境内，水源均来自于堵河，上游控制流域面积约11140 km2,占堵河流域面积的95%左右(见图7)。堵河是汉江最大的支流，流域属于亚热带半湿润季风气候，雨量丰沛，多年平均年降水量约1000 mm, 但时空分配不均匀，在空间分配上呈现出上游少下游多的态势，在时间分配上呈现出汛期(6—10月)较为集中、频繁，且暴雨持续时间较长的特点。本次研究的多维时间序列选用1956—2005年黄龙滩水库上游控制流域的洪峰以及最大7d洪量系列。
图7 黄龙滩水库控制流域位置
2.2 突变点检测
采用上述基于启发式分割算法的多维时间序列突变点检测方法对洪量-洪峰二维序列进行突变点检测。洪峰以及最大7 d洪量的时间序列如图8所示，首先采用BG算法分别检测各维的突变点，从图中可以看出，洪峰流量突变点为1990年，最大7d洪量突变点为1985年，对于描述同一洪水事件的峰量二维序列，它们的突变点却并不在同一年，那么该洪水事件究竟是在哪一年发生变异的呢?为此，计算洪峰-洪量的综合统计量，为综合考虑各维变量对整体的影响，取权重α=0.5,绘制各自及综合的统计量过程线如图9所示，图中菱形数据点标志为二维序列的综合突变点，三角形数据点标志为洪峰单维度上检测出的突变点(洪量单维度突变点与综合突变点重合),从图中可以看出，由综合统计量得出的突变点为1985年，因此判定为该洪量-洪峰二维序列在1985年发生了变异。
图8 洪峰-洪量时间序列
图9 洪峰-洪量综合统计量过程线
对这一结论进行合理性分析，从图9可以看出，洪峰与洪量的统计量具有较强相关性，Pearson相关系数为0.826,说明两者的变异特性相差不大，这是符合自然规律的。洪量-洪峰二维序列的突变点为1985年，与洪量单维度的检测结果一致，而洪峰突变点并非1985年。对洪峰序列的变异性统计量规律进行分析，可以看出，洪峰单维度的T统计量在1985年为局域极大值，即若仅从洪峰序列单维度进行分析，1985年是一个稍次于1990年的突变点，且也具有一定程度的变异显著性，这与综合变异点为1985年的结论不矛盾，因此，该综合统计量的检测结果是合理可信的。综上，对于同一水文事件，不宜仅从描述该事件的某单一要素来判别整个事件的变异特性，而应结合各要素进行综合分析，以对水文事件整体得出最终的结论。
2.3 结果讨论
本文通过综合统计检验判定黄龙滩水库的入库洪峰-洪量二维序列在1985年发生了变异，该结果与牛利强、黄萱、梁小青等对堵河流域内若干站点非一致性的研究结论基本一致。其中，梁小青等采用滑动t检验法对竹山站年平均径流序列进行检验(竹山水文站位于黄龙滩水库以上92 km),结果表明，竹山站径流量在1985年前后发生变异，整体上呈现先上升后下降的趋势。
目前黄龙滩水库所在的堵河流域的骨干水库格局已基本形成(见图7),堵河流域上游共有5座骨干水库，包括干流的龙背湾、松树岭、潘口以及一级支流泗河上的鄂坪和竹溪河水库。除了上述骨干水库外，在80年代前后流域内还兴建了众多的中小型水库。这些水库的建设运行，对位于流域最末级的黄龙滩水库的入库流量必将产生较大影响，是导致其洪水极值系列发生变异的主要原因；另外，自上世纪70年代末改革开放以来，堵河流域的土地利用也发生较大变化。这些原因均会对黄龙滩水库水文系列的非一致性产生影响。
3 主要结论多维时间序列的突变点检测与单维序列的突变点检测有所区别，且各单维度的突变点检测结果有可能存在一定的差异，难以得到统一的结论。本文将启发式分割BG算法与线性加权综合统计量结合，构建一种多维水文时序变异点的检测方法，主要结论如下：
(1)统计试验结果表明，该方法检测出的突变点位置与真实突变点位置基本相同，试验中产生的微小偏差是样本生成的随机性造成的，表明该方法是合理可行的。
(2)对黄龙滩水库入库洪水的示例应用结果表明，洪峰和最大7 d洪量系列单变量检测的突变点分别为1990年和1985年，通过变异性综合统计量判定该二维序列最终的突变点为1985年，解决了单要素突变点不统一的问题。
本文在进行统计试验时，仅考虑均值变异条件下生成序列的变异点检测问题，对于因CV、CS或多种变异相互叠加情形的变异点检测，尚待进一步研究。
水利水电技术（中英文）
水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。
