基于近似熵的突變檢測新方法*

何文平 何 濤 成海英 張 文 吳 瓊

1)(國家氣候中心，北京 100081)2)(濟(jì)南市環(huán)境保護(hù)監(jiān)測站，濟(jì)南 250014)3)(鹽城工學(xué)院基礎(chǔ)教學(xué)部，鹽城 224002)4)(南京郵電大學(xué)通達(dá)學(xué)院，南京 210003)5)(國家衛(wèi)星氣象中心，北京 100081)(2010年6月6日收到;2010年7月12日收到修改稿)

基于近似熵的突變檢測新方法*

何文平1)何 濤2)成海英3)張 文4)吳 瓊5)

1)(國家氣候中心，北京 100081)2)(濟(jì)南市環(huán)境保護(hù)監(jiān)測站，濟(jì)南 250014)3)(鹽城工學(xué)院基礎(chǔ)教學(xué)部，鹽城 224002)4)(南京郵電大學(xué)通達(dá)學(xué)院，南京 210003)5)(國家衛(wèi)星氣象中心，北京 100081)(2010年6月6日收到;2010年7月12日收到修改稿)

近似熵是一個(gè)有效的非線性動力學(xué)指數(shù)，能夠用于表征時(shí)間序列的復(fù)雜性，通過滑動窗口技術(shù)，近似熵對于一維時(shí)間序列的動力學(xué)結(jié)構(gòu)突變具有一定的識別能力，但其突變檢測結(jié)果依賴于子序列長度的選擇，且不能準(zhǔn)確定位突變點(diǎn).鑒于此，本文提出了一種新的突變檢測方法——滑動移除近似熵.測試結(jié)果表明，滑動移除近似熵具有檢測結(jié)果穩(wěn)定性好、準(zhǔn)確性高等特點(diǎn)，明顯優(yōu)于滑動近似熵和Mann-Kendall方法，其在實(shí)際觀測資料中的應(yīng)用進(jìn)一步證實(shí)了新方法的可靠性.

近似熵，滑動移除近似熵，突變檢測

PACS:92.60.Wc

1.引 言

近年來，隨著非線性科學(xué)的不斷發(fā)展，非線性動力學(xué)領(lǐng)域中的諸多指數(shù)被用于描述系統(tǒng)的動力學(xué)特征，如關(guān)聯(lián)維數(shù)、李雅普諾夫指數(shù)、柯爾莫哥羅夫熵等［1—8］，其中關(guān)聯(lián)維數(shù)描述了一個(gè)復(fù)雜系統(tǒng)所需要的最少的有效狀態(tài)變量;李雅普諾夫指數(shù)則被用于度量非線性耗散系統(tǒng)的兩個(gè)運(yùn)動特性:收縮與發(fā)散，表征了相空間中在不同方向上軌道的收縮或發(fā)散的速率，是判斷系統(tǒng)動力學(xué)行為是否為混沌的一個(gè)很好的物理指標(biāo);由于每一個(gè)正的李雅普諾夫指數(shù)表示該方向上長度是伸長的，因此該方向上初始信息不斷損失.所有正的李雅普諾夫指數(shù)之和代表了總的信息損失率，此即為柯爾莫哥羅夫熵的定義，它是表征系統(tǒng)產(chǎn)生信息量的多少、快慢的物理量.

然而，要直接根據(jù)這些指數(shù)的定義來計(jì)算分維和李雅普諾夫指數(shù)并不容易，實(shí)際中往往會面臨單變量的一維時(shí)間序列，很顯然，系統(tǒng)的非線性特征信息蘊(yùn)藏在這些序列中，如何從中提取這類信息成為一個(gè)非常困難的問題［9］.鑒于此，Packara提出了通過時(shí)間序列相空間重構(gòu)來計(jì)算李雅普諾夫指數(shù)和關(guān)聯(lián)維數(shù)，但是為了得到一個(gè)好的最大李雅普諾夫指數(shù)估計(jì)值，必須要求:1)所給定的時(shí)間序列有足夠的長度，Wolf認(rèn)為，這個(gè)長度的估計(jì)范圍為10D—30D(D為吸引子的分維);2)延遲時(shí)間的選取;3)臨近相空間點(diǎn)的選取.但目前，相空間重構(gòu)中關(guān)于延遲時(shí)間的選取是一個(gè)難題，而對初始相點(diǎn)的最鄰近點(diǎn)的選擇更為困難，尤其是需要花費(fèi)大量的機(jī)時(shí)不斷增加嵌入空間的維數(shù)以便獲得穩(wěn)健的李雅普諾夫指數(shù).因此，這些非線性動力學(xué)指數(shù)在實(shí)際應(yīng)用中顯得并不方便，需要尋找更為有效、穩(wěn)健且對數(shù)據(jù)長度要求較短的新方法.

近似熵［10—12］(approximate entropy，ApEn)方法無疑是一個(gè)理想的非線性動力學(xué)指數(shù)，它具有所需數(shù)據(jù)點(diǎn)較少、有較好的抗噪和抗干擾能力且適用于隨機(jī)信號、確定性信號或兩種混合信號等特點(diǎn).因此，ApEn方法在醫(yī)學(xué)、生物學(xué)、機(jī)械設(shè)備故障診斷、語音信號端點(diǎn)檢測等方面得到了廣泛應(yīng)用［13—17］.

氣候系統(tǒng)是一個(gè)復(fù)雜的非線性動力學(xué)系統(tǒng)，如何從過去的氣候記錄中檢測氣候突變，對于預(yù)測未來的氣候演變趨勢有著重要的現(xiàn)實(shí)意義［18—26］.已有學(xué)者對此做了大量卓有成效的工作［27—35］，例如基于標(biāo)度理論先后發(fā)展了滑動去趨勢波動分析方法、滑動移除去趨勢波動分析方法、滑動移除重標(biāo)極差方法［25—28］，這些方法的出現(xiàn)在很大程度上豐富了突變檢測的手段，并能夠很好的對大量實(shí)際觀測資料進(jìn)行動力學(xué)結(jié)構(gòu)突變檢測，但是必須指出的是，這些方法也存在一個(gè)共同的缺陷，即應(yīng)用這三種方法的前提是待分析時(shí)間序列具有分形特征，但并不是所有的時(shí)間序列都能夠滿足這一條件.ApEn對于各種信號的廣泛適用性，使得其成為一種檢測突變的另一途徑.類似于滑動-t方法，王啟光等［36］將 ApEn方法引入到氣候系統(tǒng)的動力學(xué)結(jié)構(gòu)突變檢測中，利用滑動技術(shù)通過理想試驗(yàn)證實(shí)了ApEn具有區(qū)分不同動力學(xué)結(jié)構(gòu)的能力，表明其能夠用于動力學(xué)結(jié)構(gòu)突變檢測中，并將其應(yīng)用于實(shí)際氣候資料的突變檢測，得到了較好的結(jié)果.但是這種滑動技術(shù)本身存在著眾所周知的先天不足:其一，檢測結(jié)果依賴于子序列長度;其二，不能夠準(zhǔn)確的定位突變點(diǎn)的位置，只能夠給出一個(gè)大致的突變區(qū)域，這顯然不能夠滿足實(shí)際應(yīng)用的需求.

借鑒于滑動移除去趨勢波動分析方法、滑動移除重標(biāo)極差方法中所采用滑動移除數(shù)據(jù)的方式，本文對ApEn在突變檢測中的應(yīng)用提出了一種改進(jìn)的新方法——滑動移除近似熵(moving cut dataapproximate entropy，MC-ApEn).試驗(yàn)結(jié)果驗(yàn)證了新方法的可靠性和準(zhǔn)確性，其檢測結(jié)果明顯優(yōu)于簡單的滑動 ApEn方法和 Mann-Kendall方法，其在逐日降水資料中的應(yīng)用進(jìn)一步證實(shí)了新方法的可靠性，展示了其潛在的應(yīng)用前景.

2.方 法

2.1.ApEn方法

ApEn是一種基于邊緣概率分布統(tǒng)計(jì)量化時(shí)間序列復(fù)雜程度的方法，其算法如下［10—12］:

1.對時(shí)間序列{u(i)，i=1，2，…，N}進(jìn)行相空間重構(gòu)，重構(gòu)維數(shù)為 m，據(jù)此可構(gòu)造一組維數(shù)為 m的新向量 X(1)，X(2)，…，X(N － m+1)，其中，

2.計(jì)算任意向量X(i)與其余向量X(j)之間的相對歐式距離 d［X(i)，X(j)］;

3.設(shè)定容許偏差為 r，統(tǒng)計(jì)每個(gè)向量 X(i)的d［X(i)，X(j)］小于r的數(shù)目，求出該數(shù)目與向量總數(shù)N－m+1的比Cmi(r);4.將Cmi(r)取對數(shù)，再求其對所有i的平均值，記作m(r);

顯然，ApEn值與m，r的取值有關(guān).Pincu建議取m=2，r=0.1—0.2σ(σ是原始數(shù)據(jù)的標(biāo)準(zhǔn)偏差).在本文中，取參數(shù)m=2，r=0.15σ.

近似熵方法反映了時(shí)間序列在模式上的自相似程度，以及當(dāng)維數(shù)變化時(shí)，產(chǎn)生新模式的可能性大小.ApEn值愈大，說明產(chǎn)生新模式的概率愈大，序列越復(fù)雜，系統(tǒng)可預(yù)測性越差.它給出新模式發(fā)生率隨維數(shù)而增減的情況，因而反映數(shù)據(jù)在結(jié)構(gòu)上的復(fù)雜性.同時(shí)它并不是企圖完全重構(gòu)吸引子，所以對各種非線性時(shí)間序列有很好的適用性［36］.

2.2.MC-ApEn方法

類似于滑動移除去趨勢波動分析方法、滑動移除重標(biāo)極差方法，MC-ApEn方法的計(jì)算過程如下:

1.首先選擇滑動移除數(shù)據(jù)的窗口長度L;

2.然后從待分析時(shí)間序列的第 i(i=1，2，…，N－L+1，N為時(shí)間序列中的總記錄個(gè)數(shù))個(gè)數(shù)據(jù)開始連續(xù)移除L個(gè)數(shù)據(jù)，再將剩余N－L個(gè)數(shù)據(jù)直接連在一起得到一個(gè)新的時(shí)間序列;

3.利用近似熵方法計(jì)算新序列的ApEn值;

4.保持移除數(shù)據(jù)的窗口長度不變，以步長L逐步移動窗口，重復(fù) 2，3步操作，直到原序列結(jié)束為止;

5.通過1—4步操作可得到一個(gè)隨著窗口移動的ApEn序列;

6.基于不同動力學(xué)性質(zhì)的數(shù)據(jù)其復(fù)雜性大小不相同，而具有相同動力學(xué)性質(zhì)的數(shù)據(jù)的復(fù)雜性差異不大這一特點(diǎn)，結(jié)合步驟5中得到的ApEn序列判斷突變點(diǎn)或突變區(qū)間.

3.MC-ApEn在突變檢測中的性能測試

3.1.MC-ApEn在線性時(shí)間序列中的性能測試

為便于對新方法與滑動近似熵的檢測結(jié)果進(jìn)行比較，線性理想時(shí)間序列仍然由王啟光等［36］所構(gòu)造的方程獲取，

圖1(a)給出該理想時(shí)間序列的演化曲線，由(2)式可知，在 t=1001時(shí)，理想序列 IS1發(fā)生了動力學(xué)結(jié)構(gòu)突變，其控制方程由簡單的正弦函數(shù)突變?yōu)檎液瘮?shù)和余弦函數(shù)所組成的新方程.

圖1 線性時(shí)間序列IS1及其突變檢測 (a)理想時(shí)間序列 IS1;(b)滑動 ApEn檢測，滑動窗口長度為S=200，滑動步長為1;(c)Mann-Kendall檢測，UF為時(shí)間序列順序統(tǒng)計(jì)量，UB為逆序統(tǒng)計(jì)量，Uα線為置信區(qū)間，顯著性水平α=0.01

滑動ApEn對IS1的分析結(jié)果表明，ApEn能夠根據(jù)時(shí)間序列的復(fù)雜性來很好的區(qū)分線性動力系統(tǒng)的異同，當(dāng)子序列來源于同一個(gè)動力學(xué)系統(tǒng)時(shí)，隨著子序列的滑動其ApEn值的變化十分穩(wěn)定，而當(dāng)子序列中的數(shù)據(jù)來自于不同的動力學(xué)系統(tǒng)時(shí)，序列的復(fù)雜性顯然較單一動力學(xué)系統(tǒng)更高，因而其ApEn值也隨之增大(圖1(b)).基于此，可以很容易地判斷出原時(shí)間序列約在t=1000附近發(fā)生了動力學(xué)結(jié)構(gòu)突變，顯然若想進(jìn)一步準(zhǔn)確定位突變點(diǎn)則存在一定的困難.

從Mann-kendall方法對IS1的檢測結(jié)果中可知(圖1(c))，自t=995開始，理想序列呈現(xiàn)出明顯的下降趨勢突變，且突變點(diǎn)通過了α=0.01的顯著性檢驗(yàn).Mann-kendall檢測到的突變點(diǎn)(t=995)與真實(shí)突變點(diǎn) t=1001非常接近，這表明 Mann-kendall在某些情況下具備一定的檢測系統(tǒng)動力學(xué)結(jié)構(gòu)突變的能力.但必須指出的是，這并不意味著 Mannkendall方法是一種檢測時(shí)間序列動力學(xué)結(jié)構(gòu)突變的好方法.從本質(zhì)上而言，Mann-kendall方法僅僅是一種檢測趨勢突變的方法，若時(shí)間序列中不存在明顯的趨勢突變，則其無法檢測出時(shí)間序列中存在的動力學(xué)結(jié)構(gòu)突變，不再贅述.

從圖2中可以發(fā)現(xiàn)，四種不同窗口長度情形下的MC-ApEn檢測結(jié)果有一個(gè)共同特點(diǎn)，即在 t﹤1000時(shí)，移除數(shù)據(jù)后得到的 ApEn值明顯大于 t﹥1000時(shí)的情形，這表明以 t=1000為界，可以將ApEn序列分為兩個(gè)不同的演化階段.根據(jù)ApEn的物理含義，值越大意味著序列的復(fù)雜性越大.因此，移除數(shù)據(jù)后相對較小的ApEn值意味著所移除的數(shù)據(jù)相對而言復(fù)雜性更大，進(jìn)而可以判斷出t﹥1000后數(shù)據(jù)的復(fù)雜性強(qiáng)于 t﹤1000前的數(shù)據(jù).IS1的MC-ApEn檢測結(jié)果表明，移除等長度的數(shù)據(jù)，若具有相同動力學(xué)屬性，則移除數(shù)據(jù)對于時(shí)間序列的ApEn估算的影響大致相同，盡管在移除數(shù)據(jù)尺度較小時(shí)，得到的ApEn值會存在局部波動相對較大的現(xiàn)象(如圖2(a))，但就統(tǒng)計(jì)意義而言，不論在滑動移除數(shù)據(jù)的窗口長度較小時(shí)(如L=10和20)，還是窗口長度較大時(shí)(如L=50和100)，無需借助其他統(tǒng)計(jì)工具就能夠十分容易的將MC-ApEn得到的近似熵值分為明顯變化的兩個(gè)時(shí)段，進(jìn)而能夠大致的判斷出約在t=1000時(shí)發(fā)生了動力學(xué)結(jié)構(gòu)突變.由于圖2中所選擇的滑動移除窗口長度均屬于一個(gè)特殊的長度，因此，得到的突變點(diǎn)完全一致，即在t=1001處發(fā)生了突變.由此引出如下問題:對于任意窗口尺度結(jié)果是否仍然完全一致?結(jié)果顯然是否定的，但不同窗口尺度檢測到的突變點(diǎn)變化不大，即窗口尺度對檢測結(jié)果影響較小.

圖2 線性時(shí)間序列IS1的MC-ApEn檢測 (a)滑動移除窗口長度L=10;(b)L=20;(c)L=50;(d)L=100

為了更為準(zhǔn)確地提取IS1中突變點(diǎn)的位置，在大致掌握 IS1中突變點(diǎn)的位置后，在本節(jié)中利用MC-ApEn對IS1中部分時(shí)間序列進(jìn)行突變檢測，這里任取了IS1中的前1200個(gè)數(shù)據(jù)，圖3(a)和(b)給出了滑動移除窗口長度L=10和20時(shí)的MC-ApEn檢測結(jié)果，相比于圖2，這樣處理的優(yōu)勢在于可以降低整個(gè)待檢測時(shí)間序列的復(fù)雜性，并將更易于區(qū)分不同復(fù)雜度的時(shí)間序列，有利于準(zhǔn)確的判斷突變點(diǎn).正如前文所言，窗口長度 L=10和20是兩個(gè)比較特殊的尺度.鑒于此，分別在滑動移除窗口尺度L為1和2時(shí)再次進(jìn)行了MC-ApEn分析，結(jié)果表明，可以非常容易發(fā)現(xiàn)在t<1001時(shí)，滑動移除數(shù)據(jù)對于所分析時(shí)間序列的ApEn的影響非常有規(guī)律，即ApEn隨著滑動窗口的移動展現(xiàn)為一種準(zhǔn)周期性的波動，且其波動幅度也大致相同(圖3(c)和(d))，但自t=1001后，這種波動規(guī)律明顯不同于ApEn在此之前的演變，這意味著前后兩段時(shí)間序列的復(fù)雜性不盡相同，因而數(shù)據(jù)的移除對于ApEn計(jì)算的影響各異.基于此，可以準(zhǔn)確地判斷IS1在t=1001處發(fā)生了動力學(xué)結(jié)構(gòu)突變.

3.2.MC-ApEn在非線性時(shí)間序列中的性能測試

為了測試MC-ApEn在非線性時(shí)間序列中的突變檢測性能，構(gòu)造了理想時(shí)間序列IS2，序列總長為2000，前 1000 個(gè)數(shù)據(jù)由 Logistic映射［37］產(chǎn)生，其方程如下:

上式中初值x0=0.8，參數(shù)u=3.8.IS2的后1000個(gè)數(shù)據(jù)為滿足正態(tài)分布的隨機(jī)數(shù)，其隨時(shí)間的演化曲線已由圖4(a)給出.從其構(gòu)造可知，IS2在t=1001時(shí)，系統(tǒng)的演化突然由確定性方程變?yōu)殡S機(jī)行為，由ApEn的物理意義可知，IS2中后1000個(gè)數(shù)據(jù)的近似熵值要大于前1000個(gè)數(shù)據(jù)，即t=1001后系統(tǒng)的復(fù)雜性突然增大了，可預(yù)測性變小了.

圖3 在t<1201時(shí)，對線性時(shí)間序列IS1中前1200個(gè)樣本的MC-ApEn檢測 (a)滑動移除窗口長度L=10;(b)L=20;(c)L=1;(d)L=2

圖4 非線性時(shí)間序列IS2及其MC-ApEn檢測 (a)理想序列IS2，長度為2000;(b)L=10;(c)L=5;(d)L=2

圖4為理想時(shí)間序列 IS2的MC-ApEn檢測結(jié)果，從圖4(b)和(c)中不難發(fā)現(xiàn)，在滑動移除窗口尺度L=10和5時(shí)，ApEn的演變趨勢非常類似，在t=1000前后，數(shù)據(jù)的移除對于ApEn值的計(jì)算影響各不相同，由于確定性方程Logistic模型產(chǎn)生的數(shù)據(jù)其復(fù)雜性要明顯低于隨機(jī)數(shù)，因此，容易理解對于移除等長度的數(shù)據(jù)，移除隨機(jī)數(shù)后所獲取的 ApEn值更小一些，即部分隨機(jī)數(shù)的移除導(dǎo)致剩余數(shù)據(jù)的復(fù)雜性減小.而且隨著滑動移除窗口長度的增大，移除具有不同動力學(xué)性質(zhì)的數(shù)據(jù)后計(jì)算得到的ApEn值之間的差異明顯增大.隨后考察了更小尺度滑動移除窗口的情形，圖4(d)給出了IS2在移除窗口長度L=2時(shí)的 MC-ApEn檢測結(jié)果，與 L=10和5時(shí)的結(jié)果一致.這表明，MC-ApEn的檢測結(jié)果對于窗口尺度的依賴性較小，該方法不僅適用于線性時(shí)間序列的動力學(xué)結(jié)構(gòu)突變檢測，而且對于非線性時(shí)間序列的動力學(xué)結(jié)構(gòu)突變也同樣具有很強(qiáng)的識別能力.

圖5 非線性時(shí)間序列IS2的滑動ApEn檢測 (a)滑動窗口長度S=100;(b)S=200;(c)S=300;(d)S=400

作為比較，圖5給出了IS2的滑動ApEn檢測結(jié)果，從圖5(a)中可以發(fā)現(xiàn)，在滑動窗口尺度較小時(shí)，ApEn的演變隨著滑動窗口的移動其值不穩(wěn)定，變化較大，這可能是ApEn算法自身的缺陷(即數(shù)據(jù)過短，得到的近似熵值不穩(wěn)健)所造成.因此，難以據(jù)此準(zhǔn)確檢測IS2中的動力學(xué)結(jié)構(gòu)突變.逐漸增大滑動窗口長度后發(fā)現(xiàn)，所得到的ApEn值漸漸變得比較平穩(wěn)，而且對于IS2中的動力學(xué)結(jié)構(gòu)突變也能夠較好地反映出來.從圖5(b)，(c)和(d)中可以發(fā)現(xiàn)三種不同窗口尺度下，均可以根據(jù)ApEn值的演變將其大致分為三個(gè)不同階段，以滑動窗口長度S=300時(shí)為例，可以粗略地分為:t∈［300，1000］，［1001，1300］，［1301，2000］等三個(gè)演變階段，在第一和第三區(qū)間ApEn值變化相對較為穩(wěn)定，而且窗口長度越大，這種穩(wěn)定性越好，而第二個(gè)區(qū)間為一個(gè)過渡階段，這一特征主要是由于在滑動的過程中，窗口序列中所包含的數(shù)據(jù)由完全是確定性方程產(chǎn)生的數(shù)據(jù)變?yōu)殡S機(jī)信號與確定信號并存，這一過程直至窗口序列中的數(shù)據(jù)全部變?yōu)殡S機(jī)信號，因此，窗口序列的ApEn值也隨之增大.基于ApEn演變的這一特征，可以大致確定突變開始于t=1000附近，但要想精確定位突變點(diǎn)還存在一定的困難.王啟光等［36］提出采用滑動-t檢驗(yàn)法來進(jìn)一步定位突變點(diǎn)，但由于滑動-t檢驗(yàn)的結(jié)果嚴(yán)重依賴于子序列長度的選取，故而未能從根本上解決此問題.因此，這是影響滑動ApEn方法在突變檢測中廣泛應(yīng)用的瓶頸之一.

當(dāng)采用 Mann-Kendall方法對非線性時(shí)間序列IS2進(jìn)行突變檢測時(shí)，發(fā)現(xiàn)順序統(tǒng)計(jì)量UF和逆序統(tǒng)計(jì)量UB雖然在t=993處存在一個(gè)交叉點(diǎn)，但其并未位于兩條信度線之間，這意味著IS2中并不存在明顯的趨勢突變.隨后，對 IS2中 t∈［501，1500］區(qū)間內(nèi)的數(shù)據(jù)進(jìn)行了Mann-Kendall檢測，發(fā)現(xiàn)統(tǒng)計(jì)量UF和UB之間的交叉點(diǎn)為t=1005，且位于兩條信度線之間，而UF在此處正好呈現(xiàn)顯著的下降趨勢，據(jù)此可以判斷序列在t=1005處存在一個(gè)下降的趨勢突變.對同一時(shí)間序列不同長度數(shù)據(jù)的 Mann-Kendall分析卻得到了完全不同的結(jié)論，表明了Mann-Kendall的分析結(jié)果嚴(yán)重依賴于時(shí)間序列的長度，顯然該方法不適合用于檢測時(shí)間序列的動力學(xué)結(jié)構(gòu)突變.

圖6 非線性時(shí)間序列IS2的Mann-Kendall檢測 (a)IS2的Mann-Kendall檢測;(b)對 IS2中 t∈［501，1500］區(qū)間內(nèi)數(shù)據(jù)的Mann-Kendall檢測

3.3.應(yīng)用實(shí)例

前文中所給出的兩個(gè)例子均為理想化序列，而實(shí)際觀測資料要比理想序列復(fù)雜得多.鑒于此，本文測試了MC-ApEn在降水觀測資料中的突變檢測性能.采用的降水資料來源于中國氣象局國家氣象信息中心資料室，對于資料中少數(shù)缺測記錄視為無降水.圖7給出了呼瑪氣象觀測站1960—2006年逐日降水資料的MC-ApEn檢測結(jié)果，與前文理想試驗(yàn)中類似的情況是，逐日降水資料的ApEn序列可以分為明顯不同的兩個(gè)穩(wěn)定變化階段，1979年以前的 ApEn值普遍大于1980年以后的情形，這意味著自1980年后至2006年期間內(nèi)，降水的復(fù)雜性突然增大了，其可預(yù)測性相對于1960—1979年而言減小了.基于此，可以判斷1980年呼瑪氣象觀測站發(fā)生了氣候突變，這與已有的研究結(jié)果基本一致［36，38］.而圖 7中滑動移除窗口長度為1年和2年時(shí)檢測到的突變點(diǎn)完全一樣，則進(jìn)一步表明本文提出的新方法其檢測結(jié)果對于窗口長度的依賴性較小.

圖7 呼瑪氣象觀測站逐日降水資料的MC-ApEn檢測，資料記錄了1960年1月1日至2006年12月31日呼瑪站的逐日降雨量 (a)滑動移除窗口長度為1年，這里取L=365天;(b)L=730天

4.結(jié)果與討論

本文基于一種非線性動力學(xué)指數(shù)——ApEn能夠反映不同動力學(xué)系統(tǒng)的復(fù)雜性特征這一特性，借鑒于滑動移除技術(shù)，提出了一種新的突變檢測方法——MC-ApEn方法.理想試驗(yàn)結(jié)果表明，滑動ApEn方法對線性和非線性時(shí)間序列的動力學(xué)結(jié)構(gòu)突變具有一定的識別能力，但是存在不能夠準(zhǔn)確定位突變點(diǎn)的問題，甚至?xí)霈F(xiàn)在滑動窗口長度較小時(shí)難以檢測突變的現(xiàn)象，因此這一缺陷制約了該方法在實(shí)際中的應(yīng)用;而MC-ApEn適合線性和非線性時(shí)間序列的動力學(xué)結(jié)構(gòu)突變檢測，其檢測結(jié)果對于滑動移除窗口尺度的依賴性較小，并能夠更為精確的檢測突變開始的時(shí)間，新方法檢測突變的性能明顯優(yōu)于滑動近似熵方法.Mann-Kendall方法雖然能夠?qū)Ρ疚闹械木€性理想時(shí)間序列中的動力學(xué)結(jié)構(gòu)突變進(jìn)行檢測，但其本質(zhì)上是一種針對時(shí)間序列趨勢突變的檢測算法，如果時(shí)間序列中不存在明顯的趨勢變化，則不可能有效的檢測序列中的動力學(xué)結(jié)構(gòu)突變.而且Mann-Kendall方法的檢測結(jié)果依賴于待分析時(shí)間序列的長度，因而它不是一種有效的動力學(xué)結(jié)構(gòu)突變檢測方法.MC-ApEn在實(shí)際降水觀測資料中的應(yīng)用則進(jìn)一步證實(shí)了新方法的有效性，展示了其潛在的應(yīng)用前景.

［1］ Wang G R，Cheng S G，Hao B L 1984 Acta Phys.Sin.33 1246(in Chinese)［王光瑞、陳式剛、郝柏林 1984物理學(xué)報(bào) 33 1246］

［2］ Cheng S G 1991 Acta Phys.Sin.40 584(in Chinese)［陳式剛1991物理學(xué)報(bào)40 584］

［3］ Wang A L，Yang C X 2002 Acta Phys.Sin.51 2719(in Chinese)［王安良、楊春信2002物理學(xué)報(bào) 51 2719］

［4］ Wan S Q，F(xiàn)eng G L，Dong W J，Li J P 2005 Acta Phys.Sin.54 5487(in Chinese)［萬仕全、封國林、董文杰、李建平 2005物理學(xué)報(bào) 54 5487］

［5］ Gong Z Q，F(xiàn)eng G L，Dong W J，Li J P 2006 Acta Phys.Sin.55 3180(in Chinese)［龔志強(qiáng)、封國林、董文杰、李建平 2006物理學(xué)報(bào) 55 3180］

［6］ Song A J，Han L 2007 Acta Phys.Sin.56 3820(in Chinese)［宋愛軍、韓 雷2007物理學(xué)報(bào) 56 3820］

［7］ Yu S M 2008 Acta Phys.Sin.57 6859(in Chinese)［禹思敏2008物理學(xué)報(bào)57 6859］

［8］ ShiN 2009 Meteorology statistical forecasting (Beijing:Metrological Press)p116(in Chinese)［施能2009氣象統(tǒng)計(jì)預(yù)報(bào)(北京:氣象出版社)第116頁］

［9］ Lin Z S 2003 Nonlinear sience and its application in Geosciences(Beijing:Metrological Press)p113(in Chinese)［林振山 2003非線性科學(xué)及其在地學(xué)中的應(yīng)用(北京:氣象出版社)第113頁］

［10］ Pincus S M 1991 Proc.Natl.Acad.Sci.USA 88 2297

［11］ Pincus S M，Goldberger A L 1994 Am.J.Physiol.266 H1643

［12］ Pincus SM 1995 Chaos 5 110

［13］ Lei G X，Zeng Y C，Li L 2007 Technical Acoustics 26 121(in Chinese)［雷國雄、曾以成、李 凌 2007聲學(xué)技術(shù) 26 121］

［14］ Li Q，Wang T Y，Leng Y G，He G Y，He H L 2007 Acta Phys.Sin.56 6803(in Chinese)［李 強(qiáng)、王太勇、冷永剛、何改云、何慧龍2007物理學(xué)報(bào)56 6803］

［15］ Li X F，Yu W L 2008 Acta Phys.Sin.57 2587(in Chinese)［劉小峰、俞文莉2007物理學(xué)報(bào) 57 2587］

［16］ Xie Y，Xu J X，Yang H J，Hu S J 2002 Acta Phys.Sin.51 205(in Chinese)［謝 勇、徐健學(xué)、楊紅軍、胡三覺 2002物理學(xué)報(bào)51 205］

［17］ Cao B，Lv X Q，Zeng M，Wang Z M，Huang S S 2006 Acta Phys.Sin.55 1696(in Chinese)［曹 彪、呂小青、曾 敏、王振民、黃石生 2006物理學(xué)報(bào) 55 1696］

［18］ Feng G L，Dong W J，Li J P 2004 Chin.Phys.13 1582

［19］ Zeng H L，Gao X Q，Zhang W 2005 Chin.Phys.14 1265

［20］ Dai X G，F(xiàn)u C B，Wang P 2005 Chin.Phys.14 850

［21］ Gao X Q，Zhang W 2005 Chin.Phys.14 2370

［22］ Shi N，Gu J Q，Yi Y M，Lin Z N 2005 Chin.Phys.14 844

［23］ Shi N，Yi Y M，Gu J Q，Xia D D 2006 Chin.Phys.15 2180

［24］ Feng G L，He W P 2007 Chin.Phys.16 2825

［25］ Zhang W，Wan S Q 2008 Chin.Phys.B 17 2311

［26］ Feng G L，Gong Z Q，Zhi R，Zhang D Q 2008 Chin.Phys.B 17 2745

［27］ Feng G L，Dai X G，Wang A H 2001 Acta Phys.Sin.50 606(in Chinese)［封國林、戴新剛、王愛慧 2001物理學(xué)報(bào) 50 606］

［28］ Feng G L，Dong W J，Jia X J 2002 Acta Phys.Sin.51 1181(in Chinese)［封國林、董文杰、賈曉靜 2002物理學(xué)報(bào) 51 1181］

［29］ Feng G L，Dong W J，Gong Z Q，Hou W，Wan S Q，Zhi R 2006 Nonlinear theories and methods on spatial-temporal distribution of the obserbational data (Beijing:Metrological Press)p27—56(in Chinese)［封國林、董文杰、龔志強(qiáng)、侯威、萬仕全、支蓉 2006觀測數(shù)據(jù)非線性時(shí)空分布理論和方法(北京:氣象出版社)第27—56頁］

［30］ Zhi R，Lian Y，F(xiàn)eng G L 2007 Acta Phys.Sin.56 1837(in Chinese)［支 蓉、廉 毅、封國林 2007物理學(xué)報(bào) 56 1837］

［31］ Feng G L，Wang Q G，Hou W，Gong Z Q，Zhi R 2009 Acta.Phys.Sin 58 2853(in Chinese)［封國林、王啟光、侯 威、龔志強(qiáng)、支 蓉 2009物理學(xué)報(bào) 58 2853］

［32］ He W P2008 Ph. D.Dissertation (Lanzhou:Lanzhou University)p23(in Chinese)［何文平 2008博士學(xué)位論文(蘭州:蘭州大學(xué))第23頁］

［33］ He W P，F(xiàn)eng G L，Wu Q，Wan S Q，Chou J F 2008 Non.Proc.Geophys.15 601

［34］ He W P，Wu Q，Zheng W，Wang Q G，Zhang Y 2009 Acta Phys.Sin.58 2862(in Chinese)［何文平、吳 瓊、張 文、王啟光、張 勇 2009物理學(xué)報(bào) 58 2862］

［35］ He W P，Deng B S，Wu Q，Zheng W，Cheng H Y 2010 Acta Phys.Sin.59 2864(in Chinese)［何文平、鄧北勝、吳 瓊、張文、成海英 2010物理學(xué)報(bào) 59 2864］

［36］ Wang Q G，Zhang Z P 2008 Acta Phys.Sin.57 1996(in Chinese)［王啟光、張?jiān)銎?2008物理學(xué)報(bào) 57 1996］

［37］ May R 1976 Nature 261 459

［38］ Zhi R，Gong Z Q 2008 Acta Phys.Sin.57 4629(in Chinese)［支 蓉 、龔志強(qiáng) 2008物理學(xué)報(bào) 57 4629］

PACS:92.60.Wc

A new method to detect abrupt change based on approximate entropy*

He Wen-Ping1)He Tao2)Cheng Hai-Ying3)Zhang Wen4)Wu Qiong5)
1)(National Climate Center，China Meteorological Administration，Beijing 100081，China)2)(JiNan Environment Protection Monitoring Station，Jinan 250014，China)3)(Department of Basic Course，Yancheng Institute of Technology，Yancheng 224002，China)4)(TongDa College，Nanjing University of post＆ telecommunications，Nanjing 210003，China)5)(National Satellite Meteorological Center，China Meteorological Administration，Beijing 100081，China)(Received 6 June 2010;revised manuscript received 12 July 2010)

Approximate entropy(ApEn)is valid index which can be used to quantitatively reflect dynamic characteristics and complexity of a time series.The ApEn has been developed to detect an abrupt change in one-dimension time series by sliding a fixed widow，which can be identified with an abrupt dynamic change to some extent，but the sliding ApEn results depend on the window scale，and cannot accurately position the time-instant of an abrupt change.Based on this，a new method is proposed in the present paper，i.e.，moving cut data-approximate entropy(MC-ApEn)，which can be used to detect an abrupt dynamic change in time series.Tests on model time series indicate that the detection results from the present method show relatively good stability and high accuracy，obviously better than those from the sliding ApEn method and the Mann-Kendall method.The applications in daily precipitation records further verify the validity of the present method.

approximate entropy，moving cut data-approximate entropy，abrupt change detection

*國家自然科學(xué)基金(批準(zhǔn)號:40905034，40875040，40930952)，公益性行業(yè)(氣象)科研專項(xiàng)基金(批準(zhǔn)號:GYHY200806005，GYHY200906019)資助的課題.

E-mail:wenping_he@163.com

*Project supported by the National Natural Science Foundation of China(Grant Nos.40905034，40875040，40930952)，the Special Scientific Research Fund of Meteorological Public Welfare Profession of China(Grant Nos.GYHY200806005，GYHY200906019).

E-mail:wenping_he@163.com.