楊 宇，王歡歡，喻鎮(zhèn)濤，袁 輝

(1.湖南大學(xué) 汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410082;2.東風(fēng)汽車公司 技術(shù)中心試驗(yàn)部，襄陽(yáng) 441004)

轉(zhuǎn)子系統(tǒng)振動(dòng)信號(hào)往往會(huì)呈現(xiàn)一定程度的非線性特征，當(dāng)轉(zhuǎn)子發(fā)生故障時(shí)，其振動(dòng)信號(hào)會(huì)混有能量較大的與轉(zhuǎn)速有關(guān)的背景信號(hào)和噪聲，傳統(tǒng)的時(shí)域或頻域方法很難從這樣的轉(zhuǎn)子振動(dòng)信號(hào)中有效提取故障信息[1]，而分形理論中的分形維數(shù)作為定量刻畫混沌吸引子的一個(gè)重要參數(shù)，在研究非線性系統(tǒng)方面具有獨(dú)特的優(yōu)勢(shì)。其中關(guān)聯(lián)維數(shù)計(jì)算簡(jiǎn)單，便于從數(shù)據(jù)序列中直接測(cè)定，因而應(yīng)用廣泛，并已應(yīng)用于機(jī)械故障診斷之中[2]，故可以用來(lái)提取轉(zhuǎn)子系統(tǒng)的故障特征。相關(guān)研究表明，關(guān)聯(lián)維數(shù)在較大噪聲水平下難以用關(guān)聯(lián)維數(shù)來(lái)準(zhǔn)確刻畫系統(tǒng)的工作狀態(tài)和故障類型[3]，因此要提高計(jì)算精度就必須對(duì)信號(hào)進(jìn)行降噪處理;而當(dāng)信號(hào)的非線性特性突出時(shí)，利用線性濾波等傳統(tǒng)方法進(jìn)行降噪的效果并不明顯[4]。采用 EMD(Empirical Mode Decomposition) 方 法[5]和 LMD (Local Mean Decomposition)方法[6]都可以實(shí)現(xiàn)對(duì)振動(dòng)信號(hào)的自適應(yīng)濾波，但這兩種方法在理論上都還存在一些問題，如EMD的過包絡(luò)、欠包絡(luò)和頻率混淆問題，LMD也存在由于解調(diào)而引起的信號(hào)突變的問題。目前這些問題仍在研究當(dāng)中[7]。為了降低轉(zhuǎn)子振動(dòng)信號(hào)的噪聲，本文采用內(nèi)稟時(shí)間尺度分解(ITD)改進(jìn)算法對(duì)轉(zhuǎn)子振動(dòng)信號(hào)進(jìn)行分解以達(dá)到降噪的目的。

ITD方法是由Frei提出的另一種自適應(yīng)時(shí)頻分析方法[8]，該方法自適應(yīng)地將非平穩(wěn)信號(hào)分解成為若干瞬時(shí)頻率具有物理意義的單分量信號(hào)(固有旋轉(zhuǎn)分量)。文獻(xiàn)[8]雖然將ITD方法和EMD方法進(jìn)行了對(duì)比，分析了ITD方法在端點(diǎn)效應(yīng)和計(jì)算速度上的優(yōu)勢(shì)，但是并沒有對(duì)ITD方法及固有旋轉(zhuǎn)分量的物理意義進(jìn)行闡述。因此本文在闡述ITD方法及其固有旋轉(zhuǎn)分量的物理意義的基礎(chǔ)上，提出了基于ITD的改進(jìn)算法，并利用ITD改進(jìn)方法對(duì)轉(zhuǎn)子振動(dòng)信號(hào)進(jìn)行分解，將包含主要故障信息的若干分量進(jìn)行組合作為降噪后的轉(zhuǎn)子振動(dòng)信號(hào)，并對(duì)降噪前后轉(zhuǎn)子系統(tǒng)5種狀態(tài)下的信號(hào)分別計(jì)算了關(guān)聯(lián)維數(shù)，對(duì)比發(fā)現(xiàn)，降噪后的關(guān)聯(lián)維數(shù)區(qū)分性明顯增強(qiáng)，因此可以用該方法識(shí)別轉(zhuǎn)子的工作狀態(tài)和故障類型。

1 ITD改進(jìn)方法

1.1 ITD方法

ITD方法將待分析的非平穩(wěn)信號(hào)分解成一系列的固有旋轉(zhuǎn)分量和一個(gè)單調(diào)的趨勢(shì)信號(hào)。設(shè)Xt是待分析的原信號(hào)，分解前先定義一個(gè)基線提取算子L，使得從原始信號(hào)中去掉該基線后剩下的余量信號(hào)成為一個(gè)固有旋轉(zhuǎn)分量。一次分解的表達(dá)式為[8]:

式中假設(shè)在[0，τk]上定義了 Lt和 Ht，而 Xt在[0，τk+2]上有定義。在連續(xù)極值點(diǎn)[τk，τk+2]區(qū)間上定義 Xt的基線提取因子L，即:

式中:a是分解時(shí)的增益控制參數(shù)，0＜a＜1;Lt保留了信號(hào)在各個(gè)極值點(diǎn)的單調(diào)性，Ht提取了各個(gè)極值點(diǎn)之間疊加的局部高頻分量信號(hào)——固有旋轉(zhuǎn)分量。重復(fù)分解過程就可得到一系列固有旋轉(zhuǎn)分量和單調(diào)趨勢(shì)信號(hào)。

給出如式(4)所示的仿真信號(hào)，它由一個(gè)調(diào)幅調(diào)頻信號(hào)和一個(gè)余弦信號(hào)組成:

圖1是該仿真信號(hào)的波形，圖2是運(yùn)用ITD算法對(duì)信號(hào)進(jìn)行分解得到的結(jié)果。

圖1 仿真信號(hào)x(t)Fig.1 Simulation signal x(t)

圖2 仿真信號(hào)x(t)的ITD分解結(jié)果Fig.2 Decomposition result of simulation signal x(t)by ITD

從圖2可以看出，分解得到的第1個(gè)分量效果較好，但是第2個(gè)分量出現(xiàn)了失真，這是因?yàn)镮TD方法采用線性變換的方法獲得基線信號(hào)，使得波形出現(xiàn)了毛刺而失真。因此，為了克服該缺陷，本文提出了改進(jìn)的ITD方法，即采用三次樣條插值代替了ITD方法中的線性變換。雖然采用了三次樣條插值，但它又不同于EMD中基于局部極值的包絡(luò)均值，因?yàn)镮TD改進(jìn)方法每次分解只需要進(jìn)行一次三次樣條插值。

1.2 ITD 改進(jìn)方法

文獻(xiàn)[8]中雖然給出了ITD方法的詳細(xì)分解過程，但是并沒有對(duì)該方法及分解得到的結(jié)果的物理意義進(jìn)行闡述;同時(shí)經(jīng)過以上的研究發(fā)現(xiàn)，由于ITD方法采用線性變換的方法實(shí)現(xiàn)對(duì)信號(hào)的分解，從第二個(gè)分量開始，出現(xiàn)了較明顯的信號(hào)失真。因此，本節(jié)首先闡述了ITD分解方法的物理意義，接著提出了ITD改進(jìn)算法。

EMD與LMD方法在提供有效的自適應(yīng)時(shí)頻分析方法的同時(shí)，提供了一條很好的自適應(yīng)時(shí)頻分析思路，即首先假設(shè)任意一個(gè)復(fù)雜信號(hào)由若干個(gè)瞬時(shí)頻率具有物理意義的單分量信號(hào)組成，并給出瞬時(shí)頻率具有物理意義的單分量信號(hào)所需要滿足的條件，然后據(jù)此條件對(duì)原始信號(hào)進(jìn)行自適應(yīng)地分解。

在ITD方法中，對(duì)于理想的固有旋轉(zhuǎn)分量，其如式(2)所示的基線信號(hào)在所有的區(qū)間[τk，τk+2]上都應(yīng)該等于零，因此如式(3)所示的Lk+1應(yīng)該等于零。理論上ITD方法對(duì)信號(hào)進(jìn)行分解得到的固有旋轉(zhuǎn)分量應(yīng)該滿足其基線信號(hào)控制點(diǎn)Lk+1等于零。本文在此基礎(chǔ)上定義了瞬時(shí)頻率具有物理意義的內(nèi)稟尺度分量(Intrinsic Scale Component，ISC)。所定義的ISC分量滿足以下條件:

(1)在整個(gè)數(shù)據(jù)段內(nèi)，任意兩個(gè)極值點(diǎn)之間呈單調(diào)性;

圖3 ISC需要滿足的條件Fig.3 The needed condition of ISC

(2) 在整個(gè)數(shù)據(jù)段內(nèi)，其極值點(diǎn)為 Xk，k=1，2，…，M，各個(gè)極值點(diǎn)相對(duì)應(yīng)的時(shí)刻為 τk，k=1，2，…，M，由任意兩個(gè)極大(小)值點(diǎn)(τk，Xk)、(τk+2，Xk+2)連接形成的線段在其中間極小(大)值點(diǎn)(τk+1，Xk+1)相對(duì)應(yīng)時(shí)刻 τk+1的函數(shù)值Xk)與該極小(大)值 Xk+1的比值關(guān)系不變，即，即滿足，如圖3 所示。其中，a∈(0，1)為一常量，典型地，a=1/2，例如正弦或余弦信號(hào)、調(diào)幅信號(hào)、調(diào)頻信號(hào)、調(diào)幅－調(diào)頻信號(hào)等。

以上兩個(gè)條件保證了ISC分量任意兩個(gè)相鄰極值點(diǎn)之間具有單一的模態(tài)，而且在局部(極值點(diǎn)與相鄰的零交叉點(diǎn)之間)近似吻合標(biāo)準(zhǔn)正弦曲線，因此瞬時(shí)頻率具有物理意義。

改進(jìn)的ITD算法假設(shè)任何復(fù)雜信號(hào)是由不同ISC分量組成，任何兩個(gè)ISC分量之間相互獨(dú)立，這樣就可以將復(fù)雜信號(hào)分解為若干個(gè)ISC分量之和，具體步驟如下:

(1)確定 X(t)的所有極值點(diǎn) Xk，k=1，2，…，M 及其相對(duì)應(yīng)的時(shí)刻 τk，k=1，2，…，M，并設(shè)置參數(shù) a，計(jì)算各基線控制點(diǎn)Lk，即式(3);然后對(duì)所有Lk進(jìn)行三次樣條插值，得到基信號(hào)線L1。

(2)將L1從原始信號(hào)中分離出來(lái)，得到P1。理想地，P1為一個(gè)ISC分量，則P1為信號(hào)x(t)的第1個(gè)分量。理想的ISC分量應(yīng)滿足Lk+1等于零，實(shí)際中可以設(shè)定變動(dòng)量Δ，當(dāng)時(shí)迭代結(jié)束。

(3)如P1不滿足ISC的條件，則將P1作為原始信號(hào)重復(fù)步驟(1)和(2)，循環(huán)k次，直到得到內(nèi)稟尺度分量Pk，Pk即為信號(hào)x(t)的第1個(gè)分量ISC1。

(4)將ISC1從x(t)中分離出來(lái)，得到一個(gè)新的信號(hào)r1，將r1作為原始信號(hào)重復(fù)步驟(1)、(2)和(3)，得到x(t)的第二個(gè)滿足ISC條件的分量ISC2。重復(fù)循環(huán)n次，得到信號(hào)x(t)的n個(gè)滿足ISC條件的分量，直到rn為一單調(diào)函數(shù)為止。這樣便可以將x(t)分解為n個(gè)內(nèi)稟尺度分量ISC和一個(gè)單調(diào)函數(shù)rn之和，即:

ISC迭代終止判據(jù)Δ的選擇對(duì)信號(hào)的分解和迭代次數(shù)都有重要的影響，類似于EMD方法內(nèi)稟模態(tài)函數(shù)(IMF)的判據(jù)。常用的判據(jù)方法有，Huang在文獻(xiàn)[9]中介紹的標(biāo)準(zhǔn)偏差法，F(xiàn)landrin在文獻(xiàn)[10]中介紹的三參數(shù)法等。其中標(biāo)準(zhǔn)偏差法與IMF分量的定義無(wú)關(guān)。因此，本文在考慮ISC分量定義的基礎(chǔ)上，結(jié)合三參數(shù)法，來(lái)確定變動(dòng)量Δ的取值。

經(jīng)過仿真信號(hào)分解實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)基線控制點(diǎn)Lk+1與信號(hào)的對(duì)應(yīng)極值點(diǎn)Xk+1之比，小于一定取值范圍，即Δ取值范圍介于0.005時(shí)，則得到的分量可以近似為一個(gè)ISC分量，此時(shí)迭代次數(shù)一般為4～10次。

為了減少ITD改進(jìn)算法中邊界效應(yīng)的影響，在處理數(shù)據(jù)之前，需對(duì)原始數(shù)據(jù)進(jìn)行延拓處理，處理方法可采用G.Rilling程序的鏡像對(duì)稱延拓方法[11]。

在處理了端點(diǎn)效應(yīng)之后，對(duì)式(4)所示的仿真信號(hào)采用改進(jìn)ITD算法進(jìn)行分解，分解時(shí)取a=0.5，迭代終止判據(jù)，分解后得到如圖4所示的兩個(gè)分量和一個(gè)余量，兩次的迭代次數(shù)分別為5次和4次。采用EMD方法對(duì)該仿真信號(hào)進(jìn)行分解，在延拓方法和迭代終止判據(jù)相同的情況下，分解結(jié)果如圖5所示，兩次的迭代次數(shù)分別為17次和105次。顯然，ITD改進(jìn)算法的迭代次數(shù)要少得多，因而分解速率更快。

圖4 仿真信號(hào)的改進(jìn)ITD分解結(jié)果(已處理端點(diǎn)效應(yīng))Fig.4 Decomposition result of simulation signal by improved ITD(with boundary processing)

圖5 仿真信號(hào)的EMD分解結(jié)果(已處理端點(diǎn)效應(yīng))Fig.5 Decomposition result of simulation signal by EMD(with boundary processing)

圖6 是由電渦流位移傳感器測(cè)得的不平衡故障狀態(tài)下的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)，對(duì)其采用ITD改進(jìn)算法進(jìn)行分解，共得到3個(gè)基本分量和一個(gè)余量，如圖7所示，從圖中可以看出，第1個(gè)分量頻率成分最高，第2個(gè)分量主要為與轉(zhuǎn)速有關(guān)的信號(hào)，因此故障信息主要包含在第1和第2個(gè)ISC分量中，可以選擇這兩個(gè)分量進(jìn)行組合，得到降噪后的轉(zhuǎn)子信號(hào)。

圖6 不平衡故障狀態(tài)下的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)Fig.6 Rotor radial displacement vibration signal with unbalance fault

圖7 不平衡故障狀態(tài)下的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)的改進(jìn)ITD分解結(jié)果Fig.7 Decomposition result of rotor vibration signal under state of unbalance by improved ITD

2 實(shí)驗(yàn)分析

本文采用南京東大測(cè)振儀器廠生產(chǎn)的ZT－3型轉(zhuǎn)子振動(dòng)模擬實(shí)驗(yàn)臺(tái)進(jìn)行轉(zhuǎn)子系統(tǒng)故障實(shí)驗(yàn)，實(shí)驗(yàn)裝置如圖8所示。實(shí)驗(yàn)時(shí)鍵相傳感器采用電渦流傳感器，可以提供相位和轉(zhuǎn)速信號(hào);轉(zhuǎn)子徑向位移振動(dòng)信號(hào)由垂直和水平安裝的電渦流傳感器拾取，信號(hào)經(jīng)過信號(hào)調(diào)理箱處理后，送入數(shù)據(jù)采集系統(tǒng)。實(shí)驗(yàn)?zāi)M了不平衡、不對(duì)中、碰摩和油膜渦動(dòng)四種故障狀態(tài)[12]，信號(hào)采集時(shí)對(duì)應(yīng)轉(zhuǎn)速為3 000 r/min，采樣率為4 096 Hz。

圖9～12分別是轉(zhuǎn)子在正常、不對(duì)中、碰摩、油膜渦動(dòng)狀態(tài)下的徑向位移振動(dòng)信號(hào)(不平衡狀態(tài)下的振動(dòng)信號(hào)見圖6)。

圖8 轉(zhuǎn)子試驗(yàn)臺(tái)裝置Fig.8 Rotor test rig

圖9 正常狀態(tài)下的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)時(shí)域波形Fig.9 Rotor radial displacement vibration signal under normal state

圖10 不對(duì)中故障的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)時(shí)域波形Fig.10 Rotor radial displacement vibration signal with misalignment fault

圖11 碰摩故障的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)時(shí)域波形Fig.11 Rotor radial displacement vibration signal with rub-impact fault

圖12 油膜渦動(dòng)故障的轉(zhuǎn)子徑向位移振動(dòng)信號(hào)時(shí)域波形Fig.12 Rotor radial displacement vibration signal with oil whirl fault

首先，原始振動(dòng)信號(hào)不進(jìn)行降噪處理，直接進(jìn)行分形維數(shù)的計(jì)算，計(jì)算結(jié)果如圖13所示?？紤]到不同嵌入維數(shù)下計(jì)算得到的關(guān)聯(lián)維數(shù)也不同，因此選取了一定范圍的嵌入維數(shù)，從圖中可以看出，各狀態(tài)對(duì)應(yīng)的分形維數(shù)相互交叉，因此對(duì)原始信號(hào)直接計(jì)算關(guān)聯(lián)維數(shù)并不能清楚地分辨出轉(zhuǎn)子的工作狀態(tài)和故障類型。

圖13 不同狀態(tài)轉(zhuǎn)子原始信號(hào)在不同嵌入維數(shù)下的關(guān)聯(lián)維數(shù)Fig.13 The correlation dimension of vibration signals under different embedding dimension

然后，采用ITD改進(jìn)算法對(duì)轉(zhuǎn)子各個(gè)狀態(tài)的原始信號(hào)進(jìn)行降噪處理。即對(duì)各信號(hào)進(jìn)行分解，得到若干分量。研究發(fā)現(xiàn)，各狀態(tài)下轉(zhuǎn)子振動(dòng)信號(hào)的故障信息主要集中在前兩個(gè)分量中，選取了不同狀態(tài)下轉(zhuǎn)子振動(dòng)信號(hào)的前兩個(gè)分量重構(gòu)信號(hào)。進(jìn)一步對(duì)這些重構(gòu)信號(hào)計(jì)算不同嵌入維數(shù)下的關(guān)聯(lián)維數(shù)，結(jié)果如圖14所示。從圖中可以看出，經(jīng)過ITD改進(jìn)算法降噪后的各狀態(tài)轉(zhuǎn)子信號(hào)在一定嵌入維數(shù)范圍內(nèi)的關(guān)聯(lián)維數(shù)明確分離，可以明顯地區(qū)別各工作狀態(tài)和故障類型，可見轉(zhuǎn)子信號(hào)的改進(jìn)ITD降噪結(jié)合關(guān)聯(lián)維數(shù)可以準(zhǔn)確識(shí)別轉(zhuǎn)子的工作狀態(tài)和故障類型。嵌入維數(shù)的選取恰當(dāng)與否對(duì)相空間的重構(gòu)質(zhì)量有重要的影響。如果選取過小，將無(wú)法容納系統(tǒng)的吸引子，會(huì)出現(xiàn)偽鄰近點(diǎn);反之太大，不僅造成計(jì)算工作量大，還減少可使用數(shù)據(jù)長(zhǎng)度，使所建相空間中的相點(diǎn)過于稀疏。關(guān)于嵌入維數(shù)的選取，實(shí)際工程應(yīng)用中，一般情況，缺少對(duì)系統(tǒng)分形維數(shù)的先驗(yàn)知識(shí)。對(duì)于實(shí)際系統(tǒng)嵌入維數(shù)m的選擇比較困難，文獻(xiàn)[13]中給出了經(jīng)驗(yàn)公式m≥D，本文中嵌入維數(shù)的范圍在5～35時(shí)各狀態(tài)關(guān)聯(lián)維數(shù)的區(qū)分效果明顯。

圖14 轉(zhuǎn)子信號(hào)經(jīng)改進(jìn)ITD分解剔除噪聲和背景信號(hào)后在不同嵌入維數(shù)下的關(guān)聯(lián)維數(shù)Fig.14 The correlation dimension of vibration signals denoised by improved ITD under different embedding dimension

3 結(jié)論

本文將改進(jìn)ITD算法和關(guān)聯(lián)維數(shù)相結(jié)合應(yīng)用于轉(zhuǎn)子的故障診斷中，先通過改進(jìn)算法對(duì)轉(zhuǎn)子信號(hào)進(jìn)行分解降噪，再計(jì)算關(guān)聯(lián)維數(shù)來(lái)提取故障特征，以識(shí)別轉(zhuǎn)子的不同工作狀態(tài)和故障類型。由于關(guān)聯(lián)維數(shù)受噪聲的影響較大，而轉(zhuǎn)子原始信號(hào)中往往含有噪聲信號(hào)和與轉(zhuǎn)速有關(guān)的背景信號(hào)，在求取關(guān)聯(lián)維數(shù)之前，通過ITD改進(jìn)算法，可將轉(zhuǎn)子原始信號(hào)的故障主要信息提取出來(lái)，以剔除噪聲和背景信號(hào)，再計(jì)算關(guān)聯(lián)維數(shù)，可獲得較好的效果，從而為轉(zhuǎn)子系統(tǒng)故障診斷提供了一種新的方法。值得提出的是，ITD方法剛剛提出，還有許多理論問題需要研究和完善，如迭代終止條件、端點(diǎn)效應(yīng)等問題，因此對(duì)ITD方法的研究還不是很深入。盡管如此，本文所采用的先定義瞬時(shí)頻率具有物理意義的單分量信號(hào)，然后再據(jù)此進(jìn)行自適應(yīng)信號(hào)分解的方法，開拓了時(shí)頻分析的新思路，對(duì)于信號(hào)分析和故障診斷都有較大的借鑒意義。

[1]McFadden P D.Window functions for calculating of the time domain averages of the vibration of the individual planet gears and sun gear in an epicyclical gearbox[J].Journal of Vibration and Acoustics，1994，116(1):179 －187.

[2]Alberto R，Maria M.A method for the correlation dimension estimation for on-line condition monitoring of large rotating machinery[J].Mechanical Systems and Signal Processing.2005，19:939 －954

[3]Argyris J.The influence of noise on the correlation dimension of chaotic attractors[J].Chaos Solitions ＆ Fractals，1998，9(3):343－361.

[4]徐美芳，韓 焱，董劍龍.數(shù)字射線圖像降噪中的線性濾波方法[J].無(wú)損檢測(cè)，2005，27(5):259 －263.XU Mei-fang，HAN Yan，DONG Jian-long.Linear filter methods for image denoising I digital radiography[J].Nondestructive Inspection，2005，27(5):259 －263.

[5]Huang N E，Shen Z，Long S R.A new view of nonlinear water waves:the Hilbert spectrum[J].Annu.Rev.Fluid Mech，1999，31:417－457.

[6]Smith J S.The local mean decomposition and its application to EEG perception data[J].Journal of the Royal Society Interface，2005，2(5):443－454.

[7]Lei Y G，He Z J，Zi Y Y.Application of the EEMD method to rotor fault diagnosis of rotating machinery[J].Mechanical Systems and Signal Processing，2009，23:1327－1338.

[8]Frei M G，Osorio I.Intrinsic time-scale decomposition:timefreqeucny-energy analysis and real-time filtering of nonstationary signals[J].Proceedings of the Royal Society，A，2007，463:321－342.

[9]Huang N E，Shen Z，Long S R.The empirical mode decomposition and the Hilbert spectrum for nonlinear nonstationary time Series Analysis[J]. Proc. Roy. Soc.London，1998，454:903 －995.

[10]Rilling G，F(xiàn)landrin P.On empirical mode decomposition and its algorithms[J].IEEE Sig.Pro.Lett，2003:326 －329.

[11]Rilling G，F(xiàn)landrin P，Goncalves P.On empirical mode decomposition and its algorithms[J].In IEEE-EURASIP Workshop on Nonlinear Signaland Image Processing，Grado，Italy，June 8－11，2003.

[12]趙利華.旋轉(zhuǎn)機(jī)械故障機(jī)理與軸心軌跡識(shí)別方法研究[D].大連:大連理工大學(xué)，2010.

[13]Maestrello L.Chaotic response of panel vibrations forced by turbulent boundary layer and sound[J].AIAA Journal，1999，37(3):289 －294.