基于卡爾曼濾波和主元分析的發(fā)動(dòng)機(jī)性能衰退

李 冬，李本威，趙鵬飛，趙 凱，朱飛翔

（1.海軍航空工程學(xué)院a.研究生管理大隊(duì)；b.飛行器工程系，山東煙臺(tái)264001；2.92407部隊(duì)，遼寧綏中125200）

研究部件性能參數(shù)衰退是開展發(fā)動(dòng)機(jī)整機(jī)性能研究的基礎(chǔ)。部件性能參數(shù)包括流量、效率等。在不采取任何維修措施的前提下，發(fā)動(dòng)機(jī)部件性能參數(shù)隨使用時(shí)間延長會(huì)逐漸衰退，使發(fā)動(dòng)機(jī)整機(jī)性能下降，如不及時(shí)采取必要的維修措施，衰退程度會(huì)進(jìn)一步加重，產(chǎn)生更為嚴(yán)重的影響[1-3]。只有準(zhǔn)確地估計(jì)部件性能參數(shù)，才能準(zhǔn)確地了解發(fā)動(dòng)機(jī)狀況，才有可能制定更為科學(xué)的維修計(jì)劃，延長發(fā)動(dòng)機(jī)使用壽命，降低維修成本。近年來，不斷出現(xiàn)研究發(fā)動(dòng)機(jī)性能變化的監(jiān)控算法，例如加權(quán)最小二乘法[4]、專家系統(tǒng)[5]、卡爾曼濾波[6]和遺傳算法[7]。

本文首先通過非線性參數(shù)估計(jì)，在給定測量參數(shù)基礎(chǔ)上得到對應(yīng)離散時(shí)刻的部件性能參數(shù)，并對性能參數(shù)增加參數(shù)變化約束范圍加快收斂。當(dāng)測量參數(shù)變化比較集中時(shí)，回代到穩(wěn)態(tài)模型利用非線性技術(shù)求解時(shí)，由于穩(wěn)態(tài)模型精度和求解效率的原因，共同方程組不一定收斂到最佳，因而部件性能參數(shù)衰退趨勢不明顯。利用卡爾曼濾波器估計(jì)參數(shù)變化，將性能參數(shù)擴(kuò)充為狀態(tài)參數(shù)，在某些離散時(shí)間點(diǎn)設(shè)置卡爾曼濾波器的狀態(tài)參數(shù)，得到部件性能參數(shù)完整衰退過程。在此性能參數(shù)基礎(chǔ)上，將其代入到穩(wěn)態(tài)模型中，得到對應(yīng)的測量參數(shù)。利用主元分析方法對測量參數(shù)在2個(gè)空間（主成分空間PCA和殘差空間RS）進(jìn)行分解，提取反映部件性能參數(shù)衰退的成分。在此基礎(chǔ)上構(gòu)造反映發(fā)動(dòng)機(jī)整機(jī)性能衰退的統(tǒng)計(jì)量，進(jìn)一步分析發(fā)動(dòng)機(jī)性能變化。為及時(shí)做好維護(hù)工作、延長發(fā)動(dòng)機(jī)使用壽命提供理論指導(dǎo)。

1 性能衰退參數(shù)特性表示及優(yōu)化選取

對于發(fā)動(dòng)機(jī)性能衰退的模擬，將發(fā)動(dòng)機(jī)部件性能參數(shù)（流量，效率）乘以相應(yīng)的衰退系數(shù)（一般介于0.95～1）。在穩(wěn)態(tài)模型的基礎(chǔ)上，得到對應(yīng)的發(fā)動(dòng)機(jī)測量參數(shù)。圖1為壓氣機(jī)性能衰退前后的特性圖。

圖1 壓氣機(jī)性能衰退前后的特性曲線

由圖1可知，當(dāng)壓氣機(jī)性能出現(xiàn)衰退時(shí)，其壓比和效率特性曲線在特性圖上表現(xiàn)為退化狀態(tài)曲線較健康部件曲線向下偏移。

研究發(fā)動(dòng)機(jī)性能參數(shù)的衰退情況，沒有必要分析所有性能參數(shù)。這樣一來，會(huì)帶來計(jì)算負(fù)擔(dān)，而且有些參數(shù)變化對發(fā)動(dòng)機(jī)整體性能變化并不起主要作用，所以要對表征部件性能衰退的參數(shù)進(jìn)行優(yōu)選。

首先，定性來看，由于葉片在使用過程中更容易發(fā)生形變，與風(fēng)扇、壓氣機(jī)和渦輪部件相比，燃燒室等部件性能參數(shù)在使用過程中變化緩慢，因而表征發(fā)動(dòng)機(jī)性能的參數(shù)初步選擇風(fēng)扇、壓氣機(jī)、高壓渦輪和低壓渦輪的流量和效率；其次，采用擾動(dòng)法，通過部件性能參數(shù)與測量參數(shù)相關(guān)性分析[8]；最終，確定的部件性能參數(shù)為：風(fēng)扇流量、壓氣機(jī)流量、高壓渦輪效率、低壓渦輪效率。定義以下4個(gè)部件性能參數(shù)變化量表征發(fā)動(dòng)機(jī)性能的變化。

風(fēng)扇流量退化參數(shù)：

壓氣機(jī)流量退化參數(shù)：

高壓渦輪效率退化參數(shù)：

低壓渦輪效率退化參數(shù)：

式（1）～（4）中：(?)af、(?)ac、(?)th、(?)tl、(?)act、(?)ref分別表示風(fēng)扇、壓氣機(jī)、高壓渦輪、低壓渦輪、實(shí)際值以及參考值；M、η、T、P分別表示流量、效率、總溫和總壓；Δ為相對變化量。

2 發(fā)動(dòng)機(jī)性能參數(shù)非線性估計(jì)

航空發(fā)動(dòng)機(jī)部件性能參數(shù)變化估計(jì)，從本質(zhì)上講是利用實(shí)際測量值與發(fā)動(dòng)機(jī)模型計(jì)算值之間的偏差來調(diào)整部件性能參數(shù)，從而使模型計(jì)算值與實(shí)際值相一致[9]。由發(fā)動(dòng)機(jī)穩(wěn)態(tài)模型得到測量參數(shù)：

式（5）中：u為模型控制量，如供油量、尾噴口面積；Δp=[ΔMaf,ΔMac,Δηth,Δηtl]，ΔMaf、ΔMac、Δηth、Δηtl分別表示風(fēng)扇流量、壓氣機(jī)流量、高壓渦輪效率和低壓渦輪效率相對變化量；y=[T3,P3,T5,P5]T為輸出的測量參數(shù)，T3、P3、T5、P5分別表示壓氣機(jī)出口總溫、總壓、低壓渦輪出口總溫、總壓。

設(shè)發(fā)動(dòng)機(jī)實(shí)際測量值為yreal，則對Δp的估計(jì)問題轉(zhuǎn)化為求解非線性方程的問題：

由于本文設(shè)定的部件性能參數(shù)個(gè)數(shù)與實(shí)際測量參數(shù)個(gè)數(shù)相等，式（6）有唯一解。本文采用在共同工作方程基礎(chǔ)上增加4個(gè)平衡方程的方法求解部件性能參數(shù)，迭代計(jì)算采用粒子群算法。由于部件性能參數(shù)不斷惡化，因而在粒子群算法中進(jìn)一步增加參數(shù)變化范圍作為約束條件加快收斂。第k時(shí)刻的性能退化量Δp(k)的約束表達(dá)式描述如下：

式（7）中：Δp為性能變化量；k為參數(shù)估計(jì)的時(shí)刻；r+和r-均為非負(fù)因子，用于保證狀態(tài)參數(shù)估計(jì)在一定范圍內(nèi)變化。

3 基于卡爾曼濾波器的性能衰退模型

當(dāng)發(fā)動(dòng)機(jī)測量參數(shù)比較集中，回代到穩(wěn)態(tài)模型利用非線性技術(shù)求解時(shí)，由于模型精度的原因，得到的性能參數(shù)可能不能如實(shí)反映發(fā)動(dòng)機(jī)性能衰退規(guī)律?？紤]到卡爾曼濾波器的參數(shù)估計(jì)能力，通過含有噪聲的發(fā)動(dòng)機(jī)，可測輸出偏差量Δy估計(jì)部件性能退化量Δp，將性能退化量并入到狀態(tài)空間參數(shù)中，在某些時(shí)刻設(shè)置卡爾曼濾波器的狀態(tài)參數(shù)。建立描述性能退化參數(shù)的發(fā)動(dòng)機(jī)線性化狀態(tài)空間模型[10]：

式（8）中：Δx˙表示狀態(tài)變量的變化率；Δx=[Δnl,Δnh]T為狀態(tài)變量，Δnl、Δnh分別表示低、高壓轉(zhuǎn)子轉(zhuǎn)速相對變化量；Δy=[Δnl,Δnh,ΔT3,ΔP3,ΔT5,ΔP5]T；Δu=[ΔMf,ΔA8]T為控制量，ΔMf、ΔA8分別表示供油量、尾噴口面積相對變化量；w1、w2分別為系統(tǒng)噪聲與測量噪聲，它們相互獨(dú)立；系數(shù)矩陣A、B、C、D、L、M采用抽功法[8]求取。

由于Δp是被估計(jì)向量，將Δp與Δx合并為增廣的狀態(tài)變量，考慮系統(tǒng)噪聲w1和測量噪聲w2，增廣狀態(tài)空間模型（ASVM）表示如下：

式（9）中：K為卡爾曼濾波增益矩陣；Δx?、Δy?、Δp?為估計(jì)值。

結(jié)合第2 節(jié)的部件性能參數(shù)估計(jì)方法，不同時(shí)刻的性能參數(shù)作為卡爾曼濾波器中對應(yīng)的狀態(tài)參數(shù)。由于得到的測量參數(shù)都是間斷的，可用插值的方法擴(kuò)充數(shù)據(jù)。具體算法如圖2所示。

1）給定發(fā)動(dòng)機(jī)某一狀態(tài)控制量，以及已知某些離散時(shí)刻對應(yīng)的測量參數(shù)；

2）通過非線性參數(shù)估計(jì)方法，估計(jì)已知離散時(shí)刻測量參數(shù)對應(yīng)的性能參數(shù)，計(jì)算要經(jīng)過充分的迭代運(yùn)算保證較高的精度；

3）采用插值方法增加測量參數(shù)數(shù)據(jù)，將2）得到的性能參數(shù)擴(kuò)展為濾波器中的狀態(tài)參數(shù)；

4）結(jié)合擴(kuò)展卡爾曼濾波器，估計(jì)性能參數(shù)完整的衰退過程。

圖2 性能參數(shù)估計(jì)

一般來講，在沒有得到維護(hù)的情況下，發(fā)動(dòng)機(jī)性能是不斷惡化的。在卡爾曼濾波參數(shù)估計(jì)中，由于在間斷時(shí)刻設(shè)置狀態(tài)參數(shù)，會(huì)在其附近產(chǎn)生性能參數(shù)突變情況?？紤]性能變化的漸變性和突變2 種情況[11]；性能漸變性變化，即性能參數(shù)變化Δp˙緩慢，在考查時(shí)間段內(nèi)，滿足Δp˙=0；突變主要表現(xiàn)為性能在較短時(shí)間內(nèi)迅速變化，本質(zhì)為部件性能參數(shù)迅速變化。突變發(fā)生時(shí)Δp˙=0不再成立，但突變時(shí)間較短，并且突變后再次穩(wěn)定工作后參數(shù)變化平穩(wěn)，滿足Δp˙=0。

4 基于PCA方法的性能監(jiān)控模型

主元分析（PCA）作為一種受歡迎的方法在工程領(lǐng)域備受關(guān)注，特別是化學(xué)工業(yè)和結(jié)構(gòu)振動(dòng)，但在發(fā)動(dòng)機(jī)性能監(jiān)控和故障診斷方面得到的關(guān)注卻很少[12]。如果某些變量或者過程是相關(guān)的，則會(huì)存在一定程度的冗余量，這樣就可以在保證信息量不變的情況下化簡消元。主元分析方法的基本思想是在最大程度保持原始數(shù)據(jù)信息的基礎(chǔ)上，盡量減少數(shù)據(jù)的維數(shù)，達(dá)到化簡的目的。圖3用一個(gè)簡單的2-D的例子說明主元分析的含義[13]。

圖3 所示的是在給定狀態(tài)下，包含z1和z22 類數(shù)據(jù)。例如發(fā)動(dòng)機(jī)的燃油流量和轉(zhuǎn)速，2 個(gè)量在一定程度上是相關(guān)的，即當(dāng)z1增加時(shí)，z2也增加。當(dāng)涉及到多個(gè)變量的時(shí)候，問題變得復(fù)雜了。還以2個(gè)變量為例，要使待分析數(shù)據(jù)在坐標(biāo)系內(nèi)得以準(zhǔn)確表達(dá)，最佳選擇就是建立直角坐標(biāo)系，其中大多數(shù)據(jù)集中在一個(gè)坐標(biāo)附近（記作PCS，主成分空間Pˉ），包含主要的數(shù)據(jù)信息；另一個(gè)坐標(biāo)（記作RS，殘差空間P?）垂直于它，攜帶少量殘余數(shù)據(jù)信息，表示噪聲以及性能退化帶來的誤差。噪聲可以通過采取濾波等方法去除，分析時(shí)只考慮性能退化的影響。因此，可以由主元分析方法完成對發(fā)動(dòng)機(jī)參數(shù)中反應(yīng)性能衰退成分的提取。

圖3 主元分析示意圖

引入如下計(jì)算步驟：

1）利用圖2的卡爾曼濾波器得到完整的部件性能參數(shù)變化規(guī)律。在此基礎(chǔ)上代入到穩(wěn)態(tài)模型得到對應(yīng)的測量參數(shù)；

2）引入R=PΛPT；其中，R是上述1）中測量參數(shù)協(xié)方差矩陣，P為特征向量，Λ為包含特征值的對角矩陣。

3）給定測量參數(shù)y，經(jīng)過標(biāo)準(zhǔn)化處理作如下變換，令其在時(shí)刻k時(shí)，PCS空間上的投影yˉk和RS空間上的投影為y?k，則：

4）發(fā)動(dòng)機(jī)性能出現(xiàn)衰退后，各部件將在新的平衡點(diǎn)重新達(dá)到穩(wěn)定關(guān)系。換成PCA語言來解釋，意味著各個(gè)變量之間的關(guān)系重新修正，代表性能參數(shù)的主要信息yˉk發(fā)生改變，但是殘差部分y?k相對無衰退情況明顯增大。因此，構(gòu)造反映發(fā)動(dòng)機(jī)性能衰退的統(tǒng)計(jì)量——圖3中的SPE，即是待測數(shù)據(jù)到PCS空間的距離：

因此，當(dāng)SPE超過一定程度（閾值）時(shí)，性能出現(xiàn)明顯衰退，甚至發(fā)生故障。

5）引入長度為l的滑動(dòng)窗口向量[SPE[k-(l-1)],SPE[k-(l-2)],…,SPE[k]]，l指滑動(dòng)窗口長度?？紤]在k時(shí)刻，統(tǒng)計(jì)量設(shè)計(jì)如下[14]：

H0和H1分別表示零假設(shè)、被擇假設(shè)。設(shè)發(fā)動(dòng)機(jī)當(dāng)前（未知）狀態(tài)和健康狀態(tài)剩余殘差SPE(k)序列的方差分別為當(dāng)前狀態(tài)（未知狀態(tài)）的方差估計(jì)值利用l個(gè)滑動(dòng)窗口樣本計(jì)算；考慮發(fā)動(dòng)機(jī)初始狀態(tài)為健康狀態(tài)，健康狀態(tài)方差利用初始l0個(gè)SPE(k)樣本計(jì)算。進(jìn)而假定統(tǒng)計(jì)量服從F分布[15]，以風(fēng)險(xiǎn)水平α設(shè)計(jì)（即錯(cuò)誤報(bào)警率為α）測試，設(shè)Fu=F(1-α)(l-1,l0-1)是分布概率為(1-α)的關(guān)鍵點(diǎn)。

圖4 基于主元分析和統(tǒng)計(jì)分析的性能監(jiān)控

F[k]≤Fu?H0接受，在時(shí)刻k性能良好；

Else?H1接受，在時(shí)刻k性能出現(xiàn)明顯衰退。

通過上述步驟的分析，得到基于主元分析和統(tǒng)計(jì)分析的性能監(jiān)控算法如圖4所示。

1）模塊1為獲取測量參數(shù)。利用第3節(jié)的算法估計(jì)出完整的性能參數(shù)變化過程，代入到穩(wěn)態(tài)模型中，得到對應(yīng)的測量參數(shù)；

2）模塊2為構(gòu)造性能參數(shù)統(tǒng)計(jì)量。進(jìn)一步判斷待分析測量參數(shù)的時(shí)間序號(hào)是否大于待分析總時(shí)間序號(hào)。若大于，表示分析完成，則退出；否則繼續(xù)向下分析；

3）利用主元分析和統(tǒng)計(jì)分析的方法對發(fā)動(dòng)機(jī)性能狀況進(jìn)行監(jiān)控。首先，對第i個(gè)時(shí)刻的測量參數(shù)進(jìn)行PCA分解，在PS空間中獲得表達(dá)性能變化的成分，構(gòu)造表達(dá)性能變化的統(tǒng)計(jì)量；

4）模塊3 為判斷性能趨勢。首先假定統(tǒng)計(jì)量SPE服從F分布，在一定置信水平下設(shè)置性能衰退報(bào)警閾值，跟蹤性能參數(shù)變化情況，判斷性能衰退是否超限。

5 仿真結(jié)果與分析

采用圖2 的算法，通過發(fā)動(dòng)機(jī)性能參數(shù)非線性估計(jì)，在已知測量參數(shù)基礎(chǔ)上估計(jì)出部件性能參數(shù)，將部件性能參數(shù)擴(kuò)充為狀態(tài)參數(shù)，結(jié)合擴(kuò)展卡爾曼濾波器，在某些離散時(shí)間點(diǎn)設(shè)置狀態(tài)參數(shù)，估計(jì)部件性能參數(shù)的變化。結(jié)果見圖5～8。圖中的非線性參數(shù)估計(jì)結(jié)果是在已知測量參數(shù)的基礎(chǔ)上，基于非線性模型得到的部件（風(fēng)扇、壓氣機(jī)、高壓渦輪和低壓渦輪）的性能參數(shù)，表征發(fā)動(dòng)機(jī)性能參數(shù)的變化趨勢。

圖5 風(fēng)扇流量衰退估計(jì)結(jié)果

圖6 壓氣機(jī)流量衰退估計(jì)結(jié)果

圖7 高壓渦輪效率衰退估計(jì)結(jié)果

圖8 低壓渦輪效率衰退估計(jì)結(jié)果

分析圖5～8 可知，卡爾曼濾波器能夠很好地估計(jì)風(fēng)扇流量、壓氣機(jī)流量、高壓渦輪效率、低壓渦輪效率的衰退過程，反映了部件性能參數(shù)逐漸衰退的過程。但在個(gè)別點(diǎn)的卡爾曼濾波的參數(shù)估計(jì)結(jié)果與非線性估計(jì)表征的性能參數(shù)變化趨勢存在較大差異，因?yàn)樵诳柭鼮V波方法中，在此間斷點(diǎn)設(shè)置了部件性能參數(shù)。圖9和圖10是圖8中時(shí)間樣本2 000～2 800點(diǎn)的性能變化情況。包括了性能變化的2種情況。

圖9 低壓渦輪效率漸變衰退

圖10 低壓渦輪效率突變衰退

分析圖9和圖10可知，低壓渦輪效率在時(shí)間樣本2 000～2 500時(shí)變化比較平穩(wěn)，此時(shí)滿足假設(shè)的Δp˙=0的條件。圖10 則說明了在時(shí)間樣本2 670點(diǎn)附近，低壓渦輪效率出現(xiàn)突變，此時(shí)已不滿足Δp˙=0，但在時(shí)間樣本2 670點(diǎn)以后，效率重新平穩(wěn)下來，滿足Δp˙=0的條件。

圖11 性能衰退統(tǒng)計(jì)量

6 結(jié)論

本文結(jié)合卡爾曼濾波算法和非線性估計(jì)理論，在給定離散點(diǎn)測量參數(shù)基礎(chǔ)上估計(jì)部件性能參數(shù)衰退過程。以部件性能參數(shù)為基礎(chǔ)，擴(kuò)展為卡爾曼濾波器的狀態(tài)參數(shù)，結(jié)合主元分析方法，提取反映部件性能衰退的成分，構(gòu)造了反映性能衰退的統(tǒng)計(jì)量，假定該統(tǒng)計(jì)量服從F分布，對性能衰退作進(jìn)一步分析，得到在10 000點(diǎn)附近性能衰退最為嚴(yán)重，已超過閾值。統(tǒng)計(jì)量分析只能說明整體性能下降，并不能具體定位是哪個(gè)部件性能下降。通過辨別統(tǒng)計(jì)量曲線來確定部件性能衰退將是以后研究的方向。

[1] 李本威，李冬，李姜華，等.單級(jí)壓氣機(jī)性能衰退定量研究[J].航空動(dòng)力學(xué)報(bào)，2010，25（7）：1588-1594.

LI BENWEI，LI DONG，LI JIANGHUA，et al. Quantative research of single stage compressor performance deterioration[J]. Journal of Aerospace Power，2010，25（7）：1588-1594.（in Chinese）[2] MATHIOUDAKIS K，KAMBOUKOS PH，STAMATIS A. Turbofan performance deterioration tracking using nonlinear models and optimization techniques[J]. ASME Laboratory of Thermal Turbo-machines，2002，124（3）：580-587.

[3] ZAITA A V，BALEY，G，KARLSONS G. Performance deterioration modeling in aircraft gas turbine engines[J].Gas Turbines Power，1998，120：344-349.

[4] DOEL D. An assessment of weight-least-squares-based gas path analysis[J]. Engineering for Gas Turbine and Power，1994，90（2）：366-373.

[5] DEPOLD H，GASS F. The application of expert systems and neural networks to gas turbine prognotics and diagnostics[J].Engineering for Gas Turbine and Power，1999，112（5）：607-612.

[6] VOLPONI A，DEPOLD H，GANGULI R，et al. The use of kalman filter and neural network methodologies in gas turbine performance diagnostics∶a comparative study[J].Engineering for Gas Turbine and Power，1999，110（4）：602-606.

[7] KOBAYASHI T，SIMON D L. A hybrid neural networkgenetic algorithm technique for aircraft engine performance diagnostics[C]//37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.2001：141-149.

[8] 尹大偉. 某型渦扇發(fā)動(dòng)機(jī)建模及模型修正研究[D]. 煙臺(tái)：海軍航空工程學(xué)院，2007.

YIN DAWEI. Research of X turbofan engine modeling and model modification[D]. Yantai：Naval Aeronautical and Astronautical University，2007.（in Chinese）

[9] 付堯明，傅強(qiáng).非線性航空發(fā)動(dòng)機(jī)性能參數(shù)估計(jì)方法研究[J].機(jī)械設(shè)計(jì)與制造，2006，32（8）：140-141.

FU RAOMING，F(xiàn)U QIANG. Research of aero-engine nonlinear performance parameter estimation method[J].Mechanical Design and Manufacture，2006，32（8）：140-141.（in Chinese）

[10]吳丹.航空推進(jìn)系統(tǒng)非線性性能尋優(yōu)控制[D].西安：西北工業(yè)大學(xué)，2004.

WU DAN. Optimal control of aerospace propulsion system nonlinear performance[D]. Xi’an：Northwestern Polytechnical University，2004.（in Chinese）

[11]張鵬.基于卡爾曼濾波的航空發(fā)動(dòng)機(jī)故障診斷技術(shù)研究[D].南京：南京航空航天大學(xué)，2008.

ZHANG PENG. Research of aeroengine fault diagnosis technology based on kalman filter[D]. Nanjing：Nanjing University of Aeronauticas and Astronauticas，2008.（in Chinese）

[12] GOMEZ G，LENDASSE A. Statistical fault isolation with principal component analysis[M]. California：IFAC Safe Process，2000：18-78.

[13] BORGUET S，LéONARD Q. Coupling principle component analysis and kalman filter algorithms for on-line aircraft engine diagnostics[J].Control Engineering Practice，2009，17（2）：494-502.

[14] FASSOIS S，SAKELLARIOU J. Time series methods for fault detection and identification in vibrating structures[J]. Philosophical Transactions of the Royal Society∶Mathematical，Physical and Engineering Sciences，2007，18（3）：121-128.

[15] DIMOGIANOPOULOS D，HIOS J，F(xiàn)ASSOIS S. Aircraft engine health management via stochastic modeling of flight data interrelation[J]. Aerospace and Technology，2011，12（6）：1-12.

[16]吳翊，李永樂，胡慶軍.應(yīng)用數(shù)理統(tǒng)計(jì)[M].長沙：國防科技大學(xué)出版社，1994：81-84.

WU YI，LI YONGLE，HU QINGJUN. Application and Statistics[M]. Changsha：National Defense Technology University Press，1994：81-84.（in Chinese）