基于Fisher信息理論的突變檢測新方法*

蔡舒平 戴理 倪淏

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

2)(江蘇大學，機械工業(yè)設施農(nóng)業(yè)測控技術與裝備重點實驗室，鎮(zhèn)江 212013)

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

1 引言

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

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

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

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

2 方法

2.1 Fisher信息

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

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

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

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

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

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

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

2.2 Fisher信息的計算

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

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

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

在把時間序列轉變?yōu)橐粋€重疊窗口序列以后，接下來是將每個滑動窗口劃分成如下I個區(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ù)個數(shù)w之比值.

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

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

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

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

3 結果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

進一步考察當子序列的長度取不同值時，滑動t檢驗和Yamamoto法的檢測結果是否仍然保持一致.圖5給出了當子序列長度取20，滑動步長取1時，滑動t檢驗和Yamamoto法對此的檢測結果.

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

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

3.2 FI方法在實際氣象觀察資料中的應用

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

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

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

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

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

4 討論與結論

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

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

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

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

[1]Feng G L，Gong Z Q，Zhi R 2008Acta Metrol.Sin.66 892(in Chinese)[封國林，龔志強，支蓉2008氣象學報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)[龔志強，封國林，萬仕全，李建平2006物理學報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)[封國林，龔志強，支蓉2006觀察數(shù)據(jù)時空分布的非線性理論和方法(北京：中國氣象出版社)第1—100頁]

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

[5]He W P，Wu Q，Zhang W，Wang Q G，Zhang Y 2009Acta Phys.Sin.58 2862(in Chinese)[何文平，吳瓊，張文，王啟光，張勇2009物理學報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物理學報55 2663]

[7]Wang Q G，Zhang Z P 2008Acta Phys.Sin.57 1976(in Chinese)[王啟光，張增平2008物理學報57 1976]

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

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

[10]Cheng H Y，He W P，He T，Wu Q 2012Acta Phys.Sin.61 039201(in Chinese)[成海英，何文平，何濤，吳瓊2012物理學報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)[龔志強，鄒明瑋，高新全，董文杰2005物理學報54 3947]

[14]Feng G L，Hou W，Dong W J 2006Acta Phys.Sin.55 962(in Chinese)[封國林，侯威，董文杰2006物理學報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物理學報54 5487]

[16]Feng G L，Wang Q G，Hou W，Gong Z Q，Zhi R 2009Acta Phys.Sin.58 2853(in Chinese)[封國林，王啟光，侯威，龔志強，支蓉2009物理學報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大氣科學19 140]

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