蔡舒平 戴理 倪淏

1)(江蘇大學(xué)電氣信息工程學(xué)院，鎮(zhèn)江 212013)

2)(江蘇大學(xué)，機(jī)械工業(yè)設(shè)施農(nóng)業(yè)測控技術(shù)與裝備重點(diǎn)實(shí)驗(yàn)室，鎮(zhèn)江 212013)

(2013年1月8日收到;2013年6月14日收到修改稿)

1 引言

由于自然界本身和人類活動的加劇，許多復(fù)雜系統(tǒng)如生態(tài)系統(tǒng)、氣象系統(tǒng)等常常表現(xiàn)為從一個(gè)相對穩(wěn)定狀態(tài)翻轉(zhuǎn)到另一個(gè)相對穩(wěn)定狀態(tài)的現(xiàn)象.生態(tài)系統(tǒng)中一個(gè)典型的例子是在一個(gè)淺水湖泊中，由于磷的流入，生態(tài)系統(tǒng)有可能從一個(gè)寡營養(yǎng)狀態(tài)翻轉(zhuǎn)到一個(gè)富營養(yǎng)狀態(tài)而導(dǎo)致藻類過度生長，進(jìn)而由于缺乏魚類生存所需要的氧氣而使生物的多樣性減少和水質(zhì)量惡化，隨之而來的是系統(tǒng)的大部分功能和效用喪失.氣候系統(tǒng)中如全球氣候在溫暖而潮濕的間冰期和涼爽而干燥的冰川期之間快速地翻轉(zhuǎn)等等.近年來，對復(fù)雜系統(tǒng)狀態(tài)翻轉(zhuǎn)的研究因其對現(xiàn)代社會的影響而引起了人們廣泛的關(guān)注.科學(xué)工作者通過深入研究后發(fā)現(xiàn)，所有這些系統(tǒng)的狀態(tài)翻轉(zhuǎn)具有一個(gè)共同的特征，那就是其內(nèi)在的動力學(xué)結(jié)構(gòu)隨著外強(qiáng)迫作用的變化而發(fā)生了突變.因此，有關(guān)動力學(xué)結(jié)構(gòu)突變檢測的研究近年來受到了普遍重視，其中氣候系統(tǒng)因其本身的復(fù)雜性及與人類日常生活的密切相關(guān)性，使得在這方面的研究顯得尤為迫切.

封國林等[1]從氣候突變的檢測和研究、氣候資料信息的提取和分離、氣候系統(tǒng)復(fù)雜性的研究、氣候系統(tǒng)動力學(xué)結(jié)構(gòu)特征及極端氣候事件的檢測等5個(gè)方面闡述了近年來在我國氣候突變檢測技術(shù)研究方面所取得的新進(jìn)展.傳統(tǒng)的檢測方法如滑動t-檢驗(yàn)，Cramer法，Yamamoto信噪比法以及Mann-Kendall(M-K)法等，在突變檢測的早期得到了廣泛應(yīng)用，然而由于這些傳統(tǒng)方法自身固有的一些缺陷使得在實(shí)際檢測中所得結(jié)果均存在不同程度的突變點(diǎn)漂移，在此情況下，一些新的檢測方法如分割方法(即BG算法)應(yīng)運(yùn)而生，龔志強(qiáng)等[2]對此進(jìn)行了實(shí)際檢驗(yàn)和應(yīng)用，并得到了一些有益的結(jié)果.近年來，國內(nèi)學(xué)者如封國林、龔志強(qiáng)、何文平、王啟光和侯威等人將非線性科學(xué)的一些研究成果應(yīng)用于動力學(xué)結(jié)構(gòu)突變檢測中，分別提出了動力學(xué)相關(guān)因子指數(shù)法[3]、動力學(xué)指數(shù)分割算法[4]、滑動去趨勢波動算法[5]、排列熵算法[6]、近似熵算法[7]、滑動移除近似熵算法[8]、重標(biāo)極差法[9]等檢測方法.成海英等[10]將偏度系數(shù)和峰度系數(shù)引入時(shí)間序列的突變檢測中，提出了基于概率密度分布型變化的突變檢測方法.此外，還有學(xué)者將小波分析、粒子濾波等方法用于突變檢測中[11,12].在與氣候突變有關(guān)的其他研究領(lǐng)域，龔志強(qiáng)等[13]把希爾伯特-黃變換(Hilbert-Huang transform，HHT)應(yīng)用于氣候資料信息的分離和提取中，且與小波變換進(jìn)行了對比分析，并針對它們各自的缺點(diǎn)提出了可能改進(jìn)的措施.與此同時(shí)，封國林等[14]介紹了一種能判斷不同時(shí)間序列之間動力學(xué)異同性的條件熵算法，討論了其物理意義，并以長三角地區(qū)為例，驗(yàn)證了這一算法的有效性和可靠性.萬仕全等[15]利用Q指數(shù)來識別時(shí)間序列之間動力結(jié)構(gòu)的異同性，并利用該方法分析了幾種典型代用資料的動力學(xué)結(jié)構(gòu)，得到了一些有益的結(jié)果.而在極端氣候事件的檢測方面，封國林等[16]利用固定閾值法研究了Lorenz系統(tǒng)極端事件序列的長程相關(guān)性特征，并在此基礎(chǔ)上分析了實(shí)際的氣象觀察資料，得到了一些有參考價(jià)值的結(jié)論.所有這些研究工作極大地豐富了突變檢測方法，拓寬了突變檢測的途徑.然而，由于所研究問題的復(fù)雜性，不可能存在一種萬能的方法一勞永逸地解決所有問題，因此，人們從未停止過對突變檢測新思想、新方法的探尋，尤其是分析那些已被噪聲污染且稀疏的實(shí)際野外觀察數(shù)據(jù)方法的研究.

實(shí)際中往往會面臨單變量的一維時(shí)間序列，很顯然，系統(tǒng)的非線性突變特征信息就蘊(yùn)藏在這些序列中，如何有效地從中提取這類信息成為一個(gè)非常關(guān)鍵的問題.

Fisher信息(FI)恰好為我們提供了表征這樣一類突變特征信息的可能性.原因是任何類型的數(shù)據(jù)和模型本質(zhì)上都可以轉(zhuǎn)換為信息而不管最初的學(xué)科是什么[17].不像系統(tǒng)信息的其他測量方法，F(xiàn)isher信息提供了一種通過監(jiān)測系統(tǒng)變量從而監(jiān)測系統(tǒng)的狀態(tài)和狀態(tài)突變的方法[18]，這種檢測狀態(tài)和狀態(tài)突變的能力使得識別發(fā)生在系統(tǒng)中的最根本的變化成為可能.事實(shí)上，F(xiàn)rieden[19,20]使用該理論推導(dǎo)出了許多基本的物理學(xué)、熱力學(xué)、人口遺傳學(xué)方程.近年來，生態(tài)學(xué)家廣泛使用Fisher信息理論來研究生態(tài)領(lǐng)域的問題，F(xiàn)ath等[18]和Mayer等[21]使用Fisher信息來作為一種復(fù)雜生態(tài)系統(tǒng)動態(tài)秩序的度量.Mayer等[21]和Karunanithi等[22]提議用來作為一種量化的指標(biāo)來檢測和評估生態(tài)系統(tǒng)狀態(tài)的翻轉(zhuǎn)以及作為一種可持續(xù)性標(biāo)準(zhǔn)[23].Fisher信息還被用于研究包含一個(gè)多間隔食物鏈穩(wěn)定性問題的模型系統(tǒng)[24-26]以及用于可持續(xù)環(huán)境管理中的動態(tài)模型系統(tǒng)的優(yōu)化控制問題[26,27].Rico-Ramirez等[28]把Fisher信息應(yīng)用于醫(yī)療化工領(lǐng)域的優(yōu)化控制問題.在了解了以上為研究生態(tài)系統(tǒng)嘗試使用Fisher信息理論所取得的成果后，本文試圖用它來解決動力學(xué)結(jié)構(gòu)的突變檢測問題.

2 方法

2.1 Fisher信息

統(tǒng)計(jì)學(xué)家Ronald Fisher在1922年提出了一種度量不確定性的方法，即現(xiàn)在稱之為費(fèi)歇信息[29].然而，由于信息只能從代表現(xiàn)行動力系統(tǒng)穩(wěn)定性的數(shù)據(jù)中獲得，使得Fisher信息尤其被認(rèn)為是一種對現(xiàn)行動力系統(tǒng)秩序的度量.許多文獻(xiàn)已經(jīng)展示了使用Fisher信息來度量生態(tài)系統(tǒng)的秩序的能力[21,22,30].一個(gè)變量單次測量的Fisher信息I計(jì)算如下：

這里P(s)為概率密度函數(shù)(PDF)，s是一個(gè)狀態(tài)變量.

對于一個(gè)穩(wěn)定的動力系統(tǒng)而言，系統(tǒng)變量的概率密度分布也較為穩(wěn)定，而當(dāng)動力學(xué)結(jié)構(gòu)發(fā)生變化后，系統(tǒng)變量的概率密度分布也會發(fā)生不同程度的變化.由于Fisher信息是基于概率密度函數(shù)的導(dǎo)數(shù)，因此它反映的是一種局部特性，這使得它能敏銳捕捉系統(tǒng)變量的概率密度分布發(fā)生微小的變化.鑒于此，本文用它檢測和識別蘊(yùn)藏在時(shí)間序列中的突變特征，為突變檢測提供一種新的途徑.

實(shí)際計(jì)算中，為了避免因除以較小的P(s)值而帶來的計(jì)算誤差，我們用一個(gè)變量的強(qiáng)度來代替(1)式中概率密度函數(shù)P(s)，即令q2(s)=P(s)，于是(1)式變?yōu)?/p>

為了使(2)式適合于分析和檢測實(shí)際系統(tǒng)的觀察序列，我們用求和來代替積分以及用差分來代替微分，即用ΔS=Si-Si+1和Δq=qi-qi+1來分別代替ds和dq，從而Fisher信息能被近似為

(3)式中，Si代表系統(tǒng)的某一個(gè)特定狀態(tài)，即S1代表狀態(tài)1，S2代表狀態(tài)2等等.于是：Si-Si+1=1，則計(jì)算Fisher信息的最終表達(dá)式為

表達(dá)式(4)將會用于我們以下所有的計(jì)算.

2.2 Fisher信息的計(jì)算

檢測系統(tǒng)動力學(xué)結(jié)構(gòu)突變的第一步是獲取描述系統(tǒng)狀態(tài)隨時(shí)間變化的觀察數(shù)據(jù)序列集，然后在數(shù)據(jù)集上定義一個(gè)滑動時(shí)間窗口，窗口寬度取決于可得到的數(shù)據(jù)量和系統(tǒng)的行為.按照經(jīng)驗(yàn)，窗口寬度至少應(yīng)該包含8個(gè)數(shù)據(jù)點(diǎn)以確保窗口中的每一點(diǎn)不會過度影響整個(gè)計(jì)算.窗口的滑動因子應(yīng)小于窗口寬度，這樣相鄰的窗口間就會出現(xiàn)一個(gè)重疊，這樣做的好處是能夠捕獲有可能已經(jīng)延伸到該窗口邊界之外的動力學(xué)結(jié)構(gòu)的突變.

即若設(shè)觀察數(shù)據(jù)序列集為D={d(k),k=1,···,N}，其中N為序列總長度，窗寬為w∈N，滑動因子為δ∈N，用數(shù)學(xué)公式可把滑動窗W描述如下：

式中m=1，2，···，M，M為窗口個(gè)數(shù)，M=(N-w)/δ.

在把時(shí)間序列轉(zhuǎn)變?yōu)橐粋€(gè)重疊窗口序列以后，接下來是將每個(gè)滑動窗口劃分成如下I個(gè)區(qū)間：

式中{Zi=[Si-1,Si),i=1,2,···,I}，互不相交.

則數(shù)據(jù)d(k)∈W(m,w,δ)落入?yún)^(qū)間Zi的概率P(Zi)等于d(k)∈W(m,w,δ)落入?yún)^(qū)間Zi的數(shù)目與W(m,w,δ)中總數(shù)據(jù)個(gè)數(shù)w之比值.

即使用如下公式計(jì)算窗口中相應(yīng)于Zi區(qū)間上的概率密度分布P(Zi)：

需要說明的是，在上述計(jì)算中，窗口中區(qū)間個(gè)數(shù)I的取值越大，區(qū)間越窄，相同條件下，所計(jì)算出的窗口概率密度分布比較平緩，從而Fisher信息值較小，這使得其對突變的反應(yīng)不夠敏感.反之，I取值越小，區(qū)間越寬，相同條件下，所計(jì)算出的窗口概率密度分布比較陡峭，從而Fisher信息值較大，而這又容易受到噪聲的干擾而引起誤判.一般情況下，I的取值應(yīng)根據(jù)待分析數(shù)據(jù)的最大波動范圍而定.

接下來，為每個(gè)狀態(tài)計(jì)算q(Zi)(q(Zi)=以及使用(4)式為每個(gè)時(shí)間窗口估計(jì)FI.值得注意的是計(jì)算的FI值被分配給每個(gè)時(shí)間窗口的末端位置.

歸納起來，計(jì)算FI的基本步驟如下：1)把時(shí)間序列數(shù)據(jù)劃分為一連串重疊的時(shí)間窗口，2)把每個(gè)時(shí)間窗口劃分成長度相等的若干個(gè)區(qū)間，清點(diǎn)落入每個(gè)區(qū)間的數(shù)據(jù)個(gè)數(shù);3)計(jì)算相應(yīng)于該窗口的概率密度分布;4)對每個(gè)時(shí)間窗口使用Fisher信息實(shí)用計(jì)算公式從概率密度函數(shù)中計(jì)算Fisher信息.

3 結(jié)果

3.1 FI方法在理想時(shí)間序列中的性能測試

在所有的突變測試函數(shù)中，正余玄復(fù)合函數(shù)因其構(gòu)造簡單，表現(xiàn)形式多樣，物理意義明確而得到了廣泛的應(yīng)用.為了便于對照，本文亦選擇它為測試函數(shù)來檢驗(yàn)新方法的有效性.

(11)式表明，上述時(shí)間序列來自兩種不同的動力學(xué)結(jié)構(gòu)，在區(qū)間[1，1000]上，由單一的正玄函數(shù)所組成，在區(qū)間(1000，2000]上，由正余玄函數(shù)復(fù)合而成.很顯然，在t=1000時(shí)，系統(tǒng)從一種形式的動力學(xué)結(jié)構(gòu)突變?yōu)榱硪环N形式的動力學(xué)結(jié)構(gòu).圖1(a)給出了該理想信號序列波形.為了進(jìn)一步測試新方法的抗噪能力，在原理想信號的基礎(chǔ)上，同時(shí)構(gòu)造如下兩種形式的加噪信號.圖1(b)為加入信噪比為30 dB高斯白噪聲后信號波形，圖1(c)為隨機(jī)加入12個(gè)尖峰噪聲后的信號波形，其大小在2—5之間變化.

圖2為采用Fisher信息方法在滑動窗口寬度w=20，滑動步長δ=1情況下，對3種信號進(jìn)行突變檢測所得結(jié)果.

圖1 測試信號 (a)理想信號;(b)加白噪聲的理想信號;(c)加尖峰噪聲的理想信號

從圖2(a)所呈現(xiàn)的檢測結(jié)果可以看出，在區(qū)間[1，1000]上，F(xiàn)I值表現(xiàn)出疏密有致的有規(guī)律波動，而在區(qū)間(1000，2000]上，F(xiàn)I值表現(xiàn)出疏密不等的無規(guī)律變化，在n=1000處，出現(xiàn)了一個(gè)幅度為4(最大的FI值)的尖峰.而由Fisher信息的定義可知，F(xiàn)isher信息是基于序列概率密度分布的導(dǎo)數(shù)，這意味著在上述兩個(gè)區(qū)間上，序列的概率密度分布肯定不同，這種不同意味著在整個(gè)區(qū)間[1，2000]上，序列的概率密度分布在n=1000處會出現(xiàn)一次“躍變”.我們知道，躍變函數(shù)的一階導(dǎo)數(shù)剛好是沖擊函數(shù)，這就不難理解為什么在n=1000處會出現(xiàn)一個(gè)尖峰沖擊.由此看來，基于Fisher信息的突變檢測方法深刻地揭示了突變發(fā)生的機(jī)理，如實(shí)地反映了突變發(fā)生的整個(gè)物理進(jìn)程.另一方面，F(xiàn)isher信息又是一種對信號穩(wěn)定程度的一種度量，從這一點(diǎn)來講，突變前序列的穩(wěn)定性要優(yōu)于突變后.

從圖2(b)的測試結(jié)果來看，該方法具有一定的抗噪能力.和圖2(a)相比較，突變前后因其噪聲的加入使得在一定程度上會影響序列的概率密度分布，又因?yàn)镕isher信息是基于概率密度分布的導(dǎo)數(shù)，反映的是一種局部特性，因此，它對任何影響序列概率密度分布發(fā)生微小變化的因素較敏感，其結(jié)果是突變前后FI值較圖2(a)出現(xiàn)了一定程度的波動，尤其是突變前，F(xiàn)I值呈現(xiàn)出一種毫無規(guī)律的小幅尖峰震蕩.但該方法對噪聲的這種敏感性僅限于此，絲毫不會影響突變檢測的最終結(jié)果，這一點(diǎn)從圖2(b)的FI分布趨勢與圖2(a)中完全一致就說明了這一點(diǎn).

圖2 三種信號的FI突變檢測 (a)理想信號;(b)加白噪聲信號;(c)加尖峰噪聲信號

圖2(c)展示的該突變檢測方法抗尖峰噪聲的能力.從圖中可以看到，信號中的尖峰噪聲僅僅會對突變前后FI值造成輕微的影響，并使突變點(diǎn)n=1000處FI值的躍變幅度有所減弱，盡管如此，從圖中仍可以清楚地識別出突變點(diǎn).有趣的是尖峰噪聲對突變前后FI值的影響恰好與白噪聲相反，這是由FI值的計(jì)算方法所決定的.FI是以窗口為單位進(jìn)行計(jì)算的，單個(gè)的尖峰噪聲對整個(gè)窗口的影響被平均掉了，這樣看來，窗口中的單一信號比復(fù)合信號能更有效地抑制尖峰噪聲的影響，這也就是圖2(c)中所看到的突變前FI的變化要比突變后FI的變化更平穩(wěn)的原因.

在以上測試過程中，子序列長度即窗口寬度固定為一個(gè)特定的值(w=20)，得出了完全相同的突變檢測結(jié)果.由此引出了另外一個(gè)問題，那就是對于任意的子序列長度檢測結(jié)果是否仍然完全相同呢?為此，需要考察不同子序列長度下該方法對同一信號的突變檢測結(jié)果.

圖3給出了子序列長度分別取10，50，100和200時(shí)，對理想信號的突變檢測結(jié)果.

從圖3不難看出，四種不同情況下，F(xiàn)I的演化規(guī)律完全類似，即以突變點(diǎn)為界，F(xiàn)I劃分為兩個(gè)明顯不同的演化階段，在兩個(gè)不同階段的銜接處即為突變點(diǎn)，此處的FI值表現(xiàn)為一個(gè)急促向上的跳變，跳變幅度及突變前后FI波動幅度隨著子序列長度的增加有所減弱，這是因?yàn)樽有蛄腥〉迷酱螅有蛄兄械膯蝹€(gè)數(shù)據(jù)點(diǎn)對由整個(gè)子序列計(jì)算出來的FI值的影響愈小，F(xiàn)I值愈穩(wěn)定，但這絲毫不會影響最終的檢測結(jié)果.不同子序列長度下檢測結(jié)果的一致性，表明了該方法的檢測結(jié)果不依賴于子序列長度.

圖3 理想序列FI突變檢測 (a)w=10;(b)w=50;(c)w=100;(d)w=200

一種方法之所以被提出必有其不同于傳統(tǒng)方法的獨(dú)特之處.為此，有必要把該方法與一些傳統(tǒng)方法進(jìn)行比較，以檢驗(yàn)其有效性.仍以上述線性理想序列為例，分別基于滑動t檢驗(yàn)和Yamamoto法對該序列進(jìn)行檢測，其中子序列長度取200，滑動步長取1，圖4所示為這兩種方法的檢測結(jié)果.

在圖4(b)和(c)中，對應(yīng)突變點(diǎn)t=1000，滑動t檢驗(yàn)和Yamamoto法均給出了一個(gè)明顯的突變區(qū)域(800—1200之間)，在突變點(diǎn)附近，無論t值還是信噪比均出現(xiàn)了最大值，且超過了0.01顯著性水平.然而，若想進(jìn)一步準(zhǔn)確定位突變點(diǎn)則存在一定的困難.而由FI法在此種情況下所給出的檢測結(jié)果中(見圖3(d))，對突變點(diǎn)的檢測能夠做到準(zhǔn)確定位.

進(jìn)一步考察當(dāng)子序列的長度取不同值時(shí)，滑動t檢驗(yàn)和Yamamoto法的檢測結(jié)果是否仍然保持一致.圖5給出了當(dāng)子序列長度取20，滑動步長取1時(shí)，滑動t檢驗(yàn)和Yamamoto法對此的檢測結(jié)果.

圖4 子序列長度取200時(shí)檢測結(jié)果 (a)理想序列;(b)滑動t檢驗(yàn);(c)Yamamoto法

從圖5中可以看出，盡管滑動t檢驗(yàn)和Yamamoto法在此情況下仍能對突變點(diǎn)做出一定的反應(yīng)，但要對檢測結(jié)果做出合理的解釋已變得不可能;而FI法在此種情況下對突變點(diǎn)的檢測結(jié)果(見圖2(a))未發(fā)生任何改變.改變子序列的長度(如取50，100)會得出類似的結(jié)果(圖略).同時(shí)對在本質(zhì)上也是基于均值思想的傳統(tǒng)檢測方法——Crammer法也進(jìn)行了類似的實(shí)驗(yàn)，發(fā)現(xiàn)存在同樣的問題.而另一種傳統(tǒng)檢測方法——M-K法常常被用來檢測序列中的趨勢突變，若時(shí)間序列中不存在明顯的趨勢突變，則對此序列中動力學(xué)結(jié)構(gòu)突變的檢測就無從談起.綜合對比以上檢測實(shí)驗(yàn)結(jié)果和分析可知，F(xiàn)I法無論在對動力學(xué)結(jié)構(gòu)突變檢測的準(zhǔn)確性還是穩(wěn)定性方面較傳統(tǒng)方法均具有明顯的優(yōu)勢.

3.2 FI方法在實(shí)際氣象觀察資料中的應(yīng)用

圖5 子序列長度取20時(shí)檢測結(jié)果 (a)理想序列;(b)滑動t檢驗(yàn);(c)Yamamoto法

氣候系統(tǒng)具有典型的時(shí)空多尺度、結(jié)構(gòu)多層次、本質(zhì)非線性特征，不同層次之間關(guān)系及其相互作用十分復(fù)雜[31].自從Lorenz[32,33]和Charney等[34]從理論上論證了氣候突變存在的可能性，有關(guān)氣候突變的研究得以廣泛開展.按照控制氣候系統(tǒng)的動力學(xué)方程是否發(fā)生變化將氣候突變分為動力學(xué)結(jié)構(gòu)突變和氣候狀態(tài)變量在統(tǒng)計(jì)意義上發(fā)生了顯著變化的突變，如均值突變、趨勢突變等.傳統(tǒng)的突變檢測方法如滑動t-檢驗(yàn)，Cramer法，Yamamoto法以及M-K法等均是基于狀態(tài)變量在統(tǒng)計(jì)意義上的顯著變化來判斷突變的發(fā)生，其檢測結(jié)果均具有多尺度特征.然而系統(tǒng)的動力學(xué)結(jié)構(gòu)突變發(fā)生的時(shí)間與特定的時(shí)間尺度沒有必然的關(guān)聯(lián)，因此基于系統(tǒng)狀態(tài)變量在統(tǒng)計(jì)意義上的顯著變化來判斷突變的傳統(tǒng)檢測方法不能夠?qū)ο到y(tǒng)的動力學(xué)結(jié)構(gòu)突變進(jìn)行有效識別[5].近年來，一些新的突變檢測方法如近似熵(ApEn)算法，Q算法等用于檢測氣候系統(tǒng)的動力學(xué)結(jié)構(gòu)突變均取得了一定的效果.然而隨著人們對復(fù)雜氣候系統(tǒng)愈來愈多的了解和掌握，一些原先未曾認(rèn)識到的問題逐漸暴露了出來，人們迫切需要探尋新的突變檢測方法來應(yīng)對不斷涌現(xiàn)的新情況、新問題.因此，有關(guān)氣候系統(tǒng)的動力學(xué)結(jié)構(gòu)突變檢測問題一直處于持續(xù)的研究和探索中.

氣候系統(tǒng)中的兩大氣象要素溫度和降水是兩個(gè)最能體現(xiàn)氣候動力學(xué)結(jié)構(gòu)發(fā)生突變的量，因此，本文采用兩大氣象要素之一的溫度來檢驗(yàn)本文所提方法在檢測氣候系統(tǒng)動力學(xué)結(jié)構(gòu)突變中的有效性.文中所用數(shù)據(jù)來自中國氣象局國家氣象信息中心“西北地區(qū)地面氣候資料日值數(shù)據(jù)集”蘭州站1960—2008年日平均氣溫資料.圖6給出了該觀察站1960—2008年逐日平均氣溫的FI檢測結(jié)果.

圖6 蘭州站逐日平均溫度序列FI變化情況 (a)子序列長度取365 d(1年);(b)子序列長度取730天(2年)

從圖中可以看出，20世紀(jì)90年代初該觀察站日均溫度發(fā)生了一次明顯的突變，這與已有該站所在地區(qū)歷史上日均溫度突變年份的研究結(jié)果相一致[35,36].另一方面，從圖6中子序列長度取1年和2年時(shí)，所檢測到的突變點(diǎn)相一致來看，進(jìn)一步表明該方法所得檢測結(jié)果不依賴于子序列的長度.

4 討論與結(jié)論

本文引入常被用來檢測生態(tài)系統(tǒng)穩(wěn)定性和生態(tài)狀態(tài)翻轉(zhuǎn)的指標(biāo)——Fisher信息來識別和檢測系統(tǒng)動力學(xué)結(jié)構(gòu)突變，是因?yàn)楸疚木o扣Fisher信息是基于概率密度分布的導(dǎo)數(shù)這一基本原則.基于不同動力學(xué)性質(zhì)的數(shù)據(jù)其概率密度分布大小不同，而具有相同動力學(xué)性質(zhì)的數(shù)據(jù)其概率密度分布差異不大.因此，當(dāng)系統(tǒng)的動力學(xué)結(jié)構(gòu)發(fā)生突變時(shí)，基于突變前后兩種不同動力學(xué)結(jié)構(gòu)的數(shù)據(jù)其概率密度分布肯定會發(fā)生不同程度的改變.這種細(xì)微的改變恰好通過Fisher信息加以定格和放大.因此，理論上Fisher信息是對動力學(xué)結(jié)構(gòu)突變的一種很好的測量.

本文首先將它用于分析理想時(shí)間序列，結(jié)果表明該方法能夠準(zhǔn)確定位突變點(diǎn)，并且具有良好的抗噪性能，同時(shí)，檢測結(jié)果不依賴于子序列長度的選擇.緊接著，比較了新方法與傳統(tǒng)檢測方法的檢測效果，彰顯了新方法的優(yōu)越性.最后，應(yīng)用于“西北地區(qū)地面氣候資料日值數(shù)據(jù)集”蘭州站1960—2008年(49年)日平均氣溫資料，發(fā)現(xiàn)檢測結(jié)果與已有該站所在地區(qū)歷史上日均溫度突變年份的研究結(jié)果完全符合，進(jìn)一步展示了該方法的有效性和實(shí)用性，表明該方法在動力學(xué)結(jié)構(gòu)突變檢測中的應(yīng)用前景和實(shí)用價(jià)值.

不像在研究氣候時(shí)空分布時(shí)那樣需要做時(shí)空插值，由于受觀察資料的限制，邊界效應(yīng)是一個(gè)不容忽視的問題.而本方法是按照所給時(shí)段內(nèi)的觀察數(shù)據(jù)，采用滑動窗技術(shù)逐點(diǎn)計(jì)算FI值，并把計(jì)算結(jié)果分配給窗口的末端，這就意味著對突變點(diǎn)的真正檢測是從觀察數(shù)據(jù)的第L點(diǎn)開始的(L為窗口寬度即子序列長度)，對前L-1點(diǎn)上的突變檢測無法進(jìn)行，這是基于滑動窗技術(shù)的此類檢測方法的一個(gè)共性問題(包括滑動t檢驗(yàn)，Yamamoto法和Crammer法等基于子序列均值的傳統(tǒng)檢測方法).

本方法可用于檢測任何時(shí)間尺度的觀察數(shù)據(jù)序列，觀察數(shù)據(jù)的時(shí)間尺度僅僅影響檢測的精度即分辨率.

[1]Feng G L，Gong Z Q，Zhi R 2008Acta Metrol.Sin.66 892(in Chinese)[封國林，龔志強(qiáng)，支蓉2008氣象學(xué)報(bào)66 892]

[2]Gong Z Q，F(xiàn)eng G L，Wan S Q，Li J P 2006Acta Phys.Sin.55 477(in Chinese)[龔志強(qiáng)，封國林，萬仕全，李建平2006物理學(xué)報(bào)55 477]

[3]Feng G L，Gong Z Q，Zhi R 2006Nonlinear Theories and Methods on Spatial Temporal Distribution of the Observational Data(Beijing：China Meteorological Press)pp1—100(in Chinese)[封國林，龔志強(qiáng)，支蓉2006觀察數(shù)據(jù)時(shí)空分布的非線性理論和方法(北京：中國氣象出版社)第1—100頁]

[4]Gong Z Q，F(xiàn)eng G L，Dong W J 2006Acta Phys.Sin.55 3180(in Chinese)[龔志強(qiáng)，封國林，董文杰2006物理學(xué)報(bào)55 3180]

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

[6]Hou W，F(xiàn)eng G L，Dong W J，Li J P 2006Acta Phys.Sin.55 2663(in Chinese)[侯威，封國林，董文杰，李建平2006物理學(xué)報(bào)55 2663]

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

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

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

[10]Cheng H Y，He W P，He T，Wu Q 2012Acta Phys.Sin.61 039201(in Chinese)[成海英，何文平，何濤，吳瓊2012物理學(xué)報(bào)61 039201]

[11]Babak A S，Krishnaprasad P S 2004EURASIP J.Appl.Signal Processing15 2295

[12]Abhisek U，Rastko Z 2006Electric Power System Research76 815

[13]Gong Z Q，Zou M W，Gao X Q，Dong W J 2005Acta Phys.Sin.54 3947(in Chinese)[龔志強(qiáng)，鄒明瑋，高新全，董文杰2005物理學(xué)報(bào)54 3947]

[14]Feng G L，Hou W，Dong W J 2006Acta Phys.Sin.55 962(in Chinese)[封國林，侯威，董文杰2006物理學(xué)報(bào)55 962]

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

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

[17]Cabezas H，F(xiàn)ath B D 2002Fluid Phase Equilibria194–197 3

[18]Fath B，Cabezas H，Pawlowski C W 2003J.Theor.Biol.222–224 517

[19]Frieden B R 1998Physics from Fisher Information：A Unification(Cambridge，UK：Cambridge University Press)

[20]Frieden B R 2001Probability,Statistical Optics and Data Testing(3th Ed.)(New York：Springer-Verlag)

[21]Mayer A L，Pawlowski C W，Cabezas H 2006Ecol.Model.195 72

[22]Karunanithi A T，Cabezas H，F(xiàn)rieden B R，Pawlowski C W 2008Ecol.Soc.13 22

[23]Fath B D，Cabezas H 2004Ecol.Model.174 25

[24]Cabezas H，Pawlowski C W，Mayer A L，Hoagland N T 2005J.Clean.Product.13 455

[25]Cabezas H，Pawlowski C W，Mayer A L，Hoagland N T 2005Resources,Conservation and Recycling44 279

[26]Shastri Y，Diwekar U，Cabezas H 2008Environ.Sci.Technol.42 5322

[27]Shastri Y，Diwekar U，Cabezas H，Williamson J 2008Environ.Sci.Technol.42 6710

[28]Rico-Ramirez V R，Qunitana-Hernandez P A，Ortiz-Cruz J A，Hernandez-Castro S 2008Computer-Aided Chem.Engin.25 1155

[29]Fisher R A 1922Philosoph.Trans.Roy.Soc.London222 309

[30]Eason T，Cabezas H 2012J.Environ.Manage.94 41

[31]Li Y Q 2001Plateau Meteorol.20 88(in Chinese)[李躍清2001高原氣象20 88]

[32]Lorenz E N 1963J.Atmos.20 130

[33]Lorenz E N 1976Quat.Res.6 495

[34]Charney J G，Devore J G 1979J.Atmos.Sci.36 1205

[35]Wei F Y，Cao H X 1995Sci.Atmos.Sin.19 140(in Chinese)[魏鳳英，曹鴻興1995大氣科學(xué)19 140]

[36]Zhao F F，Xu Z X 2006Sci.Atmos.Sin.64 246(in Chinese)[趙芳芳，徐宗學(xué)2006氣象學(xué)報(bào)64 246]