含裂紋故障的齒輪系統(tǒng)動(dòng)力學(xué)特性研究及其故障特征分析

王 旭, 伍 星, 肖正明, 曹永立， 劉 韜

(1.昆明理工大學(xué) 機(jī)電工程學(xué)院，昆明 650500； 2.昆明云內(nèi)動(dòng)力股份有限公司，昆明 650500)

含裂紋故障的齒輪系統(tǒng)動(dòng)力學(xué)特性研究及其故障特征分析

王 旭1, 伍 星1, 肖正明1, 曹永立2， 劉 韜1

(1.昆明理工大學(xué) 機(jī)電工程學(xué)院，昆明 650500； 2.昆明云內(nèi)動(dòng)力股份有限公司，昆明 650500)

考慮齒輪的時(shí)變嚙合剛度、傳動(dòng)誤差和軸承支撐剛度的影響，建立含齒根裂紋故障的齒輪系統(tǒng)多自由度力學(xué)模型，基于動(dòng)力學(xué)方法對(duì)其故障機(jī)理進(jìn)行研究。通過材料力學(xué)的方法計(jì)算齒輪在正常和含裂紋兩種情況下的嚙合剛度，對(duì)比兩種剛度曲線的變化趨勢(shì)，便于進(jìn)行精確的動(dòng)力學(xué)特性分析;對(duì)建立的模型求解系統(tǒng)的動(dòng)態(tài)響應(yīng)，結(jié)果表明當(dāng)齒根存在裂紋時(shí)，其時(shí)域波形中會(huì)出現(xiàn)周期性的沖擊現(xiàn)象，頻譜中在嚙合頻率的基頻及其倍頻等地方形成一系列等間隔的邊頻譜線，其間隔大小等于故障齒輪的轉(zhuǎn)頻;這些邊頻成分幅值較低，能量分散且分布不均勻，在不同頻帶的幅值大小存在差異。針對(duì)上述特點(diǎn)，通過正交小波包方法對(duì)信號(hào)的頻帶進(jìn)行分解，應(yīng)用倒頻譜分析各子頻帶信號(hào)的邊頻成分;結(jié)果表明，該方法能夠有效的提高信號(hào)的信噪比，有助于識(shí)別和提取信號(hào)中由裂紋故障引起的邊頻成分。

齒輪系統(tǒng);裂紋故障;嚙合剛度;邊頻譜線;正交小波包;倒頻譜

隨著工業(yè)領(lǐng)域中生產(chǎn)技術(shù)的進(jìn)步和自動(dòng)化程度的提高，針對(duì)設(shè)備的在線監(jiān)測(cè)和故障診斷技術(shù)正在受到更多的重視。齒輪傳動(dòng)系統(tǒng)作為機(jī)械設(shè)備中重要的運(yùn)動(dòng)和動(dòng)力傳遞部件在高速和重載等惡劣的工作環(huán)境中容易發(fā)生故障而影響到設(shè)備的可靠運(yùn)行，安全生產(chǎn)和經(jīng)濟(jì)效益[1]。因此，針對(duì)齒輪傳動(dòng)裝置的設(shè)備監(jiān)測(cè)和故障診斷技術(shù)的發(fā)展和應(yīng)用也顯得越來越重要。在實(shí)際的振動(dòng)監(jiān)測(cè)過程中齒輪的故障特征信號(hào)多具有信號(hào)微弱，噪聲干擾強(qiáng)，故障種類復(fù)雜和不易識(shí)別等特點(diǎn)，目前針對(duì)齒輪傳動(dòng)裝置的故障診斷技術(shù)還不夠成熟，有待進(jìn)一步研究和完善。

針對(duì)齒輪故障的研究，國內(nèi)外學(xué)者分別從齒輪的動(dòng)力學(xué)機(jī)理研究和信號(hào)特征提取等方面做了大量的工作。Wu等[2]基于齒輪系統(tǒng)的6自由度動(dòng)力學(xué)模型分析了齒輪剛度隨裂紋擴(kuò)展的變化趨勢(shì)。Ding等[3]研究了齒輪磨損故障的動(dòng)力學(xué)行為。馬銳等[4]建立了單級(jí)齒輪的扭轉(zhuǎn)動(dòng)力學(xué)模型，對(duì)齒輪在含裂紋情況下的非線性動(dòng)力學(xué)特性進(jìn)行研究。Mohammed等[5]使用短時(shí)傅里葉變換方法對(duì)齒輪裂紋故障進(jìn)行了特診提取。萬志國等[6]提出了齒輪嚙合剛度的一種修正計(jì)算方法，并采用小波方法對(duì)含裂紋故障的齒輪振動(dòng)信號(hào)進(jìn)行特征提取。李云等[7]提出一種基于包絡(luò)角域同步平均的齒輪故障診斷方法。

本文在總結(jié)前人研究成果的基礎(chǔ)上結(jié)合動(dòng)力學(xué)仿真和故障特征分析方法開展齒根裂紋故障的研究。首先通過動(dòng)力學(xué)的方法揭示仿真信號(hào)的邊頻特征，然后采用正交小波包和倒頻譜相結(jié)合的方法對(duì)齒輪裂紋故障特征進(jìn)行分析。

1 齒輪的嚙合剛度計(jì)算

1.1 正常齒輪的嚙合剛度計(jì)算

根據(jù)材料力學(xué)方法可將輪齒簡(jiǎn)化為變截面的懸臂梁(見圖1)，此時(shí)嚙合輪齒的綜合彈性剛度由懸臂梁的彎曲剛度kb，剪切剛度ks和徑向壓縮剛度kc三部分組成，對(duì)應(yīng)的公式表達(dá)如下

圖1 輪齒受力示意圖

(1)

(2)

(3)

(4)

(5)

(6)

Ai=2Whi

(7)

(8)

式中:E為彈性模量;W為齒寬；Hp為節(jié)圓齒厚;v為泊松比。

實(shí)際輪齒的根部是彈性的輪緣，因此需要考慮齒輪基體的柔性變形，齒輪基體柔性變形所對(duì)應(yīng)的剛度可表示為

(9)

式中,各參數(shù)詳見文獻(xiàn)[8]。

輪齒嚙合中的赫茲剛度表示為

(10)

綜上所述，一對(duì)齒輪副總的嚙合剛度可表示為

(11)

1.2 含裂紋齒輪的嚙合剛度計(jì)算

圖2為輪齒裂紋的擴(kuò)展模型，圖2中裂紋存在于齒根處，假設(shè)裂紋和齒輪中心線的夾角為αc，裂紋深度為q。

圖2 裂紋的擴(kuò)展模型

當(dāng)齒根產(chǎn)生裂紋時(shí)，輪齒的受力區(qū)域?qū)l(fā)生變化，剛度計(jì)算中的截面模量Ii別和截面面積Ai將變?yōu)?/p>

(12)

(13)

將式(12)和式(13)代入齒輪剛度計(jì)算公式中可得到對(duì)應(yīng)的含裂紋情況下的輪齒嚙合剛度。齒輪參數(shù)如表1所示，假設(shè)在從動(dòng)輪的齒根處產(chǎn)生裂紋，剛度計(jì)算結(jié)果如圖3所示，圖中分別包含了正常齒輪，裂紋深度為1 mm和裂紋深度為2 mm的嚙合剛度曲線。

表1 齒輪參數(shù)

圖3 嚙合剛度曲線

為進(jìn)一步揭示嚙合剛度的頻域特征，本文擬對(duì)剛度曲線進(jìn)行階次譜分析。值得強(qiáng)調(diào)的是，在轉(zhuǎn)速恒定且未知的情況下，采用階次譜比頻譜更為方便，階次譜不需要考慮時(shí)間信息，而得到的譜線結(jié)構(gòu)與頻譜一致，只是橫坐標(biāo)的表示有所區(qū)別。如圖4所示，圖4(a)為正常齒輪嚙合剛度曲線對(duì)應(yīng)的階次譜，齒輪的1階代表轉(zhuǎn)頻，圖4中的譜線成分包括75階(代表齒輪的嚙合頻率)及其倍階次成分；圖4(b)為齒根裂紋為2 mm時(shí)對(duì)應(yīng)嚙合剛度曲線的階次譜，對(duì)比兩圖可以發(fā)現(xiàn)，齒根產(chǎn)生裂紋時(shí)會(huì)在低階部分出現(xiàn)一系列間隔為1階的譜線，這些譜線占據(jù)的譜帶較寬，說明其在對(duì)應(yīng)的頻譜中的頻帶也較寬。這些低頻譜線在動(dòng)力學(xué)響應(yīng)中與嚙合頻率及其倍頻成分形成頻率耦合而出現(xiàn)一系列的組合頻率，這些組合頻率圍繞在嚙合頻率及其倍頻附近，構(gòu)成其邊頻成分。具體內(nèi)容將在第3節(jié)(動(dòng)力學(xué)仿真)作進(jìn)一步說明。

2 齒輪系統(tǒng)動(dòng)力學(xué)模型

根據(jù)考慮的因素的不同，齒輪動(dòng)力學(xué)模型的建立包括多種方法，為了精確的描述齒輪系統(tǒng)的動(dòng)態(tài)行為，本文考慮了軸承的支撐剛度，建立齒輪傳動(dòng)系統(tǒng)的6自由度動(dòng)力學(xué)模型，如圖5所示。圖5中，Gear 1和Gear 2分別代表主動(dòng)齒輪和從動(dòng)齒輪。

(a) 正常齒輪嚙合階次譜

(b) 齒根裂紋2 mm的嚙合階次譜

圖5 齒輪6自由度動(dòng)力學(xué)模型

根據(jù)牛頓力學(xué)的方法建立方程組如下

(14)

其中，

式中:p12(t)為主動(dòng)輪相對(duì)于從動(dòng)輪在輪齒接觸點(diǎn)公法線方向上的相對(duì)位移；m1、m2為主動(dòng)輪和從動(dòng)輪的質(zhì)量；r1、r2為主動(dòng)輪和被動(dòng)輪的基圓半徑；J1、J2為主動(dòng)輪和被動(dòng)輪的轉(zhuǎn)動(dòng)慣量；k12(t) 、c12(t)為齒輪的嚙合剛度和阻尼；k1bx、k1by、k2bx、k2by和c1bx、c1by、c2bx、c2by分別為主動(dòng)輪和從動(dòng)輪的軸承剛度和阻尼；e12(t)為齒輪副的嚙合誤差，其公式為

e12(t)=Afcos(wft+φf)

(15)

式中，wf為動(dòng)輪的嚙合頻率。

方程組中存在剛體位移，消除剛體位移后方程組可表示為

(16)

3 動(dòng)力學(xué)仿真

相較于正常齒輪，輪齒裂紋的存在會(huì)導(dǎo)致齒輪系統(tǒng)在運(yùn)動(dòng)過程中產(chǎn)生周期性的沖擊響應(yīng)。在齒輪傳動(dòng)過程中，當(dāng)含有齒根裂紋的輪齒參與嚙合時(shí)，其嚙合剛度減小，對(duì)應(yīng)的剛度激勵(lì)改變，進(jìn)而引起系統(tǒng)的沖擊響應(yīng)。這種沖擊具有周期性，沖擊頻率為從動(dòng)輪的轉(zhuǎn)頻。由于裂紋的存在，齒輪嚙合剛度的頻譜會(huì)在低頻處形成一系列等間隔的譜線，這些譜線占據(jù)的頻帶較寬(如圖4所示)，在齒輪系統(tǒng)的動(dòng)態(tài)響應(yīng)中這些低頻成分會(huì)與齒輪的嚙合頻率及其倍頻成分等進(jìn)行頻率耦合而形成組合頻率，在圖譜上表現(xiàn)為齒輪的嚙合頻率及其倍頻等附近的邊頻成分，對(duì)應(yīng)的邊頻范圍較寬，邊頻間隔等于從動(dòng)輪的轉(zhuǎn)頻。

(a) 時(shí)域波形

(b) 幅值譜

Fig.6 Time domain waveforms and amplitude spectrum(normal)

正常齒輪的頻譜主要包含嚙合頻率及其倍頻信息，其基頻為1 125 Hz，等于齒輪的嚙合頻率,如圖6所示。而在圖7～圖8中還出現(xiàn)了明顯的邊頻帶，邊頻的間隔為15 Hz，等于從動(dòng)輪的轉(zhuǎn)頻。隨著裂紋的擴(kuò)展，時(shí)域信號(hào)的沖擊會(huì)不斷增強(qiáng)，相較于圖7，圖8中的時(shí)域波形的沖擊明顯增強(qiáng)，邊頻幅值明顯增大，從幅值譜圖的放大顯示(見圖8(c))可以看到，當(dāng)齒根裂紋為1 mm時(shí)，頻譜中出現(xiàn)的邊頻幅值最高為0.21，而齒根裂紋為2 mm時(shí)，頻譜中出現(xiàn)的邊頻幅值最高為0.52，邊頻成分明顯增強(qiáng)。在實(shí)際的齒輪故障診斷中，當(dāng)齒輪產(chǎn)生齒根裂紋時(shí)，其振動(dòng)信號(hào)會(huì)在時(shí)域波形中出現(xiàn)周期性的沖擊，沖擊頻率為裂紋齒輪所在軸的轉(zhuǎn)頻，對(duì)應(yīng)的頻譜在其嚙合頻率及倍頻等位置附近出現(xiàn)邊頻帶，邊頻間隔等于沖擊頻率。對(duì)比試驗(yàn)分析結(jié)果可知，該仿真信號(hào)符合齒輪箱故障診斷中實(shí)測(cè)信號(hào)的規(guī)律，驗(yàn)證了該仿真模型的可靠性，同時(shí)能夠?yàn)楣こ虒?shí)際信號(hào)的邊頻現(xiàn)象提供一種基于動(dòng)力學(xué)原理的解釋。

這些由齒根裂紋故障所形成的邊頻成分具有幅值較低，能量分散且分布不均勻的特點(diǎn)，在不同頻帶的幅值大小存在差異，如圖7(c)和圖8(c)所示。在嚙合頻率的基頻或不同倍頻附近出現(xiàn)的邊頻譜線的幅值具有明顯的差異，其中在嚙合頻率的5倍頻附近的邊頻幅值最高。在實(shí)際的測(cè)試中，振動(dòng)信號(hào)包含有復(fù)雜的噪聲，由裂紋故障引起的邊頻分量往往會(huì)被噪聲信息淹沒而難以分析。因此，針對(duì)邊頻能量在不同頻帶分布不均勻的特點(diǎn)，可以通過帶通濾波將頻帶分解為多個(gè)子頻段，然后對(duì)每個(gè)子頻段進(jìn)行分析。該方法一方面可以有效濾除寬頻噪聲的干擾，另一方面能夠篩選出邊頻幅值較高的頻段，突出邊頻信息，提高信噪比，能夠更為有效的發(fā)現(xiàn)和提取信號(hào)的邊頻特征。

(a) 時(shí)域波形

(b) 幅值譜

(c) 幅值譜放大圖

Fig.7 The time domain waveforms and amplitude spectrums(with crack fault:1 mm)

(a) 時(shí)域波形

(b) 幅值譜

(c) 幅值譜放大圖

Fig.8 The time domain waveforms and amplitude spectrums(with crack fault:2 mm)

4 故障特征分析

根據(jù)上節(jié)對(duì)信號(hào)特征的分析，本節(jié)采用正交小波包與倒頻譜相結(jié)合的方法來對(duì)信號(hào)的邊頻譜線進(jìn)行分析與討論。正交小波包能夠?yàn)樾盘?hào)提供更加精細(xì)的分析方法，可以同時(shí)對(duì)信號(hào)的低頻和高頻進(jìn)行多級(jí)分解，這種分解既無冗余，也無疏漏，便于實(shí)現(xiàn)信號(hào)的時(shí)頻局部化分析[9]。而且，分解后的子頻帶有效消除了齒輪嚙合頻率及其倍頻成分的干擾，便于利用倒頻譜對(duì)邊頻成分進(jìn)行分析。正交小波包分解過程可以形象的表示為一個(gè)二進(jìn)制樹結(jié)構(gòu)圖，如圖9所示為小波包的3層分解過程。

在Daubechies小波系中濾波器長度系數(shù)為5的“db5”小波對(duì)突變信號(hào)比較靈敏且階數(shù)比較低，本文選用“db5”小波函數(shù)作小波包分解，并選取小波包的第3層分解結(jié)果進(jìn)行討論。圖9中(0,0)表示原始信號(hào)，dij表示在節(jié)點(diǎn)(i,j)處的小波包系數(shù)。對(duì)小波包進(jìn)行重構(gòu)后，用sij表示dij的重構(gòu)信號(hào)，則小波包第3層分解重構(gòu)后的信號(hào)依次為s30～s37，表2表示了這8個(gè)信號(hào)所代表的原信號(hào)的頻率范圍，表中fc代表原信號(hào)的截止頻率。

圖9 二進(jìn)制樹結(jié)構(gòu)圖

表2 8個(gè)信號(hào)所代表的原信號(hào)的頻率范圍

Tab.2 The frequency range represented by the 8 signals

信號(hào)頻帶范圍信號(hào)頻帶范圍s30[0,fc/8]s34[7fc/8,fc]s31[fc/8,2fc/8]s35[6fc/8,7fc/8]s32[3fc/8,4fc/8]s36[4fc/8,5fc/8]s33[2fc/8,3fc/8]s37[5fc/8,6fc/8]

圖10～圖11為s30～s37信號(hào)的時(shí)域波形，在進(jìn)行分析時(shí)將分別對(duì)各個(gè)信號(hào)單獨(dú)進(jìn)行分析，以便于更為有效的分析邊頻成分。

通過正交小波包方法將信號(hào)分解為多個(gè)子頻段的信號(hào)后，可以對(duì)各個(gè)信號(hào)單獨(dú)進(jìn)行分析，針對(duì)信號(hào)在嚙合頻率及其倍頻等地方出現(xiàn)的邊頻譜線具有數(shù)量較多，幅值較小，而且具有明顯的周期性，適宜采用倒頻譜進(jìn)行故障特征的提取，倒頻譜分析方法能將圖譜上同一系列的邊頻譜線簡(jiǎn)化為倒頻譜上的單根或幾根譜線，譜線的位置是原譜圖上邊頻的頻率間隔的倒數(shù)，譜線的高度反應(yīng)了這一系列邊頻成分的強(qiáng)度[10]。

(a) s30信號(hào)

(b) s31信號(hào)

(c) s32信號(hào)

(d) s33信號(hào)

(a) s34信號(hào)

(b) s35信號(hào)

(c) s36信號(hào)

(d) s37信號(hào)

圖12～圖13為s30～s37信號(hào)對(duì)應(yīng)的倒頻譜，圖譜中主要譜線的間隔在0.066～0.067，對(duì)應(yīng)的頻率近似為15 Hz，等于從動(dòng)輪的轉(zhuǎn)頻。各個(gè)倒頻率對(duì)應(yīng)的幅值大小不同，其中s36中的倒頻率幅值最高，第一條譜線的值約為0.115 8，反映出的信號(hào)的邊頻特征最為明顯。

5 結(jié) 論

本文通過材料力學(xué)的方法計(jì)算含裂紋齒輪的嚙合剛度，將其代入齒輪系統(tǒng)的動(dòng)力學(xué)方程組中求解系統(tǒng)的動(dòng)態(tài)響應(yīng)，從時(shí)域中可以發(fā)現(xiàn)波形中存在周期性的沖擊成分，在頻譜嚙合頻率的基頻及其倍頻等地方形成一系列等間隔的邊頻譜線，其間隔大小等于故障齒輪的轉(zhuǎn)頻。這些邊頻成分幅值較低，能量分散且分布不均勻，而且在不同頻帶的幅值大小存在差異。在實(shí)際的測(cè)試中，振動(dòng)信號(hào)包含有復(fù)雜的噪聲，由裂紋故障引起的邊頻分量往往會(huì)被噪聲信息淹沒而難以分析?；谏鲜鎏卣?，本文首先通過正交小波包方法對(duì)信號(hào)的頻帶進(jìn)行分解，然后將獲得的各個(gè)子頻帶信號(hào)應(yīng)用倒頻譜進(jìn)行分析，提取信號(hào)的邊頻成分。結(jié)果表明，采用正交小波包與倒頻譜相結(jié)合的方法更加有助于發(fā)現(xiàn)和提取信號(hào)中由裂紋故障引起的邊頻成分。

(a) s30信號(hào)

(b) s31信號(hào)

(c) s32信號(hào)

(d) s33信號(hào)

(a) s34信號(hào)

(b) s35信號(hào)

(c) s36信號(hào)

(d) s37信號(hào)

[1] 李潤方, 王建軍. 齒輪系統(tǒng)動(dòng)力學(xué)—振動(dòng)·沖擊·噪聲[M]. 北京: 科學(xué)出版社, 1997.

[2] WU Siyan, ZUO Mingjian, PAREY A. Simulation of spur gear dynamics and estimation of fault growth[J]. Journal of Sound and Vibration, 2008, 317(3): 608-624.

[3] DING H, KAHRAMAN A. Interaction between nonlinear spur gear dynamic and surface wear[J]. Journal of Sound and Vibration, 2007, 307(3): 662-679.

[4] 馬銳, 陳予恕. 含裂紋故障齒輪系統(tǒng)的非線性動(dòng)力學(xué)研究[J]. 機(jī)械工程學(xué)報(bào), 2011, 47(21): 84-90.

MA Rui, CHEN Yushu. Nonlinear dynamic research on gear system with cracked failure[J]. Journal of Mechanical Engineering, 2011, 47(21): 84-90.

[5] MOHAMMED O D, RANTATALO M. Dynamic response and time-frequency analysis for gear tooth crack detection[J]. Mechanical Systems and Signal Processing, 2016, 66-67: 612-624.

[6] 萬志國, 訾艷陽,曹宏瑞,等. 時(shí)變嚙合剛度算法修正與齒根裂紋動(dòng)力學(xué)建模[J]. 機(jī)械工程學(xué)報(bào), 2013, 49(11): 153-160.

WAN Zhiguo, ZI Yanyang, CAO Hongrui, et al. Time-varying mesh stiffness algorithm correction and tooth crack dynamic modeling[J]. Journal of Mechanical Engineering, 2013, 49(11): 153-160.

[7] 李云, 郭瑜, 那靖,等. 基于包絡(luò)同步平均的齒輪故障診斷[J]. 振動(dòng)與沖擊, 2013, 32(19): 17-62.

LI Yun, GUO Yu, NA Jing, et al. Gear fault diagnosis based on envelope synchronous average[J]. Journal of Vibration and Shock, 2013, 32(19): 17-62.

[8] CHAARI F, FAKHFAKH T, HADDAR M. Analytical modeling of spur gear tooth crack and influence on gear mesh stiffness[J]. European Journal of Mechanics A/Solids, 2009, 28(3): 461-468.

[9] 胡廣書. 現(xiàn)代信號(hào)處理教程[M]. 北京: 清華大學(xué)出版社, 2004.

[10] 楊國安. 齒輪故障診斷實(shí)用技術(shù)[M]. 北京: 中國石化出版社, 2012.

Dynamic characteristics of a gear system with crack fault and its fault feature analysis

WANG Xu1, WU Xing1, XIAO Zhengming1, CAO Yongli2, LIU Tao1

(1. Faculty of Mechanical and Electrical Engineering, Kunming University of Science and Technology, Kunming 650500, China; 2. Kunming Yunnei Power Co., Ltd., Kunming 650500, China)

A multi-DOF mechanical model of a gear system with crack fault was established to study its failure mechanism based on the dynamic method, considering effects of time-varying mesh stiffness, transmission error and bearing support stiffness. First of all, meshing stiffnesses of gears in two cases including normal condition and the one with crack fault were calculated for the accurate dynamic characteristics analysis. The results of the system’s dynamic responses solved with the built dynamical model showed that there are some periodic impulses in time domain waveform, and there are some side frequencies equally spaced near the mesh fundamental frequency and its harmonic components, the interval between two side-frequencies is equal to the rotating frequency of the faulty gear; the amplitudes of side frequencies are lower and their energy is dispersive and not uniform, these amplitudes within different frequency bands are different. Aiming at the above mentioned features, the signal frequency band was decomposed with the orthogonal wavelet packet method and each sub-frequency band signal was analyzed using the cepstrum analysis method. The results showed that the proposed method can effectively improve the signal-to-noise ratio of signals, and make it easier to identify and extract side frequency components caused due to crack fault in signals to be detected.

gear system; crack fault; mesh stiffness; side frequencies; orthogonal wavelet packet； cepstrum

國家自然科學(xué)基金資助項(xiàng)目(51465021)；云南省教育廳科學(xué)研究基金重大專項(xiàng)項(xiàng)目(ZD2013004)；云南省自然科學(xué)基金(2013FB014)

2015-11-18 修改稿收到日期：2016-03-10

王旭 男，博士生，1987年生

伍星 男，博士，教授，博士生導(dǎo)師，1973年生

TH113

A

10.13465/j.cnki.jvs.2017.09.012