對角切片譜提升的形態(tài)學濾波方法及其在萬向軸故障識別中的應用

李奕璠， 林建輝， 劉偉渭

(1. 西南交通大學 機械工程學院，成都 610031； 2. 西南交通大學 牽引動力國家重點實驗室，成都 610031)

對角切片譜提升的形態(tài)學濾波方法及其在萬向軸故障識別中的應用

李奕璠1， 林建輝2， 劉偉渭1

(1. 西南交通大學 機械工程學院，成都 610031； 2. 西南交通大學 牽引動力國家重點實驗室，成都 610031)

針對多尺度形態(tài)學濾波器將所有尺度的濾波結果進行算術平均這一問題，提出利用三階累積量對角切片譜對其進行改進，用于對高速列車萬向軸不平衡故障進行檢測與識別。該方法首先對安裝在齒輪箱的振動傳感器采集到的信號進行多尺度形態(tài)學濾波，得到不同尺度下的濾波結果，再計算濾波后信號的三階累積量及其對角切片譜，最后依據對角切片譜的特征頻率系數，選取出最能凸顯故障特征的切片譜，從而避免了多尺度濾波器濾波結果的算術平均問題。在萬向軸不平衡試驗臺進行了試驗，結果表明，該方法能有效地識別出萬向軸不平衡引起的基頻和倍頻故障特征，與傳統的多尺度形態(tài)濾波相比，此方法更能彰顯故障特征。

高速列車；萬向軸；多尺度形態(tài)濾波；對角切片譜

萬向軸是CRH5型列車傳動系統的重要組成部件，其兩端通過十字萬向節(jié)分別與牽引電機和齒輪箱相連[1]。由于萬向軸為細長桿狀結構，其彎曲和扭轉剛度都較小。而傳動時，萬向軸既要傳遞牽引力矩，又要適應各種運動關系[2]，在此過程中引起的附加力矩容易導致扭轉振動[3]。此外，列車長期運行也會使得萬向節(jié)軸的磨損間隙增大，導致萬向軸偏心。上述因素均可能使萬向軸產生不平衡，從而加劇傳動系統的振動，加速傳動系統萬向節(jié)和軸承的破壞，危及行車安全。因此，對萬向軸工作狀態(tài)進行監(jiān)測與識別很有必要。

目前，對萬向軸狀態(tài)進行監(jiān)測和識別的研究并不多見。文獻[4]使用第二代小波變化結合奇異值分解的方法對此問題進行了研究，但基于小波分解的方法始終存在基小波與相應參數選取主觀性的問題。文獻[5]使用經驗模態(tài)分解，根據特征模態(tài)分量的能量大小識別是否存在萬向軸不平衡，但其正常樣本與故障樣本之間的差異并不顯著。為此，本文使用基于數學形態(tài)學濾波的方法對萬向軸不平衡進行研究。相較于文獻[4-5]，形態(tài)濾波分析方法不需要人為設定參數[6]，結果具有更好的復現性，同時其計算過程只包含加、減和布爾運算，計算速度更快，實時性更好。

但是，傳統的多尺度形態(tài)濾波器的輸出是將每個尺度的濾波結果進行算術平均，這種做法并不合理，因為那些遠離理論中心尺度的濾波結果往往被噪聲污染，難以有效反映故障特征[7]。針對此問題，本文提出一種利用三階累積量對角切片譜的方法對多尺度形態(tài)濾波器進行改進，以期實現萬向軸不平衡的檢測。

1 多尺度形態(tài)濾波

形態(tài)濾波是基于數學形態(tài)學變換的非線性濾波方法。它依據待分析信號的局部形態(tài)特征，通過一定的形態(tài)學變換，將信號與噪聲分離。形態(tài)學濾波的核心思想是通過結構元素的連續(xù)移動，對待分析信號進行匹配，以達到提取信號細節(jié)、抑制噪聲的目的[8]。

膨脹和腐蝕是形態(tài)學變換中兩種最基本的運算。設原始信號為f(n) (n=1,2,…,N)，結構元素為g(m) (m=1,2,…,M)，N≥M。f(n)關于g(m)的膨脹和腐蝕分別定義為

(f⊕g)(n)=max[f(n-m)+g(m)]

(1)

(fΘg)(n)=min[f(n+m)-g(m)]

(2)

式中，符號⊕和Θ分別表示膨脹和腐蝕。

f(n)關于g(m)的開運算和閉運算可分別表示為

(f°g)(n)=(fΘg⊕g)(n)

(3)

(f·g)(n)=(f⊕gΘg)(n)

(4)

式中，符號° 和·分別表示開、閉運算。

常用的組合形態(tài)濾波器由開、閉運算組成

y(n)=[(f°g·g)(n)+(f·g°g)(n)]/2

(5)

形態(tài)學信號處理的實質是基于試探，即用不同的結構元素去匹配待分析信號[9]。傳統的形態(tài)學濾波方法中，結構元素一旦選定，其長度和高度在整個濾波過程中不會改變，這種做法對頻率結構單一的信號很有效，但是，如果待分析信號的頻率成分復雜，單一尺度的形態(tài)濾波器難以達到較好的濾波效果。

多尺度形態(tài)濾波采用不同尺度的結構元素對信號進行分析。設λ(λ=1,2,…,k)為多尺度濾波器的尺度，g(m)為單位結構元素，多尺度濾波器的結構元素可表示為

(6)

經過多尺度膨脹、腐蝕、開、閉運算后，得到的多尺度組合形態(tài)濾波器在尺度為λ的濾波結果為yλ(n)=[(f·λg°λg)(n)+(f°λg·λg)(n)]/2

(7)

多尺度濾波器的最終輸出為所有尺度濾波結果的算術平均

(8)

文獻[10]的研究結果表明，式(8)的處理方式會破壞多尺度濾波器特征提取的效果。因此，本文在第二部分引入三階累積量對角切片譜方法對多尺度濾波結果進行提升。

2 累積量對角切片譜

長期以來，一階和二階統計量被廣泛用于振動信號處理，如均值、方差、自相關函數和功率譜等，但這些方法并不能準確描述非線性和非高斯信號，這就需要從更高的階次上來反映信號特征。目前使用最為廣泛的高階統計量是三階累計量，三階累計量的傅里葉變換被稱為雙譜。由于雙譜是一個二維矩陣，計算量較大，難以實現較大數據量的實時計算[11]，而如果使用較少的數據量計算雙譜又會增加估計偏差。所以，在實際應用中，常采用三階累積量的某種一維切片來提取信號的特征，本文使用對角切片。

假設x(t) (t=0,1,2,3, …)為零均值實平穩(wěn)離散隨機過程，其三階累積量定義為

C3x(τ1,τ2)=E[x(t)x(t+τ1)x(t+τ2)]

(9)

式中：E[·]是求數學期望操作;τ1和τ2是時間滯后。三階累積量對角切片是三階累積量在(τ1,τ2)平面內一條τ1=τ2=τ的直線，經過對角切片之后，一個二維矩陣轉化為了一個數據序列

C3x(τ)=E[x(t)x(t+τ)x(t+τ)]

(10)

三階累積量對角切片C3x(τ,τ)的離散傅里葉變換稱為對角切片譜或1.5維譜

(11)

相較于雙譜，對角切片譜計算量更小，但它保留了高階統計量分析的優(yōu)點[12]：

(1) 若x(t)是高斯信號，Sx(ω)=0。這表明對角切片譜可以抑制高斯噪聲。

(2) 若x(t)是對稱分布的信號，Sx(ω)=0。這說明對角切片譜可以抑制對稱分布的噪聲。

(3) 若x(t)=p(t)+q(t)，p(t)和q(t)相互獨立，q(t)為高斯信號，那么Sx(ω)=Sp(ω)。這意味著對角切片譜可以用來分離加性非高斯信號與高斯噪聲。機械振動信號往往被視為非高斯信號，通過對角切片譜分析能夠將噪聲有效分離。

(4) 對角切片譜能夠識別二次相位耦合關系、抑制獨立頻率分量[13]。

在對萬向軸振動信號進行采集時，傳感器無法直接安裝在萬向軸上，只能退而求其次，將其安裝在臨近的齒輪箱上，這就使得采集到的信號結構較為復雜，它不僅包含了萬向軸振動信息，還有齒輪箱自身的振動信息，噪聲，甚至從其他部件傳遞來的振動信息，萬向軸故障特征往往被掩蓋，而性質1～3能有效降噪，提高信噪比，性質4能凸顯故障特征，這就是本文使用累積量對角切片譜對形態(tài)學濾波器進行改進的原因。

此外，萬向軸不平衡故障的本質為機械轉子不平衡，與正常轉子相比，不平衡轉子的故障特征主要變現為在轉動頻率及其倍頻處振動能量的顯著變化。由于萬向軸的轉動頻率較低，而齒輪嚙合的頻率較高，齒輪箱振動特征與萬向軸振動特征容易區(qū)分。

3 累積量切片譜提升的多尺度形態(tài)濾波方法

本文提出一種累積量切片譜提升的多尺度形態(tài)濾波方法，具體步驟如下：①對原始信號進行多尺度形態(tài)學濾波，得到不同尺度下的濾波結果；②計算每個尺度濾波結果的三階累積量對角切片；③得到每個三階累積量對角切片的對角切片譜；④計算每個對角切片譜的特征頻率強度系數[14]，特征頻率強度系數越大，故障特征越顯著；⑤選出最大的特征頻率強度系數對應的尺度，只有這個尺度的對角切片譜被用于故障檢測。

由于扁平結構元素結構簡單，運算速度快，在進行多尺度濾波時，本文均使用扁平結構元素。對角切片譜提升的多尺度形態(tài)濾波方法唯一需要事先確定的參數只有一個，就是多尺度濾波器的最大分解尺度。文獻[15-16]將最大的尺度設為采樣頻率與故障特征頻率的商，本文采用同樣的方法。

需要指出的是，即使在不知道故障特征頻率的情況下，本文方法依然可以進行有效分析。因為本文方法可以從多個濾波結果中選取出一個最能反映故障特征的結果，這樣就只需要將最大分解尺度設置為一個較大值就可以了。事實上，這個較大值往往不超過采樣頻率與故障特征頻率的商的0.3倍[17]。

4 試驗驗證

為了驗證本文方法的效果，在萬向軸不平衡試驗臺(圖1)進行了試驗。試驗中，萬向軸一端與齒輪箱相連，齒輪箱由電機驅動。將萬向軸另一端支撐座墊高，確保萬向軸在垂向存在傾角，以模擬列車實際運行狀態(tài)。振動傳感器安裝在齒輪箱上最靠近萬向軸的非旋轉部位，采樣頻率20 kHz。選用同一型號的兩根萬向軸進行試驗，一根是新軸，一根是存在輕微不平衡的軸，在此將其稱之為故障軸。為了將齒輪箱中齒輪嚙合頻率去除，對采集到的信號進行了1 kHz低通濾波。

圖1 萬向軸試驗臺Fig.1 The test rig of cardan shaft

4.1 正常軸的分析結果

萬向軸的轉速為880 r/min，對應的轉頻為14.7 Hz。圖2為正常萬向軸振動的時間歷程曲線和頻譜。從圖2(b)的頻譜中，無法找到萬向軸轉頻。

圖2 正常軸振動波形及頻譜Fig.2 Vibration waveform and frequency spectrum of normal cardan shaft

使用本文方法對正常軸振動信號進行處理。最小的尺度為1，由于進行了低通濾波，最大的尺度設為1 kHz與特征頻率的商。尺度的增量為1。因此，共有68個濾波尺度。圖3(a)為經多尺度形態(tài)濾波后得到的每一個尺度的濾波結果，為了方便觀察，只給出了尺度1,7,13,…,61,67等12個結果。對每個尺度的濾波結果求取對角切片譜，同樣為了便于觀察，圖3(b)也只給出了上述12個尺度對應的對角切片譜。每一個對角切片譜的特征頻率強度系數如圖3(c)所示?？梢姡叨?為最優(yōu)的濾波尺度，因為其特征頻率強度系數最大。尺度5的對角切片譜如圖3(d)所示，從圖3(d)中可以清楚地識別萬向軸轉頻和二倍頻。

使用傳統的多尺度濾波方法對圖2(a)所示的信號進行分析，相應的濾波結果及頻譜如圖4所示。為了進行公平的比較，與圖3的分析使用同樣的結構元素和尺度。從圖4(b)中可以勉強出識別萬向軸的轉頻，但頻譜中最為突出的頻率成分與轉頻及其倍頻毫無關系。造成此結果的原因在于，多尺度形態(tài)濾波器中，較小的尺度適合于提取沖擊特征，較大的尺度適合于抑制噪聲，對所有尺度的濾波結果進行加權平均會削弱沖擊特征。

(a) 每個尺度的濾波結果

(b) 每個尺度的1.5維譜

(c) 最優(yōu)尺度的選擇

(d) 最優(yōu)尺度的1.5維譜圖3 本文方法對圖2(a)信號的分析過程與結果Fig.3 Analysis processes and results for the signal presented in Fig. 2(a) by the proposed method

對比圖3(d)與圖4(b)可以看到，本文方法在抑制噪聲、提取信號特征方面的優(yōu)越性。從多個尺度的濾波結果中選擇一個最優(yōu)的濾波結果的做法是合理的，同時基于累積量對角切片譜的方法能夠進一步抑制噪聲，強化故障特征。

4.2 故障軸的分析結果

在同樣的試驗條件下，對存在輕微不平衡故障的萬向軸進行了測試，得到的時域波形及頻譜如圖5所示。從圖5(b)的頻譜中很難識別出萬向軸的轉頻和倍頻。對比圖5與圖2，無論在時域還是頻域，正常軸與故障軸的差異均不顯著。

圖4 傳統的多尺度形態(tài)濾波對圖2(a)信號的處理結果Fig.4 Analysis results of the signal presented in Fig. 2(a) by conventional multi-scale morphology filtering

圖5 故障軸振動波形及頻譜Fig.5 Vibration waveform and frequency spectrum of fault cardan shaft

使用本文提出的方法，對圖5(a)的信號進行分析。與4.1節(jié)采用同樣的結構元素，濾波的尺度同樣設為68，經多尺度濾波得到的信號見圖6(a)，與圖3一樣，為了便于觀察，只畫出了其中12個尺度的結果。圖6(b)是圖6(a)信號的對角切片譜。圖6(c)為所有尺度對角切片譜的特征頻率強度系數，尺度11的特征頻率強度系數最大，因此，確定尺度11為最優(yōu)尺度。圖6(d)為尺度11的濾波結果的對角切片譜，從圖6(d)中可以清晰地看到14.7 Hz的轉頻和二倍頻29.4 Hz。在轉頻(14.7 Hz)的幅值為正常萬向軸同一位置(圖3(d))幅值的1.45倍，同時在二倍頻位置(29.4 Hz)的振幅較圖3(d)也有所增加，而這正是典型的轉子不平衡故障特征。

此外，從圖3(c)和圖6(c)可以看到，最優(yōu)的濾波尺度分別為尺度5和11，均小于最大分析尺度的0.3倍，從而印證了文獻[17]的結論。

圖7為使用傳統的多尺度形態(tài)濾波方法對故障軸振動信號進行處理得到的時域波形和頻譜。從圖7(b)中可以檢測到萬向軸的轉頻成分，但此頻率被噪聲污染，與圖6(d)相比，故障特征并不突出，這再次說明了本文方法在故障特征提取方面更為有效。圖7(b)在轉頻(14.7 Hz)的幅值為圖4(b)相應位置的1.28倍，這雖然可以用以區(qū)分正常軸與故障軸，但差異不如本文方法顯著。

(a) 每個尺度的濾波結果

(b) 每個尺度的1.5維譜

(c) 最優(yōu)尺度的選擇

(d) 最優(yōu)尺度的1.5維譜圖6 本文方法對圖5(a)信號的處理過程與結果Fig.6 Analysis processes and results for the signal presented in Fig. 5(a) by the proposed method

圖7 傳統的多尺度形態(tài)濾波對圖5(a)信號的分析結果Fig.7 Analysis results of the signal presented in Fig. 5(a) by conventional multi-scale morphology filtering

4.3 討 論

臺架試驗中，傳感器安裝在齒輪箱上，因此，信號中齒輪嚙合的信息占有主體地位，信號中反映萬向軸振動特征的信息相對較弱。為此，在4.1節(jié)和4.2節(jié)對試驗數據進行處理時，首先進行了1 kHz低通濾波，以去除齒輪嚙合頻率，提高萬向軸故障識別能力。另一方面，如果不進行1 kHz低通濾波，直接用本文的方法對原始信號進行分析，會遇到另外一個問題，那就是多尺度形態(tài)濾波的分析尺度。一般來說，多尺度形態(tài)濾波的最大分析尺度設為采樣頻率與故障特征頻率的商，本文的采樣頻率為20 kHz，故障特征頻率為14.7 Hz，這就意味著需要使用1 360個尺度對原始信號進行分析。而多尺度形態(tài)濾波的原理決定了尺度越大，所需要的計算時間越久，文獻[18]研究發(fā)現如果尺度1的分析時間為T，那么，尺度10的分析時間約為10T，尺度1 360所需的分析時間大致為1 360T。多尺度濾波總的分析時間為每一個尺度分析時間的累積，即T+2T+3T+…+1 360T，這就使得計算代價極大。同時，反映故障沖擊特征的信息往往蘊藏在較小的分析尺度中，分析尺度設定太大會降低計算效率，對提高分析精度也沒有幫助。

5 結 論

多尺度形態(tài)濾波器中較小尺度的濾波器善于保持信號細節(jié)、提取沖擊特征，較大尺度的濾波器適用于對信號進行平滑、抑制噪聲干擾。在使用中，往往將所有尺度的濾波結果進行算術平均作為多尺度形態(tài)濾波器的輸出，這一做法會影響故障特征提取的效果。為此，本文使用一種三階累積量對角切片譜方法對傳統的多尺度形態(tài)濾波器進行改進，提出一種高速列車萬向軸故障識別的方法。本文方法的主要優(yōu)點在于，三階累積量對角切片譜具有較強的抑制隨機噪聲的能力，能夠對多尺度形態(tài)濾波的結果中存在噪聲進一步進行抑制，去除故障無關的頻率成分；最優(yōu)尺度濾波器的選擇，能進一步的凸顯故障特征。利用臺架試驗的數據對提出的方法進行了驗證，結果表明，本文方法能夠從齒輪箱測點的振動信號中有效地提取出萬向軸振動的特征，識別正常萬向軸與存在不平衡的萬向軸。相對于多尺度形態(tài)濾波器的分析結果，無論是在譜的清晰度，還是特征頻率的振動能量都得到了顯著增強。以外，本文方法計算速度快，能夠滿足實時監(jiān)測的要求。

[1] 冷揚立，李秋澤，李慶國，等. CRH5型動車組萬向軸結構及臨界轉速分析[J]. 鐵道車輛，2010，48(12)：6-11. LENG Yangli，LI Qiuze，LI Qingguo，et al. Analysis of the structure and critical rotating speed of the cardan shafts for CRH5 multiple units [J]. Railway Vehicles，2010，48(12)：6-11.

[2] 張紅軍，姚遠，羅赟，等. CRH5型動車萬向軸傳動系統技術特征分析[J]. 鐵道學報，2009，31(2)：115-119. ZHANG Hongjun，YAO Yuan，LUO Yun，et al. Analysis on technical characteristics of CRH5 cardan drive system [J]. Journal of the China Railway Society，2009，31(2)：115-119.

[3] 姚遠，張紅軍，羅赟. CRH5 型動車萬向軸扭轉振動分析[J]. 中國鐵道科學，2009，30(2)：82-86. YAO Yuan，ZHANG Hongjun，LUO Yun. Analysis on the torsional vibration of cardan shaft on CRH5 motor car [J]. China Railway Science，2009，30(2)：82-86.

[4] 丁建明，林建輝，王晗，等. 萬向軸動不平衡檢測的二代小波變換奇異值方法 [J]. 機械工程學報，2014, 50(12)：110-117. DING Jianming, LIN Jianhui, WANG Han, et al. Detection of the dynamic unbalance with cardan shaft applying the second wavelet transform and singular value decomposition [J]. Journal of Mechanical Engineering, 2014, 50(12)：110-117.

[5] 阮騰達. CRH5動車組傳動系統故障診斷方法研究 [D]. 成都：西南交通大學，2014: 51-73.

[6] CHEN Q, CHEN Z, SUN W, et al. A new structuring element for multi-scale morphology analysis and its application in rolling element bearing fault diagnosis [J]. Journal of Vibration and Control, 2015, 21(4): 765-789.

[7] LI C, LIANG M. Continuous-scale mathematical morphology-based optimal scale band demodulation of impulsive feature for bearing defect diagnosis [J]. Journal of Sound and Vibration, 2012, 331(26): 5864-5879.

[8] 胡愛軍，唐貴基，安連鎖. 基于數字形態(tài)學的旋轉機械振動信號降噪方法[J]. 機械工程學報，2006, 42(4):127-130. HU Aijun, TANG Guiji, AN Liansuo. De-noising technique for vibration signals of rotating machinery based on mathematical morphology filter [J]. Chinese Journal of Mechanical Engineering, 2006, 42(4):127-130.

[9] ZHANG L, XU J, YANG J, et al. Multiscale morphology analysis and its application to fault diagnosis [J]. Mechanical System and Signal Processing, 2008, 22(3) :597-610.

[10] 李奕璠，劉建新，林建輝，等. 基于自適應多尺度形態(tài)學分析的車輪扁疤故障診斷方法 [J]. 交通運輸工程學報，2015，15(1)：58-65. LI Yifan, LIU Jianxin, LIN Jianhui, et al. Fault diagnosis method of railway vehicle with wheel flat based on self-adaptive multi-scale morphology analysis [J]. Journal of Traffic and Transportation Engineering, 2015, 15(1):58-65.

[11] SHEN G, MCLAUGHLIN S, XU Y, et al. Theoretical and experimental analysis of bispectrum of vibration signals for fault diagnosis of gears [J]. Mechanical System and Signal Processing, 2014, 43(1):76-89.

[12] 樊養(yǎng)余，陶寶祺，熊克，等. 艦船噪聲的1(1/2)維譜特征提取 [J]. 聲學學報，2002, 27(1)：71-76. FAN Yangyu, TAO Baoqi, XIONG Ke, et al. Feature extraction of ship-radiated noise by 1(1/2)spectrum [J]. Acta Acustica, 2002, 27(1):71-76.

[13] 張嚴，王樹勛，李生紅. 二次相位耦合的1(1/2)維譜分析 [J]. 電子學報，1996, 24(4)：109-112. ZHANG Yan, WANG Shuxun, LI Shenghong. 1(1/2)dimension spectrum approach to quadratic phase coupling estimation [J]. Acta Electronica Sinica, 1996, 24(4):109-112.

[14] LI B, ZHANG P, WANG Z, et al. Gear fault detection using multi-scale morphological filters [J]. Measurement, 2011, 44(10):2078-2089.

[15] DONG Y, LIAO M, ZHANG X, et al. Faults diagnosis of rolling element bearings based on modified morphological method [J]. Mechanical System and Signal Processing, 2011, 25(3) :1276-1286.

[16] NIKOLAOU N G, ANTONIADIS I A. Application of morphological operators as envelope extractors for impulsive-type period signals [J]. Mechanical System and Signal Processing, 2003, 17(6):1147-1162.

[17] 李兵，張培林，米雙山，等. 機械故障信號的數字形態(tài)學分析與智能分類 [M]. 北京：國防工業(yè)出版社，2011.

Diagonal slice spectrum lifted morphology filtering method and its application in cardan shaft fault detection

LI Yifan1, LIN Jianhui2, LIU Weiwei1

(1. School of Mechanical Engineering, Southwest Jiaotong University, Chengdu 610031, China;2. Traction Power State Key Lab, Southwest Jiaotong University, Chengdu 610031, China)

Aiming at solving the problem that the arithmetic average of filtering results is carried out over all scales in multi-scale morphological filter, a third-order cumulant diagonal slice spectrum method was proposed to improve the conventional multi-scale morphological filter for its use in the fault detection of high speed railway cardan shaft. The multi-scale filtering of the signals collected by the vibration sensor installed on a gearbox was conducted in advance to obtain the filtering results in different scales. Then the third-order cumulants and corresponding diagonal slices of the filtered signals were calculated, and the optimal scale which can best highlight the fault characteristics was selected based on a characteristic frequency intensity coefficient. As a consequence, the arithmetic average problem of multi-scale morphological filter is avoided. A test rig experiment was conducted on an unbalance test bench of cardan shaft. The results show that the proposed method is effective to identify the fault features of fundamental frequency and multiple frequencies caused by the unbalance of cardan shaft. Comparing with the traditional multi-scale morphological filter, the method is better to manifest fault features.

high speed train, cardan shaft, multi-scale morphology filtering, diagonal slice spectrum

國家自然科學基金資助項目(61134002)；中央高?；究蒲袠I(yè)務費(SWJTU2682014BR001EM)

2015-10-08 修改稿收到日期： 2016-03-20

李奕璠 男，博士，講師，1985年4月生

U270

A

10.13465/j.cnki.jvs.2017.10.003