熊 炘，楊世錫，甘春標(biāo)，周曉峰

(浙江大學(xué) 機(jī)械工程學(xué)系 流體傳動(dòng)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室，杭州 310027)

轉(zhuǎn)子全周碰摩與局部碰摩的識(shí)別方法研究

熊 炘，楊世錫，甘春標(biāo)，周曉峰

(浙江大學(xué) 機(jī)械工程學(xué)系 流體傳動(dòng)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室，杭州 310027)

超臨界汽輪發(fā)電機(jī)組的結(jié)構(gòu)和工況復(fù)雜，容易引起轉(zhuǎn)、靜子間的碰摩。根據(jù)碰摩誘發(fā)因素的不同，可將其分為全周碰摩與局部碰摩。由于兩種碰摩故障的時(shí)、頻特征相似，傳統(tǒng)的時(shí)、頻域分析方法很難準(zhǔn)確提取它們的故障特征。針對(duì)這一不足，提出一種基于經(jīng)驗(yàn)?zāi)Ｊ椒纸猓娈愔捣纸?EMD－SVD)與支持向量機(jī)(SVM)的碰摩故障識(shí)別方法，用于對(duì)轉(zhuǎn)子全周碰摩與局部碰摩故障進(jìn)行識(shí)別。首先，通過(guò)EMD獲取碰摩信號(hào)的固有模式函數(shù)(IMF);然后，提取表征信號(hào)主要能量的前四階IMF組成特征矩陣并進(jìn)行SVD分解，得到關(guān)于原信號(hào)的一組特征值;最后，將特征值輸入SVM，對(duì)原信號(hào)進(jìn)行分類(lèi)識(shí)別。轉(zhuǎn)子試驗(yàn)臺(tái)全周碰摩與局部碰摩試驗(yàn)結(jié)果表明，該方法對(duì)轉(zhuǎn)子全周碰摩與局部碰摩故障的分類(lèi)準(zhǔn)確率高，其中以徑向基函數(shù)作為核函數(shù)的SVM分類(lèi)準(zhǔn)確率達(dá)到96.0%。

轉(zhuǎn)子;全周碰摩;局部碰摩;經(jīng)驗(yàn)?zāi)Ｊ椒纸?奇異值分解;支持向量機(jī)

隨著我國(guó)淘汰落后產(chǎn)能、實(shí)現(xiàn)節(jié)能減排等重大戰(zhàn)略的實(shí)施，大容量、高參數(shù)化的超臨界汽輪發(fā)電機(jī)組在我國(guó)火力發(fā)電裝備中所占的比例不斷增加。較之于同容量的亞臨界機(jī)組，超臨界機(jī)組軸系跨數(shù)多、跨度大，冷、熱態(tài)溫度變化劇烈。與此同時(shí)，為了提高運(yùn)行效率，這類(lèi)機(jī)組的動(dòng)靜間隙往往控制在較小范圍內(nèi)。因此，超臨界汽輪發(fā)電機(jī)組比一般亞臨界機(jī)組更易發(fā)生轉(zhuǎn)、靜子之間的碰摩現(xiàn)象，必須加強(qiáng)對(duì)其運(yùn)行狀態(tài)的監(jiān)測(cè)，及早發(fā)現(xiàn)并消除故障隱患。

轉(zhuǎn)、靜子之間的碰摩故障一般為機(jī)組的二次故障，其誘發(fā)因素多為轉(zhuǎn)子不平衡量過(guò)大、軸系不對(duì)中、流體激勵(lì)引起的軸系失穩(wěn)等一次故障［1］。當(dāng)一次故障使轉(zhuǎn)子與碰摩點(diǎn)始終保持接觸時(shí)，稱為全周碰摩。Bently等［2］對(duì)轉(zhuǎn)子與汽封的全周碰摩現(xiàn)象進(jìn)行了系統(tǒng)的實(shí)驗(yàn)研究;Choi等［3］通過(guò)構(gòu)建碰摩試驗(yàn)裝置，對(duì)全周碰摩的機(jī)理進(jìn)行了數(shù)值仿真研究;徐尉南等［4］采用帶不平衡量的柔性靜子－柔性單圓盤(pán)轉(zhuǎn)子模型對(duì)同步全周碰摩現(xiàn)象進(jìn)行了分析研究。當(dāng)一次故障只使轉(zhuǎn)子的部分弧段與碰摩點(diǎn)發(fā)生接觸時(shí)，則為局部碰摩。Pennacchi等［5］分別利用試驗(yàn)?zāi)B(tài)分析和有限元建模方法對(duì)轉(zhuǎn)子的早期局部碰摩現(xiàn)象進(jìn)行了對(duì)比研究;王正浩等［6］分析了2圓盤(pán)、8自由度彈性轉(zhuǎn)子局部碰摩時(shí)的分叉與混沌行為;Cheng等［7］提出了一種基于經(jīng)驗(yàn)?zāi)Ｊ椒纸?Empirical Mode Decomposition，EMD)的轉(zhuǎn)子局部碰摩故障診斷方法。從文獻(xiàn)的梳理和對(duì)比分析中可以發(fā)現(xiàn)，多數(shù)轉(zhuǎn)子碰摩的研究成果集中在對(duì)轉(zhuǎn)子系統(tǒng)的數(shù)學(xué)建模和振動(dòng)響應(yīng)分析兩方面，且只針對(duì)全周碰摩和局部碰摩兩者中的一種進(jìn)行單獨(dú)研究，而將兩種碰摩故障進(jìn)行對(duì)比分析的研究工作很少。由于兩種類(lèi)型碰摩故障的振動(dòng)信號(hào)極為相似，因此，如何從振動(dòng)信號(hào)入手，正確識(shí)別碰摩故障的類(lèi)型，對(duì)于確定碰摩位置，做出正確診斷決策至關(guān)重要。

Muszynska［1］指出，碰摩發(fā)生時(shí)轉(zhuǎn)、靜子之間同時(shí)存在摩擦、沖擊、扭轉(zhuǎn)和剛度耦合等四種效應(yīng)。當(dāng)發(fā)生全周碰摩時(shí)，摩擦效應(yīng)占主導(dǎo)地位;而發(fā)生局部碰摩時(shí)，沖擊效應(yīng)占主導(dǎo)地位。因此，不同類(lèi)型的碰摩故障，其能量分布也是不同的。本文通過(guò)EMD方法與奇異值分解(Singular Value Decomposition，SVD)方法［8－9］的有效結(jié)合，提取出原信號(hào)能量分布的特征值，從而獲得表征原信號(hào)能量分布規(guī)律的特征參量。然后將經(jīng)特征參量訓(xùn)練過(guò)的支持向量機(jī)［10］(Support Vector Machine，SVM)用于對(duì)不同類(lèi)型的碰摩故障進(jìn)行分類(lèi)。對(duì)轉(zhuǎn)子試驗(yàn)臺(tái)上的碰摩仿真試驗(yàn)的應(yīng)用結(jié)果表明，本方法能夠有效區(qū)分全周碰摩與局部碰摩故障，且分類(lèi)準(zhǔn)確率高。

1 碰摩信號(hào)的EMD－SVD算法

碰摩信號(hào)特征提取算法首先對(duì)原信號(hào)進(jìn)行EMD，得到一組固有模式函數(shù)(Intrinsic Mode Function，IMF)，然后提取前四階IMF組成特征矩陣，最后利用SVD獲得原信號(hào)的特征值，算法流程如圖1所示。

圖1 碰摩信號(hào)特征提取算法Fig.1 Feature extraction algorithm for rub signals

1.1 EMD基本原理

在對(duì)原信號(hào)進(jìn)行EMD時(shí)，首先需要搜索原始信號(hào)x(t)的所有極值點(diǎn)，然后應(yīng)用三次樣條曲線對(duì)所有極大值點(diǎn)和極小值點(diǎn)分別進(jìn)行插值，獲得x(t)的上、下包絡(luò)線emax(t)和emin(t)。通過(guò)計(jì)算上、下包絡(luò)線的均值：

檢查d(t)是否滿足IMF條件：①對(duì)于某段時(shí)間序列，極值點(diǎn)和過(guò)零點(diǎn)數(shù)目必須相等或至多相差一點(diǎn);②在任意點(diǎn)處由局部極大點(diǎn)和局部極小點(diǎn)構(gòu)成的包絡(luò)平均值為零。若d(t)滿足以上兩條件，則可判定d(t)為IMF;否則，需要重復(fù)計(jì)算包絡(luò)線均值m(t)并提取信號(hào)細(xì)節(jié)成分d(t)，直至d(t)滿足IMF條件。這時(shí)稱d1(t)(d1(t)=d(t))為x(t)的一階IMF。

進(jìn)一步將d1(t)從x(t)中分離出去，可得余量：

以h1(t)作為原始數(shù)據(jù)，重復(fù)執(zhí)行上述步驟，直至分離出二階IMF。以此類(lèi)推，將x(t)分解為從高頻到低頻的若干階IMF和趨勢(shì)項(xiàng)：

式中：dk(t)為第k階IMF，mK(t)為趨勢(shì)項(xiàng)，且dk(t)和mK(t)均為行向量。應(yīng)用EMD后，全周和局部碰摩信號(hào)的前四階IMF如圖2(a)、圖2(b)所示。

圖2 碰摩信號(hào)1～4階IMFFig.2 No.1～4 IMF of rub signal

1.2 SVD基本原理

利用EMD所得的前四階IMF，組成關(guān)于原信號(hào)的特征矩陣：

其中“［ ］T”為矩陣轉(zhuǎn)置符號(hào)。

對(duì)于特征矩陣A∈Rm×n(m和n分別表示矩陣的行數(shù)和列數(shù))，存在正交矩陣 U∈Rm×m和 V∈Rn×n

使得：

且∑1=diag(σ1，σ2，…，σr)，其對(duì)角元素 σi(i=1，2，…，r)為矩陣A的非零奇異值，排列順序?yàn)棣?≥σ2≥…≥σr＞0，r=rank(A);U和V為酉矩陣，U和V的每列分別稱為左、右奇異向量。σi的大小表示奇異向量在特征矩陣中貢獻(xiàn)的大小。對(duì)于原信號(hào)的特征矩陣A，其特征值σi越大，表示相應(yīng)奇異向量包含的關(guān)于原信號(hào)的信息越多;反之，則越少。∑1具有良好的穩(wěn)定性，對(duì)于原信號(hào)中噪聲的干擾免疫力強(qiáng)。因此，可將其作為原信號(hào)特征信息的重要參量，用作識(shí)別碰摩類(lèi)型所用 SVM 的訓(xùn)練樣本［8，10］。

2 碰摩類(lèi)型的SVM識(shí)別算法

要應(yīng)用SVM對(duì)轉(zhuǎn)子碰摩類(lèi)型進(jìn)行識(shí)別，必須利用樣本數(shù)據(jù)進(jìn)行訓(xùn)練。

給定訓(xùn)練樣本集合 D={(xl，yl)，l=1，…，L}(其中輸入樣本 xl∈Rd，類(lèi)別標(biāo)識(shí) yl∈{+1，－1}，“+1”表示全周碰摩，“－1”表示局部碰摩)，對(duì)于兩類(lèi)線性可分的訓(xùn)練樣本，存在一組超平面：

能將兩類(lèi)樣本點(diǎn)完全分離。其中w為權(quán)重向量，b為偏置量。在這組超平面中，通過(guò)定義最優(yōu)分類(lèi)超平面，使得目標(biāo)函數(shù)最小化：

在大多數(shù)情況下，滿足不等式(10)的最優(yōu)超平面是不存在的。因此，引入松弛變量ξl＞0，則約束條件(10)變?yōu)椋?/p>

平衡因子C的值由用戶定義，其作用是使得最近樣本點(diǎn)到超平面的距離最大化以及分類(lèi)誤差最小化。

上述最優(yōu)化問(wèn)題式(11)～(12)的解變?yōu)椋?/p>

其中：αl為拉格朗日乘子，xa和xb為對(duì)應(yīng)于αl≠0的兩類(lèi)碰摩樣本數(shù)據(jù)中存在的支持向量，且有αa，αb＞0，ya=－1，yb=1。最后得最優(yōu)分類(lèi)器為：

實(shí)際碰摩數(shù)據(jù)的分類(lèi)問(wèn)題往往屬于非線性分類(lèi)的范疇。因此，進(jìn)一步引入核函數(shù)的概念，將樣本空間映射到特征空間，然后在特征空間中進(jìn)行線性分類(lèi)，這樣處理后所得最優(yōu)分類(lèi)器為：

其中：SVs表示支持向量集合，K(x，x')為核函數(shù)。常用的核函數(shù)有線性核函數(shù)、多項(xiàng)式核函數(shù)和徑向基核函數(shù)等。

3 碰摩分類(lèi)試驗(yàn)分析

為了驗(yàn)證EMD－SVD與SVM相結(jié)合的碰摩故障識(shí)別方法對(duì)不同類(lèi)型碰摩故障識(shí)別的有效性，本節(jié)利用Bently RK4轉(zhuǎn)子試驗(yàn)臺(tái)分別模擬轉(zhuǎn)子的全周碰摩與局部碰摩故障，并使用訓(xùn)練完的SVM對(duì)碰摩故障進(jìn)行分類(lèi)測(cè)試。

3.1 試驗(yàn)臺(tái)設(shè)置及參數(shù)選擇

Bently RK4轉(zhuǎn)子試驗(yàn)臺(tái)如圖3(a)所示。轉(zhuǎn)子上安置左、右兩個(gè)滑動(dòng)軸承，轉(zhuǎn)子中部安置調(diào)心框架一個(gè)，如圖3(b)所示。調(diào)心框架左右約15mm處各固定一質(zhì)量圓盤(pán)，用于保障高速運(yùn)轉(zhuǎn)條件下轉(zhuǎn)子的穩(wěn)定運(yùn)行。調(diào)心框架上、左、右側(cè)各安裝一根彈簧，用于調(diào)節(jié)轉(zhuǎn)子的軸心位置。此外，通過(guò)調(diào)節(jié)三根彈簧的拉力，使轉(zhuǎn)子軸心與右軸瓦中心重合。

圖3 Bently RK4轉(zhuǎn)子試驗(yàn)臺(tái)Fig.3 Bently RK4 test bench

全周碰摩通過(guò)圖3(c)的方式模擬。通過(guò)控制盒調(diào)節(jié)轉(zhuǎn)子運(yùn)行狀態(tài)，使轉(zhuǎn)子在1 800 r/min的轉(zhuǎn)速下平穩(wěn)運(yùn)行。將蠟棒通過(guò)人工按壓的方式在運(yùn)行轉(zhuǎn)子上加載10 s，通過(guò)安裝在轉(zhuǎn)子試驗(yàn)臺(tái)上的電渦流傳感器以1 280 Hz的采樣頻率采集振動(dòng)數(shù)據(jù)。從所采數(shù)據(jù)中每隔1 s選取時(shí)長(zhǎng)為0.4 s的一組樣本數(shù)據(jù)，共取50組樣本數(shù)據(jù)。由于蠟棒剛度遠(yuǎn)低于轉(zhuǎn)子系統(tǒng)，可認(rèn)為轉(zhuǎn)子系統(tǒng)發(fā)生了全周碰摩故障;局部碰摩裝置如圖3(d)所示，它由固定在支座上的碰摩支架與碰摩銅棒構(gòu)成。轉(zhuǎn)子啟動(dòng)前，銅棒與轉(zhuǎn)子分離。轉(zhuǎn)子啟動(dòng)后，將銅棒緩慢旋入支架螺紋孔直到與轉(zhuǎn)子接觸并保持一定的預(yù)緊力，然后將螺母旋緊，固定銅棒與支架的相對(duì)位置。由于銅棒的剛度與轉(zhuǎn)子系統(tǒng)剛度接近，可認(rèn)為這時(shí)轉(zhuǎn)子系統(tǒng)發(fā)生局部碰摩。最后用相同方式采集局部碰摩樣本數(shù)據(jù)50組。

3.2 試驗(yàn)結(jié)果分析

從50組全周碰摩樣本數(shù)據(jù)和50組局部碰摩樣本數(shù)據(jù)中各取一組進(jìn)行EMD并畫(huà)Hilbert時(shí)頻譜，如圖4所示。圖4(a)和圖4(b)分別為轉(zhuǎn)子全周碰摩與局部碰摩振動(dòng)信號(hào)的Hilbert時(shí)頻譜，可以看到兩類(lèi)碰摩故障的時(shí)頻特征十分相似，很難通過(guò)傳統(tǒng)的時(shí)、頻域及時(shí)頻分析方法將它們區(qū)分開(kāi)。因此，接下來(lái)應(yīng)用EMDSVD與SVM相結(jié)合的方法對(duì)兩種碰摩故障數(shù)據(jù)進(jìn)行識(shí)別。

圖4 碰摩信號(hào)Hilbert時(shí)頻譜Fig.4 Hilbert spectrogram of rub signal

首先，通過(guò)隨機(jī)排序的方法，將50組全周碰摩數(shù)據(jù)和50組局部碰摩數(shù)據(jù)混合成100組轉(zhuǎn)子碰摩數(shù)據(jù)。然后，任取其中50組，通過(guò)EMD－SVD對(duì)其進(jìn)行故障特征提取，將所得特征值輸入SVM進(jìn)行訓(xùn)練。試驗(yàn)所用SVM采用線性、多項(xiàng)式和徑向基三種核函數(shù)分別對(duì)訓(xùn)練樣本進(jìn)行學(xué)習(xí)。每種核函數(shù)中所用參數(shù)的選擇依據(jù)參考文獻(xiàn)［11］。取其中最優(yōu)結(jié)果列于表1第二列。將訓(xùn)練完的SVM用于剩余50組碰摩數(shù)據(jù)的分類(lèi)測(cè)試，并以式(17)作為衡量分類(lèi)準(zhǔn)確率的指標(biāo)。

其中：err表示分類(lèi)錯(cuò)誤率，PY為碰摩數(shù)據(jù)被SVM分類(lèi)后與實(shí)際所屬碰摩類(lèi)別不符的樣本數(shù)據(jù)個(gè)數(shù)。

為了保證分類(lèi)結(jié)果穩(wěn)定性，本試驗(yàn)對(duì)100組混合后的碰摩數(shù)據(jù)再進(jìn)行4次50組訓(xùn)練樣本的隨機(jī)提取，將5次試驗(yàn)所得分類(lèi)準(zhǔn)確率求平均值，其結(jié)果見(jiàn)表1第三列。由表可知，線性函數(shù)的分類(lèi)準(zhǔn)確率為86.3%，分類(lèi)效果一般。而多項(xiàng)式函數(shù)和徑向基函數(shù)的分類(lèi)準(zhǔn)確率都在94%以上，對(duì)全周碰摩與局部碰摩的分類(lèi)效果十分理想。其中，分類(lèi)結(jié)果最為準(zhǔn)確的是以C=80，σ=10的徑向基函數(shù)為核函數(shù)的SVM，分類(lèi)準(zhǔn)確率達(dá)到96.0%。

表1 不同核函數(shù)的分類(lèi)效果Tab.1 Classification efficiency of different Kernel functions

4 結(jié)論

轉(zhuǎn)子全周碰摩與局部碰摩信號(hào)的故障特征相似，用傳統(tǒng)的時(shí)、頻域特征提取方法很難將它們區(qū)分開(kāi)。本文提出了一種基于EMD－SVD與SVM的故障識(shí)別方法，即：首先，從采集的碰摩信號(hào)中選取一部分信號(hào)，利用EMD將碰摩信號(hào)分解為IMF函數(shù)組;然后，提取表征信號(hào)主要能量分布的前4階IMF組成信號(hào)特征矩陣，進(jìn)而對(duì)該矩陣進(jìn)行SVD，得到原信號(hào)的特征值并用其訓(xùn)練SVM;最后，將訓(xùn)練完的SVM用于剩余碰摩信號(hào)的分類(lèi)試驗(yàn)。基于轉(zhuǎn)子試驗(yàn)臺(tái)全周碰摩與局部碰摩的試驗(yàn)結(jié)果表明，EMD－SVD與SVM的碰摩故障識(shí)別方法是有效的，對(duì)于兩種不同類(lèi)型碰摩故障的分類(lèi)準(zhǔn)確率高。其中，以徑向基函數(shù)為核函數(shù)的SVM分類(lèi)準(zhǔn)確率達(dá)到96.0%。

需要指出的是，由于轉(zhuǎn)子系統(tǒng)在不同工況下的剛度存在較大變化，因此，分別利用蠟棒和銅棒模擬全周碰摩與局部碰摩與實(shí)際情況存在一定偏差，今后需要專(zhuān)門(mén)建立一套精確模擬全周碰摩與局部碰摩的試驗(yàn)裝置;此外，超臨界汽輪發(fā)電機(jī)組結(jié)構(gòu)和工況都較為復(fù)雜，對(duì)機(jī)組碰摩故障的機(jī)理進(jìn)行分析對(duì)于全周碰摩與局部碰摩故障的區(qū)分意義重大。

［1］Musznska A.Rotor－to－stationary element rub－related vibration phenomena in rotating machinery：Literature survey［J］.The Shock and Vibration Digest，1989，21(3)：3 － 11.

［2］ Bently D E，Yu J J，Goldman P，et al.Full annular rub in mechanical seals， Part Ⅰ：Experimental results［J］.International Journal of Rotating Machinery，2002，8(5)：329－336.

［3］ Choi Y S.Dynamics of rotor rub in annular clearance with experimental evaluation［J］.Journal of Mechanical Science and Technology，1994，8(4)：404 －413.

［4］徐尉南，張 文，許 斌.柔性轉(zhuǎn)子－柔性靜子系統(tǒng)的同步全周碰摩分析［J］. 振動(dòng)與沖擊，2007，26(2)：1 －5，9.

［5］Pennacchi P，Bachschmid N，Tanzi E.Light and short arc rubs in rotating machines：Experimental tests and modelling［J］.Mechanical Systems and Signal Processing，2009，23(7)：2205－2227.

［6］王正浩，袁惠群，范改燕.考慮軸初彎曲時(shí)轉(zhuǎn)子系統(tǒng)局部碰摩的分叉與混沌行為［J］.機(jī)械設(shè)計(jì)，2008，25(2)：37－41.

［7］ Cheng J S，Yu D J，Tang J S，et al.Local rub－impact fault diagnosis of the rotor systems based on EMD［J］.Mechanism and Machine Theory，2009，44(4)：784 －791.

［8］ Vanlanduit S，Parloo E，Cauberghe B，et al.A robust singular value decomposition for damage detection under changing operating conditions and structural uncertainties［J］.Journal of Sound and Vibration，2005，284(3 －5)：1033 －1050.

［9］段向陽(yáng)，王永生，蘇永生.基于奇異值分解的信號(hào)特征提取方法研究［J］. 振動(dòng)與沖擊，2009，28(11)：30－33.

［10］ Gunn S R.Support vector Machines for classification and regression［R］.Southampton：Image Speech and Intelligent Systems Research Group，University of Southampton，1997.

［11］ Ai L M，Wang J，Yao R X.Classification of parkinsonian and essential tremor using empirical mode decomposition and support vector machine［J］.Digitial Signal Process，2011，21(4)：543－550.

［12］李錄平，盧緒祥.汽輪發(fā)電機(jī)組振動(dòng)與處理［M］.北京：中國(guó)電力出版社，2007.

Recognition method for rotor system under full annular rub and partial rub

XIONG Xin，YANG Shi－xi，GAN Chun－biao，ZHOU Xiao－feng
(Department of Mechanical Engineering＆The State Key Lab of Liquid Power Transmission and Control，Zhejiang University，Hangzhou 310027，China)

Supercritical steam turbosets are highly complex in their structure and always run under various complicated working conditions.They are prone to rub between static and dynamic parts.Induced by different factors，rubbing faults can be divided into full annular rub and partial rub.Since the time－frequency characteristics of both are similar to each other，discrimination of the two kinds of rub faults is hard to proceed by using the traditional spectrum analysis methods.In response to make up for this shortage，an intelligent recognition method based on the EMD－SVD and SVM was proposed.IMFs were collected through EMD and the first－four－order IMFs，which contain the main power of the original signal，were extracted to form the characteristic matrix.The SVD was applied to obtain a series of eigenvalues，which were then inputted to train the SVM in order to classify rub faults.The newly developed intelligent recogonition method was used to analyze the signals collected from rotor test－bed under both full annular rub and partial rub conditions.The experiment results show that，classification accuracy of the method is high，especially for the SVM using radial basis as kernel function，where the classification accuracy is up to 96.0%.

rotor;full annular rub;partial rub;empirical mode decomposition(EMD);singular value decomposition(SVD);support vector machine(SVM)

TN911.7;TH165.3

A

國(guó)家自然科學(xué)基金(11072214);國(guó)家“863”高技術(shù)研究發(fā)展計(jì)劃(2008AA04Z410)

2011－03－25 修改稿收到日期：2011－07－26

熊 炘 男，博士生，1983年11月生