王其昂， 吳子燕， 劉 露

(西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院 西安，710129)

?

基于Hilbert-Huang變換與理想帶通濾波器的系統(tǒng)識(shí)別*

王其昂， 吳子燕， 劉 露

(西北工業(yè)大學(xué)力學(xué)與土木建筑學(xué)院 西安，710129)

針對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解存在模態(tài)混疊現(xiàn)象，提出基于Hilbert-Huang變換與理想帶通濾波器的系統(tǒng)識(shí)別方法。該方法利用傅里葉變換得到結(jié)構(gòu)加速度響應(yīng)頻響函數(shù)，粗略估計(jì)固有頻率范圍，通過(guò)半功率帶寬法設(shè)計(jì)理想帶通濾波器，定量化確定通帶帶寬，使信號(hào)在經(jīng)過(guò)濾波器后頻域內(nèi)零相移，同時(shí)不改變其幅值譜。結(jié)構(gòu)響應(yīng)通過(guò)指定頻帶的理想帶通濾波器產(chǎn)生若干窄帶信號(hào)，利用經(jīng)驗(yàn)?zāi)B(tài)分解獲取結(jié)構(gòu)模態(tài)響應(yīng)，經(jīng)Hilbert變換構(gòu)造模態(tài)響應(yīng)解析信號(hào)，并通過(guò)線性最小二乘擬合提取結(jié)構(gòu)模態(tài)參數(shù)與物理參數(shù)。結(jié)果表明：半功率帶寬法可實(shí)現(xiàn)帶通濾波器頻帶的定量化設(shè)計(jì),理想帶通濾波器的零相移特點(diǎn)較好契合Hilbert-Huang變換用于系統(tǒng)識(shí)別的要求，兩者結(jié)合可有效地解決模態(tài)混疊現(xiàn)象，減少虛假模態(tài)，大大提高結(jié)構(gòu)系統(tǒng)識(shí)別精度。

系統(tǒng)識(shí)別; Hilbert-Huang變換; 理想帶通濾波器; 經(jīng)驗(yàn)?zāi)B(tài)分解; 半功率帶寬法; 模態(tài)響應(yīng)

引 言

工程結(jié)構(gòu)系統(tǒng)識(shí)別與損傷診斷是結(jié)構(gòu)健康檢測(cè)的核心技術(shù)，處于該學(xué)科研究前沿。其首要任務(wù)是從測(cè)得數(shù)據(jù)中確定結(jié)構(gòu)模態(tài)參數(shù)，從整體上描述系統(tǒng)動(dòng)力特性，進(jìn)一步辨識(shí)結(jié)構(gòu)物理參數(shù)，為故障診斷、可靠性評(píng)估提供理論依據(jù)。

系統(tǒng)識(shí)別和損傷檢測(cè)Benchmark模型在哥倫比亞大學(xué)UBC實(shí)驗(yàn)室建立[1]，并發(fā)展了卡爾曼濾波、小波變換、子空間識(shí)別等一系列方法[2-4]。小波變換[5-6]具有良好的時(shí)頻分辨能力，避免傅里葉變換的Gibbs效應(yīng)，利用時(shí)域脈沖響應(yīng)提取結(jié)構(gòu)模態(tài)參數(shù)，但其本質(zhì)上為線性變換，不能有效處理非線性問(wèn)題?；诮?jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition, 簡(jiǎn)稱(chēng)EMD)的希爾伯特黃變換(Hilbert-Huang transform, 簡(jiǎn)稱(chēng)HHT)[7-8]，具有完全自適應(yīng)性，不受Heisenberg測(cè)不準(zhǔn)原理制約等優(yōu)點(diǎn)，能處理非線性、非平穩(wěn)信號(hào)。當(dāng)結(jié)構(gòu)兩階固有頻率相似，基于Hilbert-Huang變換的模態(tài)參數(shù)識(shí)別較小波變換更為準(zhǔn)確。但EMD方法存在模態(tài)混疊，限制了它在實(shí)際中的應(yīng)用，其他學(xué)者提出改進(jìn)的Hilbert-Huang方法[9-10]，利用傅里葉變換，估計(jì)固有頻率范圍，然后讓信號(hào)通過(guò)指定頻帶的帶通濾波器，改善模態(tài)混疊，但未涉及濾波器帶寬定量化設(shè)計(jì)。 Yang等[11]利用Hilbert-Huang變換識(shí)別多自由度體系模態(tài)參數(shù)，并指出該方法用于系統(tǒng)識(shí)別中，帶通濾波器相移應(yīng)盡可能小。

由此，筆者提出基于Hilbert-Huang變換與理想帶通濾波器的系統(tǒng)識(shí)別方法，利用傅里葉變換得到結(jié)構(gòu)加速度頻響函數(shù)，估計(jì)固有頻率范圍。為避免模態(tài)混疊，減少虛假模態(tài)響應(yīng)，在使用EMD方法進(jìn)行模態(tài)分解前，使信號(hào)通過(guò)指定頻帶的理想帶通濾波器，將寬頻信號(hào)分解為若干窄帶信號(hào)。設(shè)計(jì)理想帶通濾波器，使結(jié)構(gòu)響應(yīng)經(jīng)濾波器后頻域內(nèi)零相移，同時(shí)不改變其幅值譜，并使用半功率帶寬法定量化設(shè)計(jì)頻帶范圍，減少模態(tài)混疊，提高結(jié)構(gòu)系統(tǒng)識(shí)別精度。

1 脈沖激勵(lì)下模態(tài)響應(yīng)

多自由度系統(tǒng)外力作用下的動(dòng)力學(xué)方程為

(1)

其中：X(t)=[x1,x2, …,xn]T為系統(tǒng)位移向量；M，C，K分別為質(zhì)量、阻尼和剛度矩陣；p(t)為外部激勵(lì)。

通過(guò)模態(tài)變換，位移、加速度響應(yīng)分解為n階實(shí)模態(tài)，考慮模態(tài)正交性可得

(2)

其中：ψj為第j階模態(tài)向量；qj(t)為模態(tài)坐標(biāo)；ωj，ζj和mj分別為第j階模態(tài)頻率、模態(tài)阻尼和模態(tài)質(zhì)量。

若有以一激勵(lì)作用在第r個(gè)自由度上，即當(dāng)j=r時(shí)，pr(t)=p0δ(t)；當(dāng)j≠r時(shí)，pj(t)=0。pj(t)為激勵(lì)向量第j個(gè)元素。由此可知，第j個(gè)模態(tài)坐標(biāo)下加速度響應(yīng)為

(3)

因此，結(jié)構(gòu)在第s(s= 1, 2, …,n)自由度上加速度脈沖響應(yīng)可以表示為

(4)

(5)

(6)

其中：φsj,r為第j階模態(tài)向量第s，r個(gè)元素相位差。

(7)

2 結(jié)構(gòu)響應(yīng)Hilbert-Huang變換

Hilbert-Huang變換[7]由EMD方法和Hilbert變換組成，其核心思想是將時(shí)間序列通過(guò)EMD分解成數(shù)個(gè)本征模函數(shù)(intrinsic mode function, 簡(jiǎn)稱(chēng)IMF)，之后利用Hilbert變換構(gòu)造解析信號(hào)，計(jì)算得到信號(hào)瞬時(shí)頻率和振幅。

2.1 EMD概述

EMD將一個(gè)復(fù)雜的信號(hào)分解為若干個(gè)IMF之和，每一個(gè)IMF均具有相同極值點(diǎn)與過(guò)零點(diǎn)數(shù)目，相鄰零點(diǎn)之間僅一個(gè)極值點(diǎn)，且上下包絡(luò)線關(guān)于時(shí)間軸對(duì)稱(chēng)。對(duì)于給定的信號(hào)x(t)(如加速度響應(yīng))，EMD基本步驟[7]如下：a.確定信號(hào)所有局部極值點(diǎn);b.用三次樣條曲線構(gòu)造信號(hào)上(極大值點(diǎn))、下(極小值點(diǎn))包絡(luò)線，分別為u(t)和l(t)；c.計(jì)算上下包絡(luò)線的平均值，記為m1(t)，求出x(t)-m1(t)=h1(t)，若h1(t)滿(mǎn)足IMF特性，即為第1個(gè)IMF分量;d. 若h1(t)不滿(mǎn)足IMF的條件，則將其視為原始數(shù)據(jù),重復(fù)步驟a～c，得到上下包絡(luò)線平均值m11(t)，判斷h11(t)=h1(t)-m11(t)是否滿(mǎn)足IMF的條件，如不滿(mǎn)足，則循環(huán)k次，得到h1k(t)=h1(k-1)(t)-m1k(t)，使得h1k(t)滿(mǎn)足IMF的條件,記c1(t)=h1k(t)，為信號(hào)x(t)的第1個(gè)IMF分量;e.將c1(t)從x(t)中分離，得到r1(t)=x(t)-c1(t)，將r1(t)作為原始數(shù)據(jù)重復(fù)以上過(guò)程，直至分離出全部IMF，將x(t)分解為n個(gè)IMF與殘余分量之和

(8)

2.2 半功率帶寬法設(shè)計(jì)理想帶通濾波器

為減少虛假模態(tài)，在使用EMD方法進(jìn)行模態(tài)分解前，首先使信號(hào)通過(guò)不同頻段的帶通濾波器，將寬頻信號(hào)分解為若干窄帶信號(hào)。筆者設(shè)計(jì)了理想帶通濾波器，使結(jié)構(gòu)響應(yīng)經(jīng)濾波器后頻域內(nèi)零相移，同時(shí)不改變其幅值譜，提高模態(tài)參數(shù)識(shí)別精度。

理想帶通濾波器具有完全平坦的通帶，通帶內(nèi)沒(méi)有放大或衰減，通帶外頻率被完全衰減，且通帶外的轉(zhuǎn)換在極小頻率范圍完成。具體設(shè)計(jì)思路如下：將加速度響應(yīng)通過(guò)傅里葉變換得頻響函數(shù)，通帶范圍內(nèi)頻響函數(shù)數(shù)據(jù)不變，通帶外變?yōu)榱?，再通過(guò)逆傅里葉變換轉(zhuǎn)化為時(shí)域數(shù)據(jù)，獲得窄帶加速度響應(yīng)。根據(jù)半功率帶寬法，保留信號(hào)峰值處多數(shù)能量，理想帶通濾波器上、下限截止頻率設(shè)計(jì)為幅值譜半功率點(diǎn)處，其帶寬的功率譜密度比峰值低3 dB。

2.3 模態(tài)響應(yīng)Hilbert變換及其解析信號(hào)

(9)

對(duì)于cj(t)的解析信號(hào)zj(t)為

(10)

其中：Aj(t)和θj(t)分別為IMF的瞬時(shí)振幅和瞬時(shí)相位。

(11)

(12)

由瞬時(shí)相位計(jì)算信號(hào)瞬時(shí)頻率

(13)

信號(hào)瞬時(shí)振幅和瞬時(shí)相位均為時(shí)間函數(shù)。測(cè)量數(shù)據(jù)x(t)首先應(yīng)用EMD分解成單一頻率或窄帶信號(hào)IMF，經(jīng)Hilbert變換可得到有意義的結(jié)果，在任意時(shí)刻t僅有單一頻率，使得瞬時(shí)頻率具有實(shí)際物理意義。對(duì)于常規(guī)信號(hào)，瞬時(shí)頻率隨時(shí)間改變而非單一頻率。

加速度模態(tài)響應(yīng)由式(5)獲得，其Hilbert變換可由Bedrosian定理計(jì)算

(14)

其中

(15)

(16)

由式(10)得模態(tài)響應(yīng)解析信號(hào)為

(17)

若阻尼較小、固有頻率較大，模態(tài)響應(yīng)Hilbert變換式(14)簡(jiǎn)化為

(18)

瞬時(shí)幅值、瞬時(shí)相位公式簡(jiǎn)化為

(19)

(20)

3 系統(tǒng)識(shí)別

3.1 固有頻率與阻尼比的識(shí)別

3.2 模態(tài)振型及物理參數(shù)識(shí)別

理論上固有頻率、阻尼比的識(shí)別，僅需獲取某一自由度上加速度響應(yīng),但對(duì)于模態(tài)振型、質(zhì)量、剛度和阻尼矩陣的識(shí)別則需測(cè)得所有自由度上的響應(yīng)。由式(6)可得到模態(tài)元素絕對(duì)值之比

(21)

兩模態(tài)元素相位差通過(guò)式(22)可得到

(22)

從式(21)和式(22)可以看出，將t0作為所取時(shí)間段的中間某點(diǎn)，得到θsj和θrj在t0處的擬合值，進(jìn)而可以確定ψsj相對(duì)于ψrj的相位。根據(jù)式(7)，確定ψsj相對(duì)于ψrj符號(hào)(正或負(fù))。由此，在j階模態(tài)向量中的所有元素相對(duì)于某一特定元素的絕對(duì)值及符號(hào)都可確定。

得到模態(tài)參數(shù)后，結(jié)構(gòu)質(zhì)量、剛度和阻尼矩陣便可得到。將式(19)、式(20)帶入式(6)中，第j階模態(tài)質(zhì)量為

(23)

其中:p0為可測(cè)沖擊力荷載。

模態(tài)剛度ki、阻尼ci為

(24)

根據(jù)模態(tài)正交性得結(jié)構(gòu)質(zhì)量、剛度及阻尼陣為

(25)

其中:Φ為模態(tài)向量矩陣。

4 算例分析

圖1 結(jié)構(gòu)系統(tǒng)識(shí)別流程圖Fig.1 Flow chart for the structural system identification

圖2 4自由度剪切梁模型Fig.2 Four degrees of freedom shear-beam structure

筆者以某4層剪切梁模型為算例，進(jìn)行實(shí)驗(yàn)?zāi)B(tài)分析，并與理論分析結(jié)果作對(duì)比，以驗(yàn)證上述系統(tǒng)識(shí)別方法的有效性及準(zhǔn)確性。具體流程見(jiàn)圖 1。4層剪切梁模型見(jiàn)圖 2，樓層質(zhì)量均為4.989 5 kg，樓層剛度k1，k2，k3與k4分別為1 401，1 576，1 226與1 051 N/m，阻尼類(lèi)型為瑞利阻尼。由質(zhì)量矩陣和剛度矩陣按比例組合構(gòu)造而成，質(zhì)量比例系數(shù)Alpha為-0.001 6，剛度比例系數(shù)Beta為3.918 8×10-4。在第4樓層施加單位脈沖激勵(lì)，構(gòu)造系統(tǒng)狀態(tài)空間函數(shù)公式。利用Matlab系統(tǒng)仿真獲取各樓層加速度，而實(shí)際上任何實(shí)測(cè)數(shù)據(jù)都不可避免受到噪聲信號(hào)的影響,因此本研究在響應(yīng)信號(hào)中添加相同長(zhǎng)度的高斯白噪聲，信噪比為26.02 dB。

該系統(tǒng)4階自振頻率及對(duì)應(yīng)阻尼比均可由任一樓層加速度脈沖響應(yīng)獲取。圖 3為頻域內(nèi)加速度脈沖響應(yīng)幅值譜，響應(yīng)中主要包含4階模態(tài)的振動(dòng)成分，固有頻率大致在0.925,2.475,3.85和4.95 Hz左右。

圖3 脈沖響應(yīng)頻域幅值譜Fig.3 Frequency spectra of impulse responses

圖4 帶通濾波器相移圖Fig.4 Phase shift of different band-pass filters

圖5 4階模態(tài)響應(yīng)Fig.5 Four modal responses

圖6 解析信號(hào)相位圖Fig.6 Phase plot of the analytical signal

圖7 解析信號(hào)幅值圖Fig.7 Amplitude plot of the analytical signal

由固有頻率和阻尼比識(shí)別步驟3得第1階模態(tài)頻率、阻尼。同時(shí)，對(duì)于其余3層響應(yīng)可得第1階模態(tài)響應(yīng)，通過(guò)Hilbert-Huang變換理論計(jì)算可得第1階模態(tài)參數(shù)，取均值得第1階固有頻率、阻尼比。重復(fù)同樣步驟得其余3階模態(tài)頻率、阻尼比，并與Butterworth濾波器識(shí)別結(jié)果相比較，結(jié)果如表 1所示，由式(7)、式(21)確定模態(tài)向量。

表1 結(jié)構(gòu)模態(tài)參數(shù)識(shí)別結(jié)果

模態(tài)頻率識(shí)別，兩種濾波器估算結(jié)果一致，最大誤差均小于0.2%。模態(tài)阻尼識(shí)別對(duì)數(shù)值誤差較敏感，尚缺乏較為理想的識(shí)別方法，理想帶通濾波器(最大誤差23%)略?xún)?yōu)于Butterworth濾波器(最大誤差28%)。模態(tài)向量識(shí)別結(jié)果與理論值相比得模態(tài)置信度(modal assurance criterion，簡(jiǎn)稱(chēng)MAC)值(見(jiàn)表 1)。4階MAC值均大于0.99(實(shí)踐中，兩模態(tài)向量MAC值大于0.9時(shí)為強(qiáng)相關(guān))，基于理想帶通濾波器的識(shí)別結(jié)果與理論振型基本重合，模態(tài)向量識(shí)別精度高。

與常規(guī)Butterworth帶通濾波器相比，理想帶通濾波略提高模態(tài)參數(shù)識(shí)別精度。由于誤差累積效應(yīng)，這一提高對(duì)于結(jié)構(gòu)物理參數(shù)的識(shí)別至關(guān)重要。1～4樓層質(zhì)量識(shí)別結(jié)果分別為4.941，5.047，4.892和4.999kg，最大誤差為1.95%，而B(niǎo)utterworth濾波器最大誤差為11.83%。1～4樓層剛度識(shí)別結(jié)果分別為1 323.1，1 654.0，1 166.4和1 037.1，最大誤差為5.56%，Butterworth濾波器最大誤差為19.31%。阻尼矩陣主對(duì)角元素識(shí)別結(jié)果分別為Cd(1,1)=1.08，Cd(2,2)=0.98，Cd(3,3)=0.84，Cd(4,4)=0.41，最大誤差為10%，常規(guī)濾波器識(shí)別結(jié)果誤差為43.1%。為驗(yàn)證該方法識(shí)別靈敏度，在原有算例基礎(chǔ)上分別修改質(zhì)量、剛度矩陣。樓層質(zhì)量修改為5.2 kg，以剛度識(shí)別結(jié)果為例，理想帶通濾波器最大誤差為3.7%，而B(niǎo)utterworth濾波器識(shí)別結(jié)果最大誤差為11.4%。修改1～4樓層剛度，理想帶通濾波器識(shí)別結(jié)果同樣優(yōu)于常規(guī)濾波器。由此，對(duì)于結(jié)構(gòu)物理參數(shù)識(shí)別，理想帶通濾波器處理后有明顯優(yōu)勢(shì)，大大提高了識(shí)別精度。

5 結(jié) 論

1) 筆者提出基于改進(jìn)的Hilbert-Huang變換的系統(tǒng)識(shí)別方法，利用傅里葉變換得到結(jié)構(gòu)加速度頻響函數(shù)，粗略估計(jì)固有頻率范圍。為減少虛假模態(tài)，使用EMD模態(tài)分解前，使信號(hào)通過(guò)指定頻帶的理想帶通濾波器，將寬頻信號(hào)分解為若干窄帶信號(hào)，構(gòu)造模態(tài)響應(yīng)解析信號(hào)，提取模態(tài)參數(shù)與物理參數(shù)。

2) 理想帶通濾波器的零相移特點(diǎn)較好契合Hilbert-Huang變換用于系統(tǒng)識(shí)別的要求，兩者結(jié)合有效地解決了EMD中模態(tài)混疊問(wèn)題。

3) 獲取加速度響應(yīng)幅值譜，由半功率帶寬法，確定理想帶通濾波器上、下限截止頻率為幅值譜半功率點(diǎn)處，實(shí)現(xiàn)頻帶定量化設(shè)計(jì)。

4) 與常規(guī)Butterworth帶通濾波器相比，理想帶通濾波略微提高模態(tài)參數(shù)識(shí)別精度。但由于誤差累積效應(yīng)，這一提高對(duì)于結(jié)構(gòu)物理參數(shù)的識(shí)別至關(guān)重要，大大提高了結(jié)構(gòu)質(zhì)量、剛度及阻尼矩陣的識(shí)別精度。

[1] Black C J, Ventura C E. Blind test on damage detection of a steel frame structure[C]∥Proceedings of SPIE. Santa Barbara CA: Society for Experimental Mechanics, 1998: 623-629.

[2] 任宜春，易偉建. 結(jié)構(gòu)物理參數(shù)識(shí)別的多尺度參數(shù)卡爾曼濾波方法[J]. 工程力學(xué)，2008，25(5)：1-5.

Ren Yichun, Yi Weijian. Identification of physical parameters by multi-scale parameter Kalman filter[J]. Engineering Mechanics, 2008, 25(5): 1-5. (in Chinese)

[3] 史治宇，沈林. 基于小波方法的時(shí)變動(dòng)力系統(tǒng)參數(shù)識(shí)別[J]. 振動(dòng)、測(cè)試與診斷，2008，28(2): 108-112.

Shi Zhiyu, Shen Lin. Parameter identification of linear time-varying dynamical system based on wavelet method[J]. Journal of Vibration, Measurement & Diagnosis, 2008, 28(2): 108-112. (in Chinese)

[4] 徐良，江見(jiàn)鯨，過(guò)靜珺. 隨機(jī)子空間識(shí)別在懸索橋?qū)嶒?yàn)?zāi)B(tài)分析中的應(yīng)用[J]. 工程力學(xué)，2002，19(4): 46-49.

Xu Liang, Jiang Jianjing, Guo Jingjun. Application of stochastic subspace method to experimental modal analysis of suspension bridges[J]. Engineering Mechanics, 2002, 19(4): 46-49. (in Chinese)

[5] Vafaei M, Alih S C, Abd Rahman A B, et al. A wavelet-based technique for damage quantification via mode shape decomposition[J]. Structure and Infrastructure Engineering, 2015, 11(7): 869-883.

[6] Dziedziech K, Staszewski W J, Uhl T. Wavelet-based modal analysis for time-variant systems[J]. Mechanical Systems and Signal Processing, 2015, 50-51: 323-337.

[7] Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Oyal Society of London Proceedings, 1998, 454(1971): 903-995.

[8] Lee J, Hussain S H, Wang S, et al. Enhanced method to reconstruct mode shapes of continuous scanning measurements using the Hilbert Huang transform and the modal analysis method[J]. Review of Scientific Instruments, 2014, 85(9): 1-11.

[9] 付春，姜紹飛. 基于改進(jìn)EMD-ICA的結(jié)構(gòu)模態(tài)參數(shù)識(shí)別研究[J]. 工程力學(xué)，2013，30(10)：199-204.

Fu Chun, Jiang Shaofei. Study on the modal parameter identification based on improved EMD and independent component analysis[J]. Engineering Mechanics, 2013, 30(10): 199-204. (in Chinese)

[10]秦毅，秦樹(shù)人，毛永芳. 改進(jìn)的Hilbert-Huang變換在信號(hào)瞬態(tài)特征提取中的應(yīng)用[J]. 振動(dòng)與沖擊， 2008，27(11)：129-133.

Qin Yi, Qin Shuren, Mao Yongfang. Application of improved Hilbert-Huang transform in transient feature extraction[J]. Journal of Vibration and Shock, 2008, 27(11): 129-133. (in Chinese)

[11]Yang J N, Lei Ying, Pan Shuwen, et al. System identification of linear structures based on Hilbert-Huang spectral analysis. part 1: normal modes[J]. Earthquake Engineering and Structural Dynamics, 2003, 32(9): 1443-1467.

[12]Shi Z Y, Law S S, Xua X. Identification of linear time-varying mdof dynamic systems from forced excitation using Hilbert transform and EMD method[J]. Journal of Sound and Vibration, 2009, 321(3-5): 572-589.

10.16450/j.cnki.issn.1004-6801.2016.06.005

*國(guó)家自然科學(xué)基金資助項(xiàng)目(51278420)；博士創(chuàng)新基金資助項(xiàng)目(CX201408)

2015-06-08；

2015-07-13

TB122; TB123

王其昂,男，1986年8月生，博士生。主要研究方向?yàn)榻Y(jié)構(gòu)健康檢測(cè)、系統(tǒng)識(shí)別及可靠性評(píng)估。曾發(fā)表《Seismic fragility analysis of highway bridges considering multi-dimensional performance limit state》(《Earthquake Engineering and Engineering Vibration》2012,Vol.11,No.2)等論文。E-mail:qawang@mail.nwpu.edu.cn