李 萍,王樹青,李華軍
(中國海洋大學(xué)工程學(xué)院,山東青島266100)
特征系統(tǒng)實現(xiàn)算法中的噪聲問題研究*
李 萍1,王樹青2,李華軍3
(中國海洋大學(xué)工程學(xué)院,山東青島266100)
研究特征系統(tǒng)實現(xiàn)算法在模態(tài)識別過程中的噪聲處理問題。特征系統(tǒng)實現(xiàn)算法為目前土木工程領(lǐng)域應(yīng)用較為廣泛的一種模態(tài)識別法。在處理實測數(shù)據(jù)時,特征系統(tǒng)實現(xiàn)算法將實測數(shù)據(jù)分解為結(jié)構(gòu)真實信號和噪聲信號,從而將噪聲消除。分塊Hankel矩陣的維數(shù)對信噪分離過程影響很大。提出當(dāng)分塊Hankel矩陣的分塊行與分塊列數(shù)接近時,真實和噪聲信號可以最好的分離。通過在實測信號里添加噪聲(產(chǎn)生噪聲-噪聲信號),將實測信號和噪聲-噪聲信號識別的結(jié)果進行比較,提出了1種驗證識別模態(tài)參數(shù)是否為結(jié)構(gòu)真實模態(tài)的檢驗方法。本文應(yīng)用導(dǎo)管架平臺實例來驗證提出方法的適用性。結(jié)果表明,當(dāng)分塊Hankel矩陣的分塊行與分塊列數(shù)接近時,真實和噪聲信號能夠有效地分離,此時添加的噪聲主要影響噪聲信號部分。由實測信號和噪聲-噪聲信號識別的模態(tài)參數(shù)非常吻合。
模態(tài)參數(shù)識別;特征系統(tǒng)實現(xiàn)算法;噪聲;模型階數(shù)
模態(tài)參數(shù)識別是通過結(jié)構(gòu)的動態(tài)特性試驗數(shù)據(jù),來估計結(jié)構(gòu)模態(tài)參數(shù)的過程[1-2]。模態(tài)參數(shù)識別結(jié)果的準(zhǔn)確性直接影響到結(jié)構(gòu)物有限元模型修正和損傷檢測工作的進行[1,3]。參數(shù)識別方法分頻域法和時域法2類。時域法是1970年代隨著計算機技術(shù)迅速發(fā)展起來的,先后經(jīng)歷了單輸入單輸出(SISO)、單輸入多輸出(SIMO)、多輸入多輸出(M IMO)3個階段。
特征系統(tǒng)實現(xiàn)算法(Eigensystem Realization A lgorithm,ERA)是美國NASA的Langley研究中心于1984年在最小實現(xiàn)理論基礎(chǔ)上提出來的[4],屬于M IMO時域整體模態(tài)參數(shù)識別方法。該方法源于控制理論中Ho-Kalman的最小實現(xiàn)理論[5],利用脈沖響應(yīng)數(shù)據(jù),采用奇異值分解(SVD),求得模態(tài)參數(shù)。目前在國外航空航天領(lǐng)域應(yīng)用廣泛。
特征系統(tǒng)實現(xiàn)算法通過M arkov參數(shù)(如:脈沖響應(yīng)數(shù)據(jù))來組建分塊Hankel矩陣,并以此矩陣作為離散時間狀態(tài)模型[6]實現(xiàn)的基礎(chǔ)。通過最小實現(xiàn)理論,ERA方法可用于實測信號的模態(tài)參數(shù)識別。ERA算法的實質(zhì)是:將實測信號分解為真實信號和噪聲信號,并借助于模態(tài)置信度(MAC)和模態(tài)奇異值(M SV)來區(qū)分真實和噪聲模態(tài)[7]。其中,分塊Hankel矩陣秩(即模型階數(shù))的確定對能否將真實信號和噪聲信號有效地分離非常關(guān)鍵。周等[8]提出了利用奇異值差值法進行模型定階,該方法具有直觀,無需進行反復(fù)試算的特點。Hu[9]等,Peeters[10]等,M isra[11]等利用了奇異值分解法進行模型定階。目前的噪聲處理方法大多對信號進行降噪預(yù)處理,得到消噪信號后再進行模態(tài)識別。林[12]等提出了利用小波去噪方法對脈沖響應(yīng)信號進行處理。Zhang[13-14]利用相關(guān)濾波對測量數(shù)據(jù)進行預(yù)處理,有效地提高了識別精度。本文研究ERA方法中的噪聲處理問題。區(qū)別于以往對實測信號進行降噪預(yù)處理的方法,本文研究模態(tài)分析過程中如何將真實信號和噪聲信號有效地分離。Hankel矩陣維數(shù)對其信噪分離過程很關(guān)鍵。如果矩陣的維數(shù)選擇不佳,可能導(dǎo)致真實和噪聲信號不能很好地分離,從而影響模態(tài)參數(shù)識別的準(zhǔn)確性。
本文提出了確定最佳分塊Hankel矩陣的方法。通過對此矩陣最佳維數(shù)的確定,可以使真實信號和噪聲信號有效分離,保證識別模態(tài)的準(zhǔn)確性。本文還提出驗證識別的模態(tài)參數(shù)魯棒性和準(zhǔn)確性的方法。此方法通過在實測信號中添加噪聲,稱之為噪聲—噪聲信號,并比較由實測信號和噪聲—噪聲信號識別的模態(tài)參數(shù)的符合性,來確定識別模態(tài)是否為真實模態(tài)。文中借助導(dǎo)管架平臺模型檢驗所提出方法的有效性。
時不變系統(tǒng)的連續(xù)時間狀態(tài)方程可表示為:

其中,x∈Rn,u∈Rr,y∈Rm分別為狀態(tài)向量,輸入向量和輸出向量;n=2M,為模型階數(shù);M為系統(tǒng)參與模態(tài)數(shù);r和m分別為輸入和輸出數(shù)量。常數(shù)矩陣Ac,Bc,C和D代表線性系統(tǒng)的內(nèi)部關(guān)系。
離散時間狀態(tài)空間模型可表示為:

其中,k代表某個時刻,常數(shù)矩陣A和B可由Ac,Bc得到。
當(dāng)輸入為脈沖激勵時,輸出脈沖響應(yīng)數(shù)據(jù)Yk∈Rm×r可表示為:

公式(3)中的常數(shù)矩陣序列又稱為M arkov參數(shù)。ERA系統(tǒng)實現(xiàn)為通過公式(3)中的M arkov參數(shù)來求解公式(2)中的矩陣組[A,B,C]的過程。同一個輸入輸出系統(tǒng)可有多個實現(xiàn)來表示。其中最小實現(xiàn)指狀態(tài)矩陣A的維數(shù)最小時對應(yīng)的模型。
假設(shè)狀態(tài)矩陣A線性獨立的特征向量為(Ψ1,Ψ2,…,Ψn),對應(yīng)的特征值為(λ1,λ2,…,λn)。定義Λ為特征值對角矩陣,Ψ為特征向量矩陣。矩陣組[A,B,C]的實現(xiàn)過程可轉(zhuǎn)換為矩陣組[Λ,Ψ-1B,CΨ]。對角矩陣Λ中包含了模態(tài)阻尼比和阻尼頻率信息;矩陣Ψ-1B包含初始模態(tài)幅值信息;矩陣CΨ包含傳感器位置處的模態(tài)振型信息。結(jié)構(gòu)的所有模態(tài)參數(shù)可由矩陣組[Λ,Ψ-1B,CΨ]來獲得。在將特征值對角陣Λ通過Λc=ln(Λ)/Δt轉(zhuǎn)換到連續(xù)時間域后,其實部和虛部分別對應(yīng)模態(tài)阻尼比和阻尼頻率。這里需要指出,從數(shù)學(xué)角度,當(dāng)狀態(tài)矩陣A的特征值存在重根時,其對應(yīng)的模態(tài)振型并不一定唯一。重根對應(yīng)的2個模態(tài)振型的任一組合均可構(gòu)成一個新的振型[6]。在工程領(lǐng)域,結(jié)構(gòu)頻率可能非常接近,出現(xiàn)重根(相同頻率)的幾率很小。
系統(tǒng)實現(xiàn)是通過公式(3)中的Markov參數(shù)來組建維數(shù)為αm×βr的分塊Hankel矩陣開始的:

其中,α為分塊行數(shù);β為分塊列數(shù)。
模型階數(shù)n的選擇對模態(tài)參數(shù)識別精度影響很大。目前,模型定階已成為模態(tài)參數(shù)識別領(lǐng)域的最為重要的問題。這里模型定階方法是通過矩陣H(0)的秩來確定n的值[9,15]。
當(dāng)k=1時,可以得到

n為模型階數(shù),當(dāng)αm≥n且βr≥n時,矩陣H(0)的秩理論上為n。
目前常用確定矩陣秩的方法為奇異值分解法。通過對矩陣H(0)進行奇異值分解(SVD),

其中,矩陣R∈Rαm×αm和S∈Rβr×βr為正交矩陣,矩陣Σ∈Rαm×βr為對角矩陣,其對角元素為降序排列的奇異值。理論上,超出矩陣秩的奇異值應(yīng)當(dāng)為0,即Σ可分解為

其中,Σn=diag[σ1,σ2,…,σn]。對實測信號,由于噪聲影響,超出矩陣秩的奇異值并不等于0,而會變得很小。這里采用的確定模型階數(shù)的方法為首先將奇異值根據(jù)第一個奇異值σ1進行歸一化處理,然后劃出歸一化奇異值圖,當(dāng)奇異值圖趨近于水平時,對應(yīng)的奇異值個數(shù)即為矩陣的秩[9]。
當(dāng)k=2時,公式(4)為

定義 Oi為階數(shù)i的空矩陣,Ii為階數(shù)i的單位陣,則以下矩陣組
為一最小實現(xiàn)。其中,Rn和Sn分別為矩陣R和S的前n列。^代表估算值。當(dāng)測量信號噪聲很小時,狀態(tài)矩陣^A的秩為n,即系統(tǒng)的階數(shù)。
M SV是量測識別的各階模態(tài)對實測響應(yīng)信號貢獻大小的一種方法。M SV的計算公式如下:

其中,


假設(shè)某實測信號有M個參與模態(tài),其對應(yīng)的分塊Hankel矩陣H可分解為2部分:

其中,S和N分別代表真實信號和噪聲信號對應(yīng)的矩陣。信噪分離過程就是基于H(0)得到S。在將矩陣H(0)逼近到S時,通常應(yīng)用最小化H(0)和S的差別的Frobenius范數(shù)的方法,即:‖H(0)-S‖2最小。矩陣S可由矩陣H(0)通過奇異值截斷得到,這種方法稱為截斷奇異值分解法(TSVD)或者Eckart-Young方法[16]。注意,矩陣S不具備原有矩陣H(0)的分塊Hankel結(jié)構(gòu)。此外,本方法不需迭代。
特征系統(tǒng)實現(xiàn)算法噪聲處理與分塊行數(shù)α和分塊列數(shù)β的選取方式關(guān)系很大。ERA法識別模態(tài)參數(shù)的準(zhǔn)確度主要依賴于能否將S和N有效地分離。1個最直接的想法是使H(0)的第n個奇異值σn與第n+1個奇異值σn+1的差別越大越好(n=2M),即由歸一化奇異值圖觀察到奇異值σn和σn+1處有一個突降。這里取α和β最接近的情況,這時矩陣H(0)的維數(shù)最大。假設(shè)每個信號只采用99個信號點,滿足α+β-1=99的情況下,則(α,β)的最優(yōu)取法為(50,50),此時矩陣H(0)維數(shù)最大,且奇異值數(shù)p也最多(p=min(αm,βr))。當(dāng)參與S的模型階數(shù)為2M時,參與N的組成個數(shù)為p-2M,這里假設(shè)所有的p-2M個組成均為噪聲。
采用某一海洋導(dǎo)管架平臺模型來驗證提出的噪聲處理方法的有效性,見圖1。平臺采用鋼材,甲板尺寸為0.6 m×0.3 m,厚度為0.01 m;平臺高度1.7 m,4條腿斜率均為10∶1,采用直徑D=0.025 m,厚度t=2.5 mm的鋼管;導(dǎo)管架設(shè)3層水平橫撐,分別位于1.35,0.9和0.5 m處,在0.9和0.5 m處設(shè)豎向斜撐。水平橫撐和斜撐的尺寸均采用直徑D=0.016 m,厚度1.5 mm的鋼管。
平臺上共布設(shè)32個傳感器(4層,每層8個傳感器),由于實際海洋平臺通常把傳感器布設(shè)在水面以上,這里只采用平臺上層的8個傳感器進行分析(見圖1),x和y向傳感器各4個。分別施加脈沖激勵于平臺甲板的中間和邊緣的x和y方向(見圖1),可獲得共32個測量信號,采樣頻率為500 Hz(Δt=0.002 s)。同時,利用ANSYS建立有限元模型,計算得到前8階模態(tài)頻率(不包括剛體模態(tài))見表1。

圖1 實驗室內(nèi)導(dǎo)管架平臺模型Fig.1 Test jacket-type p latform in the lab

表1 有限元得到的前八階模態(tài)Table 1 Modal p roperties of the first eight modes of the FEmodel fo r a jacket-type offsho re p latform
對每個實測信號,只取200個測點進行分析,見圖2。為了檢驗識別模態(tài)的真實性,對測量信號后期添加噪聲,產(chǎn)生噪聲—噪聲信號,進而比較測量信號和噪聲—噪聲信號識別的模態(tài)參數(shù)的符合性。如果他們相一致,則表明添加的噪聲沒有影響真實信號,則測量信號中的噪聲對其識別的準(zhǔn)確性也影響不大。假設(shè)如下:如果S和N能有效地分離開來,則添加的噪聲將主要影響N,因而基于S的模態(tài)識別變化會非常小。噪聲—噪聲信號的模擬是通過在測量信號基礎(chǔ)上添加高斯白噪聲得到的。添加噪聲的概率密度函數(shù)滿足正態(tài)分布,同時它的功率譜密度函數(shù)理論上是常數(shù)。噪聲水平通過一個百分比來定量描述,該百分比定義為白噪聲的標(biāo)準(zhǔn)差和測量信號的標(biāo)準(zhǔn)差之比。當(dāng)噪聲水平很低時,其對矩陣秩的影響并不明顯,只有噪聲達到一定水平時,矩陣的秩才會發(fā)生明顯變化。同時,圖2中還給出的是某信號添加的15%噪聲。圖3為圖2中給出的某15%噪聲的功率譜密度函數(shù)。由圖3可知,其功率譜密度函數(shù)大體在某常數(shù)值附近波動。

圖2 某實測信號截取的200個測點和添加的15%噪聲Fig.2 An examp le of themeasured and injected 15%noise signal

圖3 圖2中給出的某15%噪聲的功率譜密度函數(shù)Fig.3 Power spectra density of injected 15%noise in Fig.2
首先,研究噪聲對單個信號模型階數(shù)的影響。對于圖2給出的單個信號實例,其測量信號和含15%的噪聲—噪聲信號對應(yīng)的H(0)(取α=100,β=100)矩陣的歸一化奇異值曲線見圖4,每個奇異值均由第一個(最大)奇異值歸一化所得。從圖4可知,對于單一信號,添加的15%噪聲對奇異值影響很大,其影響已從第4個奇異值開始。
對實測的32個信號,采用前199個測點來組建分塊Hankel矩陣H(0)。可通過H(0)秩的大小來獲得參與模態(tài)數(shù)量,理論上H(0)的秩(即模型階數(shù))為模態(tài)數(shù)的2倍[17]。

圖4 某實測信號和含15%噪聲-噪聲信號對應(yīng)的歸一化奇異值圖Fig.4 Normalized singular values of the Hankel matrix from one signal
這里研究矩陣H(0)的2種情況:(1)α=100,β=100,矩陣H(0)∈R800×400;(2)α=3,β=197,矩陣H(0)R24×788。2種情況對應(yīng)的歸一化奇異值曲線見圖5。在這里當(dāng)歸一化奇異值曲線突降接近水平時,認(rèn)為此奇異值處對應(yīng)分塊Hankel矩陣的秩。從圖5可以觀察,對于情況一,可以明顯觀察到在第8個奇異值處曲線有明顯的突降,而對于情況二則沒有此現(xiàn)象。因此,可以從情況一中得出矩陣的秩為8。對于第二種情況,采用相同的模型階數(shù)。

圖5 2種情況下實測信號對應(yīng)的歸一化奇異值圖Fig.5 Normalized singular values of the block Hankelmatrices from 128 signals
在32個測量信號上分別添加5%和15%的白噪聲來產(chǎn)生噪聲—噪聲信號。2種情況下對應(yīng)的測量和噪聲—噪聲信號的歸一化奇異值曲線分別見圖6和圖7。
從圖6可知,對第一種情況,測量信號、含5%和15%噪聲-噪聲信號對應(yīng)矩陣的秩皆為8。且添加的噪聲絕大部分影響第8個之后的奇異值,對前8個奇異值影響非常小。而對于情況二(見圖7),添加噪聲主要影響第8個之后的奇異值,但對前8個奇異值(5~8個)也有影響,由于不能由圖7明確判定矩陣的秩,這里采用與情況一相同的模型階數(shù)。

圖6 測量和噪聲—噪聲信號對應(yīng)的歸一化奇異值圖(情況一)Fig.6 No rmalized singular valuesof the block Hankelmatrices for case 1

圖7 測量和噪聲—噪聲信號對應(yīng)的歸一化奇異值圖(情況二)Fig.7 Normalized singular values of the block Hankelmatrices for case 2

表2 測量信號和噪聲—噪聲信號識別的模態(tài)頻率(情況一)Table 2 Estimated frequencies from measured,5%-noise and 15%-noise added signals(Case 1)

表3 測量信號和噪聲—噪聲信號識別的模態(tài)頻率(情況二)Table 3 Estimated frequencies f rom measured,5%-noise and 15%-noise added signals(Case 2)
在確定模型階數(shù)之后,根據(jù)公式(9)可求得一最小實現(xiàn),進而求得模態(tài)參數(shù)。表2和表3中分別為兩種情況下由測量信號、含5%和15%噪聲-噪聲信號識別的模態(tài)頻率。對比表2和表3可知,情況一識別的模態(tài)頻率在添加噪聲前后差別很小,相對誤差均小于0.005%,而情況二識別頻率的最大相對誤差為0.632%。此外,情況一的識別頻率與有限元得出的頻率也更接近。
表4和5中分別為2種情況下由測量信號、含5%和15%噪聲-噪聲信號識別的模態(tài)阻尼比。由表4可知,情況一中識別的阻尼比相對誤差最大處為第4階識別模態(tài),為10%。對比情況二(表5),其阻尼比最大相對誤差為990.476%,為第四階模態(tài)。可知當(dāng)α和β取值接近,即當(dāng)矩陣H(0)分塊行和列的數(shù)目相近時,識別結(jié)果明顯提高,尤其是阻尼比,此時噪聲信號已被有效地分離并剔除,識別結(jié)果比較準(zhǔn)確。

表4 測量信號和噪聲—噪聲信號識別的模態(tài)阻尼比(情況一)Table 4 Estimated damping ratios from measured,5%-noise and 15%-noise added signals(Case 1)

表5 測量信號和噪聲—噪聲信號識別的模態(tài)阻尼比(情況二)Table 5 Estimated damping ratios from measured,5%-noise and 15%-noise added signals(Case 2)
表6和7分別為2種情況下由測量信號、含5%和15%噪聲-噪聲信號識別的模態(tài)奇異值(M SV)。M SV是各識別模態(tài)對測量信號貢獻大小的量度。由表可見,情況一(見表6)識別的M SV值一致性要比情況二(見表7)的好。此外,由表6和表7可知,第3階識別模態(tài)貢獻最大,第4階識別模態(tài)貢獻最小。

表6 測量信號和噪聲—噪聲信號識別的M SV值(情況一)Table 6 Estimated MSV from measured,5%-noise and 15%-noise added signals(Case 1)

表7 測量信號和噪聲—噪聲信號識別的M SV值(情況二)Table 7 Estimated MSV from measured,5%-noise and 15%-noise added signals(Case 2)
針對特征系統(tǒng)實現(xiàn)算法,本文提出了基于分塊Hankel矩陣維數(shù)選擇來進行信噪分離,從而將噪聲剔除的方法。通過采取分塊行和列的數(shù)目最接近的Hankel矩陣結(jié)構(gòu),識別的模態(tài)參數(shù)準(zhǔn)確性有明顯提高,尤其是識別阻尼比有明顯改善。通過在測量信號中添加噪聲,并比較測量信號和噪聲—噪聲信號識別模態(tài)參數(shù)的一致性,來判別識別模態(tài)是否為結(jié)構(gòu)真實模態(tài)。物理模型試驗結(jié)果表明,通過選擇最佳的分塊Hankel矩陣維數(shù),可以有效地將真實信號和噪聲信號分離,提高識別精度。
[1] Ew ins D J,Modal Testing:Theory,Practice and Applications[M],2nd Edition.Baldock,Hertfordshire,England:Research Studies Press,2000.
[2] Maia N,Silva J,He J,et al.Theoretical and Experimental Modal A-nalysis[M].Taunton,Somerset,England:Research Studies Press,1997.
[3] Friswell M I,Mottershead J E,Finite Element Model Updating in Structural Dynamics[M].1st edition,Nethelands:Kluwer Acadamic Publishers,Sp ringer,1995.
[4] Juang J N,Pappa R S,An eigensystem realization algo rithm fo r modal parameter identification and model reduction[J].Journal of Guidance,Control,and Dynamics.1985,8(5):239-260.
[5] Ho B L,Kalman R E,Effective construction of linear state variablemodels from input/output data[C].∥Proceedingsof the 3rd Annual A llerton Conference on Circuit and System Theory.U rbana:university of Illinois,1966,14:545-548.
[6] Juang J N.Applied system identification[M].1st Edition:USA:Prentice-Hall Inc.,1994.
[7] Juang J N,Pappa R.S,Effectsof noise onmodal parameters identified by the eigensystem realization algorithm[J].Journal of Guidance,Control,and Dynamics.1986,9(3):294-303.
[8] 周幫友,胡邵全,杜強,特征系統(tǒng)實現(xiàn)算法中的模型定階方法研究[J].科學(xué)技術(shù)與工程,2009,9(10):2715-2716.
[9] Hu SL J,Bao X X,Li H J.Model order determination and noise removal for modal parameter estimation[J].Mechanical Systems and Signal Processing,2010,24(6):1605-1620.
[10] Peeters B,De Roeck G.Reference-based stochastic subspace identification for output-only modal analysis[J].Mechanical Systems and Signal Processing,1999,13(6):855-878.
[11] Misra P.Nikolaou M,Input design for model order determination in subspace identification[J].AlChEJournal,2003,49(8):2124-2132.
[12] 林貴斌,陸秋海,郭鐵能.特征系統(tǒng)實現(xiàn)算法中的小波去噪方法研究[J].工程力學(xué),2004,21(6):91-96.
[13] Zhang L M.A polyreference frequency domain method formodal parameter identification[C].Ohio:ASME Design Engineering Division Conference and Exbit on Mechanical Vibration and Noise OH IO,1985,85-DET-106.
[14] Zhang L M,An imp roved time domain polyreference method for modal identification[J].Mechanical Systems and Signal Process-ing,1987,1(4):399-413.
[15] Sanliturk K Y,Cakar O,Noise elimination from measured frequency response functions[J],Mechanical Systems and Signal Processing,2005,19(3):615-631.
[16] Eckart C,Young G,The approximation of onematrix by another of lower rank[J].Psychometrika,1936,1(3):211-218.
[17] Li P,Li H J,Hu SL J.Estimating modal parameters from free vibration response of jacket-type platforms[C].China:Proceedings of the Twentieth International Offshore and Polar Engineering Conference,2010:20-25.
Noise Handling of the Eigensystem Realization Algorithm
L IPing1,WANG Shu-Qing2,L IHua-Jun3
(College of Engineering,Ocean University of China,Qingdao 266100,China)
In dealing with noisy measurement data,the Eigensystem Realization Algo rithm(ERA)partitions the realized model into signal and noise portions so that the noise portion can be disregarded.A critical issue is the determination of the dimensions of the block Hankelmatrix.We show that the signal and noise matrices can be better separated when the number of block-row sand number of block-columns of the corresponding block Hankelmatrix are chosen to be close to each other.We propose a verification procedure to justify that the estimated modal parameters are noise insensitive and thus indeed associated with the true system.The procedure involves artificially injecting random noise into the measured signals to create noisy-noisy signals,then comparing the identification results obtained respectively from the measured and noisy-noisy signals.Using experimental data collected from a jacket-type p latform,we demonstrate that signal and noise portions can be p roperly separated w hile the number of the block-rows and number of block-columns are close.
model parameters identification;eigensystem realization algorithm;noise;model order
TP273
A
1672-5174(2011)7/8-176-07
國家自然科學(xué)基金項目(51079134),重大國際合作研究項目(51010009)資助
2011-01-08;
2011-4-14
李 萍(1983-),女,博士生。E-mail:lipingljk@gmail.com
責(zé)任編輯 陳呈超