孫東永 張洪波 王義民

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

2)(西安理工大學水利水電學院,西北旱區(qū)生態(tài)水利工程國家重點實驗室培育基地,西安 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分析中的極差,具有較強的穩(wěn)定性,相關實驗證明其對于標度指數在臨界值0.5附近的的估算比R/S更加有效[19],但其計算效率偏低;由于在各階趨勢成分處理上的優(yōu)勢,DFA方法非常適用于具有各種尺度噪聲及趨勢的非平穩(wěn)序列的標度計算[20,21];這些方法對于中小數據集序列的估計具有一定的有效性,但對于超大數據集序列的分析則需要進行復雜的計算和很高的內存要求,且結果有一定的偏差[8,22].WT法是在序列尺度和時間域上進行,其多尺度特性與自相似過程的尺度不變性有著自然的聯系[23],可以快速地對數據集序列進行不同尺度的分解,通過分析不同尺度下各小波變幅的標度關系來計算標度指數,具有計算速度快、收斂性好的特點[24],節(jié)省時間和內存;其次,通過改變小波基消失矩的數目,數據集序列的多項式趨勢能夠被嚴格地剔除,而相關數值試驗模擬證明小波分析還具有很強的抗噪能力[9,25],因而適用于超大數據集的非平穩(wěn)序列的分析.

由于氣候系統是一個非線性系統,它的觀測數據量大,且常常呈現出一些非線性現象,如復雜周期、趨勢、突變等,給序列標度指數的快速、準確計算帶來一定的困難.通過WT計算標度指數可為解決這一困難提供一條思路.本文參照文獻[3]將滑動移除技術與WT相融合,發(fā)展了一種新的動力學結構突變檢測方法——滑動移除小波分析法(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)個數據開始連續(xù)移除L個數據,形成int(N/L)(int表示取整)個長度為N?L的子序列;

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

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

與MC-R/S方法相同,對于無動力學結構突變且具有相關性的時間序列,任意移除該序列的數據,對其標度指數計算的影響幾乎可以忽略.因此,可以通過1)—4)檢測時間序列不同時間段內數據對于整個序列標度指數貢獻的大小來對系統的動力學結構突變進行檢測.

3 數值實驗

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

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

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

采用文獻[29]中構造的理想時間序列IS1(圖4),序列前1000個數據由Logistic映射產生,后1000個數據由滿足正態(tài)分布的隨機數組成,序列在t=1001處發(fā)生了突變,由一種非線性狀態(tài)轉變?yōu)榱硪环N隨機狀態(tài).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處發(fā)生了突變,突變前后呈現兩種狀態(tài),表現為由Logistic映射所產生的數據的標度指數序列變化幅度相對平穩(wěn),而由隨機數據生成的標度指數序列其變化幅度相對較大,表明數據的移除對于隨機序列的影響較大.同時隨著移除窗口L的增大,其序列動力學結構的突變更加明顯,這說明MC-WT對于非線性時間序列的動力學結構突變同樣有著良好的檢測能力,且對移除窗口L的長度依賴性較小.作為比較,圖6給出了IS1的滑動t檢驗(n1=10,n2=10,n1,n2分別為基準點前后子序列的長度)和Mann-Kendall的檢測結果,由圖6(a)可以看到,曲線呈現兩個明顯的階段,約在t=1000左右發(fā)生了動力學結構突變,但很難準確定位突變點的位置;圖6(b)中UF和UB線在置信區(qū)間(α=0.05)內t=1000左右均發(fā)生了改變,但UF和UB線并沒有相交,依據Mann-Kendall定義判斷此處并沒有發(fā)生突變,與實際情況不符.

圖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在非線性時間序列區(qū)間動力學結構突變中的檢測

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

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

說明數據的移除對于該區(qū)間標度指數估計影響較大,這與IS2構造的突變區(qū)間基本符合,表明MCWT對序列的區(qū)間突變有著良好的檢測能力,且對移除窗口長度依賴較小.同時進一步證明文獻[3]提到的具有相同動力學性質的數據對于序列標度指數計算的貢獻度大致相同,而具有不同動力學屬性的數據對計算整個序列標度指數的貢獻存在著顯著的差異.為了進一步驗證突變區(qū)間的準確性,采用文獻[3]提出的方差分析方法來定量區(qū)分不同動力學特性對于標度指數估算的貢獻,即定義方差閾值為三倍平均標準方差,超過該值即認為系統發(fā)生了突變.圖9分別為滑動移除窗口L=5,10,15,30時的方差貢獻圖,可以看到,除了在區(qū)間[301,331]內標度指數計算的方差貢獻超過了三倍方差閾值,在其他區(qū)域內方差貢獻基本接近于0值,可以判定序列在區(qū)間[301,331]發(fā)生了突變,與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在處理大數據時將有明顯的優(yōu)勢.最后,為了測試高斯白噪聲對MC-WT方法檢測結果的影響程度,分別對理想時間序列IS2依次添加SNR為15,20,25,30 dB的高斯白噪聲,圖10分別給出了在滑動步長L=10情況下,加噪后IS2序列的MC-WT動力學突變檢測方差貢獻圖,可以看到方差貢獻的突變區(qū)間基本與真實區(qū)間一致(除SNR=15 dB情況下個別點超出閾值),且隨著SNR的逐漸增大突變區(qū)間愈加清晰,沒有出現虛假的突變區(qū)間,說明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在理想時間序列動力學突變中的應用,而實測資料則呈現出更加復雜的非線性動態(tài)特性.鑒于此,本文擬以實測溫度資料測試MC-WT在突變檢測中的性能.實測溫度資料采用渭河流域佛坪站1960.1.1—2012.7.31(共19207個數據)逐日最高溫度數據,資料來源于中國氣象數據網(http://data.cma.cn/),質量得到了控制.圖11給出了佛坪站日最高溫度的MC-WT檢測結果,可以看到,與理想試驗結果類似,不論滑動移除窗口L=365 d(d=1日)或L=730 d,逐日最高溫度的標度指數γ序列出現了一個基本相同的突變區(qū)間:1972(1973)—1978年,在1978年以后系統的動力學結構發(fā)生了突變,由一種狀態(tài)進入到另一種狀態(tài),系統的標度指數降低,隨機性加大,這與20世紀70年代末期全球的氣候突變相一致[32?36].作為對比,圖12給出了在滑動移除窗口L=365 d的情況下MC-WT和MC-R/S方法突變檢測的方差貢獻(MC-V/S方差貢獻超出方差閾值,故剔除).可以看到兩種方法所得到的突變區(qū)間完全一致,而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的檢測結果對滑動移除窗口的長度依賴小,對噪聲具有一定的抗干擾能力,不僅能對線性序列的動力學結構突變實現準確檢測,且對非線性序列的動力學結構變點、突變區(qū)間同樣具有很好的檢測能力.實測資料的突變檢測結果進一步印證了以上結論,并證明其在更復雜的實測序列上仍能獲得較好的檢測效果.與MC-R/S,MC-V/S相關時間序列動力學結構突變分析方法相比,MC-WT檢測不僅具有相當的精確度,且檢測速度優(yōu)勢明顯,在大量數據分析中具有一定的優(yōu)勢,可為相關時間序列的動力學結構突變分析提供一條新的途徑.同時研究中也注意到,在某些情況下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]