滑動(dòng)移除小波分析法在動(dòng)力學(xué)結(jié)構(gòu)突變檢驗(yàn)中的應(yīng)用?

孫東永 張洪波 王義民

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

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

1 引 言

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

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

2 MC-WT

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 MC-WT方法

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

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

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

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

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

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

3 數(shù)值實(shí)驗(yàn)

3.1 MC-WT在線性時(shí)間序列突變中的檢測(cè)

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

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

圖1 理想時(shí)間序列IS0Fig.1.The ideal time series IS0.

圖3 加噪后理想序列IS0的MC-WT突變檢測(cè)(L=5)Fig.3.The MC-WT mutation detection of ideal time series IS0 after adding noise(L=5).

3.2 MC-WT在非線性時(shí)間序列中的動(dòng)力學(xué)結(jié)構(gòu)突變檢測(cè)

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

圖4 理想時(shí)間序列IS1Fig.4.The ideal time series IS1.

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

圖5 理想時(shí)間序列IS1的MC-WT檢測(cè)結(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 理想時(shí)間序列IS1檢測(cè)結(jié)果 (a)滑動(dòng)t檢驗(yàn);(b)Mann-KendallFig.6.The detection result of ideal time series IS1:(a)Moving t-test;(b)Mann-Kendall.

3.3 MC-WT在非線性時(shí)間序列區(qū)間動(dòng)力學(xué)結(jié)構(gòu)突變中的檢測(cè)

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

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

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

圖7 理想時(shí)間序列IS2Fig.7.The ideal time series IS2.

圖8 理想時(shí)間序列IS2的MC-WT檢測(cè)結(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 理想時(shí)間序列IS2的MC-WT方差貢獻(xiàn) (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方差貢獻(xiàn) (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運(yùn)行時(shí)間(單位:s)Table 1.The run time of MC-WT,MC-R/S,and MCV/S under different remove windows(unit:s).

4 MC-WT在實(shí)測(cè)資料中的應(yīng)用

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

圖11 佛坪站逐日最高溫度序列MC-WT檢測(cè)結(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突變檢測(cè)方差貢獻(xiàn)圖(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運(yùn)行時(shí)間(單位:s)Table 2.The run time of MC-WT,MC-R/S under different remove windows in Foping station(unit:s).

5 結(jié) 論

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

參考文獻(xià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é)報(bào)582862]

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

[5]Sun D Y,Zhang H B,Huang Q 2014Acta Phys.Sin.63209203(in Chinese)[孫東永,張洪波,黃強(qiáng) 2014物理學(xué)報(bào)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計(jì)算機(jī)工程與應(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)[李賢彬,丁晶,李后強(qiáng) 1999水科學(xué)進(jìn)展 10144]

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

[28]Wang Q G,Zhang Z P 2008Acta Phys.Sin.571976(in Chinese)[王啟光,張?jiān)銎?2008物理學(xué)報(bào) 571976]

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

[30]Jin H M,He W P,Zhang W,Feng A X,Hou W 2012Acta Phys.Sin.61129202(in Chinese)[金紅梅,何文平,張文,馮愛(ài)霞,侯威2012物理學(xué)報(bào)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)[封國(guó)林,龔志強(qiáng),支蓉2008氣象學(xué)報(bào)66892]

[34]Shi N,Chen J Q,Tu Q P 1995Acta Meteor.Sin.53431(in Chinese)[施能,陳家其,屠其璞1995氣象學(xué)報(bào)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)[萬(wàn)玲,張揚(yáng),林君,蔣川東,林婷婷2016地球物理學(xué)報(bào)592290]