小波相關排列熵和HMM 在故障預測中的應用＊

馮輔周，司愛威，江鵬程

（裝甲兵工程學院機械工程系，北京 100072）

引 言

機械系統的大部分故障都不是瞬時發生的，而是漸進性的，其故障演變往往經歷發生、發展直至惡化的過程，而且有征兆可尋。統計表明裝備機械系統的大部分故障是具有時間依存性的、可預知的趨勢性故障，采用科學有效的故障預測方法往往能夠揭示故障的發展變化，有利于避免設備惡性事故和繼發性事故的發生。

故障預測包括3個重要環節：特征提取、退化狀態識別和剩余壽命或故障概率預測。其中，特征提取是故障預測的關鍵技術，直接關系到退化狀態識別的準確性和故障預測的精度。然而，在故障發生前系統發展變化的特征信息是一種早期故障信息，其具有的弱信息特征導致信噪比較低，在故障預測中有用信息往往被復雜系統的時變、非線性運行特征、工況和負載變化、環境干擾和測試系統噪聲等非故障信息所淹沒，如何從強噪聲背景中提取微弱特征信號，是該領域研究的一個熱點和難點。

目前，對于微弱信號的檢測還沒有一種特別有效的方法。自從Donoho等提出了小波閾值去噪方法后［1］，小波分析便廣泛應用于微弱信號檢測，但在實際應用中效果并非十分理想。而且Donoho提出的硬閾值和軟閾值方法中均有一定的缺陷：在硬閾值方法中，由于閾值函數在所取的閾值處不連續，用該方法得到的小波系數對信號進行重構時會引起波動；由軟閾值方法得到的小波系數雖然連續，但這些系數與實際信號經小波分析得到的系數存在誤差，這就直接影響著重構信號與真實信號的逼近程度。有學者利用小波多分辨率分析的特點［2，3］，對微弱信號進行處理，在強背景噪聲下顯現并增強信號的目標特性，提高信噪比。但當信號非常微弱時，信號經小波分析處理后，有用信號分量的功率仍然小于噪聲信號分量功率，甚至有可能仍然相當微弱，比噪聲小幾個數量級甚至被噪聲淹沒，或者在某些特定場合下噪聲特性與白噪聲特性不符合，使得利用小波變換降噪受到限制。同時，當噪聲頻率與信號頻率接近或重合時，抑制噪聲的同時，有用信號往往也不可避免地受到損害，這極大影響了微弱信號的檢測。

針對這一問題，本文基于小波熵理論的基本思想，引入濾波效果良好的小波相關濾波法和對信號微弱變化特征敏感的排列熵算法，定義了一種新的小波熵概念——小波相關排列熵（wavelet correlation permutation entropy，WCPE），進而提出了一種基于小波相關排列熵的特征提取新方法。試驗表明，該方法能有效地表征軸承早期故障特征。

另外，傳統的故障預測方法停留在靜態觀測的基礎上，僅取某一時刻的特征信息進行比較分析，這往往忽視了故障發生前后的關聯信息，沒有揭示系統潛在狀態所發生的變化特征，因而也難于對故障的發生和發展做出準確預測。

隱馬 爾 可 夫 模 型 （Hidden Markov Model，HMM）是一種動態模式識別工具，適用于動態過程時間序列建模并具有較強的時序模式分類能力，而且特別適用于非平穩、重復再現性不佳的信號分析。本文引入HMM，結合基于小波相關排列熵的特征提取方法，探討其應用于機械系統故障預測的技術和方法。試驗結果表明該方法能有效地識別系統的退化狀態，并能準確地預測故障發生的概率。

1 小波相關排列熵特征提取方法

1．1 小波相關濾波法原理

小波相關濾波法降噪的基本思想是基于小波分解后相鄰層系數的相關性，相鄰分解層的小波系數直接相乘后，突變信號分量會得到增強，而噪聲信號分量變得更加微弱，利用閾值檢驗，從噪聲中檢出重要的信號邊緣，并移除噪聲，使得信號的信噪比大大提高［4，5］。小波相關濾波方法具有容易實現和魯棒性好等特點，可以應用于在線監測與診斷系統進行實時信息處理。

直接將信號小波變換系數值在幾個分解層次上進行相關運算（記為Cl（m，n）），可以相當準確地確定信號邊緣或其他重要特征的位置。

式中W表示信號的小波變換，m表示分解層數，N表示離散信號的點數，n表示時間，l表示直接相乘所包含的層數。若M表示總的分解層數，則應該有m＜M－l＋1。當某個局部存在信號邊緣或其他重要特征時，這種方法可以將背景噪聲除去。直接將信號數據進行尺度空間相關，將銳化和增強信號的主要邊緣，同時噪聲和次要邊緣將被壓縮。小波相關濾波法的具體計算過程可參考文獻［6］。

1．2 排列熵理論

排列熵（Permutation Entropy，PE）是近年來新興的一門研究自然界中的不規則性以及非線性系統的算法，它可以將以前不能定量描述或是很難定量描述的復雜系統用一種較為簡便的方法描述出來。排列熵反映了一維時間序列復雜度，對信號變化具有較高的敏感性，可以很好地放大系統的微變信號，同時檢測出復雜系統的動力學突變［7］。由于系統設備在即將發生異常或故障時，其非線性因素的影響是不同的，所表現的信號復雜性也會不一樣。因此，利用排列熵進行早期故障診斷是可行的。

排列熵的大小表示時間序列｛X（i），i＝1，2，…，n｝的隨機程度：排列熵的值越小，說明時間序列越規則，反之，則時間序列越接近隨機。排列熵的變化反映并放大了時間序列的微小細節變化。排列熵的具體計算方法可參考文獻［8］。

1．3 小波相關排列熵原理

當機械系統發生早期故障時，故障信息微弱，信噪比較低，直接進行小波分解無法提取出微弱的故障信息。小波相關濾波法不僅可以對目標信號有效降噪，而且不會有分辨率的損失。因此，采用小波相關濾波法可以提取信號微弱的故障特征信息，得到信噪比較高的各層重構信號。基于上述小波熵理論和對突變信號變化敏感的排列熵算法，對小波相關濾波處理后的各頻段重構信號進行排列熵分析，即可得到各頻段信號的排列熵復雜度，由此來識別機械系統的運行狀態。

信號x（n）經過小波相關濾波降噪處理后，得到各分解層次下較高信噪比的高頻系數：Dj｛dj（k），k＝1，2，…，N，j＝1，2，…M｝和低頻系數C，它們分別包含了信號從高頻到低頻不同頻帶的信息，體現了不同分解層次下局部特征的直觀估計。為了定量描述各頻段小波系數所包含的故障特征信息，還需從得到的各層小波系數中提取反映設備運行狀態的特征參數，有利于系統運行狀態的分析和識別。在此可以把小波相關濾波降噪處理后的得到的各層小波系數Dj（j＝1，2，…，m，m＋1）看成對信號x（n）的一種劃分，定義這種劃分的測度為

式中dF（j）（k）為dj（k）的傅里葉變換。

借鑒小波熵的基本思想，定義小波相關排列熵為：將pj（k）進行相空間重構，按上述排列熵算法，依次求取經小波相關濾波后各層的排列熵值。

1．4 小波相關排列熵特征提取步驟

根據上述小波排列熵的計算方法，結合振動信號特征提取的實際需求，可總結基于小波相關排列熵的特征提取步驟如圖1所示。

圖1 小波相關排列熵特征提取步驟Fig．1 Feature extraction approach of WCPE

首先是對信號進行m層離散小波分解，得到各層的小波系數Dj＝｛dj（k），j＝1，2，…，m，m＋1｝，并對各層小波系數分層進行重構；然后，利用小波相關濾波法對各層小波重構后的信號進行逐層降噪處理，得到信噪比較高的各層分量；之后，在依據式（2）對降噪后的各層分量進行歸一化處理的基礎上，選取合適的嵌入維數和延遲時間，對其進行分層相空間重構，得到各層分量的相空間矩陣；最后，依據排列熵算法計算出各層分量的小波相關排列熵值Wp為

2 HMM的基本理論和算法改進

HMM是一個雙重隨機過程，不僅狀態之間的轉移是隨機的，而且每個狀態的觀測符號也是隨機的。模型中真實的狀態是“隱藏的”，不能直接觀測，只能通過觀測矢量感知它的存在。機械系統故障演化過程通常經歷正常、性能下降、故障、失效等幾個健康狀態衰退階段，這些健康狀態衰退過程在實際中一般不能直接觀察到，能觀測到的只是系統表現出來的性能或狀態征兆，即系統故障衰退健康狀態要通過其表現來感知，這和HMM在本質上是相通的。而且HMM可以較好地描述隱藏狀態與觀測狀態間的關聯關系。HMM可表示為

式中λ為初始概率分布矢量，A為狀態轉移概率矩陣，B為觀察值概率矩陣。圖2可以清晰地表示出一個左右型HMM的基本結構。

圖2 左右型HMM拓撲結構圖Fig．2 Framework chart of right－and－left HMM

圖2中第1行圓圈表示觀測不到的隱狀態；第2行圓圈表示可觀測到的觀測變量。其中π，A和B為 HMM 的3個基本參數，｛o1，o2，…，oT｝為觀測序列，它由狀態序列｛q1，q2，…，qT｝按觀測矩陣B產生；｛q1，q2，…，qT｝被稱為狀態路徑，其各狀態間的依賴關系是由轉移矩陣A決定。模型中的初始向量π決定初始時刻系統所處的狀態。

針對HMM的概率計算、最優狀態序列選擇和模型訓練（參數估計）等3個問題，通常對應有前向后向算法、Viterbi算法和Baum－Welch算法來解決［9］。而參數估計在HMM的應用中占有舉足輕重的地位。傳統的Baum－Welch算法（BW）雖然可以確保似然值單調上升，且該算法具有收斂性，但是該算法存在收斂速度較慢，數值運算困難等缺點，而且容易陷入局部最優。因此尋求更好的參數估計方法具有十分重要的意義。本文從優化角度出發，將遺傳算法引入HMM的參數估計中。

HMM的參數估計就是根據極大似然準則尋找合適的模型參數λ，使得輸出概率P（O｜λ）取得最大值，該過程中還存在約束條件

因而，最大化P（O｜λ）其實就是一個帶有約束的優化問題，基于此，本文采用基于懲罰策略的遺傳算法來處理這個問題。首先，定義一個模型參數λ違反約束程度的量

式中M1和M2為極大正數。

然后，可以確定適應度函數為然后即可運用遺傳算法進行HMM參數估計。

3 故障預測模型的建立

基于以上的分析，本文建立了基于HMM的預測模型框架，如圖3所示。

整個預測流程為：首先提取歷史觀測數據的特征信息輸入HMM模型，并用HMM的參數估計方法分別訓練各個退化狀態識別模型以及全壽命模型，然后就可以提取當前觀測數據的特征信息并輸入到訓練好的各個退化狀態識別模型，根據最大似然概率值即可判斷當前觀測數據代表的退化狀態，同時結合全壽命狀態轉移概率矩陣即可計算當前退化狀態下故障發生的概率。

圖3 基于WCPE和HMM的預測模型框架Fig．3 Framework of prognostic model based on WCPE and HMM

4 實例分析

4．1 數據來源及分析

為驗證該方法的有效性，選取來自美國NSFI／UCR智能維護系統中心（IMS）的軸承振動數據進行分析［10］。試驗時，在一個軸上安裝了4個軸承，軸的轉速固定在 2 000r／min，6 000 1b（約等于26 671N）徑向載荷用彈性系統加載在軸和軸承上，每個軸承上安裝2個美國PCB公司生產的353B33型高靈敏度ICP加速度傳感器（一個水平方向，一個垂直方向）。

本文的數據采用了該試驗的第2次的測試數據，其采樣頻率為20kHz，每個文件包含20 480個數據點，每10min記錄一個文件，該次試驗共采集了984個文件，即試驗持續了164h）。試驗結束之后，發現軸承1發生外圈磨損失效，這說明該試驗記錄了軸承1的全壽命數據。圖4示出了軸承1垂直方向傳感器采集的各文件數據的均方根值Xrms的趨勢圖。

圖4 軸承1垂直方向振動信號均方根值趨勢圖Fig．4 RMStrend chart of bearing 1vibration signal on vertical direction along the time

從圖4可以看出，在118．0h之前，軸承1運行比較平穩；當軸承1運行至118．0h時，其振動信號出現微小的跳變，說明軸承已經開始出現異常；118．0～160．0h，振動信號在上下波動，但波動幅度不大，說明軸承在帶故障運行，但故障不太嚴重；當軸承1運行超過160．0h，其振動信號出現劇烈變化，并在163．3h其均值達到最大值，此時軸承已出現嚴重故障，達到其壽命極限。

另外，從圖4可以看出，基于振動信號特征值的趨勢分析雖然能夠大致反映軸承運行的全壽命過程，但其發生早期輕微故障時（118．0～160．0h），信號反映的特征非常微弱；而且如果在實際應用的復雜機械系統中，特征信號更容易淹沒在各種噪聲的干擾之中，這就導致在故障發生的早期，檢測和發現故障非常困難；而當特征信號表現明顯（超過160．0 h）被發現的時候，故障已相當嚴重，而且此時留給維護人員做出維修決策的時間又非常短暫，往往會導致設備被迫停車或嚴重損壞的后果，損失將不可估量。

4．2 特征提取

為驗證小波相關排列熵特征提取算法在早期故障診斷中的有效性，用該算法對上述軸承原始振動數據進行分析。對每個文件的數據按照基于小波相關排列熵的特征提取步驟依次求取其小波相關排列熵值，其中小波分解層數為4層，由于Daubechies系列小波正交、緊支，對不規則信號較為敏感，同時依據文獻［11］的方法選用db4小波基函數進行小波分解。排列熵計算時分別利用互信息法和偽近鄰法選取嵌入維數為d＝5［12］，延遲時間為τ＝4。同時，為減少計算量，計算每個文件數據各層的排列熵時，將總長為20 480的數據分為20段，每段1 024點分別求其排列熵值，然后求這20個排列熵值的均值作為該文件數據該層的排列熵。圖5示出了軸承1振動信號第一層小波相關排列熵的變化趨勢圖。

圖5 軸承1振動信號第一層小波相關排列熵的趨勢圖Fig．5 The first layer WCPE trend chart of bearing 1vibration signal along the time

從圖5中可以看出：當軸承運行至86．8h，第一層小波相關排列熵值開始跳變，而且變化幅度較大，趨勢明顯，即振動信號的小波相關排列熵可以檢測出軸承的早期故障，而且其檢出時間比用均方根值特征檢測提前了31．2h。同時也表明了小波相關排列熵算法對微弱信號特征提取的有效性。

為了進行比較，這里直接選取軸承1垂直方向振動信號進行排列熵分析，其中數據分段方法與前文相同。軸承1時域振動信號的排列熵分析結果如圖6所示。

圖6 軸承1振動信號的排列熵趨勢圖Fig．6 PE trend chart of bearing 1vibration signal along the time

從圖6中可以看出，時域信號的排列熵特征值也能反映滾動軸承運行的全壽命過程，但對比圖5和6，可以明顯看出，經這兩種算法提取的特征向量檢測出軸承狀態開始變化的時間是不同的，即小波相關排列熵算法提取的特征向量較早地檢測到了軸承的早期故障，而且比排列熵算法提前了6．7h。值得說明的是，該軸承振動數據是在大載荷、高轉速的加速壽命試驗的條件下采集的，如果是在實際正常條件下運行，提前的時間將會更長。

4．3 退化狀態模型訓練及識別

為實現軸承的退化狀態識別并最終實現故障預測，就要分別建立各退化狀態的HMM以及整個全壽命數據的HMM，根據文獻［13］算法的分析，該軸承振動數據的最佳退化狀態數為6（其中包含1個正常狀態、4個退化狀態和1個故障狀態）。因此，確定訓練退化狀態模型數為6個，分別為HMM1（λ1），HMM2（λ2），…，HMM6（λ6）。

選用上一節提取的軸承振動數據的小波相關排列熵作為特征信息點進行退化狀態的識別和故障預測。從圖4可以看出，當軸承運行超過120h（720個文件記錄）時，已處于嚴重故障狀態，因此這里僅選取前120h的特征信息作為模型訓練和檢測的依據。將720個小波相關排列熵特征點劃分為120個樣本（每小時一個樣本），每個樣本6個特征點，即觀測值數為6。其中，前1～80樣本為正常狀態的80個樣本，81～88為退化狀態1的8個樣本，89～96為退化狀態2的8個樣本，97～104為退化狀態3的8個樣本，105～112為退化狀態4的8個樣本，113～120為故障狀態的8個樣本。同時，選取每個狀態前一半樣本作為訓練樣本，后一半作為測試樣本。

HMM的一個重要參數就是馬爾可夫鏈的狀態數，它描述了狀態的轉移，直接影響模型的性能。本文選取ADC準則來確定各退化狀態HMM的馬爾可夫鏈的狀態數。通過計算，最佳狀態數為5（詳細可參考文獻［14］）。因此采用5狀態HMM初始化模型。因初始狀態轉移矩陣和初始觀測值概率矩陣對離散HMM訓練影響較小，這里采用等概率方式產生。分別采用BW算法和GA算法對各退化狀態模型進行訓練，同時設定訓練過程中最大迭代步數和進化代數為100。用遺傳算法訓練時，種群大小設為40，最大進化代數設為100，交叉概率設為0．8，變異概率設為0．02。圖7示出了分別用BW算法和GA算法訓練正常狀態模型時，各自的相對誤差隨迭代次數（進化代數）變化的趨勢圖。

圖7 BW算法和GA算法訓練誤差變化趨勢圖Fig．7 Training error trend chart of BW and GA

從圖7可以看出，使用GA算法訓練時的誤差曲線明顯比使用GA算法訓練時的誤差曲線平滑，而且BW算法的訓練誤差曲線存在明顯的“平坦區”，這說明BW算法更容易陷入局部最優。同時使用GA算法訓練模型時，進化代數達到29時即可收斂，而BW算法需要43步，即GA算法具有更好的收斂速度。另外，使用GA算法收斂時的相對誤差要遠低于使用BW算法收斂時的相對誤差，表明GA算法具有更好的估計精度。

在各退化狀態HMM模型訓練完畢之后，選取測試樣本進行退化狀態的識別。識別時，同樣將每個測試樣本分別輸入到6個訓練好的退化狀態模型CHMM1（λ1），CHMM2（λ2），…，CHMM6（λ6）中求其各自的輸出概率，輸出概率最大的模型所對應的退化狀態即代表該測試樣本的識別結果。識別結果表明：正常狀態、退化狀態1～4和故障狀態的測試樣本識別正確率分別為90%，75%，100%，100%，100%和100%。因此，可計算6個退化狀態模型的總體識別率為93．3%，識別效果比較理想。

另外，本文選取了BP神經網絡和SVM模型，利用和HMM相同的訓練樣本和測試樣本進行訓練和識別，其中BP神經網絡的輸入層節點數為6，輸入層節點數為3，隱含層節點數為13，最大迭代次數為1 000，誤差精度為0．001；SVM 模型采用徑向基核函數，并選用“1對1”的多類分類算法。結果顯示，BP神經網絡的識別率正確率為78．3%，而SVM模型的識別率正確率為86．7%。因此，HMM的識別效果明顯優于BP神經網絡和SVM模型。

4．4 全壽命模型訓練及故障預測

因軸承的整個全壽命過程可分為6個運行狀態，在此全壽命預測模型將采用6狀態HMM進行訓練。為保證訓練數據的覆蓋性和訓練效果，按照各退化狀態數據交叉組合的原則，通過前文用到的720個小波相關排列熵特征點構建48個訓練樣本，每個樣本包含6個數據點，分別對應6種退化狀態。其中，正常狀態的480個特征點中每10點選取第1點作為一個訓練樣本的第一個觀測值，退化狀態1～4和故障狀態分別包含的48個特征點依次選取作為訓練樣本的第二、第三、第四、第五和第六個觀測值。

將48個全壽命樣本輸入到全壽命模型中，采用GA算法來訓練模型，參數設置與退化狀態模型訓練時相同，當進化代數達到37時，模型收斂。此時得到訓練后的全壽命模型的狀態轉移矩陣A＊為

在退化狀態識別和全壽命模型訓練完畢的基礎上，即可得出所處狀態的發生故障的概率。設軸承處于正常狀態和退化狀態1～4時發生故障的概率分別為P1，P2，P3，P4和P5，則可計算如下

分析式（9）可知，軸承在正常狀態發生故障概率較小，而一旦進入退化狀態，其故障概率就隨之不斷增大，這也與機械系統的退化規律相吻合。

在獲取處于各退化狀態時發生故障概率的基礎上，結合新樣本退化狀態的識別結果，即可得到新樣本所處狀態發生故障的概率。

5 結 論

（1）軸承振動信號小波相關排列熵特征能有效地表征軸承的運行狀態，并較好地檢測出軸承故障早期的微弱變化；

（2）基于小波相關排列熵和HMM的故障預測方法能準確地識別機械系統的退化狀態，并有效地預測故障發生的概率。

［1］ Donoho D L，Jonnston I M．Adapting to unknown smoothness via wavelet shrinkage［J］．J．Am．Stat．Assoc．，1995，90：1 200—1 224．

［2］ 王仲生，何紅，陳錢．小波分析在發動機早期故障識別中的應用研究［J］．西北工業大學學報，2006，24（1）：68—71．

Wang Zhongsheng，He Hong，Chen Qian．Exploring effective early identification of aero－engine rotor faults［J］．Journal of Northwestern Polytechnical University，2006，24（1）：68—71．

［3］ Duan Chendong，He Zhengjia，Jiang Hongkai．New method for weak fault feature extraction based on second generation wavelet transform and its application［J］．Chinese Journal of Mechanical Engineering，2004，17（4）：543—547．

［4］ 楊江天，周培鈺．經驗模態分解和Laplace小波在機車柴油機齒輪系故障診斷中的應用［J］．機械工程學報，2011，47（7）：109—115．

Yang Jiangtian，Zhou Peiyu．Fault diagnosis for gear train of locomotive diesel engine based on empirical mode decomposition and Laplace wavelet［J］．Chinese Journal of Mechanical Engineering，2011，47（7）：109—115．

［5］ 曾慶虎，邱靜，劉冠軍．基于小波相關特征尺度熵的HSMM設備退化狀態識別與故障預測方法［J］．儀器儀表學報，2008，29（12）：2 559—2 564．

ZENG Qinghu，QIU Jing，LIU Guanjun．Equipment degradation state recognition and fault prognosis method based on wavelet correlation feature scale entropy and HSMM［J］．Chinese Journal of Scientific Instrument，2008，29（12）：2 559—2 564．

［6］ 馮輔周，司愛威，饒國強，等．基于小波相關排列熵的軸承早期故障診斷技術［J］．機械工程學報，2012，48（13）：73—79．

Feng Fuzhou，Si Aiwei，Rao Guoqiang，et al．Early fault diagnosis technology for bearing based on wavelet correlation permutation entropy［J］．Chinese Journal of Mechanical Engineering，2012，48（13）：73—79．

［7］ 劉永斌，龍潛，馮志華，等．一種非平穩、非線性振動信號檢測方法的研究［J］．振動與沖擊，2007，26（12）：131—134．

Liu Yongbin，Long Qian，Feng Zhihua，et al．Detection method for nonlinear and non－stationary signals［J］．Journal of Vibration and Shock，2007，26（12）：131—134．

［8］ 馮輔周，饒國強，司愛威，等．排列熵算法研究及其在振動信號突變檢測中的應用［J］．振動工程學報，2012，25（2）：221—224．

Feng Fuzhou，Rao Guoqiang，Si Aiwei，et al．Research and application of the arithmetic of PE in testing the sudden change of vibration signal［J］．Journal of Vibration Engineering，2012，25（2）：221—224．

［9］ Bunks C，Mccarthy D．Condition－based maintenance of machines using hidden Markov models［J］．Mechanical Systems and Signal Processing，2000，14：597—612．

［10］Qiu Hai，Lee J，Lin Jing，et al．Wavelet filter－based weak signature detection method and its application on roller bearing prognostics［J］．Journal of Sound and Vibration，2006，289：1 066—1 090．

［11］Yang H T，Liao C C．A de－noising scheme for enhancing wavelet－based power quality monitoring systems［J］．IEEE Transactions on Power Delivery，2001，16（3）：353—359．

［12］謝忠玉，張立．相空間重構參數選擇方法的研究［J］．中國科技信息，2009，（16）：42—43．

Xie Zhong－yu，Zhang Li．Selection of embedding parameters in phase space reconstruction［J］．China Science and Technology Information，2009，（16）：42—43．

［13］呂克洪．基于時間應力分析的BIT降虛警與故障預測技術研究［D］．長沙：國防科學技術大學，2008．

Lv Kehong．Research on BIT false alarm reducing and fault prediction technologies based on time stress analysis［D］．Changsha：National University of Defense Technology，2008．

［14］舒曉娟，陳洋波，任啟偉．模型選擇準則在洪水頻率分析中的應用［J］．水利學報，2010，41（1）：80—85．

Shu Xiaojuan，Chen Yangbo，Ren Qiwei．Model selection criteria for flood frequency analysis［J］．Journal of Hydraulic Engineering，2010，41（1）：80—85．