近似熵和数据滑动技术在径流突变检测中的运用

1、近似熵和数据滑动技术在径流突变检测中的运用0 引言径流时间序列突变分析是变化环境下水文分析计算的基础,是当前一个热点前沿问题.由于全球气候变化和人类活动的影响,流域下垫面发生了变化,水文系统循环时空范围与尺度发生了改变,造成了径流时间序列的时空变异,带来了一系列诸如水文极值事件、洪涝灾害、严重干旱频繁发生等新的问题,破坏了水文序列的一致性,颠覆了工程水文学的计算基础,给水利工程规划、设计及运行带来风险.因此,研究变化环境下的径流序列变异规律,对于水资源高效持续利用、兴利减灾、经济社会的稳定发展等具有重要的现实意义.径流时间序列的变异研究是对流域水文序列的变异点的识别与检验,目前主要采用数学统计

2、的方法,主要包括 Mann-Kendall 法1、两阶段线性回归2、有序聚类法3、滑动 F 检验4、去趋势波动分析5等方法,但这些方法主要从线性统计的角度出发,对于变化环境下的小数据量数据集以及高度复杂、动态、非线性的径流时间序列检测能力不足,对变异前后序列的状态趋势不能有效识别,无法有效揭示变化环境下的水文系统演变过程和机理.近些年来,非线性科学得到了飞速的发展,大量非线性分析方法被引入河川径流时间序列的分析当中.近似熵(Approximate Entropy,简记为 ApEn)是一种有效的非线性时间序列分析方法,由 Pincus 在进行时间序列复杂度的分析中首次提出6,因其在数据信号分析所

3、需数据量小、较强的抗噪声干扰能力、适用于随机信号和确定性信号等优点被广泛应用于医学、生物科学、机械设备故障诊断、气候突变712等方面.然而近似熵在径流时间序列的突变分析中的还没有涉及.本文将 ApEn和数据滑动技术相结合引入到径流突变的检测中,分析径流突变时空的不均一性和其内部演化的复杂性,揭示和辨别径流序列突变前后的动力学变化特征.1 方法ApEn 值是一个非负标量,用来表示时间序列的复杂性,熵值越大序列复杂度越高.对于任意一时间序列样本,给定维数 m 和容许偏差 r 值,公式推导如下6:1)构建一组维数为 m 的向量:2)定义向量 X(i)和 X(j)之间的欧式距离为:3)计算 i=1,

4、2,… ,n-m+1值下 dX(i), X(j)小于 r 的数目,记为C(r)4)对 C(r)mi取对数,对所有 i 取平均值,记为Φm:5)将维数 m 加 1,重复(1)(4)步骤,计算Φm+16)计算时间序列近似熵:由(6)式可知 ApEn 的值与 m,r 的取值有关,一般取 m=2 ,r =0.1- 0.2( σ为原始序列u(i)的标准差)由步骤(1)(6)可以看出,ApEn 反映的是序列重构后 m 维情况下两点的自相似程度大小,及序列在向量重构后,维数 m 增加时产生新模式的概率大小,熵值越大,产生新模式的几率就越大,因而序列就越复杂.不同时间

5、阶段下 ApEn 值反映序列各阶段保持自相似性状态的能力.水文系统可以看做为一个动力系统,某一径流要素的变化可以看做是系统运动的一个轨迹,因此可以采用动力学的方法对其进行研究,系统的动力学结构不发生变化时,对应的 ApEn 值相对平稳,在系统的动力学结构发生改变时,对应的 ApEn 值会发生明显的跃变,动力学结构越复杂,熵值就越大.数据滑动是一种有效的数据分析方法,可以通过数据滑动的方式实现子序列的提取,分析子序列与原序列、子序列与子序列之间的动力学结构异同,进而对序列进行突变分析,将近似熵与数据滑动相结合,称之为滑动近似熵(Moving Approximate Entropy,简称为 M-A

6、pEn),具体步骤如下: 初选模型参数.根据序列的长度初步选择滑动窗口长度 h、滑动步长 L、维数 m 和容许偏差 r; 计算子序列.由序列的第 i( i=1, 2,… ,N-h+1,N 为序列总个数)个数据开始选取长度为 h 的新 利用 ApEn 计算各新子序列的熵值; 保持 h 不变,以滑动步长 L 移动窗口,重复步骤,直至原序列结束; 通过步骤,可以得到一个长度为 N-h+1 的 ApEn 序列; 对得到的 ApEn 值进行均值分析并进行显着性检验,进行突变点或区间的最终确定.2 实例分析林家村站位于渭河流域关中段,受全球变暖、降水丰枯变化、蒸发的影响,流域降水发生了明显的

7、减少,同时由于近 50 年来人类活动(大、中及小型蓄水、引水工程)对流域下垫面的影响,径流发生了明显的改变,具体表现在年际丰枯变化剧烈(Cv 值由60 年代的 0.36 变为 90 年代的 0.54),径流总量呈现明显下降趋势值(如图 1).选取渭河流域林家村站1960.1.1-2000.12.31 共 41 年的日径流时间序列来进行分析计算,其中 2 月统一按 28 天处理.2.1 林家村日径流序列的 M-ApEn 突变性分析首先进行参数的选取,取 h=365d (其中 d 为 1 天),L=365d,m=2,r=0.15 ,通过(1)(6)的步骤计算,如图 2(a)所示 ApEn 值变化呈

8、现出四个变化阶段,发生了四次均值突变,分别为 1960-1972 年、1973-1982 年、1983-1993 年、1994-2000 年,对应突变点 1972 年、1982 年、1994 年,对突变点前后序列进行 t 检验(表 1),可以看到在置信度α=0.01的情况下,突变点前后的 ApEn 序列均存在显着地差异,即发生了突变,同时可以清楚看到各阶段的复杂性变化情况:1960-1972 年和 1983-1993 年两阶段相对于1973-1982 年和 1994-2000 年熵值大,复杂性强,可预测性差.图 3 中给出了径流序列滑动 t 检验的分析结果,给定显着性水平 &alp

9、ha;=0.01,可以看到,自 1960 年以来,t 统计量有一处(1994 年)超过 0.01 显着性水平,说明林家村站发生过一次明显的突变,对比 M-ApEn 的检测结果,滑动 t 检验易造成突变点的遗漏,且不能对突变点前后的动力学特征进进行有效的识别.2.2 滑动窗口和步长对检测结果的影响M-ApEn 检测结果与滑动窗口长度的选取有关,窗口的选择过大或过小将会对检测结果造成一定的影响9.在进行滑动窗口和滑动步长的影响分析时,首先假定一方不发生变化,如进行滑动窗口 h 对检测结果影响分析时,假定滑动步长 L 不变.取 h=730、1095d,结果如图 2(b)、(c),t 检验结果如表 1

10、,可以看出 ApEn 值变化规律基本和 h=365d 结果一致,均通过了α= 0.01的 t 检验,突变检验结果较为稳定,但如果滑动窗口长度过大,如图 2(d)h=3650d 时,各阶段序列熵值的阶段性变化不明显,使得突变点无法检测.故选择合适的滑动窗口长度对于检测结果的准确性具有重要的意义.初步分析中常采用较长的滑动窗口检测序列突变点或区间的大致范围,然后逐步缩小,直至序列突变点或区域稳定为止.对于日径流尺度的突变分析来说,建议滑动窗口 h 取 365d,730d,在此滑动窗口长度下,ApEn 能够很好地识别序列的突变点.如图 4 给出了咸阳、华县站在 h 取 365d,730d

11、 时径流序列突变分析结果,均取得了很好的突变检测效果.滑动步长 L 的选取不易过长或过短,应不小于滑动窗口 h 的长度,对于径流突变来说,一般分析以年为单位,建议取 L=365d,过长容易造成突变点的淹没,过短易造成计算量增加.3 突变发生的气候因素和人类活动影响M-ApEn 分析得出林家村站径流序列发生在 1972,1982,1994 年,这是气候因素和人类活动双重因素作用的结果.1)人类活动因素人类活动对于径流突变的影响主要表现在下垫面的影响.70 年代以前,流域内人口少,水库、灌溉等水利水保措施较少,工农业和生活用水较少,径流基本上处于天然状态.70 年代初,流域内陆续新建了大量的灌溉工

12、程和相应的水利措施,加之人口的不断增加,造成了径流量的减少,如 1970 年建成的东方红抽渭灌溉工程、1971 年建成的宝鸡峡塬上干渠于、灌溉工程等水利水土保持措施造成流域源区产流效率降低,径流减小,导致渭河的径流在 1971 开始发生了显着的变化13-14,在 1972 年发生了突变.80 年代以后,经济的发展使得水利工程(水库)的修建进入了一个高峰期,大量的径流被拦蓄,造成流域下垫面因素的改变,径流继续减少.如 1982 年冯家山水库的建成,最终导致 1982 年林家村的径流时间序列发生了突变.2)气候因素气候因素的影响主要和降水、温度、蒸发等有关.而这些因素与大气环流的变化和 ENSO

13、事件的发生有着直接的关系.渭河流域处于西太平洋副热带高压的影响之下,副热带的高压的变化直接影响到流域降水的多少,1972 年西太平洋副热高压极弱,西伸脊点偏东,加之处于厄尔尼诺年,在此环流背景下,流域降水急剧较少,形成了径流突变的气候原因.1982 年西太平洋副热带高压异常强大,西伸脊点偏西,加之 20 世纪最强的一次厄尔尼诺增暖事件,使得降水主要集中于江淮一带,而渭河流域的降水剧减,发生了重大转折15,形成了径流突变的气候条件.90 年代开始全球气候经历了剧烈的变化.1993 年开始流域的气温发生了大的转折16,而 1994 年西太平洋热带高压偏强、偏西和偏北,厄尔尼诺事件的发生使得流域降水

14、量大幅度下降,比 80 年代减少了14.86%,年平均径流量减少了 42.4%,与人类活动与共同促成了径流的突变17,18.4 结论 针对径流序列的非线性特征,通过近似熵结合滑动数据的方法进行分析,可得到更为合理的结果.实际过程应用中,通过对各子序列熵值的计算,画图分析其变化规律,并结合成因分析可以明确地诊断出序列的突变点,得到合理的结果. 滑动窗口以及步长的选择对于 M-ApEn 结果的准确性有着重要的影响,经验的作法是反复进行计算比对来确定最终的长度,需特别予以注意. M-ApEn 与其他方法相比不仅能够对序列的突变进行分析,还能诊断突变前后阶段的复杂性,有效地揭示径流序列在不同阶段的动力

15、学结构特征.对于变化环境下的高度复杂、非线性径流时间序列的突变分析有着很好的实用性. 由于气候和人类活动的影响径流变异表现出非线性的特点,滑动近似熵从非线性的角度出发进行突变的分析是行之有效的,可为流域水文变异诊断提供一条思路.参考文献:1 谢平, 陈广才, 李德, 等. 水文变异综合诊断方法及其应用研究 J.水电能源科学, 2005,23(2): 11-14.XIE Ping, CHEN Guangcai, LI De, et al. Comprehensive diagnosis method of hydrologic time series change-point analysis

16、J.Hydroelectric Energy, 2005, 23(2): 11-14. (in Chinese)2 ZHANG Q, XU C, CHEN Y D, et al. Abrupt behaviors of the stream flow of the Pearl river basin and implications forhydrological alterations across the Pearl river delta, China J. Journal of Hydrology, 2009, 377(3-4): 274-283.3 丁晶. 洪水时间序列干扰点的统计推

17、估 J.武汉水利水电学院学报, 1986, (5): 35-37.DING Jing. Statistical detection for transition point in flood time sequences J. Journal of Wuhan University of Hydraulic andElectric Engineering, 1986, (5): 35-37. (in Chinese)4 陈广才, 谢平. 水文变异的滑动 F 检验与检验方法J. 水文, 2006, 26(5): 57-60.CHEN Guangcai, XIE Ping .Slide F tes

18、t of change-point analysis J. Journal of China Hydrology, 2006, 26(5): 57-60. (inChinese)5 燕爱玲, 黄强, 王义民, 等.河川径流演变的非趋势波动分析 J. 水力发电学报, 2007, 26(3): 1-4.YAN Ailing, HUANG Qiang, WANG Yimin, et al. Detrended fluctuation analysis of river runoff evolvement J. Journal ofHydroelectric Engineering, 2007, 26

19、(3): 1-4. (in Chinese)6 PINCUS Steven M. Approximate entropy as a measure of system complexity J. Proc Natl Acad sic USA, 1991, 88(6):2297-2301.7 贺书云, 胡三觉, 王贤辉, 等.心脏对异丙肾上腺素反应性与心搏间期近似熵的关系 J. 生物物理学报, 2002, 18(2):231-235.HE Shuyun, HU Sanjue, WANG Xianhui et al. Relationship between the sensitivity to isoprenaline and the approximateentropy of r-r interval J. Acta Biophysica Sinica, 2002, 18(2): 231-235. (in Chinese)8 何文平, 吴琼, 张文, 等.滑动去趋势波动分析与近似熵在动力学结构突变检测中的性能比较 J.物理学报, 2009, 38(4):2862-2871.HE Wenping, WU Qiong, ZHANG Wen, et al. Comparison of c