朱會杰，王新晴，芮 挺，李艷峰，張紅濤，趙 洋

(1.解放軍理工大學野戰(zhàn)工程學院，江蘇南京210007；2.防空兵指揮學院，河南鄭州450052)

基于移不變稀疏編碼的單通道機械信號盲源分離

朱會杰1，王新晴1，芮 挺1，李艷峰1，張紅濤2，趙 洋1

(1.解放軍理工大學野戰(zhàn)工程學院，江蘇南京210007；2.防空兵指揮學院，河南鄭州450052)

針對特征反復出現的機械信號，提出了一種使用移不變稀疏編碼的單通道盲源分離方法。移不變稀疏編碼將原始信號看成多個基與系數的卷積，能夠根據信號的統(tǒng)計分布，利用信號自身特征自適應地學習到匹配的基和稀疏的系數。在恒定工況下，不同的信號源具有不同的特征，同一信號源的特征結構相似，將學習到的不同特性的基分別重構即可得到相應的源信號。將該方案應用于仿真的齒輪故障和軸承故障振動信號盲源分離問題中，以及用來提取實測的液壓泵壓力脈動。結果顯示，這種方法較其他方法有所改進，所需人工經驗少、抗噪能力強、信號恢復精度高、魯棒性好，適用于單通道機械信號盲源分離，為單通道信號盲源分離提供了一種新思路。

信號處理；移不變稀疏編碼；盲源分離；正交匹配追蹤；字典學習

引 言

在機械工程領域，很多觀測信號都是多種源信號的組合，需要使用盲源分離(blind source separation，BSS)從若干觀測信號中恢復出各個源信號。當源信號數量多于觀測信號數量是欠定盲源分離，常規(guī)盲源分離方法無法奏效。尤其是單通道盲源分離，已被公認為盲信號處理中的挑戰(zhàn)性難題［1-2］。

單通道盲源分離已被應用于聲 音處理3］、醫(yī)學［4］、故障診斷［5］、雷達信號分析［6］等多個領域，雞尾酒會問題本質也屬于這一問題。為此，有許多學者嘗試解決該問題。針對具有稀疏性的欠定線性瞬時混合盲源分離問題，Roan［7］，Shoko［8］等采用聚類的方法進行了嘗試，但對于單通道盲源分離問題，線性稀疏性條件往往不能滿足。還有一些學者將源信號分解，增加觀測信號的個數，將欠定問題轉化為適定或超定問題。成謝鋒［2］基于圓周卷積分層將原始信號由一維向量轉化成為多維向量，也有一些學者使用經驗模態(tài)分解(Empirical modes decomposition，EMD)［9-10］、集成經驗模態(tài)分解［11］、極值域均值模式分解［12］、局域均值分解［13］以及奇異值分解［14］等方法將信號分解成多個分量，增加觀測信號的個數，再使用獨立分量分析(independent component analysis，ICA)進行盲提取，但這種方式改變了信號的結構特征，受人工經驗影響較大，仍不能達到理想的效果。近些年來，稀疏編碼被用于解決欠定盲源分離問題，提高了分離的效果，但文獻［15］假設源信號的混合矩陣已知，文獻［16］使用預先設定的字典對信號的匹配能力仍有待提高。

Smith［17］提出的移不變稀疏編碼(Shift invariant sparse coding，SISC)，將信號表示為多個基(特征或原子)與系數的卷積形式，自適應地學習基與系數達到對原信號的稀疏表達。這一理論被發(fā)表于《Nature》雜志以來，因其自適應性好、抗噪能力強、所需人工經驗少、魯棒性強等優(yōu)點，已經引起了廣泛關注［3-5，18］。常見的機械信號特征都是重復出現的，比如故障軸承的沖擊、往復式活塞的壓力變化、液壓泵的流量變化等等。在恒定工況下，由于不同的源信號具有不同的特征，同一源信號的特征結構相似，使用SISC學習到這些特征后，分別將各個基進行重構，可得到相應的源信號，實現了單通道機械信號的盲源分離。

1 理論方法

1.1 SISC的原理

稀疏分解理論認為，使用一個包含K個基的過完備字典D=［d1，d2，…，dK］∈Rn×K，單通道信號y∈Rn可以被表示成這些基的線性組合，即y=Dx+ε，其中ε為標準白噪聲，向量x∈RK為稀疏的。

SISC將信號看成基與系數的卷積［17］

式中 *表示卷積，dj∈Rm，xj∈Rl是稀疏的，m＜n，l＜n。由稀疏分解理論可知［17］，SISC是依據統(tǒng)計概率對信號進行稀疏表示，只有學習到的基與信號本質特征相匹配時才能足夠稀疏地表示原信號。在恒定工況下，不同源信號具有不同的特征，相同的源信號具有類似的特征，可以把不同的基與對應系數的卷積(dj*xj)看作一個源信號或者源信號的一部分。

d和x都是未知的，無法同時對它們進行求解，但可以利用最大后驗概率將該問題轉化為系數求解和字典學習兩個過程，反復迭代，直到收斂或達到停止條件［17］，下面依次介紹這兩個過程［18］。

1.1.1 系數求解

本文選用正交匹配追蹤(Orthogonal matching pursuit，OMP)求解稀疏表示的系數，這一過程相當于固定所有的基，求解對應的系數

其中·表示非零系數的個數，T0為稀疏度。

OMP是匹配追蹤(Matching pursuit，MP)的改進，也是一種迭代的貪婪算法［19］。OMP每次選取與信號最匹配的原子之后，都要將全部所選原子正交化處理，并將它們重新在原始信號上進行投影。這種方式保證了殘余與所有已選原子正交，并且選取的原子都互不相關，因此殘差信號的分解過程也非常迅速。OMP加快了收斂速度又保證了最優(yōu)迭代性，因此本文選用OMP。

令D={gj｝1K為用于信號稀疏分解的過完備原子庫，其中gj為dj補零后的原子，其長度與原信號相同，且‖gj‖=1。OMP過程如下：

(1)設y為待處理的信號序列，信號殘余r0=y，已選原子的集合ψ-1為空集，迭代次數n=0。

(3)更新已選原子的集合ψn=ψn-1∪ {｝。

(4)根據最小二乘法計算投影系數an=(ψnTψn)-1·ψnTy，因此重構信號y~n=ψnan，殘余rn=y-n。

(5)更新迭代次數n=n+1，判斷殘余信號與原始信號的能量比‖rn‖2/‖y‖2是否小于設定值，或者是否達到稀疏度T0，若不滿足則返回步驟(2)，若滿足則輸出重構信號~yn，殘余為rn。

在選擇最匹配的原子時，需要大量的內積運算，使用快速卷積算法可以加快運算。

1.1.2 字典學習

字典學習相當于固定所有的系數，求解與信號最匹配的基：

直接求解該問題的計算量很大，已經有學者提出多種快速優(yōu)化算法。Olshausen［20］和Blumensath［3-4］利用差分求導來求取基的數值解，Szlam［21］使用主成分法進行基的更新。Grosse［18］使用了一種巧妙的優(yōu)化方式，將這一問題轉化到頻域求解，本文也借鑒了這種思想［18］。由帕斯瓦爾定理可知，式(4)可以轉化為下面的優(yōu)化問題

式中 f對應于頻率，系數λ∈RK，單位向量l∈RK，并且

使用牛頓法對參數λ進行優(yōu)化，代入式(7)即可完成對基的更新［18］。

1.2 基于SISC的單通道信號盲源分離流程

采用本文方法對恒定工況下的單通道信號進行盲源分離的流程如圖1所示，步驟如下：

(1)首先對算法進行初始化，設定最大迭代次數、稀疏度等。

(2)對觀測信號進行移不變稀疏分解，求得對應于不同源信號的基，并根據基與對應的重構信號的特性確定其來源。

(3)將不同的基和對應的系數進行重構，得到不同的源信號。

圖1 本文方案的流程圖Fig.1 The flow chart of proposed scheme

2 實驗驗證

2.1 仿真實驗驗證

齒輪和軸承是常用的機械部件，經常將它們組合使用，因此對它們的信號進行盲分離具有實際意義。根據齒輪和軸承的故障模型［22］，仿真了一個源信號為齒輪故障信號和軸承故障信號線性疊加的單通道含噪信號。齒輪故障時會出現幅值調制和相位調制現象，它的振動信號的模型為

式中 k為諧波成分的階次；Ak(t)，φk分別為第k次諧波成分的幅值和相位；fm為嚙合頻率。

故障軸承的沖擊響應可以簡化為

式中 a( t)為調幅幅值；p為衰減系數，代表了沖擊衰減的快慢；fb為共振的固有頻率；φ為相位。

因此故障軸承的振動信號為

這里 m為沖擊的次數，T為沖擊的周期，τm≈0.01～0.02 T，為隨機滑動。

仿真信號y=ym+yb+yn，yn為高斯白噪聲，信噪比為6 dB。各參數取值為：采樣頻率fs=12 000 Hz，信號長度為4 096點，K=5，Ak(t)=2.5·［sin( 60 π t)+2］/k，fm=340 Hz，φk=kπ/4；p=1 200，φ=0，fb=3 400 Hz，a( t)=20·［sin( 60 π t)+2］，T=0.01 s。觀測信號y的波形如圖2所示(由于篇幅有限，僅取其中1 000個點，下同)，源信號yb，ym和yn的波形如圖3所示。

圖2 觀測信號Fig.2 Measured signal

圖3 源信號Fig.3 Original signals

由圖2可知，軸承和齒輪故障信號混合在一起，已經無法看出它們的原始波形。從圖3可以看出噪聲的幅值已經與部分軸承和齒輪故障的沖擊接近，噪聲污染比較嚴重。

使用SISC對觀測信號進行分解，首先要確定源信號的個數。盲源分離包括盲源個數的估計和盲源的分離兩個部分，本文重點研究盲源的分離，盲源個數的估計可以參考其他學者提出的方法，比如文獻［2，14，16］的方法。本文提出了一種使用基的頻域相關度進行源數估計的方法。對于未知類別的源信號，首先假定一個比預估源數大一些的數作為基的個數進行分解，然后計算各基之間的相關度，如果有明顯相關的情況，則減少基的個數再次對信號分解，以此類推，直到沒有任何兩個基有明顯相關。定義相關度pi，j=ρ［FFT( di)，FFT( dj ) ］，其中FFT表示傅里葉變換，ρ(x，y)表示相關系數。這里使用基的傅里葉變換而非基本身求取相關系數是為了消除基的平移所引起相關系數的變化。

圖4 學習到的4個基Fig.4 Learned four bases

大量實驗表明，當Pij＞0.7時，可以認為基di與dj相關，應減少源數的個數；當Pij＜0.5時，可以認為基di與dj相關性不強，應屬于不同的源信號。而0.5＜Pij＜0.7時則要根據實際情況進一步討論或者參考其他方法。

對于上述仿真信號y，首先根據經驗預估一個比真實源數大的值作為基的個數，進行試探性分解，這里首先選用4個基，基的長度都選為35個點，稀疏度T0=180，最大迭代次數設為30次。基的波形如圖4所示，各個基之間的相關系數如表1所示。

表1 基的相關度Tab.1 The correlations between bases

由表1和圖4可知，基1與基2，基1與基4，基2與基4相互之間都存在著高度相關，并且它們的波形也比較相似，它們應該屬于同一個信號源。為了謹慎起見，將源信號個數也即基的個數減小為3，波形如圖5所示，基之間的相關度如表2所示。

由圖5、表2可知，基1與基3相關度較高，波形相似，基2與基1或基3的相關度低，波形結構差異大，因此源數應該為2。選擇源數(基的個數)為2進行分解的波形如圖6所示，它們的相關度為0.27，相關度較低，波形結構差異大，可以推斷該源信號包含兩個源信號。

圖5 學習到的3個基Fig.5 Learned three bases

表2 基的相關度Tab.2 The correlations between bases

由圖6可知，基1是一個衰減震蕩，屬于典型的軸承故障信號，基2是一個單峰沖擊，符合齒輪故障信號的特征。分別將它們與對應的系數進行重構，得到兩個估計的源信號，如圖7(a)，(b)所示。由圖7(a)可以看出這是一個調幅的周期性沖擊振蕩，再次驗證了該信號為軸承故障信號。而圖7(b)為調幅的周期性沖擊，為典型的齒輪故障信號。信號殘余為白噪聲，如圖7(c)所示。本文方案準確地將軸承故障信號和齒輪故障信號分離和恢復出來，僅在標注處與原信號有微小誤差，保留了源信號的細節(jié)。此外由于在SISC的模型中已經將高斯白噪聲考慮在內，其分解效果受噪聲干擾很小，抗噪能力強。

為進一步驗證該方法的有效性，還與其他學者的方法進行了對比。李志農［9］用EMD將源信號分解為多個分量，之后使用ICA進行盲分離，使用他的方法所得結果如圖8所示。董紹江［16］使用數學形態(tài)法(Mathematical morphology，MM)進行濾波，之后使用基于Gabor原子的匹配追蹤將源信號進行分解，把分解信號與源信號組合成新的觀測信號，再使用快速核獨立分量分析(Fast KICA)進行盲源分離。參考文獻［16］，并經實驗對比發(fā)現，形態(tài)結構元素為扁平型、長度為3、幅值為0時效果較好，故選用該形態(tài)結構元素，盲分離的結果如圖9所示。

圖6 學習到的2個基Fig.6 Learned two bases

圖7 本文方法估計的源信號Fig.7 Estimated signals by proposed method

為了進一步驗證這幾種方法的效果，分別使用信噪比(Signal to noise ratio，SNR)和均方誤差(Mean square error，MSE)對估計的源信號進行評價，5次實驗平均的結果如表3所示。

表3 不同方法效果對比Tab.3 The comparison of different algorithms

圖8 估計的源信號(EMD+ICA)Fig.8 Estimated signals by EMD+ICA

圖9 估計的源信號(MM+MP+Fast KICA)Fig.9 Estimated signals by MM+MP+Fast KICA

由圖8，9和表3可知，EMD+ICA和MM+MP+Fast KICA估計的軸承故障信號yb和齒輪故障信號ym都發(fā)生了明顯的失真，尤其是估計的軸次試驗的平均結果如圖10所示。隨著噪聲的增加，本文方法分解的性能有所下降，但無論在SNR還是MSE的表現上，都明顯優(yōu)于EMD+ICA和MM+MP+Fast KICA。實驗結果表明本文方法具有良好的魯棒性，在不同信噪比下均具有一定優(yōu)勢。承故障信號yb都混入了大量的噪聲，而且在SNR和MSE的對比中也都不如本文方案優(yōu)異。

為了驗證本文方法在不同信噪比下的性能，對仿真信號依次添加1～8 dB的白噪聲，然后計算信號yb與ym的信噪比與均方誤差，并與EMD+ICA方法以及MM+MP+Fast KICA方法進行對比，5

圖10 不同信噪比下各種方法效果對比Fig.10 The comparison of different algorithms for different SNR

2.2 在液壓壓力信號中的應用

圖11為液壓泵壓力信號采集平臺，液壓泵轉速為900 r/min，液壓泵有7個柱塞，采樣頻率為5 000 Hz。在測試過程中，壓力信號受到了電機、變頻器的電磁脈沖干擾，信號被嚴重污染。實測壓力信號如圖11(a)所示，可見油壓脈動已經被電磁干擾和白噪聲所淹沒，為此使用SISC對其進行盲源分離，提取油壓脈動成分。首先按照2.1節(jié)所述方法確定源數為2，因此選用2個基，一個用來匹配壓力脈動，一個來匹配電磁脈沖，分解的信號殘余則為白噪聲。設定基的長度為48個點，稀疏度為T0=150，最大迭代次數為20次。

兩個基的波形如圖12所示，盲源分離結果如圖13(b)，(c)，(d)所示。由圖12中可以看出，基1具有短時脈沖特性，基2形狀是先降后升，類似于壓力脈沖中壓力的變化波形，圖13(b)和(c)分別是由基1和基2的重構圖形，由重構波形可以進一步看出基1和基2分別匹配電磁干擾和壓力脈動。結果顯示，使用本文方法能將壓力脈動成分準確地從原信號中提取出來，消除了電磁干擾和白噪聲的影響，為下一步的分析奠定了基礎。

圖11 液壓泵壓力信號采集平臺Fig.11 The test rig of pressure signal acquisition for hydraulic pump

圖12 學習到的基Fig.12 Learned bases

圖13 原始壓力信號與分解結果Fig.13 Original pressure signal and its decompositions

3 結 論

本文基于SISC提出了一種對恒定工況下單通道信號進行盲源分離的新方法，并將其用于源信號為仿真的軸承故障信號和齒輪故障信號的盲分離問題，以及提取實測壓力信號的脈動成分，并與其他學者的方法進行了對比。經實驗驗證，這種方法不僅能夠實現單通機械信號的盲源分離，為單通道盲源分離以及盲源分離提供了一種新的思路，而且信號恢復精度高，具有較強的抗噪能力和魯棒性，是一種優(yōu)異的盲分離方法。

［1］張赟，李本威，賈舒宜，等.航空發(fā)動機混疊振動信號的欠定盲源分離方法［J］.推進技術，2014，35(04)：552—558.Zhang Yun，Li Ben-wei，Jia Shu-yi，et al.Underdetermined blind source sparation of aeroengine vibration signal mixtures［J］.Journal of Propulsion Tchnology，2014，35(04)：552—558.

［2］成謝鋒，馬勇，張學軍，等.一種不用先驗知識的單路混合信號的盲源分離新方法［J］.電子學報，2011，39(10)：2 317—2 321.Cheng Xie-feng，Ma Yong，Zhang Xue-jun，et al.Asingle-channel mixed signal BSS new method without using the prior knowledge［J］.Acta Electronica Sinica，2011，39(10)：2 317—2 321.

［3］Blumensath T，Davies M.Sparse and shift-invariant representations of music［J］.Audio，Speech，and Language Processing，IEEE Transactions on，2006，14(1)：50—57.

［4］Blumensath T，Davies M.Shift-invariant sparse coding for single channel blind source separation［J］.SPARS，2005，5：75—78.

［5］Tang H，Chen J，Dong G.Sparse representation based latent components analysis for machinery weak fault detection［J］.Mechanical Systems and Signal Processing，2014，http：//dx.doi.org/10.1016/j.ymssp.2014.01.011.

［6］陳曉軍，成昊，唐斌.基于ICA的雷達信號欠定盲分離算法［J］.電子與信息學報，2010，32(4)：919—924.Chen Xiao-jun，Cheng Hao，Tang Bin.Underdetermined blind radar signal separation based on ICA［J］.Journal of Electronics&Information Technology，2010，32(4)：919—924.

［7］Roan M，Erling J，Sibul L.Anew，non-linear，adaptive，blind source separation approach to gear tooth failure detection and analysis［J］.Mechanical Systems and Signal Processing，2002，16(5)：719—740.

［8］Araki S，Sawada H，Mukai R，et al.Underdetermined blind sparse source separation for arbitrarily arranged multiple sensors［J］.Signal Processing，2007，87(8)：1 833—1 847.

［9］李志農，呂亞平，范濤，等.基于經驗模態(tài)分解的機械故障欠定盲源分離方法［J］.航空動力學報，2009，24(8)：1 886—1 892.Li Zhi-nong，Lv Ya-ping，Fan Tao，et al.Underdetermined blind source separation method of machine faults based on empirical mode decomposition［J］.Journal of Aerospace Power，2009，24(8)：1 886—1 892.

［10］孫潔娣，郝雅立，溫江濤，等.基于EMD的高壓燃氣管道泄漏信號欠定盲分離方法［J］.振動與沖擊，2013，32(18)：81—86.Sun Jie-di，Hao Ya-li，Wen Jiang-tao，et al.Underdetermined blind source separation method of pipeline leakage signals based on empirical mode decomposition［J］.Journal of Vibration and Shock，2013，32(18)：81—86.

［11］劉佳，楊士莪，樸勝春.基于EEMD的地聲信號單通道盲源分離算法［J］.哈爾濱工程大學學報，2011，32(2)：194—199.Liu Jia，Yang Shi-e，Piao Sheng-chun.The single channels eismic-acoustic signal blind-source separation method based on EEMD［J］.Journal of Harbin Engineering University，2011，32(2)：194—199.

［12］孟宗，梁智.基于EMMD和BSS的單通道旋轉機械故障診斷方法［J］.儀器儀表學報，2013，34(3)：635—642.Meng Zong，Liang Zhi.Fault diagnosis method for single channel rotating machinery based on EMMD and BSS［J］.Chinese Journal of Scientific Instrument，2013，34(3)：635—642.

［13］李志農，劉衛(wèi)兵，易小兵.基于局域均值分解的機械故障欠定盲源分離方法研究［J］.機械工程學報，2011，47(7)：97—102.Li Zhi-nong，Liu Wei-bing，Yi Xiao-bing.Underdetermined blind source separation method of machine faults based on local mean decomposition［J］.Journal of Mechanical Engineering，2011，47(7)：97—102.

［14］申永軍，楊紹普，孔德順.基于奇異值分解的欠定盲信號分離新方法及應用［J］.機械工程學報，2009，45(8)：64—70.Shen Yong-jun，Yang Shao-pu，Kong De-shun.Newmethod of blind source separation in under-determined mixtures based on singular value Decomposition and application［J］.Journal of Mechanical Engineering，2009，45(8)：64—70.

［15］余豐，奚吉，趙力，等.基于CS與K-SVD的欠定盲源分離稀疏分量分析［J］.東南大學學報：自然科學版，2012，41(6)：1 127—1 131.Yu Feng，Xi Ji，Zhao Li，et al.Sparse presentation of underdetermined blind source separation based on compressed sensing and K-SVD［J］.Journal of Southeast University(Natural Science Edition)，2012，41(6)：1 127—1 131.

［16］董紹江，湯寶平，張焱.基于最優(yōu)匹配跟蹤算法的單通道機械信號盲源分離［J］.振動工程學報，2013，25(6)：724—731.Dong Shao-jiang，Tang Bao-ping，Zhang Yan.Blind source separation of single-channel mechanical signal based on optimal matching pursuit algorithm［J］.Journalof Vibration Engineering，2013，25(6)：724—731.

［17］Smith E C，Lewicki MS.Efficient auditory coding［J］.Nature，2006，439(7079)：978—982.

［18］Grosse R，Raina R，Kwong H，et al.Shift-invariance sparse coding for audio classification［A］.Conference on Uncertainty in AI［C］.2007：149—158.

［19］Needell D，Vershynin R.Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit［J］.Foundations of Computational Mathematics，2007，9(3)：317—334

［20］Olshausen B A.Sparse coding of time-varying natural images［A］.Proceeding of Internatinal Conference of Independent Component Analysis and Blind Source Separation(ICA)［C］.2000：603—608

［21］Szlam A，Kavukcuoglu K，LeCun Y.Convolutional matching pursuit and dictionary training［OL］.http：//ar Xiv preprint ar Xiv：1010.0422，2010.

［22］Sawalhi N，Randall R.Simulating gear and bearing interactions in the presence of faults：Part I.The combined gear bearing dynamic model and the simulation of localised bearing faults［J］.Mechanical Systems and Signal Processing，2008，22(8)：1 924—1 951.

Shift invariant sparse coding for blind source separation of single channel mechanical signal

ZHU Hui-jie1，WANG Xin-qing1，RUI Ting1，LI Yan-feng1，ZHANG Hong-tao2，ZHAO Yang1

(1.College of Field Engineering，PLAUniversity of Science and Technology，Nanjing 210007，China；2.PLAAir Defense Forces Command College，Zhengzhou 450052，China)

For the single channel mechanical signal with repeated features，the method for blind source separation based on shift invariant sparse coding was proposed in this paper.In the literatures of shift invariant sparse coding，a signal is described as the convolutions of multi bases and their coefficients.According to statistical distribution of a signal，shift invariant sparse coding could adaptively learn its bases and the sparse coefficients from the structures of the signal itself.Under stable condition，different signal sources have different features，and the features from the same source are similar，thus the learned bases with different features could be used to reconstruct corresponding signal sources.This scheme was applied in the blind source separation of simulated vibration signals of faulty gear and bearing，as well as the extraction of pressure pulsation of hydraulic pump.The result showed that this algorithm has improved a lot compared to other algorithms，and this algorithm needs less expertise，has strong anti-interference ability，in addition，it is robust and could recover original signals more accurately.Therefore，

this technique is appropriate to blind source separation for single channel mechanical signal，and provides a new way for single channel blind source separation.

signal processing；shift invariant sparse coding；blind source separation；orthogonal matching pursuit；dictionary learning

TN911.7；TH165+.3

A

1004-4523(2015)04-0625-08

10.16385/j.cnki.issn.1004-4523.2015.04.016

朱會杰(1987—)，男，博士研究生。電話：(025)80821440；E-mail：zhuhuijiehao@163.com

2014-05-03；

2014-07-18

國家自然科學基金資助項目(61472444)