馮馳，胡楊，王兆豐

1.哈爾濱工程大學(xué)信息與通信工程學(xué)院，黑龍江哈爾濱150001 2.西安航空發(fā)動(dòng)機(jī)(集團(tuán))有限公司，陜西西安710021

基于分形理論的渦輪葉片特征提取

馮馳1，胡楊1，王兆豐2

1.哈爾濱工程大學(xué)信息與通信工程學(xué)院，黑龍江哈爾濱150001 2.西安航空發(fā)動(dòng)機(jī)(集團(tuán))有限公司，陜西西安710021

采用了葉片溫度這一葉片質(zhì)量的重要指標(biāo)，對(duì)數(shù)據(jù)預(yù)處理并進(jìn)行特征提取，為渦輪葉片建立起特征模型?；诜中卫碚撎崛∪~片溫度信號(hào)的3種分形維數(shù)特征，結(jié)合K-means聚類分析和ReliefF算法計(jì)算各特征值的權(quán)重，從而建立起渦輪葉片的溫度特征模型，實(shí)現(xiàn)對(duì)渦輪葉片故障的早期預(yù)警。統(tǒng)計(jì)結(jié)果表明，該特征模型能夠較好地反映出處于故障狀態(tài)的渦輪葉片的狀態(tài)。

渦輪葉片;特征提取;分形理論; K-means; ReliefF

網(wǎng)絡(luò)出版地址: http://www.cnki.net/kcms/detail/23.1191.U.20150727.0857.001.html

燃?xì)廨啓C(jī)作為國(guó)家高技術(shù)水平和科技實(shí)力的重要標(biāo)志之一被廣泛應(yīng)用于航空國(guó)防和工業(yè)發(fā)電領(lǐng)域。在它的工作過(guò)程中，發(fā)動(dòng)機(jī)噴出高速噴氣流，同時(shí)帶動(dòng)壓氣機(jī)和渦輪繼續(xù)旋轉(zhuǎn)，渦輪葉片長(zhǎng)期工作在高溫高壓的惡劣環(huán)境下，其自身溫度過(guò)熱使壽命受到極大影響，因此對(duì)渦輪葉片工作狀態(tài)的監(jiān)測(cè)具有十分重要的意義。分形理論最基本特點(diǎn)是用分形的數(shù)學(xué)工具來(lái)描述研究客觀事物。它跳出了一維的線、二維的面、三維的立體乃至四維時(shí)空的傳統(tǒng)束縛，更加趨近復(fù)雜系統(tǒng)的真實(shí)屬性與狀態(tài)的描述，更加符合客觀事物的多樣性與復(fù)雜性［1］。綜上所述，文中將分形理論應(yīng)用于渦輪葉片模型的建立，作為葉片故障診斷的理論依據(jù)。文中對(duì)渦輪葉片旋轉(zhuǎn)一個(gè)周期的溫度信號(hào)進(jìn)行預(yù)處理，以得到單個(gè)葉片的溫度信號(hào)，提取3種分形維數(shù)組成的渦輪葉片特征向量，同時(shí)引入K-means聚類分析方法對(duì)葉片進(jìn)行分類，根據(jù)分類結(jié)果應(yīng)用ReliefF算法計(jì)算各分形維數(shù)所占權(quán)重，繼而形成每個(gè)葉片的特征模型。渦輪葉片輻射測(cè)溫系統(tǒng)的輸出信號(hào)為電壓信號(hào)，電壓經(jīng)擬合后變換為溫度，為方便起見(jiàn)，文中均直接采用電壓表示葉片溫度，并未進(jìn)行電壓溫度轉(zhuǎn)換。

1　單葉片溫度數(shù)據(jù)分割

文中所采集的每組數(shù)據(jù)包括多個(gè)渦輪旋轉(zhuǎn)周期，短時(shí)間內(nèi)(100 ms)，每個(gè)葉片的溫度分布幾乎不變，因此可將多個(gè)周期的數(shù)據(jù)對(duì)齊取平均以減少噪聲帶來(lái)的影響。然而因數(shù)據(jù)的每個(gè)渦輪周期之間不可避免地存在相位偏差，為提高對(duì)齊精度，需先對(duì)數(shù)據(jù)進(jìn)行10倍插值。文中采用的三次樣條插值法具有良好的收斂性與穩(wěn)定性，又有二階光滑性。使用三次樣條法插值前后溫度數(shù)據(jù)局部對(duì)比如圖1所示。從圖中可以看出，插值后數(shù)據(jù)沒(méi)有出現(xiàn)畸變，而是變得更加平滑。

圖1　插值前后波形對(duì)比

以每周期溫度數(shù)據(jù)的最大值為基準(zhǔn)分割出17個(gè)周期，計(jì)算余下每周期與第一周期數(shù)據(jù)的相關(guān)性，并按照相關(guān)性最強(qiáng)點(diǎn)調(diào)整對(duì)齊。圖2為調(diào)整前后對(duì)比圖。

圖2　對(duì)齊前后波形對(duì)比

對(duì)齊之后取平均，對(duì)這一周期的平均溫度信號(hào)按極大值分割法(以每個(gè)葉片波形的最大的極大值點(diǎn)為基準(zhǔn)逐一分割)分出每個(gè)渦輪葉片的波形，如圖3所示。

圖3　葉片波形分割圖

文中算法流程如圖4所示。

圖4　特征提取算法流程

2　特征值計(jì)算

2.1分形維數(shù)

分形維數(shù)［2］是分形意義上由標(biāo)度關(guān)系得出的一個(gè)定量數(shù)值，可以定量地描述分形集的不規(guī)則程度和復(fù)雜程度。

2.1.1分形的標(biāo)度不變性

使用分形理論的前提是需分析的信號(hào)具有分形性質(zhì)。事實(shí)上，滿足一定條件的動(dòng)力學(xué)系統(tǒng)都會(huì)產(chǎn)生分形;反之，判別任意形體是否分形，沒(méi)有必要從這個(gè)形體的產(chǎn)生機(jī)制上入手，只要判斷它是否具有標(biāo)度不變性就可以了。

所謂標(biāo)度不變性是指無(wú)論測(cè)量尺度如何改變，所測(cè)量對(duì)象的特性(如形態(tài)特性、復(fù)雜程度、不規(guī)則性、統(tǒng)計(jì)特性等)均不發(fā)生變化。除了嚴(yán)格的數(shù)學(xué)模型以外，如Koch曲線，對(duì)于實(shí)際的分形集來(lái)說(shuō)，這種標(biāo)度不變性只在一定的范圍內(nèi)適用。通常把標(biāo)度不變性適用的空間稱之為該分形體的無(wú)標(biāo)度區(qū)間。

設(shè)平面R2內(nèi)有圖形F，在平面內(nèi)作間距為δ的方格網(wǎng)，則F與方格網(wǎng)相交的方格數(shù)Nδ(F )稱為圖形F在標(biāo)度(分辨率)δ下的盒數(shù)。如果存在一個(gè)標(biāo)度范圍(δ1，δ2)，在此范圍內(nèi)log Nδ(F )與log δ保持大約恒定的斜率，則這個(gè)范圍就稱為無(wú)標(biāo)度區(qū)。在無(wú)標(biāo)度區(qū)內(nèi)，可以將該形體視為分形。即雙對(duì)數(shù)關(guān)系曲線中存在一段直線(或近似直線)，則該段就是無(wú)標(biāo)度區(qū)。在該段直線的標(biāo)度范圍(δ1，δ2)內(nèi)，可以認(rèn)為波形F是分形，且這段直線的平均斜率就是分形的盒維數(shù)，即

2.1.2分形特性驗(yàn)證實(shí)驗(yàn)

基于上述理論，文中將標(biāo)度定位于(1，80)實(shí)驗(yàn)驗(yàn)證渦輪葉片溫度曲線是否具有分形特性。從實(shí)驗(yàn)所得的log Nδ(F )－log δ雙對(duì)數(shù)曲線來(lái)看不易準(zhǔn)確判斷定位是否為直線段，為此作出斜率與標(biāo)度的關(guān)系曲線，如圖5所示。

圖5 logNδ(F )－logδ雙對(duì)數(shù)曲線及曲線斜率

圖中可看出，δ在(1，15)的范圍內(nèi)，斜率值基本保持不變，故可認(rèn)為無(wú)標(biāo)度區(qū)為(1，15)，此處無(wú)需精確無(wú)標(biāo)度區(qū)范圍，只需證明其存在無(wú)標(biāo)度區(qū)即可。綜上，渦輪葉片的溫度信號(hào)具有分形特性，事實(shí)上，根據(jù)分形的定義，文中提取的葉片溫度信號(hào)的分形特征所反映出的是葉片在分形算法設(shè)定區(qū)域內(nèi)的波形分布特征。

2.1.3分形特征值計(jì)算

在設(shè)備的故障診斷方面，單重分形維數(shù)(例如盒維數(shù))只能從單一測(cè)度出發(fā)描述設(shè)備狀態(tài)信號(hào)的分形特征，并不能全面地反映其完整特性。因此文中選用多重分形，計(jì)算渦輪葉片溫度信號(hào)的分形特征，取權(quán)重因子q=0、1、2，所對(duì)應(yīng)的分形維數(shù)分別為容量維數(shù)(盒維數(shù))、信息維數(shù)、關(guān)聯(lián)維數(shù)［3］。部分葉片的3種分形維數(shù)計(jì)算結(jié)果如表1所示。

表1　第1～5號(hào)渦輪葉片的3種分形維數(shù)

3　分配特征權(quán)重

當(dāng)一個(gè)特征能夠清晰刻畫形體的時(shí)候，類似形體的特征應(yīng)是相近的，而非類似的形體特征相距較遠(yuǎn)。因此，具有相似特性的葉片的特征向量之間距離較近，反之較遠(yuǎn)。而好的特征應(yīng)該使屬于同一類的葉片樣本接近，使不同類的葉片樣本之間遠(yuǎn)離。鑒于這種思想，文中使用最常用的K-means聚類算法將葉片分為N類，根據(jù)結(jié)果使用ReliefF算法計(jì)算各分形維數(shù)所占權(quán)值。

3.1改進(jìn)的K-means聚類算法

聚類算法是給予數(shù)據(jù)自然上的相似劃分，要求得到的結(jié)果是每個(gè)類別內(nèi)部數(shù)據(jù)盡可能的相似而類別之間要盡可能大的存在差異。K-means算法是常用的基于劃分的聚類算法，首先隨機(jī)選擇k個(gè)對(duì)象作為初始的k個(gè)簇的質(zhì)心，然后將其余對(duì)象根據(jù)其與各個(gè)簇的質(zhì)心的距離分配到最近的簇;最后重新計(jì)算各個(gè)簇的質(zhì)心，不斷重復(fù)此過(guò)程，直到目標(biāo)函數(shù)最小為止［4］。

3.1.1改進(jìn)方法

K-means算法存在以下缺陷［5］:

1)聚類個(gè)數(shù)k需要預(yù)先給定，而k值選定難估計(jì);

2)對(duì)初始聚類中心的選取具有依賴性，算法常常因此陷入局部最優(yōu)解。

基于上述缺陷，文中采用了改進(jìn)的K-means算法:

1)最佳聚類數(shù)的選擇。根據(jù)DB Index準(zhǔn)則，DBk值越小，說(shuō)明聚類的效果越好，一般情況下，最佳聚類個(gè)數(shù)不會(huì)超過(guò)個(gè)，因而迭代算法可以在2～之間進(jìn)行，在文中選擇2～10，選擇DB指標(biāo)最小的聚類數(shù)為最佳聚類數(shù)。

2)初始聚類中心的選擇。普通K-means聚類方法使用隨機(jī)選擇的聚類中心，這使得算法極易陷入局部最優(yōu)解。文中通過(guò)選舉機(jī)制產(chǎn)生初始聚類中心。具體做法是首先計(jì)算最大最小樣本點(diǎn)之間的距離d。指定常整數(shù)p，k≤p≤N。在N個(gè)模式中，找到p個(gè)相距較遠(yuǎn)的樣本點(diǎn)。再對(duì)p個(gè)樣本點(diǎn)計(jì)算以本身為球心、以d(指定的常數(shù))為半徑的球內(nèi)所包含的樣本點(diǎn)個(gè)數(shù)。選舉個(gè)數(shù)較多的前k個(gè)樣本點(diǎn)作為初始聚類中心。

3.1.2聚類效果演示

經(jīng)計(jì)算，根據(jù)DB Index準(zhǔn)則，計(jì)算得最佳聚類數(shù)為9，圖6所示的聚類初始中心點(diǎn)與更新中心點(diǎn)對(duì)比圖，可見(jiàn)改進(jìn)的K-means聚類方法沒(méi)有使聚類結(jié)果陷入局部最優(yōu)解。

圖6　聚類初始中心點(diǎn)與更新中心點(diǎn)圖

最后的聚類結(jié)果分類如圖7所示，表2記錄了各分類所包含的渦輪葉片號(hào)，基于分類結(jié)果，文中采用ReliefF算法計(jì)算各分形維數(shù)對(duì)分類所做的貢獻(xiàn)值，以貢獻(xiàn)值大小作為該種分形維數(shù)的權(quán)重。

圖7　聚類結(jié)果分類

表2　各類別所包含的葉片代號(hào)(共86個(gè)葉片)

3.2基于聚類結(jié)果的ReliefF算法

ReliefF算法從訓(xùn)練集D中隨機(jī)選擇一個(gè)樣本R，然后從和R同類的樣本中尋找最近鄰的k個(gè)樣本Hii=1，2，…，k ()，稱為Near Hits，從和R不同類的樣本中尋找最近鄰的k個(gè)樣本Mi(C ) i=1，2，…，k ()，稱為Near Misses，然后根據(jù)式(1)更新某個(gè)特征A的權(quán)重W(A)。

式中: m為重復(fù)次數(shù)，diff(A，R，S )指的是樣本R和S在特征A上的差，p(C )指樣本落在集合C中的概率［6-8］。特征的權(quán)重越大，表示該特征的分類能力越強(qiáng)，反之，表示該特征分類能力越弱。經(jīng)計(jì)算，葉片的3種分形維數(shù)所對(duì)應(yīng)的權(quán)值如表3所示。

表3 3種分形維數(shù)所對(duì)應(yīng)的權(quán)值

將葉片各特征乘以其所對(duì)應(yīng)的權(quán)值后，再按改進(jìn)的K-means算法計(jì)算分類后的DB指標(biāo)，指標(biāo)從0.286 8降到0.155 7，可見(jiàn)ReliefF算法能有效地提高聚類性能，并且能定量分析各維特征對(duì)分類的貢獻(xiàn)程度，使葉片溫度信號(hào)特征更具有代表性。

3.3模擬故障信號(hào)及其識(shí)別實(shí)驗(yàn)

建立葉片特征模型的最終目的是為了能夠在葉片發(fā)生故障的初期進(jìn)行預(yù)警，提示操作人員停機(jī)檢查。文中主要研究葉片故障中2種情況:一類是葉片氣道大面積堵塞導(dǎo)致整體溫度過(guò)高。整體溫度過(guò)高又分為3種情況:前后氣道均堵塞，中部正常導(dǎo)致的前緣溫度過(guò)高，尾緣溫度也過(guò)高;葉片前緣氣道堵塞導(dǎo)致前緣溫度過(guò)高，使得尾緣氣流通過(guò)量增大導(dǎo)致尾緣溫度過(guò)低;前緣溫度過(guò)低，尾緣溫度過(guò)高。另一類是葉片局部毛細(xì)氣道堵塞或者隔熱涂層小面積脫落導(dǎo)致局部小范圍溫度過(guò)高。隔熱涂層的發(fā)射率大概在0.7左右，而金屬發(fā)射率在0.9左右，因此若葉片的隔熱涂層脫落會(huì)引起葉片局部的溫度偏高，通過(guò)斯特藩—玻爾茲曼定律，葉片工作溫度在650～750℃的時(shí)候，隔熱涂層脫落會(huì)使得葉片溫度上升40℃左右，轉(zhuǎn)換成電壓即上升0.14 V左右，為模擬涂層脫落故障的漸變過(guò)程，文中僅將電壓調(diào)整為上升0.8 V。整體溫度偏高的故障較為少見(jiàn)，常見(jiàn)故障以局部故障為主［9-10］。

文中使用數(shù)學(xué)軟件MATLAB R2014a版本進(jìn)行模擬計(jì)算，MATLAB具有高效的數(shù)值計(jì)算及符號(hào)計(jì)算功能，且提供了大量方便實(shí)用的處理工具箱。隨機(jī)選取正常的第11號(hào)葉片溫度數(shù)據(jù)，模擬6個(gè)模擬故障信號(hào)，其中1號(hào)和2號(hào)故障數(shù)據(jù)是通過(guò)將正常數(shù)據(jù)乘以Kaiser窗獲得的，其中一個(gè)Kaiser窗的衰減參數(shù)設(shè)為0.3，另一個(gè)設(shè)為0.4，來(lái)模擬葉片前緣、尾緣溫度都過(guò)高的情況; 3、4號(hào)故障數(shù)據(jù)是通過(guò)將正常葉片數(shù)據(jù)分別通過(guò)幅度為0.01和－0.01的正弦窗獲得的，用來(lái)模擬葉片前緣溫度高、尾緣溫度低和前緣溫度低、尾緣溫度高的情況; 5、6號(hào)故障數(shù)據(jù)是將正常數(shù)據(jù)加上一個(gè)以第50數(shù)據(jù)點(diǎn)和第100數(shù)據(jù)點(diǎn)為中心，寬度為50個(gè)數(shù)據(jù)點(diǎn)，高度為0.8的高斯窗，用來(lái)模擬第50數(shù)據(jù)點(diǎn)和第100數(shù)據(jù)點(diǎn)發(fā)生局部故障的情況。1、3、5、6號(hào)故障數(shù)據(jù)模擬前后對(duì)比如圖8所示，圖中實(shí)線為正常數(shù)據(jù)波形，而虛線則為故障模擬數(shù)據(jù)波形。

表4所示的是11號(hào)葉片溫度信號(hào)在未進(jìn)行求和取平均時(shí)在前10個(gè)周期的波形所對(duì)應(yīng)的乘權(quán)值后的分形特征值及因這10個(gè)周期11號(hào)葉片波形變化導(dǎo)致的分形特征值變化的范圍。由表4可以看出，3種分形維數(shù)的變化范圍均不超過(guò)2%，故將σ=－2%+ 2% []設(shè)為置信區(qū)間，用來(lái)判斷葉片是否處于正常狀態(tài)，此置信區(qū)間是根據(jù)文中所用數(shù)據(jù)設(shè)置的故障識(shí)別門限，因文中的數(shù)據(jù)量有限，置信區(qū)間與實(shí)際可能存在出入，具體應(yīng)用時(shí)可根據(jù)大量數(shù)據(jù)樣本對(duì)故障識(shí)別門限進(jìn)行修改。

表4 11號(hào)葉片10個(gè)周期的分形特征及其波動(dòng)范圍

與11號(hào)葉片類似，文中使用同樣的方式模擬出其他85個(gè)葉片對(duì)應(yīng)的6種故障數(shù)據(jù)，即每種故障均有86個(gè)模擬樣本。將各故障類型的分形特征與其原始正常波形的分形特征相比較，統(tǒng)計(jì)故障分形特征落在置信區(qū)間之外不同范圍的葉片個(gè)數(shù)，并計(jì)算識(shí)別率組成表5。設(shè)需判別的某模擬葉片的某分形特征為Dxi，其所對(duì)應(yīng)的正常葉片數(shù)據(jù)xi同一分形特征為dxi，則識(shí)別率η計(jì)算公式如下:

表5　各故障類型分形特征與正常相比偏離超過(guò)2%的葉片個(gè)數(shù)統(tǒng)計(jì)表

由表可以看出，模擬的各類故障波形與正常波形相比，均能在分形特征上體現(xiàn)出顯著差異，且在關(guān)聯(lián)維數(shù)上體現(xiàn)得尤其明顯，信息維數(shù)其次，容量維數(shù)最次，且當(dāng)增大故障變化的幅度時(shí)，識(shí)別率能夠達(dá)到更高。因此，分形特征能夠幫助監(jiān)測(cè)人員識(shí)別故障葉片以及預(yù)測(cè)發(fā)生何種故障。對(duì)于不同種故障，3種分形特征會(huì)體現(xiàn)出不同的識(shí)別情況，因模擬的故障類型有限，不能對(duì)各分形維數(shù)最適于識(shí)別何種故障類型進(jìn)行一一討論。值得一提的是，若使用將原正常波形平移或壓縮伸展的方式來(lái)模擬故障信號(hào)的話，提取出的故障信號(hào)的分形特征與正常信號(hào)相比，將不會(huì)發(fā)生變化，這是因?yàn)橛?jì)算分形維數(shù)時(shí)，算法會(huì)先將波形伸展壓縮在一個(gè)等邊長(zhǎng)的正方形區(qū)域里，這導(dǎo)致無(wú)論將原始波形做任何伸展平移，經(jīng)計(jì)算后都會(huì)產(chǎn)生同樣的結(jié)果，事實(shí)上，葉片發(fā)生故障時(shí)不會(huì)僅僅是正常波形的單純平移伸展，所以這種情況可以忽略。綜上，文中所研究的渦輪葉片分形特征能顯著地反映出葉片故障情況，可以較好地應(yīng)用于葉片的故障識(shí)別中。但故障識(shí)別門限以及故障識(shí)別類型需大量的具有代表性的正常葉片波形和故障葉片波形樣本，因研究條件有限故文中不做深入討論。

4　結(jié)束語(yǔ)

針對(duì)渦輪葉片溫度信號(hào)進(jìn)行了插值、對(duì)齊、分割等數(shù)據(jù)預(yù)處理，通過(guò)分形理論提取了溫度數(shù)據(jù)的3種分形特征，并引用經(jīng)改善過(guò)的K-means算法結(jié)果和應(yīng)用ReliefF算法，使得葉片特征按重要性進(jìn)行權(quán)值分配，最后進(jìn)行了常見(jiàn)葉片故障波形仿真。經(jīng)統(tǒng)計(jì)結(jié)果顯示，文中提取的特征值能夠有效地反映出葉片故障情況，可以應(yīng)用于葉片的故障識(shí)別中。但因本課題研究條件有限，不能深入探究故障識(shí)別的具體門限值，接下來(lái)的研究重點(diǎn)將放在基于大量具有代表性葉片波形樣本的故障識(shí)別方向上。

［1］趙玉春.基于混沌分形與模糊聚類的滾動(dòng)軸承故障診斷［D］.秦皇島:燕山大學(xué)，2011: 36-42.

［2］郝研.分形維數(shù)特性分析及故障診斷分形方法研究［D］.天津:天津大學(xué)，2012: 35-51.

［3］李兆飛.振動(dòng)故障分形特征提取及診斷方法研究［D］.重慶:重慶大學(xué)，2013: 36-48.

［4］王千，王成，馮振元，等.K-means聚類算法研究綜述［J］.電子設(shè)計(jì)工程，2012(7) : 21-24.

［5］袁方，周志勇，宋鑫.初始聚類中心優(yōu)化的K-means算法［J］.計(jì)算機(jī)工程，2007，33(3) : 65-66.

［6］張麗新，王家，趙雁南，等.基于Relief的組合式特征選擇［J］.復(fù)旦學(xué)報(bào)，2004，43(5) : 893-897.

［7］張勇.基于ReliefF算法的模糊聚類新算法［J］.應(yīng)用技術(shù)，2009(1) : 43-46.

［8］李曉嵐.基于Relief特征選擇算法的研究與應(yīng)用［D］.大連:大連理工大學(xué)，2013: 11-30.

［9］劉英乾.渦輪葉片故障診斷與模擬研究［D］.哈爾濱:哈爾濱工程大學(xué)，2013: 37-39.

［10］劉大響.航空發(fā)動(dòng)機(jī)葉片故障及預(yù)防研討會(huì)論文集［C］//北京:航空航天工業(yè)出版社，2005: 91-97.

Feature extraction of turbine blades based on the fractal theory

FENG Chi1，HU Yang1，WANG Zhaofeng2

1.College of Information and Communication Engineering，Harbin Engineering University，Harbin 150001，China 2.Xi’an Aero-engine (Group) Ltd.，Xi’an 710021，China

The data is preprocessed and features are extracted to establish a feature model of turbine blades by blade temperature，which is an important indicator of the turbine blade＇s quality.Based on the fractal theory，the three fractal dimension features of blade temperature signals are extracted，the weight of each feature is calculated by combining with K－means clustering analysis and ReliefF algorithm，and thereby establish the temperature feature model of turbine blades，achieving early warning of the turbine blade failure.Statistical results show that the feature model can effectively reflect the state of turbine blades with failure.

turbine blades; feature extraction; fractal theory; K-means; ReliefF

TK473

A

1009-671X(2015) 04-064-06

10.3969/j.issn.1009-671X.201411006

2014-11-14.網(wǎng)絡(luò)出版日期: 2015-07-27.

黑龍江省自然科學(xué)基金資助項(xiàng)目(F201413).

馮馳(1961-)，男，教授;胡楊(1990-)，女，碩士研究生.

胡楊，E-mail: huyang900218@163.com.