程衛東, 王天楊, Wang Jin-jiang, Robert X Gao, 溫偉剛, 李建勇
(1.北京交通大學機械與電子控制工程學院, 北京 100044;2.康涅狄格大學機械工程系 斯托斯 CT 06269 美國)
包絡變形對階比分析結果的影響及消除方法
程衛東1, 王天楊1, Wang Jin-jiang2, Robert X Gao2, 溫偉剛1, 李建勇1
(1.北京交通大學機械與電子控制工程學院, 北京 100044;2.康涅狄格大學機械工程系 斯托斯 CT 06269 美國)
轉速變化的滾動軸承故障診斷通常先采用計算階比跟蹤技術消除變轉速的影響,再進行階比分析獲得故障特征。然而,等角度重采樣會引起沖擊包絡變形進而導致階比分析結果產生偏差,在轉速變化范圍較大時這個偏差必須要受到重視。介紹了等角度重采樣方法,仿真了包絡變形對階比分析結果的影響,給出了峰值點間隔變化量詳細的理論推導,提出了通過修正角度變換值來消除包絡變形影響的方法,并進行了實驗驗證。研究結果顯示,包絡在寬窄方向的變形對階比分析結果沒有影響;相鄰包絡峰值點間隔的變化對階比分析結果產生影響;在角度序列逐點減去峰值點間隔變化量可消除包絡變形對階比分析結果的影響。
故障診斷; 滾動軸承; 計算階比跟蹤; 角域重采樣; 包絡變形
在機械、能源、石化、運載和國防等行業,間歇性工作的機械設備普遍存在、量大面廣。例如礦用電鏟(一種大型挖掘機),它用鏟斗以間歇重復的工作循環進行工作。一個工作循環由挖掘、滿斗回轉至卸載地點、卸載、空斗回轉至挖掘地點4個工序構成。工作節拍是提升機構配合推壓機構完成鏟裝動作,回轉機構配合行走機構完成卸料動作,每個工作節拍之間有一定的時間間隔。各工作部的傳動機構間歇性工作,齒輪箱中的軸及滾動軸承轉速變化不定,以靜止-升速-轉速波動持續-降速-靜止的變化狀態工作。該類設備上的軸承幾乎沒有恒定的轉速,并且轉速變化范圍大。對這類設備上的滾動軸承進行故障診斷,必須要從轉速大范圍變化并且頻繁啟停工況下的振動信號中提取故障特征[1-4]。
滾動軸承故障產生的振動信號已經被廣泛研究[5-8],共振解調技術[9]及其改進技術[10-11]由于實用有效而被廣泛使用。包絡分析方法是40多年前提出的滾動軸承故障診斷的基礎方法,滾動軸承的故障診斷技術發展至今,很多技術或方法都是圍繞包絡分析方法的改進和提高。包絡分析法的核心是先對振動信號進行帶通濾波,帶通濾波的頻帶選擇采用小波變換[12]和譜峭度[13]等技術。濾波后進行幅值解調獲得包絡信號,再對其進行頻譜分析得到故障特征信息。然而,這些技術僅適用于軸承在恒定轉速的工況。在變轉速工況下,需要使用計算階比跟蹤技術來去除轉速變化的影響[14-18]。計算階比跟蹤和包絡分析法結合的集成方法被提出來用于變轉速工況下滾動軸承的故障診斷[3,19]。集成方法的流程如圖1所示,轉速的變化引起故障沖擊包絡疏密變化,不能直接對其進行頻譜分析。使用計算階比跟蹤技術,即利用轉速信息和插值算法對包絡信號進行等角度重采樣得到同步包絡信號。對同步包絡信號的譜分析(稱為階比分析)才可以得到變轉速軸承的故障特征。

圖1 計算階比跟蹤與包絡分析集成方法Fig.1 Flow chart of integrative envelope analysis and computed order tracking
然而,針對滾動軸承的故障信號使用計算階比跟蹤技術有一定的特殊性。在實現等角度重采樣時,首先需要設定一個轉速參考基準。如果參考基準設定為轉速變化范圍的最小值,那么等角度重采樣時要增加采樣點;如果參考基準設定為轉速變化范圍的最大值,則等度角重采樣時會減少采樣點。對于實際診斷系統來說,選擇轉速的最小值作為參考基準更切實可行,這樣可以避免因減少采樣點帶來的包絡信息的丟失風險,又不至于實際轉速意外超出由轉速最大值確定的參考基準而發生錯誤。
采用增加采樣點的方法(即設定轉速的最小值作為參考基準),不同轉速時所插數據點數不同會導致包絡發生變形。也就是說,沖擊包絡的峰值點與該包絡的起始點之間被拉長或被壓縮,引起相鄰包絡峰值點之間的間隔變化,進而影響階比分析的結果[20]。
針對這個問題,本文的主要工作是研究包絡變形對后續階比分析結果的影響因素及規律,并進一步提出消除方法。
對包絡進行頻譜分析時,包絡峰值之間的間隔決定著包絡周期性出現的頻率,包絡形狀及峰值的大小影響周期分量的振幅值。對于變化轉速滾動軸承的故障產生的包絡,需要知道故障包絡與變化的轉速之間的關系。
2.1 故障沖擊包絡峰值與轉速的關系
圖2顯示的沖擊響應信號是在Machinery Fault Simulator實驗平臺上獲得,是滾動軸承在900 r/min轉速時外圈單點故障產生的沖擊響應。定義峰值Amax所在時刻與包絡開始時刻之間的時間為包絡上升時間tr;峰值與包絡衰減為零時刻之間的時間為衰減持續時間tf。
包絡峰值Amax體現了沖擊能量的大小,而沖擊能量大小受故障點通過軸承載荷區速度的影響[5],即同一個故障點轉速快時激發的包絡峰值要大于轉速慢時激發的包絡峰值。

圖2 滾動軸承外圈單點故障的沖擊響應Fig.2 Waveform of fault impulse response
2.2 包絡上升時間與轉速的關系
滾動軸承變轉速運行時,同一個故障點引發的兩次相鄰沖擊處于不同的轉速工況,決定了這兩次沖擊具有不同的沖擊能量。利用Machinery Fault Simulator實驗平臺,獲得軸承外圈單點故障轉速在600~2 400 r/min之間的7種穩定轉速下的響應信號(采樣頻率24 kHz),統計各種轉速時的包絡上升時間。由于轉速的快慢影響本底噪聲水平,統計時以0.05 V做為門檻電壓,該門檻與包絡的交叉點記為沖擊響應的開始時刻。每一種轉速選擇10個相鄰的沖擊包絡來計算平均的包絡上升時間。統計結果見表1所示。
表1 不同轉速時包絡上升時間統計
Tab.1 Statistical analysis of peak time under different rotational speed

軸承轉速/(r·min-1)上升時間/s均值/s標準差/s6000.087089000.0859612000.0875815000.085790.086510.0006318000.0865421000.0862924000.08629
由各轉速下上升時間的統計結果看,滾動軸承在600~2 400 r/min轉速范圍產生的故障沖擊包絡的上升時間雖然各不相同,但是沒有明顯的變化趨勢,各轉速下的均值為0.086 51 s,標準差為0.000 63 s。由沖擊能量與響應上升時間的理論分析可以得到同樣結論(本文略)。故在后文推導過程假設不同轉速下沖擊響應的包絡上升時間不變。
3.1 等角度重采樣方法
等角度重采樣之前,時間軸是描述時域信號的橫軸,單位為s。時域振動信號數據點之間是等時間間隔Δt為
(1)
式中fst為時域采樣頻率,單位為Hz。
等角度重采樣之后,角度軸是重采樣之后用于描述角域信號的橫軸,單位設為r。時間到角度的轉換采用下式
(2)
式中φ表示t時間內軸承所轉過的角度,通過軸承實時轉速對時間的積分求得。n(t)為軸承的實時轉速,單位為r/min。
利用角域數據點構建三次樣條曲線,再對該曲線進行等角間隔重采樣。在等角度重采樣處理過程中,設定最低轉速nmin為參考基準。在最低轉速時軸承每轉過一周的采樣點數等同于時域采樣頻率,就可保證重采樣滿足采樣定理,角域重采樣頻率fsφ由下式確定
(3)
式中fsφ為角域采樣頻率,單位為“1/r”。重采樣之后角度軸上數據點之間的角度間隔Δφ為
(4)
從上述描述可見,由于實際轉速都高于參考基準,所以重采樣后角域數據點數多于原時域數據點數,并且轉速高的數據段重采樣之后的數據點數多于轉速低的數據段的數據點數,導致了不同轉速時不同程度的包絡變形。
3.2 等角度重采樣仿真
如圖1所示,故障引起的振動沖擊最終體現為包絡解調之后的上半周包絡信號。為了說明對包絡進行等角度重采樣產生的包絡變形,采用三角波模擬沖擊包絡進行仿真。仿真時設軸承旋轉一周均勻出現5個沖擊包絡,軸承轉速均勻上升時,起始轉速60 r/min,加速度為40 r/s2;軸承轉速均勻下降時,起始轉速600 r/min,加速度-40 r/s2。因為沖擊包絡的發生只和軸承轉過的角度有關,所以包絡的起點在旋轉一周內均布。圖3(a)是轉速由低向高變化的時域包絡圖,包絡的起點在時間軸上出現的位置由疏到密變化,包絡的寬度不變;圖3(c)是轉速由高向低變化的時域包絡圖,包絡起點位置由密到疏變化,包絡的寬度不變。圖3(b)是圖3(a)經過等角度重采樣后的角域包絡圖,圖3(d) 是圖3(c)經過等角度重采樣后的角域包絡圖。
從圖3(b)可見,在加速工況下等角度重采樣后包絡起點均勻分布在角度軸的一周之內,包絡的波形隨著轉速升高被逐漸拉寬,并且隨著包絡波形變寬峰值點也相對右移,相鄰包絡峰值點的間隔也逐漸變大。圖3(d)所示的減速工況,包絡的波形隨著轉速降低逐漸被壓窄,并且相鄰峰值點的間隔逐漸變小。

圖3 轉速上升、下降時等角度重采樣仿真Fig.3 Simulation in different domain during speed-up and speed-down
角域重采樣時引起的包絡變形體現在兩方面:(1)包絡橫向被拉寬或壓窄;(2)峰值點間隔逐漸變大或變小。仿真時采用三角波模擬故障引起的沖擊包絡,轉速均勻上升時,加速度為40 r/s2,起始轉速60 r/min,等角度重采樣后的角域信號包含了包絡橫向被拉寬并且峰值點間隔逐漸變大的因素,如圖4(c)所示。轉速下降時,減速度為-40 r/s2,起始轉速600 r/min,等角度重采樣后的角域信號包含了包絡橫向被壓窄并且峰值點間隔逐漸變小的因素,如圖4(d)所示。把圖4(c)和(d)兩圖中包絡的峰值點移成等間隔分布,對應成圖4(a)和(b),分別只含有包絡橫向被拉寬或壓窄的因素,峰值點間隔相同。

圖4 不同的包絡變形Fig.4 The different envelope deformation
對圖4的4種包絡做FFT計算,即階比分析,結果見圖5所示。觀察峰值點間隔相同的(a)和(b)兩條包絡的階比分析結果,故障特征峰都出現在5階比處,橫向變形沒有影響特征頻率。看(c)和(d)兩條包絡的階比分析結果,故障特征峰出現在4.82和5.13階比處,說明包絡峰值點間隔逐漸變大使得階比分析結果變小,而峰值點間隔逐漸變小導致階比分析結果變大。仿真研究結果顯示,包絡峰值點間隔逐漸變大或變小影響特征在頻率軸上的位置,而包絡拉寬或壓窄只影響特征幅值大小。

圖5 針對圖4的不同變形包絡的階比分析結果Fig.5 Order analysis on different envelope deformation
考慮到變化轉速滾動軸承在實際診斷中只取一小段信號來分析,故設軸承轉速以恒加速度均勻變化,并設相鄰的兩個包絡上升時間tr相同。如圖6所示,第1個沖擊起始點所在角度為φ1,對應的時刻為t1,瞬時轉速為n1;峰值點所在角度為φ3,對應的時刻為t3,轉速為n3。峰值點與起始點之間的關系表示為:
(5)
(6)
(7)

圖6 相鄰包絡起點和峰值點相互關系Fig.6 Relationship between adjacent impulses
隨后到達的沖擊包絡起始點所在角度為φ2,對應的時刻為t2,實時轉速為n2。由于包絡起始點的出現只與軸承轉過的角度有關,所以,相鄰沖擊起點位置間隔就是一個不變的角度φc(見圖6所示)。第2個沖擊與第1個沖擊的相互關系:
(8)
(9)
相鄰沖擊間隔角為
(10)
由公式(9)和(10),可以得出新的關系
(11)
第2個沖擊包絡的峰值點所在角度為φ4,對應的時刻為t4,實時轉速為n4,峰值點和起點之間的關系為
t4=t2+tr
(12)
參照公式(7)和(11),t4時刻對應的角度為
(13)
根據公式(7)和(13),可以求出兩個相鄰包絡峰值點的間隔
(14)
根據公式(14),可以得出包絡峰值點的間隔與包絡起點的間隔存在差異
(15)
即,由于包絡變形引起了相鄰峰值間隔的角度變化量為Δφc。相對變化量表示為
(16)
由上述理論分析可知,如果包絡的起點等間隔分布,則第2個包絡與第1個包絡的峰值點間隔增大Δφc,第3個包絡與第1個包絡的峰值點間隔增大量的計算要將2φc代入式(15),即
可見峰值點間隔的增大呈非線性。公式(15)顯示了引起包絡變形的具體因素,即軸承的瞬時轉速、加速度、相鄰包絡的間隔和包絡的上升時間。通過公式(15)還可以看到,如果上升時間tr為零,包絡峰值點間隔的變化量就為零,這也正是其他的研究認為故障沖擊時間很短可以忽略,當成一個瞬間脈沖來處理的結果。公式(15)顯示當轉速變化范圍大并且變化劇烈時,加速度及瞬時轉速對包絡變形有很大影響,所以包絡上升時間tr不能忽略。
由公式(16)可以得出結論:加速度大(轉速變化劇烈),包絡峰值點間隔變化大;加速度小(轉速變化緩慢),包絡峰值點間隔變化小。加速度為正時(升速階段)包絡峰值點間隔逐漸變大;加速度為負時(降速階段)包絡峰值點間隔逐漸變小。瞬時轉速小(低轉速段)包絡峰值點間隔變化大,瞬時轉速大(高轉速段)包絡峰值點間隔變化小。
由上文第4小節的分析可知,包絡峰值點間隔的變化是影響階比分析結果的直接原因。本文的第5小節的理論分析結論是包絡峰值點間隔的變化量與瞬時轉速、加速度、相鄰故障沖擊的間隔和沖擊包絡的上升時間這4個因素有關。在研究方法的改進時,要依據理論分析結果去除包絡峰值點間隔的變化量。
在圖7中,與時軸上的時刻t1~t6對應的轉速分別為n1~n6。按照上文分析的公式(13)有
(17)
(18)
(19)
比較上述3個公式發現,可以把各式減去峰值點間隔的變化量

圖7 時域角域變換對照Fig.7 Adjacent impulses between time-domain and angle-domain
(20)
式中φc=0。

(21)
(22)
根據公式(20),(21)和(22),可以得到各相鄰峰值之間的間隔
(23)
即,每個數據點的角度減去峰值點間隔變化量之后,可以保證重采樣后包絡峰值點間隔等同于起點間隔,即如果包絡起點間隔是固定值,那么包絡峰值點間隔也是固定值。
改進后計算階比跟蹤方法的實現步驟為:1)采用公式(2)作時間到角度的轉換;2)每個數據點的角度減去公式(15)描述的變化量;3)利用角域數據點構建三次樣條曲線,對該曲線進行等角間隔重采樣。
為了驗證上述的理論分析,在實驗臺上測得滾動軸承的振動信號。軸承的轉速由靜止上升到2 400 r/min然后再下降到靜止狀態,其中的加速度和減速度的平均值為3.021 58 r/s2。實驗平臺(如圖8所示)是SpectraQuest公司生產的Machinery Fault Simulator實驗臺。轉速信號通過電機輸出軸上的鑒相裝置(旋轉一圈產生一個光電脈沖)獲得。軸承故障類型是外滾道單點損傷,軸承型號為 MB ER-10K,故障特征頻率fBPFO=3.028×轉軸頻率。使用NI-PCI 6259數據采集系統采集振動信號和轉速信號,采樣頻率24 kHz。

圖8 實驗平臺Fig.8 Experimental setup
對振動數據去均值處理,并進行6~12 kHz的帶通濾波,采用希爾伯特法進行包絡解調獲得包絡信號。對包絡信號進行等角度重采樣,設nmin=30 r/min為參考基準,故角域采樣頻率為4 800 (1/r)。對重采樣后的角域包絡做短時傅里葉變換計算(STFT),窗口寬度為7 200數據點。
在轉速上升時計算了4段數據,4段數據對應的轉速分別是984,1 830,2 028,2 226 r/min。從圖9(a)可以看到,在轉速由低到高變化過程中,故障特征的階比值從2.86逐漸接近理論值3.028。可以證明加速上升時(加速度為正),階比結果小于勻速時的階比結果。
在轉速下降時也計算了4段數據,4段數據對應的轉速分別是2 331,2 163,1 951和1 145 r/min。從圖9(b)可以看到,在轉速由高到低變化過程中,故障特征的階比值從接近理論值3.028逐漸增加到3.24。可以證明減速下降時(加速度為負),階比結果大于勻速時的階比結果。并且,圖9(a)和(b)兩圖還可以證明低轉速時階比結果的偏差大,高轉速時階比結果的偏差小。
與圖9(a)的數據處理過程相同,根據本文第6小節的修正方法,設定包絡上升時間為0.086 51 s(依據前文的實驗結果),加速度為3.021 58 r/s2,相鄰包絡的間隔為1/3.028 r,在角度轉換時減去變化量Δφc,再進行等角度重采樣,STFT計算結果(階比分析結果)見圖8(c)所示。984, 1 830, 2 028和2 226 r/min 4種轉速時對應數據段的階比分析結果都落在了理論值3.028處。減速下降工況的修正效果同圖9(c)(略),驗證了修正方法的有效。

圖9 實驗數據的階比分析結果Fig.9 The order analysis results of the tested vibration signal
針對等角度重采樣過程中的包絡變形問題,本文分析了包絡變形對后續階比分析結果的影響因素及影響規律,并給出了影響的消除方法。
影響因素有4個,分別是加速度、瞬時轉速、包絡的上升時間和相鄰包絡的間隔。與滾動軸承工況相關的加速度及瞬時轉速兩個因素對階比分析結果的影響規律是:(1)與恒定轉速下的階比分析結果相比,加速度大(轉速變化劇烈),階比分析結果的偏差大;加速度小(轉速變化緩慢),階比分析結果偏差小。加速度為正時(升速階段)分析結果小于恒定轉速時的階比分析結果;加速度為負時(降速階段)分析結果大于恒定轉速時的階比分析結果。(2)瞬時轉速對階比分析結果的影響是在低轉速段上階比分析結果偏差大,在高轉速段上階比分析結果偏差小。
消除方法的原理是調整變形后包絡峰值點的間隔,在數據點從時域轉換到角域時引入一個去除項,仿真及實驗結果驗證了該方法的有效性。
對轉速大范圍變化并且頻繁啟停工況下的滾動軸承進行故障診斷時,通常使用計算階比跟蹤技術來去除轉速變化的影響,然后再對等角度采樣信號做進一步的頻譜分析。本文研究成果的意義在于,描述了階比跟蹤時包絡變形對后續的階比分析結果的影響因素及規律,并且給出了消除方法,使得計算階比跟蹤和包絡分析相結合的集成方法更為準確,對于這類工況下滾動軸承的故障診斷具有實際意義。
致謝:感謝美國康涅狄格大學機械工程系EMS實驗室對本文實驗的支持,特別感謝Robert X Gao教授及WANG Jinjiang博士對于本文成文的指導和修改。
[1] Amirata Y, Benbouzida M E H, Al-Ahmar E, et al. A brief status on condition monitoring and fault diagnosis in wind energy conversion systems[J]. Renewable and Sustainable Energy Reviews, 2009, 13:2 629—2 636.
[2] Cheng W, Wang T, Wen W, et al. Anomaly detection for equipment condition via frequency spectrum entropy[J], Advanced Materials Research, 2012, 433-440: 3 753—3 758.
[3] Wang J, Gao R, Yan R. A hybrid approach to bearing defect diagnosis in rotary machines[J]. CIRP Journal of Manufacturing Science and Technology, 2012,5: 357—365.
[4] Tadina M, Boltezar M. Improved model of a ball bearing for the simulation of vibration signals due to faults during run-up[J]. Journal of Sound of Vibration, 2011, 330(17):4 287—4 301.
[5] Randall R, Antoni J. Rolling element bearing diagnostics-A tutorial[J]. Mechanical Systems and Signal Processing, 2011,25: 485—520.
[6] Yan R, Gao R. Approximate entropy as a diagnostic tool for machine health monitoring[J]. Mechanical Systems and Signal Processing, 2007,21: 824—839.
[7] Sinha J K, Elbhbah K, A future possibility of vibration based condition monitoring of rotating machines[J]. Mechanical Systems and Signal Processing, 2013,34: 231—240.
[8] Jamaludin N, Mba D. Monitoring extremely slow rolling element bearings: part I[J]. NDT & E International, 2002,35: 349—358.
[9] McFadden P D, Smith D D. Vibration monitoring of rolling element bearings by the high-frequency resonance technique-a review[J]. Tribology International,1984,17: 3—10.
[10]Tandon N, Choudhury A. A review of vibration and acoustic measurement methods for the detection of defects in rolling element bearings[J]. Tribology International, 1999,32: 469—480.
[11]Sheen Y T. A complex filter for vibration signal demodulation in bearing defect diagnosis[J].Journal of Sound and Vibration,2004,276: 105—119.
[12]Yan R. Gao R, Multi-scale enveloping spectrogram for vibration analysis in bearing defect diagnosis[J]. Tribology International, 2009, 42: 293—302.
[13]Antoni J, Randall R B. The spectral kurtosis: application to vibratory surveillance and diagnostics of rotating machines[J]. Mechanical System and Signal Processing, 2006, 20: 308—331.
[14]Borghesani P, Pennacchi P, Randall R B, et al. Order tracking for discrete-random separation in variable speed conditions[J]. Mechanical Systems and Signal Processing, 2012,30: 1—22.
[15]Wang K, Heyns P. Application of computed order tracking Vold-Kalman filtering and EMD in rotating machine vibration[J]. Mechanical Systems and Signal Processing, 2011,25: 416—430.
[16]Fyfe K, Munck E. Analysis of computed order tracking[J]. Mechanical Systems and Signal Processing, 1997,11: 187—205.
[17]Brandt A, Lag? A, Ahlin K, et al. Main principles and limitations of current order tracking methods[J]. Sound and Vibration, 2005,39(3): 19—22.
[18]Saavedra P N, Rodriguez C G. Accurate assessment of computed order tracking[J]. Shock and Vibration, 2006,13: 13—32.
[19]Luo H G, Qiu H, George G, et al. Synthesized synchronous sampling technique for differential bearing damage detection[J]. ASME Journal of Engineering for Gas Turbines and Power, 2010,132: 1—16.
[20]Cheng W, Gao R, Wang J, et al. Envelope deformation in computed order tracking and error in order analysis[J]. Mechanical System and Signal Processing, 2014,48:92—102.
Effects and its elimination method of envelope deformation on order analysis
CHENGWei-dong1,WANGTian-yang1,WANGJin-jiang2,RobertXGao2,WENWei-gang1,LIJian-yong1
(1.School of Mechanical, Electronic and Control Engineering, Beijing Jiaotong University,Beijing 100044, China; 2.Department of Mechanical Engineering, University of Connecticut, Storrs, CT 06269, USA)
Defects diagnosis of rolling element bearings operating under time-varying rotational speeds entails order tracking and analysis techniques that convert a vibration signal from the time domain to the angle domain to eliminate the effect of speed variations. Angular domain re-sampling will induce shock envelope deformation leading to errors of analysis, and the effect is particularly noteworthy when the varying rate of speed is significant. This paper presents a quantitative analysis of key factors affecting the accuracy of order analysis on rotating machines with varying speeds. An analytical model is established, simulated, and experimentally evaluated. Elimination methods are proposed for which the angle value should be modified by removing the peak interval variations. The research findings show that deformation of the envelope in the width direction has no effect on results of the order analysis; interval variations of adjacent envelope peaks have influence on the results of the order analysis; removing the peak interval variations in angle sequence can eliminate the effects of envelope deformation.
fault diagnosis; rolling bearing; computed order tracking; angular resampling; deformation of envelope
2014-02-09;
2014-05-06
國家自然科學基金資助項目(51275030); 中央高校基本業務專項資金資助項目(2011BM093)
TH165+.3; TP277
A
1004-4523(2015)03-0470-08
10.16385/j.cnki.issn.1004-4523.2015.03.018
程衛東(1967—),男,副教授。電話: (010)51683874;E-mail:wdcheng@bjtu.edu.cn