基于改進HHT和馬氏距離的齒輪故障診斷

周小龍， 劉薇娜， 姜振海； 馬風雷

(1. 長春理工大學 機電工程學院， 長春 130022； 2. 長春工業大學 機電工程學院， 長春 130012)

基于改進HHT和馬氏距離的齒輪故障診斷

周小龍1， 劉薇娜1， 姜振海2； 馬風雷2

(1. 長春理工大學 機電工程學院， 長春 130022； 2. 長春工業大學 機電工程學院， 長春 130012)

針對齒輪振動信號非線性和非平穩的特點，提出一種基于改進希爾伯特-黃變換與馬氏距離相結合的故障診斷方法。利用自適應白噪聲的完備經驗模態分解將齒輪振動信號分解成一系列固有模態函數，并采用敏感固有模態函數判別算法判斷出對故障信息敏感的模態函數；通過對敏感固有模態分量的局部希爾伯特瞬時能量譜的分析，得出信號能量隨時間變化的精確表達；以不同故障信號局部希爾伯特瞬時能量譜的最大峰值作為特征向量，采用馬氏距離對齒輪故障進行狀態識別。試驗結果表明，該方法可有效識提取齒輪故障特征，實現不同故障狀態識別。

齒輪；自適應白噪聲完備經驗模態分解；瞬時能量譜；馬氏距離；故障診斷

齒輪是機械設備中的常用部件，而齒輪傳動也是機械設備中常見的傳動方式。通常，齒輪故障是導致機械設備失效的主要原因，因此，對齒輪進行故障診斷具有重要意義[1]。

當齒輪發生故障時，其振動信號往往是非線性、非平穩的時變信號，對于該類信號時頻分析是有效的方法[2]。以經驗模態分解(Empirical Mode Decomposition，EMD)為核心的希爾伯特-黃變換(Hilbert-Huang Transform，HHT)方法是近年來由Huang等[3]提出的一種處理非線性、非平穩信號的有效方法。然而在EMD分解過程中易產生模態混疊現象，這將影響EMD分解效果。為解決該問題，Wu等[4]提出了平均總體經驗模態分解(Ensemble Empirical Mode Decomposition，EEMD)方法，并被廣泛應用于旋轉機械的故障診斷中。李輝等[5]提出一種基于EEMD和Teager-Huang變換的齒輪箱故障診斷方法；張超等[6]采用EEMD和支持向量機相結合的方法對齒輪故障進行了有效診斷；吳小濤等[7]針對軸承振動信號的特點，提出了基于峭度準則EEMD和改進形態濾波方法的軸承故障診斷方法。雖然基于EEMD的HHT方法在機械故障診斷中取得了一定成果，但該方法仍然存在不足，并限制該方法在故障診斷領域的發展。其主要表現為：①EEMD分解后的固有模態函數(Intrinsic Mode Function，IMF)內仍有一定幅值的噪聲殘留，運算過程中選擇加入白噪聲的幅值和分解次數將直接影響著分解結果，如果參數選擇不合理，將會降低分解效率，并使分解結果中會產生較多虛假IMF[8-9]。②EEMD分解出的IMF中，只有部分IMF能夠反映故障特征，而其它的為干擾成分，如果干擾成分不剔除，將會影響HHT的故障診斷精度。

針對上述問題，提出一種基于自適應白噪聲的完備經驗模態分解(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise，CEEMDAN)局部Hilbert瞬時能量譜和馬氏距離相結合的故障診斷方法。首先，利用CEEMDAN將非平穩的齒輪振動信號分解成一系列IMF，并采用敏感IMF判別算法判斷出對故障信息敏感的IMF；其次，根據HHT方法的基本原理構建出敏感IMF分量的局部Hilbert瞬時能量譜，用于提取故障特征；最后，利用馬氏距離對齒輪故障狀態進行識別。試驗結果表明：該方法可有效識別齒輪故障狀態。

1 改進的HHT

1.1 CEEMDAN

相較于EEMD方法，CEEMDAN在分解過程中，采用對加入高斯白噪聲進行EMD分解的方法代替EEMD分解過程中每次加入的高斯白噪聲，可實現自適應加噪[10]，并能在較少平均次數下，實現重構誤差幾乎為零[11]。因此，CEEMDAN不僅能夠克服EMD的模態混疊問題，同時能夠提高分解效率，減少虛假IMF的產生。

步驟1與EEMD方法相同，在待分解信號x(t)中加入不同幅值的高斯白噪聲ωi(t)，并進行I次分解。由EEMD對信號x(t)+ε0ωi(t)分解出的第1個IMF分量即為CEEMDAN的第1個模態分量。

步驟2在第一階段(k=1)，計算第1個殘余信號

(1)

步驟3對信號r1(t)+ε1E1(ωi(t)) (i=1,2,…,I)進行i次分解，直到得到第1個EMD模態分量為止，并計算第2個模態分量

(2)

步驟4對其余每個階段(k=2,…,K)(i=1, 2,…,I)計算第k個殘余信號

(3)

步驟5分解x(t)+εkEk(ωi(t))(i=1,2,…,I)～第1個EMD模態分量，并定義第k+1個模態分量為

(4)

步驟6執行步驟4，并將k+1。當殘余信號的極值點個數不超過2個的時候分解結束。

分解結束時，設所有模態分量的數量為K，則最終的殘差信號為

(5)

因此，原信號經CEEMDAN分解后可表示為

(6)

由CEEMDAN分解過程可知，其分解過程相對于EEMD更加完整，能夠對待分解進行精確重構。在每一個模態分解階段，根據系數εk可實現信噪比的自主選擇。

為驗證CEEMDAN的有效性，仿真一個由正弦波和沖擊成分組成的信號。分別使用EEMD和CEEMDAN對仿真信號進行分解，分解過程中加入的白噪聲標準差為0.2，分解次數500次，分解結果如圖1和圖2所示。由圖1和圖2可知，EEMD和CEEMDAN都將正弦信號和沖擊信號分解出來，但由于參數選擇問題，EEMD分解結果中出現了較多虛假IMF分量。而CEEMDAN有效減少了虛假IMF分量的產生，提高了分解準確性。

(a)仿真信號 (b)IMF1 (c)IMF2

(d)IMF3 (e)IMF4 (f)IMF5圖1 EEMD分解結果Fig.1 EEMD decomposition results

(a)仿真信號 (b)IMF1 (c)IMF2 (d) IMF3圖2 CEEMDAN分解結果Fig.2 CEEMDAN decomposition results

圖3為EEMD和CEEMDAN兩種分解方法的重構誤差。由圖3可知，CEEMDAN的重構誤差小于10-15，而EEMD由于受到參數選擇的影響，出現較大重構誤差。因此，EEMD在一定程度上影響了信號分解的完整性，若想減小重構誤差，需要提高分解次數，但這將導致計算規模的增加并降低計算速率。

(a)EEMD重構誤差 (b) CEEMDAN重構誤差圖3 不同分解方法的重構誤差Fig.3 Reconstruction error of different decomposition methods

1.2 敏感IMF判別算法

信號經CEEMDAN分解最終得到一系列IMF，雖然CEEMDAN能夠有效減少虛假IMF分量的產生，但對于機械故障診斷而言，虛假IMF會影響故障診斷的準確性，同時，真實IMF中也只有部分IMF是包含故障特征或對故障特征敏感的。因此，在對IMF進行Hilbert變換前，需判別各IMF分量與所診斷故障的相關程度，保證故障特征提取的精度與診斷的有效性。在此采用一種基于信號相關性的敏感IMF判斷法[12]。具體算法如下:

(1)計算故障信號x(t)的各IMF分量c1(t),…，cn(t)與原信號x(t)的相關系數αi；

(2)計算故障信號x(t)的各IMF分量c1(t),…，cn(t)與正常狀態信號y(t)的相關系數βi；

(3)由兩個相關系數計算得到各IMF分量包含故障信息的敏感相關系數γi，記作

γi=|αi-βi| (i=1,2,…,n)

(7)

(4)計算各IMF的故障敏感系數

(8)

(9)

(6)判斷出敏感系數差值最大值所對應的序號k，那么排序后的前k個IMF即為故障敏感的IMF。

對于該方法，其主要原理是：與故障信號x(t)相關系數越大的IMF分量，所包含信號x(t)的故障信息越多；計算與正常狀態信號y(t)的相關系數βi，可知每個IMF分量中與故障信息無關的正常信息。敏感IMF判別方法將各IMF與故障信號的相關性同正常信號的相關性相結合，配合敏感系數，能突顯出故障信息，削弱正常信息的影響。

對于旋轉機械系統而言，當系統構件出現故障時，會產生瞬態沖擊。圖1中的仿真信號可以看成一故障信號，經CEEMDAN分解后，IMF1為故障特征成分，而IMF2為虛假分量，IMF3為正常狀態信息。采用敏感IMF判別算法計算各相關系數和敏感系數，結果如圖4所示。

圖4 仿真信號的敏感IMF判別結果Fig.4 The sensitive correlation coefficient results of simulation signal

由圖4可知，故障特征成分IMF1被準確判別出來，而虛假成分IMF2和含有正常狀態信息成分的IMF3被剔除，表明該算法的準確性。

1.3 局部Hilbert瞬時能量譜

對CEEMDAD分解后的各敏感IMF分量進行Hilbert變換

(10)

構造解析信號

(11)

進一步可求出瞬時頻率

(12)

將上述敏感IMF的瞬時頻率和瞬時幅值組合，可得信號的局部Hilbert譜

(13)

(14)

由此可以定義

(15)

式(15)為局部Hilbert瞬時能量譜(簡稱：局部瞬時能量譜)，它反映了信號能量隨時間分布的情況。

常用的齒輪故障診斷方法一般基于振動能量的變化，但由于齒輪故障多為發生在各別齒上，所引起的振動能量變化在整個齒輪振動中所占比例較小，因此對于這種變化檢測比較困難。當齒輪存在局部故障時，故障齒輪會產生較大的動載荷，并會影響到相鄰齒，當故障程度增加時，動載荷增加，并會使振動能量發生較大變化。由敏感IMF求得的局部瞬時能量譜是對故障能量隨時間變化的精確表示，可有效提取相關故障信息。

由之前分析可知，仿真信號的IMF1是敏感IMF，圖5為其局部瞬時能量譜。在發生故障時，信號中的動載荷增加、振動加強，使相關時間段內信號的振動能量出現較大變化，結合時域信號特征，在振動開始和結束的階段沖擊最強。在圖5中這些現象得到準確反映，顯示出局部瞬時能量譜在故障特征提取方面的有效性。

圖5 仿真信號的局部瞬時能量譜Fig.5 The local instantaneous energy spectrum of simulation signal

2 試驗研究

為驗證所提方法在故障特征提取方面的有效性，在齒輪試驗臺上進行了齒輪具有斷齒、磨損和點蝕等三種故障的試驗，其中斷齒和點蝕故障發生在從動齒輪，磨損故障發生在驅動齒輪。試驗所用的齒輪試驗臺型號為PQZZ-II型，圖6為試驗臺裝置圖。

在齒輪故障模擬時，使用鐵錘將齒輪故障實驗臺的從動齒輪的一個齒敲掉，折斷線位于分度圓附近，用于模擬斷齒故障；對另一從動齒輪的一個齒面進行點焊，并采用油石磨蝕使焊接殘留的焊渣剝落，在其齒面上形成坑點，點蝕面積率為2%，用來模擬中等程度的點蝕故障；為模擬齒輪磨損故障，設置驅動齒輪的磨損為單齒齒面嚴重磨損，齒面磨損原因為人為磨銼，磨損程度為齒厚的15%。

圖6 齒輪試驗臺Fig.6 Gear tester

試驗過程中，電動機轉頻為50 Hz，對旋轉機械而言，軸承座是設備易產生劣化現象的易損點，為保證所測數據的準確性，同時為全面監測齒輪箱的振動狀態，加速度傳感器安裝在齒輪箱軸承座上。數據采集卡采用ADA16-8/2(LPCI)型高速多功能采集卡，設置采樣頻率為5 120 Hz，采樣時間為1 s，采樣精度為16 bt。平行軸齒輪箱中各齒輪均為標準直齒輪，其主要參數如表1所示。

表1 平行軸齒輪箱參數

試驗采集到的各種狀態條件下齒輪加速度信號如圖7所示。從圖7可知，在信號采集過程中，背景噪聲較強，而由于沒有消噪室，使信號的信噪比降低，這將影響改進HHT分析的準確性，同時增加故障檢測難度。為此，采用快速獨立分量分析(Fast Independent Component Analysis，FastICA)法對時域信號進行降噪處理，降噪后信號如圖8所示。

(a)齒輪正常狀態信號 (b)齒輪斷齒故障信號 (c)齒輪磨損故障信號 (d)齒輪點蝕故障信號圖7 齒輪振動信號Fig.7 Gear vibration signals

(a)齒輪正常狀態降噪信號 (b)齒輪斷齒降噪信號 (c)齒輪磨損降噪信號 (d)齒輪點蝕降噪信號圖8 降噪后齒輪振動信號Fig.8 Noise reduction gear vibration signals

從圖8可知，環境噪聲等因素對于信號的干擾已得到明顯改善。大部分的無用高頻分量被濾去，降噪信號的波形更加突顯出原信號的信息。各種齒輪故障信號在時域中均表現為脈沖沖擊，同高斯幅值調制信號相類似。雖然不同故障的時域信號存在差異，但是以此為依據不能準確區分故障類型。采用基于CEEMDAN和敏感IMF的改進HHT方法對降噪后齒輪故障信號進行分析，3類故障信號經CEEMDAN分解后，斷齒信號得到6個IMF分量，齒輪磨損和點蝕信號分別得到5個IMF分量。根據敏感IMF判別算法，計算各故障類型IMF分量的相關系數(αn，βn)和敏感系數(λn),如圖9所示，通過敏感IMF判別算法，選擇斷齒信號的IMF1、磨損信號的IMF1和點蝕信號的IMF1 ～ IMF3為敏感IMF分量。各故障信號中其它IMF分量中所包含的是正常信號與故障信號的共同信息，所以將它們剔除。

(a)齒輪斷齒信號的敏感IMF判別結果 (b)齒輪磨損信號的敏感IMF判別結果 (c)齒輪點蝕信號的敏感IMF判別結果圖9 故障信號的敏感IMF判別結果Fig.9 The sensitive correlation coefficient results of fault signals

對齒輪各狀態下的敏感IMF分量進行Hilbert變換，相應的局部瞬時能量譜,如圖10所示。從圖10可知，相較于齒輪正常嚙合信號，各故障信號的局部瞬時能量譜具有明顯的沖擊特征，由于齒輪存在故障，這將導致齒輪在嚙合過程中產生瞬態沖擊的較大能量，不同故障類型所產生的瞬態沖擊能量的大小也不相同。對各狀態下信號的局部瞬時能量譜分析可知，當齒輪正常嚙合時，局部瞬時能量譜有4個主要峰值，峰值均值為55.75，其中最大峰值為70；齒輪斷齒狀態下有4個主要峰值，其平均值為568.25，其中最大峰值為762； 齒輪磨損狀態下有3個主要峰值，其平均值為2 430.47，其中最大峰值為3 752；齒輪點蝕狀態下有4個主要峰值，其平均值為1 024.21，其中最大峰值為1 315。由上述分析可知，對于齒輪各種狀態信號而言，瞬時能量譜峰值的最大值和均值存在明顯差別，可將其作為故障特征進行診斷研究。

(a)齒輪正常信號局部 (b)齒輪斷齒信號局部 (c)齒輪磨損信號局部 (d)齒輪點蝕信號局部 瞬時能量譜 瞬時能量譜 瞬時能量譜 瞬時能量譜圖10 齒輪不同狀態的局部瞬時能量譜Fig.10 Local instantaneous energy spectrum of different gear conditions

3 馬氏距離的故障診斷分類

馬氏距離是多元數理統計理論中常用的判別方法之一，用于表示數據的協方差距離。該方法的計算是建立在總體樣本的基礎上，排除了變量之間相關性的干擾，是一種衡量兩個未知樣本集相似程度的有效方法[13]。由于算法簡單，并適用于小樣本問題的處理，因此，被廣泛應用于故障診斷[14]領域。

根據上節分析，選擇不同狀態下齒輪瞬時能量譜的最大值作為訓練樣本S進行馬氏距離分析，具體步驟如下：

步驟1在齒輪故障綜合診斷實驗臺上分別獲取四種工況各30組數據；

步驟2為消除背景噪聲對于故障診斷準確性的影響，采用FastICA降噪；

步驟3對降噪后各個樣本信號進行改進HHT分析，得到瞬時能量譜；

步驟4以瞬時能量譜的最大值作為訓練樣本S；

步驟6對待測齒輪故障信號x(t)按步驟2和步驟3求出其待測特征值Sx;

(16)

式中,k=1,2,3,4；分別代表齒輪正常、磨損、斷齒和點蝕四種工作狀態。

步驟8比較d1、d2、d3和d4的大小，并將最小判別距離所對應的狀態作為待測信號的故障類型。

基于局部瞬時能量譜和馬氏距離的齒輪故障診斷法流程圖，如圖11所示。

圖11 齒輪故障診斷方法流程Fig.11 The technology process of gear fault diagnosis method

試驗表明：通過小樣本即可實現對齒輪故障的識別，對于未知樣本的分類速度和識別準確率都較高，40組數據中出現了2組誤判，總體識別正確率為95%。而誤判數據的產生可能是由于測量過程中的測量誤差和各種干擾因素的影響。部分識別結果如表2所示。表2部分識別結果

表2 部分識別結果

按上述方法，分別對制動力矩的輸出電流在0.1 A和0.2 A工況條件下的60組數據進行識別，識別結果與實際故障類型基本一致，其中正常狀態全部識別正確，磨損和斷齒故障各有1組識別錯誤，齒輪點蝕有2組識別錯誤，總體診斷正確率為93.33%，進一步表明該方法的有效性。

4 結 論

(1)利用CEEMDAN分析信號，不但能夠克服EMD分解過程中出現的模態混疊現象，同時可以解決EEMD分解效率低的問題。

(2)通過敏感IMF算法可以有效剔除與故障信息無關的IMF，增加故障特征提取的準確性。

(3)Hilbert瞬時能量譜能夠準確反映出齒輪故障振動信號的能量隨時間的分布情況，利用Hilbert瞬時能量譜可以較好的提取出齒輪故障信號的特征。

(4)基于改進HHT瞬時能量譜和馬氏距離的故障診斷方法可以有效識別齒輪的故障狀態，滿足齒輪故障診斷要求，同時，實際應用時應盡量排除各種干擾因素，保證數據樣本的數量。

[ 1 ] 于德介，程軍圣，楊宇. Hilbert-Huang 變換在齒輪故障診斷中的應用[J]. 機械工程學報， 2005，41(6)： 102-107.

YU Dejie，CHENG Junsheng，YANG Yu. Application of Hilbert-Huang transform method to gear fault diagnosis[J]. Chinese Journal of Mechanical Engineering, 2005，41(6)： 102-107.

[ 2 ] 程軍圣，楊宇，于德介. 局部均值分解方法及其在齒輪故障診斷中的應用[J]. 振動工程學報，2009，22(1)： 76-84.

CHENG Junsheng，YANG Yu，YU Dejie. The local mean decomposition method and its application to gear fault diagnosis[J]. Journal of Vibration Engineering, 2009，22(1)： 76-84.

[ 3 ] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [J]. Proceeding of the Royal Society of London-Series A: Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903-995.

[ 4 ] WU Z, HUANG N E. Ensemble empirical mode decomposition: a noise-assisted data analysis method [J]. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41.

[ 5 ] 李輝，鄭海起，唐力偉. 基于EEMD 和THT 的齒輪故障診斷方法[J]. 振動、測試與診斷，2011，31(4)： 496-500.

LI Hui，ZHENG Haiqi，TANG Liwei. Gear fault diagnosis method based on EEMD and THT [J]. Journal of Vibration，Measurement amp; Diagnosis, 2011，31(4)： 496-500.

[ 6 ] 張超，陳建軍，郭迅. 基于EEMD能量熵和支持向量機的齒輪故障診斷方法[J]. 中南大學學報(自然科學版)，2012，43(3)： 932-939.

ZHANG Chao，CHEN Jianjun，GUO Xun. Gear fault diagnosis method based on ensemble empirical mode decomposition energy entropy and support vector machine[J]. Journal of Central South University (Science and Technology), 2012，43(3)： 932-939.

[ 7 ] 吳小濤，楊錳，袁曉輝，等. 基于峭度準則EEMD及改進形態濾波方法的軸承故障診斷[J]. 振動與沖擊，2015，34(2)： 38-44.

WU Xiaotao, YANG Meng, YUAN Xiaohui, et al. Bearing fault diagnosis using EEMD and improved morphological filtering method based on kurtosis criterion [J]. Journal of Vibration and Shock, 2015，34(2)： 38-44.

[ 8 ] TORRES M E, COLOMINAS M A, SCHLOTTHAUER G, et al. A complete ensemble empirical mode decomposition with adaptive noise [C] // IEEE International Conference on Speech and Signal Processing. Prague: ICASSP, 2011: 4144-4147.

[ 9 ] COLOMINAS M A, SCHLOTTHAUER G, TORRES M E. Improved complete ensemble EMD: a suitable tool for biomedical signal processing [J]. Biomedical Signal Processing amp; Control, 2014, 14(1): 19-29.

[10] 韓慶陽，孫強，王曉東，等. CEEMDAN去噪在拉曼光譜中的應用研究[J]. 激光與光電子學進展，2015，52(11)： 267-274.

HAN Qingyang，SUN Qiang，WANG Xiaodong，et al. Application of CEEMDAN in Raman spectroscopy denoising[J]. Laser amp; Optoelectronics Progress, 2015，52(11)： 268-274.

[11] 李軍，李青. 基于CEEMDAN-排列熵和泄漏積分ESN的中期電力負荷預測研究[J]. 電機與控制學報，2015，19(8)： 70-80.

LI Jun，LI Qing. Medium term electricity load forecasting based on CEEMDAN-permutation entropy and ESN with leaky integrator neurons[J]. Electric Machines and Control, 2015，19(8)： 70-80.

[12] 雷亞國. 基于Hilbert-Huang變換的機械故障診斷[J]. 機械工程學報，2011，47(5)： 71-77.

LEI Yaguo. Machinery fault diagnosis based on improved Hilbert-Huang transform[J]. Journal of Mechanical Engineering, 2011，47(5)： 71-77.

[13] 劉昱，張俊紅，畢鳳榮，等. 基于LMD邊際譜的柴油機氣門故障診斷[J]. 內燃機學報，2014，35(6)： 96-100.

LIU Yu，ZHANG Junhong，BI Fengrong，et al. Diesel engine valve fault diagnosis based on LMD marginal spectrum[J]. Chinese Internal Combustion Engine Engineering, 2014，35(6)： 96-100.

[14] 駱志高，李旭東，趙俊麗，等. 利用馬氏距離判別法準確實現對裂紋的識別[J]. 振動與沖擊，2013，32(21)：186-188.

LUO Zhigao, LI Xudong, ZHAO Junli, et al. Crack identification with mahalanobis distance discrimination method [J]. Journal of Vibration and Shock, 2013，32(21)：186-188.

GearfaultdiagnosisbasedonimprovedHHTandMahalanobisdistance

ZHOU Xiaolong1， LIU Weina1， JIANG Zhenhai2， MA Fenglei2

(1. College of Mechanical and Electric Engineering, Changchun University of Science and Technology, Changchun 130022, China; 2.School of Mechatronic Engineering, Changchun University of Technology, Changchun 130012, China)

In view of nonlinear and non-stationary characteristics of gear vibration signals, a fault diagnosis method based on improved Hilbert-Huang transform and Mahalanobis distance was proposed. The gear vibration signals were decomposed by complete ensemble empirical mode decomposition with adaptive noise, the intrinsic mode functions were obtained and sensitive intrinsic mode functions were selected by the sensitivity evaluation method. Then, the local Hilbert instantaneous energy spectrum of the sensitive intrinsic mode components was analyzed, and the fault information can be extracted from the distribution of the energy of the gear vibration signal with the change of time. Finally, the maximum peak value of the local Hilbert instantaneous energy spectrum was treated as the fault features and the Mahalanobis distance method was used for judging the gear fault. Experimental results show that the method can effectively extract gear fault features and apply for different fault identification.

gear; complete ensemble empirical mode decomposition with adaptive noise; instantaneous energy spectrum; Mahalanobis distance; fault diagnosis

國家自然科學基金資助項目(51075041)；吉林省教育廳科技發展項目(2014124)

2016-05-11 修改稿收到日期： 2016-09-20

周小龍 男，博士生，1987年生

劉薇娜 女，博士，教授，1956年生

TH17

A

10.13465/j.cnki.jvs.2017.22.034