孙东永 张洪波 王义民

1)(长安大学环境科学与工程学院,旱区地下水文与生态效应教育部重点实验室,西安 710054)

2)(西安理工大学水利水电学院,西北旱区生态水利工程国家重点实验室培育基地,西安 710048)

1 引 言

在地球物理学时间序列相关动力学进程的研究过程中,分形维数和标度指数一直是目前广被认可的方法[1],尤其是标度指数法.如在系统动力学结构突变检测过程中将滑动数据移除和标度指数相结合提出的滑动移除去趋势波动分析(moving cut detrended fluctuation analysis,MC-DFA)、滑动移除重标极差分析(moving cut data-rescaled range analysis,MC-R/S)、滑动移除重标方差分析(moving cut data-rescaled variance analysis,MCV/S)等[2−5],理想时间序列数值试验和实测资料分析验证结果都表明,相对于传统统计方法Mann-Kendall、滑动t检验等,这些方法不仅能够有效地检测系统动力学结构的突变点,而且能够表征系统动力学结构突变前后的性质变化,极大地丰富了系统动力学结构突变检测理论与方法体系.其中,标度指数的快速、准确的计算是这些系统动力学结构突变方法的核心.标度指数的计算方法主要包括重标极差分析法(rescaled range analysis,R/S)[6]、去趋势波动分析(detrended fluctuation analysis,DFA)[7]、小波分析法(wavelet transformation,WT)[8−15]、重标方差(rescaled variance analysis,V/S)法等[16].R/S分析法是最常用的非参数标度指数计算方法,但在长序列分析中易受短期相关性和周期的影响,结果会出现一定的偏差[17,18];V/S分析法通过方差代替R/S分析中的极差,具有较强的稳定性,相关实验证明其对于标度指数在临界值0.5附近的的估算比R/S更加有效[19],但其计算效率偏低;由于在各阶趋势成分处理上的优势,DFA方法非常适用于具有各种尺度噪声及趋势的非平稳序列的标度计算[20,21];这些方法对于中小数据集序列的估计具有一定的有效性,但对于超大数据集序列的分析则需要进行复杂的计算和很高的内存要求,且结果有一定的偏差[8,22].WT法是在序列尺度和时间域上进行,其多尺度特性与自相似过程的尺度不变性有着自然的联系[23],可以快速地对数据集序列进行不同尺度的分解,通过分析不同尺度下各小波变幅的标度关系来计算标度指数,具有计算速度快、收敛性好的特点[24],节省时间和内存;其次,通过改变小波基消失矩的数目,数据集序列的多项式趋势能够被严格地剔除,而相关数值试验模拟证明小波分析还具有很强的抗噪能力[9,25],因而适用于超大数据集的非平稳序列的分析.

由于气候系统是一个非线性系统,它的观测数据量大,且常常呈现出一些非线性现象,如复杂周期、趋势、突变等,给序列标度指数的快速、准确计算带来一定的困难.通过WT计算标度指数可为解决这一困难提供一条思路.本文参照文献[3]将滑动移除技术与WT相融合,发展了一种新的动力学结构突变检测方法——滑动移除小波分析法(moving cut data wavelet transformation,MC-WT).该方法与MC-R/S类似,是基于数据的移除对于具有相同动力学属性的相关序列标度指数的估算几乎没有影响的这一特征而提出.为了全方面检验MC-WT方法在动力学结构突变检测中的性能,文中首先通过构造线性和非线性两种理想时间序列,分别检测MC-WT方法的有效性,再以佛坪站日最高温度实测资料对方法进行验证.

2 MC-WT

2.1 标度指数γ的小波估计原理

目前通过小波分析来估计标度指数的方法主要有wavelet-based analysis,averaged wavelet coefficient,wavelet transform modulus maxima等,本文采用文献[9]提出的wavelet-based analysis方法,该方法计算的标度指数在高斯假设条件下是一个无偏估计量,且概念简单、能够对超大数据集的进行直接有效的分析.原理如下:

对于任一时间序列x(t)(t=1,2,···,N,N为序列长度),其能量谱ΓX(ω)满足

式中,ω为能量谱的频率;cf=cπΛ(2γ−1)sin(π−πγ),c为一正常数;Λ为Gamma函数;γ为标度指数;

通过小波变换得到小波变换系数{dj,k}(j=1,2,···,M;k=1,2,···,2−j/N),其中,j为尺度参数,k为位置参数,M为分解尺度;该系数可以度量x(t)在时间2jk,频率为2−jω0处的能量,ω0为小波函数ψ(t)的参考频率;对给定的尺度j,为能量谱估计量,即

式中,Nj=2−jN为尺度j小波变换系数个数;又因

E[ΓX(2−jω0)]表示ΓX(2−jω0)的数学期望, 将(1)式代入(3)式,可得

式中,cg为与cf有关的一常数,对(4)式两边取2为底的对数,

通过该式log2与j之间的线性回归可得标度指数的估计量γ.

2.2 标度指数γ的小波估计的计算步骤[26,27]

1)依据时间序列x(t)的长度N选取分解尺度M,对其进行Mallat一维小波分解,计算小波系数dj,k(j=1,2,···,M;k=1,2,···,2−j/N).

2)由小波系数计算中间参量ηj,sj(j=1,2,···,M)

3)计算标度指数γ的小波估计值γ(j1,j2),其中1≤j1≤j2≤M,

2.3 MC-WT方法

[3]中MC-R/S方法,本文给出MCWT分析方法的具体步骤:

1)依据序列长度N选择移除窗口长度L;

2)取滑动步长为L,从序列x(t)的第t(t=1,2,···,N−L+1)个数据开始连续移除L个数据,形成int(N/L)(int表示取整)个长度为N−L的子序列;

3)通过小波估计各子序列的标度指数γ,可以得到一个长度为int(N/L)的标度指数序列;

4)对标度指数序列进行方差分析,根据方差贡献大小确定原序列的突变点或突变区间.

与MC-R/S方法相同,对于无动力学结构突变且具有相关性的时间序列,任意移除该序列的数据,对其标度指数计算的影响几乎可以忽略.因此,可以通过1)—4)检测时间序列不同时间段内数据对于整个序列标度指数贡献的大小来对系统的动力学结构突变进行检测.

3 数值实验

3.1 MC-WT在线性时间序列突变中的检测

为了对MC-WT方法的性能进行全面了解,首先进行线性序列的动力学突变检测试验.理想序列IS0采用如下方程构建[28]:

由(8)式可知,序列y(t)在t=1000处发生了动力学结构突变,序列由初始的正弦函数方程突变为正余弦函数控制的方程,如图1所示.选取滑动移除窗口L=2,采用Mallat离散小波变换算法计算各子序列标度指数γ,其中滤波器组选用sym8,根据序列长度N=2000,取M=9,j1=1,j2=M.图2(a)给出了在滑动移除窗口L=2情况下理想序列IS0的MC-WT检测结果,容易看到,在t=1001处,标度指数γ发生了一次显著的均值突变,突变前后呈现明显的两种动力学结构特征,准确地刻画了原系统的动力学结构变化,因此可通过标度指数的变化确定原序列的动力学突变.图2(b)—(d)分别给出了序列在滑动移除窗口L=5,10,50情况下的检测结果,可以看到,不论滑动移除窗口L如何变化,序列均在t=1001处发生了突变,说明MC-WT方法在对线性序列的突变检测中受移除窗口长度的影响较小,能够准确地检测系统的动力学突变.相关研究表明,在信号处理领域,由于电子设备或通信系统内部缺陷(如电路电流突变、元件静电感应、磁感应等)和外部电磁干扰(如太阳辐射电磁波、信号发射基站信号等),信号从输入端开始不可避免地叠加了不同程度的噪声,使得信道中的模拟信号受到干扰,输出信号可能出现失真、误码等情况,因此在进行数据分析时必须考虑强噪声对检测结果造成的影响[29,30].为了测试噪声对MC-WT方法检测结果的影响程度,分别对理想时间序列IS0依次添加信噪比(signal-noise ratio,SNR)为20,25,30 dB的高斯白噪声,图3可以看到,在SNR=20,25,30 dB情况下加噪后的理想时间序列(滑动移除窗口L=5)标度指数γ均在t=1001处发生了突变,说明MC-WT具有较强的抗噪能力,其他滑动步长结果类似.以上分析表明,对于线性序列的动力学结构突变MC-WT方法有着很强的检测能力,然而在自然界中,系统的演化呈现出复杂、动态的非线性特征,MC-WT方法的适用性如何,需要进一步检测.

图1 理想时间序列IS0Fig.1.The ideal time series IS0.

图3 加噪后理想序列IS0的MC-WT突变检测(L=5)Fig.3.The MC-WT mutation detection of ideal time series IS0 after adding noise(L=5).

3.2 MC-WT在非线性时间序列中的动力学结构突变检测

采用文献[29]中构造的理想时间序列IS1(图4),序列前1000个数据由Logistic映射产生,后1000个数据由满足正态分布的随机数组成,序列在t=1001处发生了突变,由一种非线性状态转变为另一种随机状态.Logistic映射方程如(9)式,其中初值x0=0.8,参数µ=3.8.

图4 理想时间序列IS1Fig.4.The ideal time series IS1.

图5为IS1在不同滑动移除窗口L下的MCWT检测结果,滤波器组选用sym8,取分解尺度M=9(j1=1,j2=M).从图5(a)—(d)可以看到,不论是滑动移除窗口L=10,还是L=20,25,50,其标度指数γ的演变趋势非常类似,均在t=1001处发生了突变,突变前后呈现两种状态,表现为由Logistic映射所产生的数据的标度指数序列变化幅度相对平稳,而由随机数据生成的标度指数序列其变化幅度相对较大,表明数据的移除对于随机序列的影响较大.同时随着移除窗口L的增大,其序列动力学结构的突变更加明显,这说明MC-WT对于非线性时间序列的动力学结构突变同样有着良好的检测能力,且对移除窗口L的长度依赖性较小.作为比较,图6给出了IS1的滑动t检验(n1=10,n2=10,n1,n2分别为基准点前后子序列的长度)和Mann-Kendall的检测结果,由图6(a)可以看到,曲线呈现两个明显的阶段,约在t=1000左右发生了动力学结构突变,但很难准确定位突变点的位置;图6(b)中UF和UB线在置信区间(α=0.05)内t=1000左右均发生了改变,但UF和UB线并没有相交,依据Mann-Kendall定义判断此处并没有发生突变,与实际情况不符.

图5 理想时间序列IS1的MC-WT检测结果 (a)L=10;(b)L=20;(c)L=25;(d)L=50Fig.5.The MC-WT detection result of ideal time series IS1:(a)L=10;(b)L=20;(c)L=25;(d)L=50.

图6 理想时间序列IS1检测结果 (a)滑动t检验;(b)Mann-KendallFig.6.The detection result of ideal time series IS1:(a)Moving t-test;(b)Mann-Kendall.

3.3 MC-WT在非线性时间序列区间动力学结构突变中的检测

以上所考虑的是单点突变的情况,即系统突然由一种状态过渡到另一种状态,而实际情况中还可能发生区间突变的情况,即系统在演变过程中某一时间段发生了动力学结构突变之后又恢复到原来的状态.依据文献[3]构造理想时间序列IS2,即在Logistic映射产生一条长度为1000的理想演化序列中,预想时间序列在区间[301,330]由确定性方程转变为随机状态.故IS2在区间[301,330]发生了一次动力学结构突变(图7).Logistic映射方程见(9)式.

图8给出了IS2序列在不同滑动移除窗口L下MC-WT方法的检测结果,滤波器组选用sym8,取分解尺度M=9(j1=1,j2=M).可以看到,滑动移除窗口长度L=5,10,15,30,在区间[301,331]内,其标度指数γ的变化明显大于其他区域,

说明数据的移除对于该区间标度指数估计影响较大,这与IS2构造的突变区间基本符合,表明MCWT对序列的区间突变有着良好的检测能力,且对移除窗口长度依赖较小.同时进一步证明文献[3]提到的具有相同动力学性质的数据对于序列标度指数计算的贡献度大致相同,而具有不同动力学属性的数据对计算整个序列标度指数的贡献存在着显著的差异.为了进一步验证突变区间的准确性,采用文献[3]提出的方差分析方法来定量区分不同动力学特性对于标度指数估算的贡献,即定义方差阈值为三倍平均标准方差,超过该值即认为系统发生了突变.图9分别为滑动移除窗口L=5,10,15,30时的方差贡献图,可以看到,除了在区间[301,331]内标度指数计算的方差贡献超过了三倍方差阈值,在其他区域内方差贡献基本接近于0值,可以判定序列在区间[301,331]发生了突变,与MC-R/S和MC-V/S分析结果一致[31],说明MCWT方法具有良好动力学突变检测能力.同时也注意到,在滑动移除窗口L=5结尾附近和L=10开始端,有个别方差贡献也超过了方差阈值,可能与算法本身有关,实验结果表明加大移除窗口的长度L可以消除该影响,如图9(c)和图9(d).其次,为了分析MC-WT的运行效率,表1给出了不同滑动移除窗口下MC-WT,MC-R/S和MC-V/S在同一电脑下(Inter Core(TM)i7-4510,2.4 GHz,4 GB,Win7)Matlab 2014 b平台的运行时间,MCWT花费时间大约为MC-R/S的1/6和MC-V/S的1/22,因而MC-WT在处理大数据时将有明显的优势.最后,为了测试高斯白噪声对MC-WT方法检测结果的影响程度,分别对理想时间序列IS2依次添加SNR为15,20,25,30 dB的高斯白噪声,图10分别给出了在滑动步长L=10情况下,加噪后IS2序列的MC-WT动力学突变检测方差贡献图,可以看到方差贡献的突变区间基本与真实区间一致(除SNR=15 dB情况下个别点超出阈值),且随着SNR的逐渐增大突变区间愈加清晰,没有出现虚假的突变区间,说明MC-WT具有很好的抗噪能力.

图7 理想时间序列IS2Fig.7.The ideal time series IS2.

图8 理想时间序列IS2的MC-WT检测结果 (a)L=5;(b)L=10;(c)L=15;(d)L=30Fig.8.The MC-WT detection result of ideal time series IS2:(a)L=5;(b)L=10;(c)L=15;(d)L=30.

图9 理想时间序列IS2的MC-WT方差贡献 (a)L=5;(b)L=10;(c)L=15;(d)L=30Fig.9.The variance contribution of MC-WT detection result for IS2:(a)L=5;(b)L=10;(c)L=15;(d)L=30.

图10 加噪后IS2序列MC-WT方差贡献 (a)SNR=15 dB;(b)SNR=20 dB;(c)SNR=25 dB;(d)SNR=30 dBFig.10.The variance contribution of MC-WT detection result for the IS2 after adding noise:(a)SNR=15 dB;(b)SNR=20 dB;(c)SNR=25 dB;(d)SNR=30 dB.

表1 不同移除窗口下MC-WT,MC-R/S和MC-V/S运行时间(单位:s)Table 1.The run time of MC-WT,MC-R/S,and MCV/S under different remove windows(unit:s).

4 MC-WT在实测资料中的应用

前文分析了MC-WT在理想时间序列动力学突变中的应用,而实测资料则呈现出更加复杂的非线性动态特性.鉴于此,本文拟以实测温度资料测试MC-WT在突变检测中的性能.实测温度资料采用渭河流域佛坪站1960.1.1—2012.7.31(共19207个数据)逐日最高温度数据,资料来源于中国气象数据网(http://data.cma.cn/),质量得到了控制.图11给出了佛坪站日最高温度的MC-WT检测结果,可以看到,与理想试验结果类似,不论滑动移除窗口L=365 d(d=1日)或L=730 d,逐日最高温度的标度指数γ序列出现了一个基本相同的突变区间:1972(1973)—1978年,在1978年以后系统的动力学结构发生了突变,由一种状态进入到另一种状态,系统的标度指数降低,随机性加大,这与20世纪70年代末期全球的气候突变相一致[32−36].作为对比,图12给出了在滑动移除窗口L=365 d的情况下MC-WT和MC-R/S方法突变检测的方差贡献(MC-V/S方差贡献超出方差阈值,故剔除).可以看到两种方法所得到的突变区间完全一致,而MC-WT方法所花费的时间是MC-R/S方法的1/25左右,说明在进行大数据分析中,MC-WT具有更高的效率.

图11 佛坪站逐日最高温度序列MC-WT检测结果(a)L=365 d;(b)L=730 dFig.11.The MC-WT detection results of daily maximum temperature sequence in Foping station:(a)L=365 d;(b)L=730 d.

图12 MC-WT和MC-R/S突变检测方差贡献图(L=365 d)Fig.12.The variance contribution of MC-WT and MC-R/S detection result(L=365 d).

表2 佛坪站极端温度不同移除窗口下MC-WT,MCR/S运行时间(单位:s)Table 2.The run time of MC-WT,MC-R/S under different remove windows in Foping station(unit:s).

5 结 论

本文通过融合小波标度指数与数据移除技术,提出一种新的动力学结构突变检测方法—MCWT.理想时间序列的试验结果表明,MC-WT的检测结果对滑动移除窗口的长度依赖小,对噪声具有一定的抗干扰能力,不仅能对线性序列的动力学结构突变实现准确检测,且对非线性序列的动力学结构变点、突变区间同样具有很好的检测能力.实测资料的突变检测结果进一步印证了以上结论,并证明其在更复杂的实测序列上仍能获得较好的检测效果.与MC-R/S,MC-V/S相关时间序列动力学结构突变分析方法相比,MC-WT检测不仅具有相当的精确度,且检测速度优势明显,在大量数据分析中具有一定的优势,可为相关时间序列的动力学结构突变分析提供一条新的途径.同时研究中也注意到,在某些情况下MC-WT在检测开始时会出现1—2个虚假的突变点,这可能与小波分解算法的选取有关,可以通过对比不同滑动窗口下检测结果予以剔除;其次,对于强噪声对信号序列的影响,文中只考虑了高斯白噪声的情况,实际情况中各种噪声(如尖峰噪声)对信号序列的影响不同[37,38],因此,下一步将展开相关研究.

参考文献

[1]Rehman S,Siddiqi A H 2009Chaos Soliton.Fract.401081

[2]He W P,Feng G L,Wu Q,He T,Wan S Q,Chou J F 2012Int.J.Climatol.321604

[3]He W P,Wu Q,Zhang W,Wang Q G,Zhang Y 2009Acta Phys.Sin.582862(in Chinese)[何文平,吴琼,张文,王启光,张勇2009物理学报582862]

[4]He W P,Deng B S,Wu Q,Zhang W,Cheng H Y 2010Acta Phys.Sin.598264(in Chinese)[何文平,邓北胜,吴琼,张文,成海英2010物理学报598264]

[5]Sun D Y,Zhang H B,Huang Q 2014Acta Phys.Sin.63209203(in Chinese)[孙东永,张洪波,黄强 2014物理学报63209203]

[6]Hurst H E 1951Trans.Am.Soc.Civ.Eng.116770

[7]Peng C K,Buldyrev S V,Havlin S,Simons M,Stanley H E,Goldberger A L 1994Phys.Rev.E491685

[8]Simonsen I,Hansen A,Nes O M 1998Phys.Rev.E582779

[9]Veitch D,Abry P C 1999IEEE Trans.Inf.Theory45878

[10]Jones C L,Lonergan G T,Mainwaring D E 1996J.Phys.A:Math.Gen.292509

[11]Hu K,Ivanov P C,Chen Z,Carpena P,Stanley H E 2001Phys.Rev.E64011114

[12]Gloter A,Hoff mann M 2007Ann.Stat.351947

[13]Manimaran P,Panigrahi P K,Parikh J C 2005Phys.Rev.E72046120

[14]Ciftlikli C,Gezer A 2010Turk.J.Elec.Eng.&Comp.Sci.18117

[15]Wu L,Ding Y M 2015Int.J.Wavelets Multiresolut Inf.Process.131550044

[16]Giraitis L,Kokoszka P,Leipus R,Teyssiere G 2003J.Econom.112265

[17]Clausel M,Roue ffF,Taqqu M S,Tudor C 2014Esaim Probab.Stat.1842

[18]Taqqu M S,Teverovsky V 1997Comm.Stat.Stoch.Model13723

[19]Cajueiro D O,Tabak B M 2005Math.Comp.Sim.70172

[20]Kantelhardt J W,Koscielny-Bunde E,Rego H H A,Havlin S,Bunde A 2001Physica A295441

[21]Matos J A O,Gama S M A,Ruskin H J,Sharkasi A A,Crane M 2008Physica A3873910

[22]Mielniczuk J,Wojdyllo P 2007Comput.Stat.Data Anal.514510

[23]Zhao Y Z,Wu L W 2014Comput.Eng.Appl.50154(in Chinese)[赵彦仲,吴立文2014计算机工程与应用50154]

[24]Dang T D,Molnar S 1999Period.Polytech,Electr.Eng.43227

[25]Giordano S,Miduri S,Pagano M,Russo F,Tartarelli S 1997Proceedings of 13th International Conference on Digital Signal ProcessingSantorini,Greece,July 2–4,1997 p479

[26]Li X B,Ding J,Li H Q 1999Adv.Water Sci.10144(in Chinese)[李贤彬,丁晶,李后强 1999水科学进展 10144]

[27]Li X B,Ding J,Li H Q 1998J.Hydraul.Eng.821(in Chinese)[李贤彬,丁晶,李后强 1998水利学报 821]

[28]Wang Q G,Zhang Z P 2008Acta Phys.Sin.571976(in Chinese)[王启光,张增平 2008物理学报 571976]

[29]He W P,He T,Cheng H Y,Zhang W,Wu Q 2011Acta Phys.Sin.60049202(in Chinese)[何文平,何涛,成海英,张文,吴琼2011物理学报60049202]

[30]Jin H M,He W P,Zhang W,Feng A X,Hou W 2012Acta Phys.Sin.61129202(in Chinese)[金红梅,何文平,张文,冯爱霞,侯威2012物理学报61129202]

[31]He W P,Liu Q Q,Jiang Y D,Lu Y 2015Chin.Phys.B24049205

[32]Powell A M,Xu J J 2011Theor.Appl.Climatol.104443

[33]Feng G L,Gong Z Q,Zhi R 2008Acta Meteor.Sin.66892(in Chinese)[封国林,龚志强,支蓉2008气象学报66892]

[34]Shi N,Chen J Q,Tu Q P 1995Acta Meteor.Sin.53431(in Chinese)[施能,陈家其,屠其璞1995气象学报53431]

[35]Tong J L,Wu H,Hou W,He W P,Zhou J 2014Chin.Phys.B23049201

[36]Wu H,Hou W,Yan P C,Zhang Z S,Wang K 2015Chin.Phys.B24089201

[37]Zhang M L,Qu H,Xie X R,Kurths J 2017Neurocomputing219333

[38]Wan L,Zhang Y,Lin J,Jiang C D,Lin T T 2016Chin.J.Geophys.592290(in Chinese)[万玲,张扬,林君,蒋川东,林婷婷2016地球物理学报592290]