孫東永 張洪波 王義民

1)(長安大學(xué)環(huán)境科學(xué)與工程學(xué)院,旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點實驗室,西安 710054)

2)(西安理工大學(xué)水利水電學(xué)院,西北旱區(qū)生態(tài)水利工程國家重點實驗室培育基地,西安 710048)

1 引 言

在地球物理學(xué)時間序列相關(guān)動力學(xué)進程的研究過程中,分形維數(shù)和標(biāo)度指數(shù)一直是目前廣被認(rèn)可的方法[1],尤其是標(biāo)度指數(shù)法.如在系統(tǒng)動力學(xué)結(jié)構(gòu)突變檢測過程中將滑動數(shù)據(jù)移除和標(biāo)度指數(shù)相結(jié)合提出的滑動移除去趨勢波動分析(moving cut detrended fluctuation analysis,MC-DFA)、滑動移除重標(biāo)極差分析(moving cut data-rescaled range analysis,MC-R/S)、滑動移除重標(biāo)方差分析(moving cut data-rescaled variance analysis,MCV/S)等[2?5],理想時間序列數(shù)值試驗和實測資料分析驗證結(jié)果都表明,相對于傳統(tǒng)統(tǒng)計方法Mann-Kendall、滑動t檢驗等,這些方法不僅能夠有效地檢測系統(tǒng)動力學(xué)結(jié)構(gòu)的突變點,而且能夠表征系統(tǒng)動力學(xué)結(jié)構(gòu)突變前后的性質(zhì)變化,極大地豐富了系統(tǒng)動力學(xué)結(jié)構(gòu)突變檢測理論與方法體系.其中,標(biāo)度指數(shù)的快速、準(zhǔn)確的計算是這些系統(tǒng)動力學(xué)結(jié)構(gòu)突變方法的核心.標(biāo)度指數(shù)的計算方法主要包括重標(biāo)極差分析法(rescaled range analysis,R/S)[6]、去趨勢波動分析(detrended fluctuation analysis,DFA)[7]、小波分析法(wavelet transformation,WT)[8?15]、重標(biāo)方差(rescaled variance analysis,V/S)法等[16].R/S分析法是最常用的非參數(shù)標(biāo)度指數(shù)計算方法,但在長序列分析中易受短期相關(guān)性和周期的影響,結(jié)果會出現(xiàn)一定的偏差[17,18];V/S分析法通過方差代替R/S分析中的極差,具有較強的穩(wěn)定性,相關(guān)實驗證明其對于標(biāo)度指數(shù)在臨界值0.5附近的的估算比R/S更加有效[19],但其計算效率偏低;由于在各階趨勢成分處理上的優(yōu)勢,DFA方法非常適用于具有各種尺度噪聲及趨勢的非平穩(wěn)序列的標(biāo)度計算[20,21];這些方法對于中小數(shù)據(jù)集序列的估計具有一定的有效性,但對于超大數(shù)據(jù)集序列的分析則需要進行復(fù)雜的計算和很高的內(nèi)存要求,且結(jié)果有一定的偏差[8,22].WT法是在序列尺度和時間域上進行,其多尺度特性與自相似過程的尺度不變性有著自然的聯(lián)系[23],可以快速地對數(shù)據(jù)集序列進行不同尺度的分解,通過分析不同尺度下各小波變幅的標(biāo)度關(guān)系來計算標(biāo)度指數(shù),具有計算速度快、收斂性好的特點[24],節(jié)省時間和內(nèi)存;其次,通過改變小波基消失矩的數(shù)目,數(shù)據(jù)集序列的多項式趨勢能夠被嚴(yán)格地剔除,而相關(guān)數(shù)值試驗?zāi)M證明小波分析還具有很強的抗噪能力[9,25],因而適用于超大數(shù)據(jù)集的非平穩(wěn)序列的分析.

由于氣候系統(tǒng)是一個非線性系統(tǒng),它的觀測數(shù)據(jù)量大,且常常呈現(xiàn)出一些非線性現(xiàn)象,如復(fù)雜周期、趨勢、突變等,給序列標(biāo)度指數(shù)的快速、準(zhǔn)確計算帶來一定的困難.通過WT計算標(biāo)度指數(shù)可為解決這一困難提供一條思路.本文參照文獻[3]將滑動移除技術(shù)與WT相融合,發(fā)展了一種新的動力學(xué)結(jié)構(gòu)突變檢測方法——滑動移除小波分析法(moving cut data wavelet transformation,MC-WT).該方法與MC-R/S類似,是基于數(shù)據(jù)的移除對于具有相同動力學(xué)屬性的相關(guān)序列標(biāo)度指數(shù)的估算幾乎沒有影響的這一特征而提出.為了全方面檢驗MC-WT方法在動力學(xué)結(jié)構(gòu)突變檢測中的性能,文中首先通過構(gòu)造線性和非線性兩種理想時間序列,分別檢測MC-WT方法的有效性,再以佛坪站日最高溫度實測資料對方法進行驗證.

2 MC-WT

2.1 標(biāo)度指數(shù)γ的小波估計原理

目前通過小波分析來估計標(biāo)度指數(shù)的方法主要有wavelet-based analysis,averaged wavelet coefficient,wavelet transform modulus maxima等,本文采用文獻[9]提出的wavelet-based analysis方法,該方法計算的標(biāo)度指數(shù)在高斯假設(shè)條件下是一個無偏估計量,且概念簡單、能夠?qū)Τ髷?shù)據(jù)集的進行直接有效的分析.原理如下:

對于任一時間序列x(t)(t=1,2,···,N,N為序列長度),其能量譜ΓX(ω)滿足

式中,ω為能量譜的頻率;cf=cπΛ(2γ?1)sin(π?πγ),c為一正常數(shù);Λ為Gamma函數(shù);γ為標(biāo)度指數(shù);

通過小波變換得到小波變換系數(shù){dj,k}(j=1,2,···,M;k=1,2,···,2?j/N),其中,j為尺度參數(shù),k為位置參數(shù),M為分解尺度;該系數(shù)可以度量x(t)在時間2jk,頻率為2?jω0處的能量,ω0為小波函數(shù)ψ(t)的參考頻率;對給定的尺度j,為能量譜估計量,即

式中,Nj=2?jN為尺度j小波變換系數(shù)個數(shù);又因

E[ΓX(2?jω0)]表示ΓX(2?jω0)的數(shù)學(xué)期望, 將(1)式代入(3)式,可得

式中,cg為與cf有關(guān)的一常數(shù),對(4)式兩邊取2為底的對數(shù),

通過該式log2與j之間的線性回歸可得標(biāo)度指數(shù)的估計量γ.

2.2 標(biāo)度指數(shù)γ的小波估計的計算步驟[26,27]

1)依據(jù)時間序列x(t)的長度N選取分解尺度M,對其進行Mallat一維小波分解,計算小波系數(shù)dj,k(j=1,2,···,M;k=1,2,···,2?j/N).

2)由小波系數(shù)計算中間參量ηj,sj(j=1,2,···,M)

3)計算標(biāo)度指數(shù)γ的小波估計值γ(j1,j2),其中1≤j1≤j2≤M,

2.3 MC-WT方法

[3]中MC-R/S方法,本文給出MCWT分析方法的具體步驟:

1)依據(jù)序列長度N選擇移除窗口長度L;

2)取滑動步長為L,從序列x(t)的第t(t=1,2,···,N?L+1)個數(shù)據(jù)開始連續(xù)移除L個數(shù)據(jù),形成int(N/L)(int表示取整)個長度為N?L的子序列;

3)通過小波估計各子序列的標(biāo)度指數(shù)γ,可以得到一個長度為int(N/L)的標(biāo)度指數(shù)序列;

4)對標(biāo)度指數(shù)序列進行方差分析,根據(jù)方差貢獻大小確定原序列的突變點或突變區(qū)間.

與MC-R/S方法相同,對于無動力學(xué)結(jié)構(gòu)突變且具有相關(guān)性的時間序列,任意移除該序列的數(shù)據(jù),對其標(biāo)度指數(shù)計算的影響幾乎可以忽略.因此,可以通過1)—4)檢測時間序列不同時間段內(nèi)數(shù)據(jù)對于整個序列標(biāo)度指數(shù)貢獻的大小來對系統(tǒng)的動力學(xué)結(jié)構(gòu)突變進行檢測.

3 數(shù)值實驗

3.1 MC-WT在線性時間序列突變中的檢測

為了對MC-WT方法的性能進行全面了解,首先進行線性序列的動力學(xué)突變檢測試驗.理想序列IS0采用如下方程構(gòu)建[28]:

由(8)式可知,序列y(t)在t=1000處發(fā)生了動力學(xué)結(jié)構(gòu)突變,序列由初始的正弦函數(shù)方程突變?yōu)檎嘞液瘮?shù)控制的方程,如圖1所示.選取滑動移除窗口L=2,采用Mallat離散小波變換算法計算各子序列標(biāo)度指數(shù)γ,其中濾波器組選用sym8,根據(jù)序列長度N=2000,取M=9,j1=1,j2=M.圖2(a)給出了在滑動移除窗口L=2情況下理想序列IS0的MC-WT檢測結(jié)果,容易看到,在t=1001處,標(biāo)度指數(shù)γ發(fā)生了一次顯著的均值突變,突變前后呈現(xiàn)明顯的兩種動力學(xué)結(jié)構(gòu)特征,準(zhǔn)確地刻畫了原系統(tǒng)的動力學(xué)結(jié)構(gòu)變化,因此可通過標(biāo)度指數(shù)的變化確定原序列的動力學(xué)突變.圖2(b)—(d)分別給出了序列在滑動移除窗口L=5,10,50情況下的檢測結(jié)果,可以看到,不論滑動移除窗口L如何變化,序列均在t=1001處發(fā)生了突變,說明MC-WT方法在對線性序列的突變檢測中受移除窗口長度的影響較小,能夠準(zhǔn)確地檢測系統(tǒng)的動力學(xué)突變.相關(guān)研究表明,在信號處理領(lǐng)域,由于電子設(shè)備或通信系統(tǒng)內(nèi)部缺陷(如電路電流突變、元件靜電感應(yīng)、磁感應(yīng)等)和外部電磁干擾(如太陽輻射電磁波、信號發(fā)射基站信號等),信號從輸入端開始不可避免地疊加了不同程度的噪聲,使得信道中的模擬信號受到干擾,輸出信號可能出現(xiàn)失真、誤碼等情況,因此在進行數(shù)據(jù)分析時必須考慮強噪聲對檢測結(jié)果造成的影響[29,30].為了測試噪聲對MC-WT方法檢測結(jié)果的影響程度,分別對理想時間序列IS0依次添加信噪比(signal-noise ratio,SNR)為20,25,30 dB的高斯白噪聲,圖3可以看到,在SNR=20,25,30 dB情況下加噪后的理想時間序列(滑動移除窗口L=5)標(biāo)度指數(shù)γ均在t=1001處發(fā)生了突變,說明MC-WT具有較強的抗噪能力,其他滑動步長結(jié)果類似.以上分析表明,對于線性序列的動力學(xué)結(jié)構(gòu)突變MC-WT方法有著很強的檢測能力,然而在自然界中,系統(tǒng)的演化呈現(xiàn)出復(fù)雜、動態(tài)的非線性特征,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在非線性時間序列中的動力學(xué)結(jié)構(gòu)突變檢測

采用文獻[29]中構(gòu)造的理想時間序列IS1(圖4),序列前1000個數(shù)據(jù)由Logistic映射產(chǎn)生,后1000個數(shù)據(jù)由滿足正態(tài)分布的隨機數(shù)組成,序列在t=1001處發(fā)生了突變,由一種非線性狀態(tài)轉(zhuǎn)變?yōu)榱硪环N隨機狀態(tài).Logistic映射方程如(9)式,其中初值x0=0.8,參數(shù)μ=3.8.

圖4 理想時間序列IS1Fig.4.The ideal time series IS1.

圖5為IS1在不同滑動移除窗口L下的MCWT檢測結(jié)果,濾波器組選用sym8,取分解尺度M=9(j1=1,j2=M).從圖5(a)—(d)可以看到,不論是滑動移除窗口L=10,還是L=20,25,50,其標(biāo)度指數(shù)γ的演變趨勢非常類似,均在t=1001處發(fā)生了突變,突變前后呈現(xiàn)兩種狀態(tài),表現(xiàn)為由Logistic映射所產(chǎn)生的數(shù)據(jù)的標(biāo)度指數(shù)序列變化幅度相對平穩(wěn),而由隨機數(shù)據(jù)生成的標(biāo)度指數(shù)序列其變化幅度相對較大,表明數(shù)據(jù)的移除對于隨機序列的影響較大.同時隨著移除窗口L的增大,其序列動力學(xué)結(jié)構(gòu)的突變更加明顯,這說明MC-WT對于非線性時間序列的動力學(xué)結(jié)構(gòu)突變同樣有著良好的檢測能力,且對移除窗口L的長度依賴性較小.作為比較,圖6給出了IS1的滑動t檢驗(n1=10,n2=10,n1,n2分別為基準(zhǔn)點前后子序列的長度)和Mann-Kendall的檢測結(jié)果,由圖6(a)可以看到,曲線呈現(xiàn)兩個明顯的階段,約在t=1000左右發(fā)生了動力學(xué)結(jié)構(gòu)突變,但很難準(zhǔn)確定位突變點的位置;圖6(b)中UF和UB線在置信區(qū)間(α=0.05)內(nèi)t=1000左右均發(fā)生了改變,但UF和UB線并沒有相交,依據(jù)Mann-Kendall定義判斷此處并沒有發(fā)生突變,與實際情況不符.

圖5 理想時間序列IS1的MC-WT檢測結(jié)果 (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檢測結(jié)果 (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在非線性時間序列區(qū)間動力學(xué)結(jié)構(gòu)突變中的檢測

以上所考慮的是單點突變的情況,即系統(tǒng)突然由一種狀態(tài)過渡到另一種狀態(tài),而實際情況中還可能發(fā)生區(qū)間突變的情況,即系統(tǒng)在演變過程中某一時間段發(fā)生了動力學(xué)結(jié)構(gòu)突變之后又恢復(fù)到原來的狀態(tài).依據(jù)文獻[3]構(gòu)造理想時間序列IS2,即在Logistic映射產(chǎn)生一條長度為1000的理想演化序列中,預(yù)想時間序列在區(qū)間[301,330]由確定性方程轉(zhuǎn)變?yōu)殡S機狀態(tài).故IS2在區(qū)間[301,330]發(fā)生了一次動力學(xué)結(jié)構(gòu)突變(圖7).Logistic映射方程見(9)式.

圖8給出了IS2序列在不同滑動移除窗口L下MC-WT方法的檢測結(jié)果,濾波器組選用sym8,取分解尺度M=9(j1=1,j2=M).可以看到,滑動移除窗口長度L=5,10,15,30,在區(qū)間[301,331]內(nèi),其標(biāo)度指數(shù)γ的變化明顯大于其他區(qū)域,

說明數(shù)據(jù)的移除對于該區(qū)間標(biāo)度指數(shù)估計影響較大,這與IS2構(gòu)造的突變區(qū)間基本符合,表明MCWT對序列的區(qū)間突變有著良好的檢測能力,且對移除窗口長度依賴較小.同時進一步證明文獻[3]提到的具有相同動力學(xué)性質(zhì)的數(shù)據(jù)對于序列標(biāo)度指數(shù)計算的貢獻度大致相同,而具有不同動力學(xué)屬性的數(shù)據(jù)對計算整個序列標(biāo)度指數(shù)的貢獻存在著顯著的差異.為了進一步驗證突變區(qū)間的準(zhǔn)確性,采用文獻[3]提出的方差分析方法來定量區(qū)分不同動力學(xué)特性對于標(biāo)度指數(shù)估算的貢獻,即定義方差閾值為三倍平均標(biāo)準(zhǔn)方差,超過該值即認(rèn)為系統(tǒng)發(fā)生了突變.圖9分別為滑動移除窗口L=5,10,15,30時的方差貢獻圖,可以看到,除了在區(qū)間[301,331]內(nèi)標(biāo)度指數(shù)計算的方差貢獻超過了三倍方差閾值,在其他區(qū)域內(nèi)方差貢獻基本接近于0值,可以判定序列在區(qū)間[301,331]發(fā)生了突變,與MC-R/S和MC-V/S分析結(jié)果一致[31],說明MCWT方法具有良好動力學(xué)突變檢測能力.同時也注意到,在滑動移除窗口L=5結(jié)尾附近和L=10開始端,有個別方差貢獻也超過了方差閾值,可能與算法本身有關(guān),實驗結(jié)果表明加大移除窗口的長度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在處理大數(shù)據(jù)時將有明顯的優(yōu)勢.最后,為了測試高斯白噪聲對MC-WT方法檢測結(jié)果的影響程度,分別對理想時間序列IS2依次添加SNR為15,20,25,30 dB的高斯白噪聲,圖10分別給出了在滑動步長L=10情況下,加噪后IS2序列的MC-WT動力學(xué)突變檢測方差貢獻圖,可以看到方差貢獻的突變區(qū)間基本與真實區(qū)間一致(除SNR=15 dB情況下個別點超出閾值),且隨著SNR的逐漸增大突變區(qū)間愈加清晰,沒有出現(xiàn)虛假的突變區(qū)間,說明MC-WT具有很好的抗噪能力.

圖7 理想時間序列IS2Fig.7.The ideal time series IS2.

圖8 理想時間序列IS2的MC-WT檢測結(jié)果 (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在實測資料中的應(yīng)用

前文分析了MC-WT在理想時間序列動力學(xué)突變中的應(yīng)用,而實測資料則呈現(xiàn)出更加復(fù)雜的非線性動態(tài)特性.鑒于此,本文擬以實測溫度資料測試MC-WT在突變檢測中的性能.實測溫度資料采用渭河流域佛坪站1960.1.1—2012.7.31(共19207個數(shù)據(jù))逐日最高溫度數(shù)據(jù),資料來源于中國氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn/),質(zhì)量得到了控制.圖11給出了佛坪站日最高溫度的MC-WT檢測結(jié)果,可以看到,與理想試驗結(jié)果類似,不論滑動移除窗口L=365 d(d=1日)或L=730 d,逐日最高溫度的標(biāo)度指數(shù)γ序列出現(xiàn)了一個基本相同的突變區(qū)間:1972(1973)—1978年,在1978年以后系統(tǒng)的動力學(xué)結(jié)構(gòu)發(fā)生了突變,由一種狀態(tài)進入到另一種狀態(tài),系統(tǒng)的標(biāo)度指數(shù)降低,隨機性加大,這與20世紀(jì)70年代末期全球的氣候突變相一致[32?36].作為對比,圖12給出了在滑動移除窗口L=365 d的情況下MC-WT和MC-R/S方法突變檢測的方差貢獻(MC-V/S方差貢獻超出方差閾值,故剔除).可以看到兩種方法所得到的突變區(qū)間完全一致,而MC-WT方法所花費的時間是MC-R/S方法的1/25左右,說明在進行大數(shù)據(jù)分析中,MC-WT具有更高的效率.

圖11 佛坪站逐日最高溫度序列MC-WT檢測結(jié)果(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 結(jié) 論

本文通過融合小波標(biāo)度指數(shù)與數(shù)據(jù)移除技術(shù),提出一種新的動力學(xué)結(jié)構(gòu)突變檢測方法—MCWT.理想時間序列的試驗結(jié)果表明,MC-WT的檢測結(jié)果對滑動移除窗口的長度依賴小,對噪聲具有一定的抗干擾能力,不僅能對線性序列的動力學(xué)結(jié)構(gòu)突變實現(xiàn)準(zhǔn)確檢測,且對非線性序列的動力學(xué)結(jié)構(gòu)變點、突變區(qū)間同樣具有很好的檢測能力.實測資料的突變檢測結(jié)果進一步印證了以上結(jié)論,并證明其在更復(fù)雜的實測序列上仍能獲得較好的檢測效果.與MC-R/S,MC-V/S相關(guān)時間序列動力學(xué)結(jié)構(gòu)突變分析方法相比,MC-WT檢測不僅具有相當(dāng)?shù)木_度,且檢測速度優(yōu)勢明顯,在大量數(shù)據(jù)分析中具有一定的優(yōu)勢,可為相關(guān)時間序列的動力學(xué)結(jié)構(gòu)突變分析提供一條新的途徑.同時研究中也注意到,在某些情況下MC-WT在檢測開始時會出現(xiàn)1—2個虛假的突變點,這可能與小波分解算法的選取有關(guān),可以通過對比不同滑動窗口下檢測結(jié)果予以剔除;其次,對于強噪聲對信號序列的影響,文中只考慮了高斯白噪聲的情況,實際情況中各種噪聲(如尖峰噪聲)對信號序列的影響不同[37,38],因此,下一步將展開相關(guān)研究.

參考文獻

[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物理學(xué)報582862]

[4]He W P,Deng B S,Wu Q,Zhang W,Cheng H Y 2010Acta Phys.Sin.598264(in Chinese)[何文平,鄧北勝,吳瓊,張文,成海英2010物理學(xué)報598264]

[5]Sun D Y,Zhang H B,Huang Q 2014Acta Phys.Sin.63209203(in Chinese)[孫東永,張洪波,黃強 2014物理學(xué)報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計算機工程與應(yīng)用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水科學(xué)進展 10144]

[27]Li X B,Ding J,Li H Q 1998J.Hydraul.Eng.821(in Chinese)[李賢彬,丁晶,李后強 1998水利學(xué)報 821]

[28]Wang Q G,Zhang Z P 2008Acta Phys.Sin.571976(in Chinese)[王啟光,張增平 2008物理學(xué)報 571976]

[29]He W P,He T,Cheng H Y,Zhang W,Wu Q 2011Acta Phys.Sin.60049202(in Chinese)[何文平,何濤,成海英,張文,吳瓊2011物理學(xué)報60049202]

[30]Jin H M,He W P,Zhang W,Feng A X,Hou W 2012Acta Phys.Sin.61129202(in Chinese)[金紅梅,何文平,張文,馮愛霞,侯威2012物理學(xué)報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氣象學(xué)報66892]

[34]Shi N,Chen J Q,Tu Q P 1995Acta Meteor.Sin.53431(in Chinese)[施能,陳家其,屠其璞1995氣象學(xué)報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地球物理學(xué)報592290]