小波包能量譜-稀疏核主元在故障檢測(cè)中的應(yīng)用

朱丹，范玉剛，鄒金慧，吳建德，黃國(guó)勇

1.昆明理工大學(xué)信息工程與自動(dòng)化學(xué)院，昆明 650500

2.云南省礦物管道輸送工程技術(shù)研究中心，昆明 650500

小波包能量譜-稀疏核主元在故障檢測(cè)中的應(yīng)用

朱丹1，2，范玉剛1，2，鄒金慧1，2，吳建德1，2，黃國(guó)勇1，2

1.昆明理工大學(xué)信息工程與自動(dòng)化學(xué)院，昆明 650500

2.云南省礦物管道輸送工程技術(shù)研究中心，昆明 650500

針對(duì)滾動(dòng)軸承故障檢測(cè)的問(wèn)題，提出一種基于小波包能量譜-稀疏核主元的滾動(dòng)軸承故障檢測(cè)方法。對(duì)振動(dòng)信號(hào)進(jìn)行小波包分解，提取信號(hào)的能量頻譜，用增量式樣本基構(gòu)造方法，提取能量頻譜的樣本基，以此樣本基建立核主元模型，來(lái)分析軸承振動(dòng)信號(hào)能量頻譜的變化。通過(guò)實(shí)驗(yàn)仿真來(lái)說(shuō)明此算法的有效性。

滾動(dòng)軸承；小波包；稀疏核主元；故障檢測(cè)

1 引言

滾動(dòng)軸承廣泛應(yīng)用于各種機(jī)械設(shè)備，其運(yùn)行狀態(tài)對(duì)整臺(tái)機(jī)器的安全運(yùn)行影響極大，一旦滾動(dòng)軸承發(fā)生故障或者失效，將對(duì)機(jī)器運(yùn)行帶來(lái)安全隱患。因此，對(duì)滾動(dòng)軸承進(jìn)行監(jiān)測(cè)、診斷具有重要意義。滾動(dòng)軸承在不同的工況下運(yùn)行，產(chǎn)生不同的振動(dòng)信號(hào)，因此可以通過(guò)采集滾動(dòng)軸承工作時(shí)的振動(dòng)信號(hào)，并對(duì)該信號(hào)進(jìn)行分析和研究，來(lái)診斷滾動(dòng)軸承的運(yùn)行狀況。對(duì)振動(dòng)信號(hào)進(jìn)行小波包分解，得到反映信號(hào)特征的頻帶能量譜。當(dāng)滾動(dòng)軸承開(kāi)始發(fā)生故障時(shí)，信號(hào)的各頻帶能量中包含了相應(yīng)的故障信息，即振動(dòng)信號(hào)各頻帶的能量比例會(huì)發(fā)生變化，因而可利用滾動(dòng)軸承振動(dòng)信號(hào)各頻段的能量比例來(lái)分析軸承的運(yùn)行情況，以此作為故障診斷的依據(jù)。在對(duì)滾動(dòng)軸承進(jìn)行監(jiān)測(cè)和診斷的過(guò)程中，故障信號(hào)的分析處理是其中的核心內(nèi)容，它在提高診斷的準(zhǔn)確性和可靠性方面處于非常重要的位置。本文首先對(duì)振動(dòng)信號(hào)提取能量頻譜，并建立能量頻譜的樣本基，以此樣本基建立多元統(tǒng)計(jì)模型，用于分析軸承振動(dòng)信號(hào)能量頻譜的變化。

對(duì)得到反映信號(hào)特征的頻帶能量譜進(jìn)行分析，是故障檢測(cè)的重要內(nèi)容。主元分析法（PCA）是一種基于多元統(tǒng)計(jì)的故障檢測(cè)方法。然而在現(xiàn)代工業(yè)系統(tǒng)中都存在著非線性問(wèn)題，因此各種非線性PCA方法相繼被提出，例如，核主元分析方法（KPCA）[1]。KPCA是一種非線性特征提取方法，得到了廣泛的應(yīng)用[2-5]。然而，核主元分析存在不適合大樣本分析、計(jì)算復(fù)雜的缺點(diǎn)，因此本文提出一種稀疏核主元故障檢測(cè)方法，對(duì)振動(dòng)信號(hào)的小波包能量譜進(jìn)行分析。首先，將滾動(dòng)軸承振動(dòng)信號(hào)進(jìn)行小波包分解，計(jì)算各子頻帶信號(hào)的能量值，并構(gòu)建能量頻譜特征向量；其次，對(duì)得到的特征向量采用增量式樣本基構(gòu)造方法提取樣本基，以此樣本基建立正常的核主元模型，用于軸承振動(dòng)信號(hào)故障檢測(cè)；最后，通過(guò)仿真分析，說(shuō)明了此算法的有效性。

2 小波包的能量譜特征量提取

2.1 小波包分解

小波包分解是在多分辨分解的基礎(chǔ)上將各尺度下的細(xì)節(jié)分量作進(jìn)一步分解，從而實(shí)現(xiàn)對(duì)隨尺度變小而變寬的頻率窗口再劃分，提高信號(hào)高頻部分頻率的分辨率，使故障特征提取能夠在更加細(xì)化的頻帶內(nèi)進(jìn)行。

其中，t=1，2，…，N/2i（N為原始信號(hào)的采集點(diǎn)數(shù)）；i= 2j，j=1，2，…，J（J為最大分解層數(shù)），原始信號(hào)經(jīng)過(guò)j層小波分解得到2j個(gè)小波包，每個(gè)小波包的采樣點(diǎn)數(shù)為N/2j。

在滾動(dòng)軸承故障診斷中，將采集到的信號(hào)進(jìn)行小波包分解后，得到2j（j表示分解層數(shù)）個(gè)分解向量，它們包含著各個(gè)頻段的信息，而故障信息包含在其中的某些頻段內(nèi)。

2.2 能量譜特征量提取

按照能量方式表示的小波包分解結(jié)果稱為小波包能量譜，在小波包能量譜中，可以選取各頻帶內(nèi)信號(hào)的平方和作為能量的標(biāo)志，為進(jìn)一步的故障識(shí)別提供了診斷的前提條件。以3層小波包分解為例，說(shuō)明滾動(dòng)軸承故障信號(hào)小波包能量譜特征量提取方法[8-9]：

（1）首先對(duì)滾動(dòng)軸承振動(dòng)信號(hào)進(jìn)行3層小波包的能量譜特征量提取。

（4）構(gòu)造軸承故障信號(hào)的特征向量。由于系統(tǒng)出現(xiàn)故障時(shí)會(huì)對(duì)各頻帶內(nèi)信號(hào)能量有較大的影響，因此，以能量為元素可以構(gòu)造一個(gè)特征向量。

（5）特征向量T構(gòu)造如下：

3 核主元分析

KPCA的基本思想是首先通過(guò)非線性映射將原輸入空間映射到一個(gè)高維的特征空間，然后在高維的特征空間上進(jìn)行主元分析，從而把輸入空間中的非線性問(wèn)題轉(zhuǎn)化為特征空間中的線性問(wèn)題。

KPCA是一種非線性方法[10-12]，它首先通過(guò)非線性映射φ將原輸入空間（x1，x2，…，xN∈Rm)（其中N為采樣樣本的數(shù)目，m為測(cè)量變量的維數(shù)）映射到一個(gè)高維的特征空間F中（即：φ:Rm→F），然后在這個(gè)高維的特征空間F內(nèi)，進(jìn)行主元分析，從而把輸入空間中的非線性問(wèn)題轉(zhuǎn)化為特征空間中的線性問(wèn)題。把xi的映射記為φ(xi)=φi，則特征空間F的協(xié)方差矩陣可以表示為：

設(shè)矩陣CF所對(duì)應(yīng)的特征值為λ，特征向量為v，則存在系數(shù)αi(i=1，2，…，N)，使得矩陣CF的特征向量v可以用φi表示為：

在特征空間F上進(jìn)行主元分析之前，應(yīng)先作標(biāo)準(zhǔn)化處理，即用下式代替矩陣K：

其中，IN等于1/N與一個(gè)N×N的單位矩陣E∈RN×N相乘。所以，在特征空間中進(jìn)行主元分析，就等價(jià)于對(duì)式（10）求解特征值問(wèn)題。結(jié)合式（10）和式（6），由矩陣K的特征向量α可以求出矩陣CF的特征向量v，且滿足：

其中p為主元個(gè)數(shù)。這樣，可以通過(guò)計(jì)算映射數(shù)據(jù)在特征向量vk上的投影來(lái)計(jì)算主元，即

為了解決式（10）的特征值問(wèn)題，并利用式（13）直接從輸入空間計(jì)算特征空間的主元向量，在特征空間中引入點(diǎn)積形式的核函數(shù)k(x，y)=φ(x)，φ(y)，避免直接計(jì)算非線性映射。核函數(shù)的選擇完全決定映射φ和特征空間F，常用的核函數(shù)有：高斯核函數(shù)（也叫徑向基核函數(shù)）：k(x，y)=exp(-‖x-y‖2/σ)；多項(xiàng)式核函數(shù)：k(x，y)= x，yd；符號(hào)核函數(shù)：k(x，y)=tanh(β0x，y+β1)。其中，σ，d，β0和β1需提前指定。

然而，對(duì)于那些大型的、復(fù)雜的工業(yè)過(guò)程來(lái)說(shuō)，樣本數(shù)據(jù)一般會(huì)很多，這就導(dǎo)致核矩陣的計(jì)算復(fù)雜，因此KPCA不適合大樣本分析。針對(duì)這一不足，本文提出稀疏核主元的故障檢測(cè)方法。

4 基于小波包能量譜的稀疏核主元分析

基于小波包能量譜的稀疏核主元分析方法，將信號(hào)各頻率段的能量值構(gòu)造特征向量T，并提取T的樣本基，以此樣本基建立核主元分析模型。因此，該算法主要包括樣本基提取、KPCA模型建立兩大部分。

4.1 樣本基提取

樣本基是稀疏核主元算法的基礎(chǔ)，樣本基的構(gòu)造實(shí)質(zhì)是從輸入樣本集X中找到一個(gè)樣本子集，滿足條件：（1）由該樣本子集映射到特征空間的核矩陣滿秩；（2）該樣本子集是滿足條件（1）的最大樣本集。

該定理的用途是：當(dāng)Kn降秩時(shí)，第n個(gè)樣本對(duì)應(yīng)的核變換可用前n-1個(gè)樣本的核變換線性表示。

由此得到如下的樣本基構(gòu)造算法：

（1）令初始樣本基僅包含任意一個(gè)樣本，d=1，計(jì)算相應(yīng)的核矩陣。

（2）逐個(gè)檢驗(yàn)樣本，計(jì)算δ，如果δ=0，則該樣本不加入樣本基S；否則將加入樣本S，并用遞推算法修改相應(yīng)的核矩陣和核矩陣的逆，令d=d+1。

（3）算法結(jié)束并得到d和樣本基S={s1，s2，…，sd}。

（4）以此樣本基取代全部訓(xùn)練樣本，建立KPCA模型，稱為稀疏核主元。

4.2 基于小波包能量譜的稀疏核主元故障檢測(cè)

基于小波包能量譜的稀疏核主元故障檢測(cè)步驟如下：

（1）首先對(duì)滾動(dòng)軸承振動(dòng)信號(hào)進(jìn)行3層小波包分解，分解算法如式（1），得到第3層從低頻到高頻的8個(gè)頻帶成分。

（2）求各子頻帶信號(hào)的能量值，計(jì)算方法如式（3）。

（3）對(duì)能量頻譜構(gòu)成的特征向量T，采用增量式樣本基的構(gòu)造方法構(gòu)造樣本基，以此樣本基作為訓(xùn)練樣本，建立正常的核主元模型。

（4）將軸承振動(dòng)信號(hào)（本文選取軸承滾動(dòng)鐵和內(nèi)圈這兩種故障）重復(fù)步驟（1）、（2），得到振動(dòng)信號(hào)分解后各頻段的能量值，以此作為測(cè)試樣本。

（5）對(duì)測(cè)試樣本進(jìn)行核主元分析，求取T2和SPE統(tǒng)計(jì)量進(jìn)行故障檢測(cè)。

5 仿真實(shí)驗(yàn)研究

本文采用來(lái)自美國(guó)西儲(chǔ)大學(xué)（Case Western Reserve）滾動(dòng)軸承實(shí)驗(yàn)中心的實(shí)驗(yàn)數(shù)據(jù)作為故障檢測(cè)的依據(jù)，從中選取的數(shù)據(jù)是6205-2RSJEMSKF型的滾動(dòng)軸承數(shù)據(jù)，實(shí)驗(yàn)裝置和軸承尺寸見(jiàn)http：//www.eecs.cwru.edu/laboratory/bearing/。其加速度傳感器安裝在電機(jī)驅(qū)動(dòng)端，空載轉(zhuǎn)速為1 797 r/min（也就是轉(zhuǎn)頻為1 797/60 r/s= 29.95 Hz），采樣頻率為48 kHz，模擬軸承滾動(dòng)鐵和內(nèi)圈正常與故障的[故障直徑為0.014 inch（0.035 56 cm）]振動(dòng)信號(hào)。利用小波包能量譜變化的主元分析方法、標(biāo)準(zhǔn)核主元分析方法和稀疏核主元分析方法，分別對(duì)軸承數(shù)據(jù)進(jìn)行仿真分析，并比較所得結(jié)果，最終驗(yàn)證了稀疏核主元分析方法的有效性。

5.1 頻譜分析

選取上面所述條件下的滾動(dòng)軸承數(shù)據(jù)，外加一組在此條件下的滾動(dòng)鐵和內(nèi)圈的故障數(shù)據(jù)（故障直徑為0.007 inch）進(jìn)行頻譜分析，得到的頻譜圖如圖1～圖5所示。

圖1 正常軸承數(shù)據(jù)頻譜圖

圖2 滾動(dòng)鐵故障頻譜圖（故障直徑為0.014 inch）

圖3 滾動(dòng)鐵故障頻譜圖（故障直徑為0.007 inch）

圖4 內(nèi)圈故障頻譜圖（故障直徑為0.014 inch）

圖5 內(nèi)圈故障頻譜圖（故障直徑為0.007 inch）

由圖可看出，正常頻譜圖大概在0.4×104Hz附近有波峰，而軸承故障頻譜圖大概在0.3×104Hz附近有波峰，直接觀察故障特征不是很明顯，所以提出主元分析方法、核主元方法和稀疏核主元方法進(jìn)一步對(duì)軸承數(shù)據(jù)進(jìn)行故障檢測(cè)。

5.2 主元分析

滾動(dòng)軸承實(shí)驗(yàn)數(shù)據(jù)故障檢測(cè)的主要步驟如下：

（1）首先選取軸承滾動(dòng)鐵和內(nèi)圈正常與故障中的400組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行三層小波包變換，由于采樣頻率為48 kHz，故分析范圍為0～24 kHz，將每組信號(hào)分解成8個(gè)頻段，各節(jié)點(diǎn)頻率范圍分別為：0～3 kHz，3～6 kHz，6～ 9 kHz，9～12 kHz，12～15 kHz，15～18 kHz，18～21 kHz，21～ 24 kHz。計(jì)算各節(jié)點(diǎn)的小波系數(shù)，并求出其能量值。

（2）分別將其作為輸入變量進(jìn)行主元分析，求取主元變量的T2和SPE統(tǒng)計(jì)量[14]，檢測(cè)故障是否發(fā)生，如圖6～圖9所示。

圖6 滾動(dòng)鐵故障的T2統(tǒng)計(jì)量故障檢測(cè)圖（PCA）

圖7 滾動(dòng)鐵故障的SPE統(tǒng)計(jì)量故障檢測(cè)圖（PCA）

圖8 內(nèi)圈故障的T2統(tǒng)計(jì)量故障檢測(cè)圖（PCA）

本文采用的小波包是db3（Daubechies）正交小波包，經(jīng)3層小波包分解后，求取各組樣本中各頻帶的頻譜能量值。

圖9 內(nèi)圈故障的SPE統(tǒng)計(jì)量故障檢測(cè)圖（PCA）

由圖6～圖9可看出，數(shù)據(jù)在中間段以后有變化，但變化不是很明顯，檢測(cè)效果不好，因此提出下面的稀疏核主元故障檢測(cè)方法。

5.3 稀疏核主元分析

滾動(dòng)軸承故障檢測(cè)步驟同上，核函數(shù)選用普遍應(yīng)用的高斯核，即k(x，y)=exp(-‖x-y‖2/σ)，這里取σ=35。

經(jīng)3層小波包分解后，求取各組樣本中各頻帶的頻譜能量值。先將軸承正常400組數(shù)據(jù)訓(xùn)練，求得軸承運(yùn)行的核主元模型，故障控制限為95%。最后分別將滾動(dòng)鐵故障和內(nèi)圈故障400組測(cè)試數(shù)據(jù)處理后輸入核主元模型，求取T2和SPE統(tǒng)計(jì)量，如圖10～圖13所示。

利用4.1節(jié)的方法得d=31及相應(yīng)的樣本基S。95%控制限對(duì)應(yīng)滾動(dòng)鐵故障、內(nèi)圈故障的統(tǒng)計(jì)量為T2(95%)= 7.403 7，SPE(95%)=0.247 9，結(jié)果見(jiàn)圖14～圖17。

圖10 滾動(dòng)鐵故障的T2統(tǒng)計(jì)量故障檢測(cè)圖（標(biāo)準(zhǔn)KPCA）

圖11 滾動(dòng)鐵故障的SPE統(tǒng)計(jì)量故障檢測(cè)圖（標(biāo)準(zhǔn)KPCA）

圖12 內(nèi)圈故障的T2統(tǒng)計(jì)量故障檢測(cè)圖（標(biāo)準(zhǔn)KPCA）

圖13 內(nèi)圈故障的SPE統(tǒng)計(jì)量故障檢測(cè)圖（標(biāo)準(zhǔn)KPCA）

圖14 滾動(dòng)鐵故障的T2統(tǒng)計(jì)量故障檢測(cè)圖（稀疏KPCA）

圖15 滾動(dòng)鐵故障的SPE統(tǒng)計(jì)量故障檢測(cè)圖（稀疏KPCA）

圖16 內(nèi)圈故障的T2統(tǒng)計(jì)量故障檢測(cè)圖（稀疏KPCA）

由圖10～圖17可看出，數(shù)據(jù)在中間段以后有明顯變化，超過(guò)限值并持續(xù)，可以認(rèn)為是軸承發(fā)生故障，應(yīng)該予以報(bào)警。但在圖14～圖17（稀疏核主元的T2和SPE統(tǒng)計(jì)量故障檢測(cè)）中，樣本數(shù)200以后超過(guò)限值部分明顯，檢測(cè)出故障信號(hào)。從400個(gè)測(cè)試樣本中提取31個(gè)樣本構(gòu)成樣本基，因此核矩陣從400×400變?yōu)?1×31，極大減少了KPCA的計(jì)算量。與標(biāo)準(zhǔn)KPCA相比，計(jì)算效率得到提高，驗(yàn)證了稀疏核主元分析方法在滾動(dòng)軸承故障檢測(cè)中的有效性。

6 結(jié)論

本文提出的基于小波包能量譜-稀疏核主元的滾動(dòng)軸承故障檢測(cè)方法，首先將正常振動(dòng)信號(hào)進(jìn)行小波包分解，建立能量頻譜，并對(duì)能量頻譜構(gòu)建特征向量T，對(duì)特征向量T采用增量式樣本基構(gòu)造方法提取樣本基，以此樣本基作為訓(xùn)練樣本，建立正常KPCA模型；其次，將軸承信號(hào)進(jìn)行核主元分析，求取T2和SPE統(tǒng)計(jì)量，進(jìn)行故障檢測(cè)。稀疏核主元方法基于樣本基，從而提高了計(jì)算效率。理論分析和實(shí)驗(yàn)結(jié)果表明，利用該方法進(jìn)行滾動(dòng)軸承的故障檢測(cè)是有效的。

圖17 內(nèi)圈故障的SPE統(tǒng)計(jì)量故障檢測(cè)圖（稀疏KPCA）

[1]Schlkopf B，Smola A，Muller K R.Nonlinear component analysis as a kernel eigenvalue problem[J].Neural Computation，1998，10（5）：1299-1319.

[2]Li Junbao，Gao Huijun.Sparse data-dependent kernel principal component analysis based on least squares support vector machine for feature extraction and recognition[J]. Neural Comput&Applic，2012，21：1971-1980.

[3]Zhang Yingwei，Li Shuai，Hu Zhiyong.Improved multiscale kernel principal component analysis and its application for fault detection[J].Chemical Engineering Research and Design，2012，90：1271-1280.

[4]Nguyen V H，Golinval J C.Fault detection based on Kernel Principal Component Analysis[J].Engineering Structures，2010，32：3683-3691.

[5]趙小強(qiáng)，王新明，王迎.基于PCA與KPCA的TE過(guò)程故障檢測(cè)應(yīng)用研究[J].自動(dòng)化儀表，2011，32（1）：8-12.

[6]張億雄，顧海明，周勇軍.小波包分析在滾動(dòng)軸承故障診斷中的應(yīng)用[J].煤礦機(jī)械，2009，10（30）：240-242.

[7]孫潔娣，靳世久.基于小波包能量及高階譜的特征提取方法[J].天津大學(xué)學(xué)報(bào)，2010，43（6）：563-566.

[8]王冬云，張文志.小波包能量譜在滾動(dòng)軸承故障診斷中的應(yīng)用[J].測(cè)量與儀器，2010（11）：32-36.

[9]姜娜，紀(jì)建偉.基于小波包分析的滾動(dòng)軸承的故障診斷方法研究[J].浙江農(nóng)業(yè)學(xué)報(bào)，2012，24（2）：310-313.

[10]Xiao Yingqun，F(xiàn)eng Lianggui.A novel neural-network approach of analog fault diagnosis based on kernel discriminant analysis and partical swarm optimization[J].App-Lied Soft Computing，2012，12：904-920.

[11]彭濤，楊慧斌，李健寶，等.基于核主元分析的滾動(dòng)軸承故障混合域特征提取方法[J].中南大學(xué)學(xué)報(bào)，2011，42（11）：3385-3391.

[12]Jin Xiangyang，Lin Lin，Zhong Shisheng，et al.Rotor fault analysis of classification accuracy optimization based on kernel principal component analysis and SVM[J].Procedia Engineering，2011，15：5279-5283.

[13]付克昌，吳鐵軍.基于特征子空間的KPCA及其在故障檢測(cè)與診斷中的應(yīng)用[J].化工學(xué)報(bào)，2006，57（17）：2265-2269.

[14]許麗，張進(jìn)明.基于PCA的滾動(dòng)軸承故障檢測(cè)方法[J].計(jì)算機(jī)仿真，2010，27（6）：325-329.

ZHU Dan1，2,FAN Yugang1，2,ZOU Jinhui1，2,WU Jiande1，2,HUANG Guoyong1，2

1.Faculty of Information Engineering andAutomation,Kunming University of Science and Technology,Kunming 650500,China
2.Engineering Research Center for Mineral Pipeline Transportation,Kunming 650500,China

For the problem of rolling bearing fault detection,a method of rolling bearing fault detection is proposed,which is based on wavelet packet energy spectrum and sparse kernel principal component.The vibration signal is decomposed by wavelet packet,in order to extract the energy spectrum of the signal.Then the sample base of energy spectrum is extracted through the method of incremental sample base.A kernel principal component model is built by the sample base for the analysis of the energy spectrum of the bearing vibration signal.The experimental simulation is presented to illustrate the effectiveness of the algorithm.

rolling bearing;wavelet packet;sparse kernel principal component;fault detection

A

TP277

10.3778/j.issn.1002-8331.1301-0027

ZHU Dan,FAN Yugang,ZOU Jinhui,et al.Application of wavelet packet energy spectrum and sparse kernel principal component in fault detection.Computer Engineering and Applications,2014,50（21）：224-229.

國(guó)家自然科學(xué)基金（No.51169007）；云南省科技計(jì)劃項(xiàng)目（No.2010DH004，No.2011DA005，No.2012CA022）；云南省中青年學(xué)術(shù)和技術(shù)帶頭人后備人才培養(yǎng)計(jì)劃項(xiàng)目（No.2011CI017）。

朱丹（1987—），女，碩士，研究領(lǐng)域?yàn)榛诼暟l(fā)射的故障檢測(cè)；范玉剛（1973—），通訊作者，男，博士，副教授，研究領(lǐng)域?yàn)榛跈C(jī)器學(xué)習(xí)的智能信息處理、數(shù)據(jù)挖掘技術(shù)等；鄒金慧，女，副教授；吳建德，男，博士，教授；黃國(guó)勇，男，博士，副教授。E-mail：406742487@qq.com

2013-01-06

2013-04-17

1002-8331（2014）21-0224-06

CNKI出版日期：2013-04-26，http://www.cnki.net/kcms/detail/11.2127.TP.20130426.1018.001.html