1 引言

突变理论最早可追溯至20世纪60年代末期,法国数学家R.Thom创立了以常微分方程为数学基础的突变理论^[1],其要点在于考察某种系统或过程从一种稳定状态到另一种稳定状态的飞跃^[2]。随后,突变理论在数学、生物、天文、气象、社会科学等领域得到了广泛的应用^[3,4,5,6,7]。从统计学的角度,可以把突变现象定义为从一个统计特性到另一个统计特性的急剧变化,如考察均值、方差等特征量的急剧变化^[8]。目前,常见的气候突变是指气候从一个平均值到另一个平均值的急剧变化,称为均值突变^[9]。迄今为止,针对均值突变的情形,已经发展出很多突变检测的方法,包括较为传统的低通滤波法、滑动t检验法、Cramer法、Yamamot法,以及后来被提出并广泛使用的Mann-Kendall突变检验法^{[10,11,12,13,14]}、Pettitt法、有序聚类分析（OC）法^[15]和BG（Bernaola-Galvan）分割法等^[16,17]。传统的4种方法都因检测结果带有主观性,所以缺乏可信度。而后4种方法虽解决了传统检测方法突变点漂移和受主观因素影响的缺点^[3],但统计量结构的合理性有待讨论,另外,这些方法都是基于序列只有一个突变点的假设条件,但是气候序列是多尺度的,大小涨落（距平）都有,且具有周期性变化的特点,因此,这些方法很难找出不同尺度和层次上的所有突变点,进而影响对突变成因的分析。因此,本文提出一种新的时间序列突变点的检测方法,并与现有常用的Mann-Kendall突变检验法、Pettitt法、OC法和BG分割法4种方法相比较,证明了该方法有诸多优点且行之有效。利用这种方法分析黄河中上游输沙量序列资料,探讨黄河泥沙的突变时间及可能的影响因素,为黄河水沙治理提供支撑。

2 滑动平均差检测法

时间序列突变点的检测,均以样本观测序列为基础,根据某种统计假设,计算其检测统计量,进而确定其突变点。这过程都会受到3个关键环节因素的影响：① 样本观测序列获得过程产生的误差,这包括基本资料观测误差、样本选择误差等。这类误差随机过程一般都具有零均值和常方差的平稳过程特征。② 时间序列本身所具有的周期性变化。这变化受时间变量本身的常规物理属性驱使,在一定的周期内有确定性变化规律,而周期与周期之间的均值接近常数。当时间变量的某个物理属性驱动因子出现异常,导致这周期均值出现显著改变时,就是时间序列要检测的突变点。③ 当检测时间序列具有反复变化（正向和反向交替出现或连续多个突变）的多个突变点时,如何能有效地检测出这些突变点,是突变检测方法的难点所在。

为解决第一个因素的影响,样本序列误差消除的有效方法是期望均值计算,而解决第二和第三个环节问题的有效手段是把周期均值作为考察变量,三个环节因素综合兼顾考虑,选择以时间变量常规物理周期为滑动期的滑动平均变量作为突变点的新考察序列是较为理想的。如时间序列X_i（i=1, 2, …, n）变量的常规物理周期为p,则可以构建正向滑动平均序列：

MUi=1k∑j=1kXi-ji=2,3,⋯,n;k=min(p,i-1)（1）

类似地可以构建逆向滑动平均序列：

MDi=1k∑j=1kXi+j-1i=2,3,⋯,n;k=min(p,n-i+1)（2）

由式（1）和式（2）可知,MU_i和MD_i即第i个样本点前后两个子序列的均值。考察一个突变点的最简单理想时间序列（没有观测误差、突变点前后都为常数）,如：

Xi=10,1≤i≤250,26≤i≤50（3）

计算其正逆向滑动平均序列得：

MUi=10,2≤i≤2636-i,27≤i≤360,37≤i≤50MDi=10,2≤i≤1626-i,17≤i≤260,27≤i≤50（4）

其正逆向滑动平均序列差为：

ΔMi=MUi-MDi=0,2≤i≤26i-16,17≤i≤2636-i,27≤i≤360,37≤i≤50（5）

把序列（3）和（5）作图,如图1所示。发现正逆向滑动平均序列的最大差值正对应的点就是突变点I_cr,其含义为该点前后均值差值达到最大,而且其突变强度为其最大差值。因此,统计量ΔM_i的极大值点表示该样本点前后子序列均值发生了急剧变化,即极大值点有可能是突变点,由此提出时间序列突变点的滑动平均差位置检测指标为：

Icr=maxiΔMi（6）

和强度检测指标M_cr：

Mcr=maxΔMΔMiMUIcr（7）

3 不同检测方法理想时间序列突变点检验比较

由上节介绍可知,滑动平均差检测方法的结构很简单,其几何解释也非常直观易理解,但方法的合理性和准确性还有待检验。为此,这里把滑动平均差检测方法与目前较为常用的BG分割法、Mann-Kendall突变检验法、Pettitt法和OC法4种方法应用于理想序列进行突变点检测比较,分别从无误差单突变点、无误差多突变点和有误差多突变点3种情况进行检验分析。在理想序列的检验中,考虑到所构造的理想序列的特点,滑动平均法的滑动周期取p = 10。

3.1 无误差单突变序列检验

无误差单突变序列如式（3）所列样本,分别用5种检测方法进行检验,结果如表1、图2和图3所示。为了把不同检测指标序列显示在同一张图中,图中的样本序列和前4种方法的检测指标都进行了标准化转换。即：

xi'=xi-minximaxxi-minxi（8）

从表1和图2、图3结果可知,5种方法除Mann-Kendall突变检验法外,都准确地找出了突变点发生的位置,滑动平均差方法还准确地算出了突变强度。而Mann-Kendall突变检验法的两个序列根本就没有交点,所以没有探测出序列的突变点。

3.2 无误差多突变序列检验

无误差多突变点的序列,设计有5个突变点,其突变点发生的位置和各点的突变强度如表2所示,5种方法检测结果如表2和图4所示。

图4中的样本序列和前4种方法的检测指标都进行了标准化转换。由表2和图4结果可知,滑动平均差方法都一次检测出了所有突变点和其突变强度,而且位置和强度都没有误差。BG法和OC法检测出1和5两个突变点位置,Pettitt方法从统计量曲线上看有1、3和5三个局部极值点,但3突变点的信号较弱,没有通过其显著性水平检验,而对于2与4两个突变点BG法、OC法和Pettitt方法均没有检测出。按照这些方法设计,依次去除已检测出突变点的序列,再重复检测剩余序列,能检测出剩余的突变点。具体结果类似于3.1节的单突变点检测结果。同样的也只有Mann-Kendall突变检验法,对所有5个突变点都没有能检测出。

把表2各突变点突变强度与各方法突变点检测值进行比较,发现各方法突变点检测值都与其突变强度真值成正比关系（图5）。为把不同方法的检测值在同一图中方便表达,图5中进行了检测值数值转换,其转换公式为：

MB=6×BG/60.3（9）

MO=6×68.6/OC-400（10）

MP=6×Pettitt-400/637（11）

各方法的检测值与实际突变强度比较,发现滑动平均差法的检测值等于其实际突变强度,方法构造信号检测值过程中没有损失其有效信息,所以所有突变点都准确无误的检测到了。而其他3种方法,构造信号检测值过程中都损失了有效信息,而且突变强度越低相对损失越大,所以对2和4两个突变点的突变强度低些,其有效信息就被全部损失而无法检测到了。这是这3种方法存在的共同缺点。

分析BG法、OC法和Pettitt法3种方法有效信息损失的原因,笔者认为主要是方法本身检测指标结构的不合理引起的。3种方法都有一个暗指的假设,即只有一个突变点,且突变点前和后为两个均值不同的平稳序列,所以当有多突变点时,突变强度低的点其信息被强度高的点所“掩盖”而不起作用。特别如BG分割法,检测指标式为

MO=μleft-μrightsd（12）

sd=sleft2+sright2nleft+nright-212×1nleft+1nright12（13）

式中：μ_left和μ_right分别为每个点左边子序列的平均值和右边子序列的平均值^{[3, 16]}。

其中,结构 1nleft+1nright12的最大值发生在两端（n_left×n_right为最小值时）、最小值发生在中间（n_left×n_right为最大值时）,与突变点毫无关系,所以不仅不能挖掘突变点的有效信息,而且还会提供与突变点无关的有害信息。

3.3 有误差多突变序列检验

有误差多突变序列,其真值序列和突变点同上,加入随机误差序列后,其突变点位置与突变强度如表3所示。

滑动平均差检测法、BG分割法、OC法、和Pettitt法4种方法检测结果与无误差多突变序列相差无几,依旧只有滑动平均差检测出了所有突变点且同时准确计算出了突变强度,而其他3种方法没有检测出突变强度较弱的3个突变点。Mann-Kendall突变检验法检测结果差别较大,该方法检测出序列的突变点在85~86位置,但这个突变点显然不是检测序列实际的突变点。

4 不同检测方法的实际时间序列检验结果分析

4.1 研究区及资料概况

20世纪70年代以来,由于气候变化和人类活动的双重影响,黄河水沙条件发生了巨大变化,黄河水沙特别是泥沙量发生趋势性锐减。潼关站平均输沙量从天然情况的16亿t减至不足1亿t^{[18,19,20,21,22,23]}。依地理位置和河流特征,黄河可分为上、中、下游三段,其中由内蒙古托克托县的河口镇至龙门为中游段,流经黄土高原和丘陵区,水土流失尤为严重,是黄河泥沙的主要来源区,因此,本文选取黄河中游干流控制站头道拐、龙门和潼关,窟野河控制站温家川、汾河控制站河津以及渭河控制站华县6个水文站点的年输沙量序列,检验不同方法在实际序列中的检测效果。站点基本情况介绍如表4所示,研究区和站点位置如图6所示。

4.2 潼关断面年输沙量突变点检测

年输沙量时间序列通常都包含年输沙量观测过程产生的误差和样本序列的年际周期变化。降雨等气候因子、植被变化和人类活动等是引起水文序列变化的主要因素,为了研究人类活动和植被变化对水文序列的影响,准确检测出水文序列的突变点,应尽可能减小或抵消降雨等气候因子的周期震荡对水文序列突变点检测的影响。因此,滑动平均差检测法的滑动周期应取气候周期,以消除气候震荡带来的误差。滑动周期太短抵消不了样本序列的观测误差和年际波动,容易造成突变点的漂移,而滑动周期过长,会导致序列过渡平滑,淹没突变点。结合黄土高原地区地形及气候变化周期特点,本研究中滑动平均差检测法的滑动周期取p = 16 a,潼关断面1950-2015年输沙量序列和不同方法的计算结果如图7所示。

从图7结果看,滑动平均差法检测出两个显著的突变点,分别在1979年和2003年,对应的突变强度分别为6.528和5.467。OC法和Pettitt法检测出突变点在1979年,BG法方法检测出的突变点在1996年,Mann-Kendall法UF和UB曲线在2000年有一个交点,但该交点不在0.1置信区间内,因此该交点不能算显著突变点。但从BG法、OC法和Pettitt法4种方法统计量曲线来看,曲线都有两个明显的极值点,第一个极值点在1979年,据年输沙量序列和滑动平均差检测结果可知,1979年前后的突变强度最大,所以该突变点的检测结果较为一致和可信;第二个极值点在1996年,从滑动平均差法曲线可以看出1996年也是一个局部极值点,突变强度为4.868,相较于1979年和2003年,该点的突变强度较弱。BG法检测出的最大突变点不是发生在均值变化最大的1979年,而是发生在突变强度很弱的1996年,且2003年应该是突变强度仅次于1979年的突变点,但BG法、OC法和Pettitt法均没有检测出。这均是由“序列只有一个突变点”的暗指假设与实际不符所导致的。综上,潼关站输沙量序列有两个显著突变点,分别发生在1979年和2003年。

4.3 其他序列突变点检测检验

头道拐、龙门、温家川、河津和华县五个断面年输沙量序列各方法的检测结果如表5所示。分析表5可得,① 当时间序列有且仅有一个突变点,且为显著突变（如头道拐年输沙量序列）时,5种方法均可检测出该突变点,且结果基本一致;② 当突变点为两个或两个以上时,只有滑动平均差检测法能同时一次性准确检测出所有突变点及相应的突变强度,BG法、OC法和Pettitt法检测出的突变点也都是滑动平均差法检测出的突变点之一,但不一定是最大突变强度的突变点,甚至出现检测出的突变点是几个突变点中突变强度最弱的情况,即3种方法都存在着突变点漂移的问题;③ Mann-Kendall方法在序列存在多个突变点的时候,统计曲线在其置信区间内无交点,找不到突变点,方法几乎崩溃。

黄河流域上游输沙量在1986年发生显著突变,20世纪80年代,黄河上游生产生活用水大大增加,比50年代增长约52%^[24],且20世纪60年代后,黄河上游干流水利工程的修建极大地改变了上游水沙的年际和年内分配,如1986年龙羊峡水库建成后,与刘家峡水库联合调度,拦蓄了大量的径流和泥沙^[25],使上游输沙量大大降低;黄河中游段输沙量主要有两个突变点,分别在1979年和1996年,黄河中游段流经黄土高原,20世纪70年代后黄土高原大规模的水土保持措施对水沙的减少有很大的促进作用^[26],90年代初期,黄土高原地区进一步展开水土保持二期治理工程（包括退耕还林、封山禁牧和大量的植树造林）,使得植被覆盖度进一步增加,水沙进一步锐减^[27]。渭河流域在1979年、1996年和2003年发生突变,且2003年和1979年为主要突变点,这一结果与潼关站检验结果相一致。潼关站在黄河流域干流的下游,来水来沙由龙门、汾河、北洛河和渭河共同组成,各种人类活动的综合作用是黄河流域输沙量减少的主要原因^[28]。

5 结论

（1）本文提出了一种新的突变点检测方法——滑动平均差检测方法,该方法能同时检测出时间序列的多个突变点及相应的突变强度。将该方法运用到理想时间序列和实际观测时间序列中进行检验,并与现有常用的4种突变检验方法的检测结果作比较,发现传统4种方法存在如下不足：① 4种方法都有“序列只有一个突变点”的默认假设,当实际情况与假设条件不符时,则统计量计算会损失有效信息而使突变强度较低的突变点无法检测出;② BG法的指标结构不太合理。而滑动平均差检测方法比较这4种检测方法具有3个明显的优势：一是结构简单、直观易理解;二是检测突变点更精确;三是能同时检测出所有突变点的突变位置和突变强度,这为流域水沙突变点检测提供了较为科学、准确的方法,值得推广使用。

（2）对黄河流域几大重点水文控制站输沙量资料进行突变分析,结果显示：黄河流域上游输沙量在1986年发生显著突变;中游段输沙量在1979年和1996年发生突变,其中1979年为主要突变点;渭河流域和潼关站输沙量有两个主要突变点,分别在1979年和2003年,大型水利工程建设和大规模水土保持措施的实施等人类活动是导致黄河流域泥沙突变的主要因素。

The authors have declared that no competing interests exist.

---

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子