室內(nèi)聲場重現(xiàn)中的時域去卷積方法
2016-04-15 11:26:02鄭四發(fā)廖祥凝連小珉
振動與沖擊 2016年4期
彭 博， 鄭四發(fā)， 廖祥凝， 連小珉
(清華大學(xué) 汽車安全與節(jié)能國家重點實驗室，北京 100084)
?
室內(nèi)聲場重現(xiàn)中的時域去卷積方法
彭博， 鄭四發(fā)， 廖祥凝， 連小珉
(清華大學(xué) 汽車安全與節(jié)能國家重點實驗室，北京100084)
摘要：為了在帶反射邊界的普通房間內(nèi)重現(xiàn)聲場，需要對原始聲場信號進(jìn)行去卷積處理，以獲得揚聲器驅(qū)動信號。建立了房間內(nèi)聲場重現(xiàn)的時域模型，在模型中采用長時長的去卷積網(wǎng)絡(luò)，以獲得高精度同時抑制了時域混疊效應(yīng)，再針對卷積濾波的海量計算問題使用了快速濾波算法，提高了計算效率。汽車車內(nèi)加速噪聲的重現(xiàn)實驗結(jié)果表明：所提出的時域模型可以高效準(zhǔn)確地完成聲壓信號的去卷積化，實現(xiàn)普通房間內(nèi)聲場的定量重現(xiàn)。
關(guān)鍵詞：聲場重現(xiàn)；非自由場；去卷積；濾波器設(shè)計；快速卷積
聲場重現(xiàn)技術(shù)通過準(zhǔn)確地恢復(fù)原始聲場的聲學(xué)參數(shù)，讓聽音者可以擁有身臨其境的感受?；趽P聲器的聲場定量重現(xiàn)系統(tǒng)為機(jī)械噪聲的聲品質(zhì)主觀評價提供了一種新的手段。
在聲場重現(xiàn)方法中，Ambisonics方法和Wave Field Synthesis(WFS)方法受到了廣泛地研究。Ambisonics方法產(chǎn)生于20世紀(jì)70年代，主要用在圓形布置或球形布置揚聲器陣列的聲場重現(xiàn)中[1-2]。該方法的核心思想是利用重現(xiàn)聲場與原始聲場模態(tài)匹配的方法來實現(xiàn)聲場重現(xiàn)，因此往往利用球諧函數(shù)或柱諧函數(shù)來表達(dá)聲場[3-4]。WFS方法最早由Berkhout等[5-6]提出，該方法通過控制連續(xù)布置于無源區(qū)域邊界的單極子和偶極子聲源實現(xiàn)區(qū)域內(nèi)聲場的控制。由于水平面內(nèi)的聲場重現(xiàn)比豎直面內(nèi)重要[7]，該方法主要研究利用平面布置的揚聲器陣列對水平面內(nèi)聲場進(jìn)行重現(xiàn)的問題。
以上兩種理論方法都受到儀器設(shè)備及實際場景帶來的限制。Ambisonics方法在原始聲場采樣時需要用到大量的傳聲器，且不規(guī)則的揚聲器陣列布局會帶來數(shù)值求解不穩(wěn)定和次優(yōu)解的問題[8]。WFS方法需要考慮實際揚聲器的特性，也需要考慮揚聲器的布置。另外，以上兩種方法的實驗大都在消聲室內(nèi)進(jìn)行，未考慮普通房間的反射對重現(xiàn)精度的影響。
為了避免上述兩種方法在應(yīng)用中的限制，通過匹配聲場中控制點聲壓以直接在頻域內(nèi)計算揚聲器驅(qū)動信號的聲場重現(xiàn)方法近期得到了研究。Kolundzija等[9]通過自由場仿真說明了此類方法可以比波束合成方法具有更高的精度。Gauthier等[10]通過此類方法在機(jī)艙模型內(nèi)完成了飛機(jī)噪聲的重現(xiàn)工作。將此類方法應(yīng)用于非自由場時，濾波器的設(shè)計方法對精度有重要影響，濾波器長度一般需要足夠長，以減小時域混疊現(xiàn)象。針對實際房間中的聲場重現(xiàn)，本文提出了聲場重現(xiàn)的時域模型，該模型采用了長時長的去卷積網(wǎng)絡(luò)，在獲得高精度的同時抑制了時域混疊效應(yīng)，接著針對運算量過大的問題采用了快速濾波算法，實現(xiàn)了聲壓信號準(zhǔn)確高效的時域去卷積化。
1非自由場內(nèi)的聲場重現(xiàn)
1.1聲場重現(xiàn)的時域模型
在聲場重現(xiàn)時，對原始聲場采樣得到的時域聲壓信號通過去卷積網(wǎng)絡(luò)的處理后獲得揚聲器驅(qū)動信號，再通過房間內(nèi)多個揚聲器聲場的合成實現(xiàn)原始聲場的重現(xiàn)，該過程如圖1所示。d(t)=[d1(t),d2(t),…,dM(t)]T為原始聲場中采集到的M個聲壓信號，[·]T表示轉(zhuǎn)置。對該信號通過時域去卷積網(wǎng)絡(luò)r(t)處理(r(t)為L×M維矩陣)，獲得揚聲器的驅(qū)動信號g(t)=[g1(t),g2(t),…,gL(t)]T。
圖1 聲場重現(xiàn)的時域模型Fig.1 Time-domain model of sound field reproduction
為檢測重現(xiàn)效果，在目標(biāo)區(qū)域內(nèi)布置M個傳聲器，其布置方式與原始聲場采樣時相同。傳聲器之間的最小距離決定了能夠有效重現(xiàn)的最高頻率[9]。為考慮聽音者對聲場的影響，在對原始聲場和重現(xiàn)聲場采樣時將人工頭置于傳聲器陣列的中心。記傳聲器采集到的重現(xiàn)聲壓信號為p(t)=[p1(t),p2(t),…,pM(t)]T。從獲得揚聲器驅(qū)動信號g(t)到采集到重現(xiàn)聲壓信號p(t)的這個過程，由一個典型的線性多輸入多輸出系統(tǒng)(即放音系統(tǒng))來保證，這個過程在時域可用以下卷積關(guān)系表示:
p(t)=h(t)*g(t)
(1)
式中:h(t)為M×L維矩陣，是放音系統(tǒng)的無限沖激響應(yīng)函數(shù)矩陣，由放音系統(tǒng)、聽音者和房間的特性決定；*為卷積運算符。從d(t)到g(t)的去卷積處理也可用卷積關(guān)系表示：
g(t)=r(t)*d(t)
(2)
為使p(t)為d(t)的重現(xiàn)信號，即p(t)=d(t)，應(yīng)滿足以下關(guān)系：
h(t)*r(t)=diag[δ1(t),δ2(t),…,δM(t)]
(3)
式中:diag[·]為由各對角線元素表示的對角矩陣，δm(t)為無限沖激函數(shù)(m=1,2,…,M)。式(3)為去卷積方程，通過求得r(t)使式(3)成立。
1.2時域去卷積網(wǎng)絡(luò)的結(jié)構(gòu)
聲場重現(xiàn)的時域模型中最關(guān)鍵的環(huán)節(jié)是去卷積網(wǎng)絡(luò)r(t)，d(t)經(jīng)過去卷積網(wǎng)絡(luò)的處理以獲得揚聲器的驅(qū)動信號，r(t)與d(t)的卷積過程可用以下框圖表示：
圖2 時域去卷積網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.2 The structure of time-domain deconvolution network
對于揚聲器l而言，其驅(qū)動信號gl(t)是d(t)與M個去卷積濾波器卷積結(jié)果的和，可表示為：
(4)
式中，rlm(τ)的時間范圍為-T≤τ≤T，長度為2T。計算時直接對采集到的離散聲壓信號進(jìn)行處理，計算公式可表示為：
(5)
式中，tk、τj為離散時間，離散時間間隔皆由采樣頻率決定，k和j為離散時間序數(shù)。Δt為離散時間間隔。rlm(τj)的離散數(shù)據(jù)點數(shù)為2N。
1.3重現(xiàn)的精度和穩(wěn)定性
本文利用平均聲壓級誤差來衡量聲場重現(xiàn)的精度，該指標(biāo)定義為重現(xiàn)聲場與原始聲場聲壓級差的絕度均值，如下式所示：
(6)
式中，SPL[·]表示求取時域聲壓級。另外，本文利用平均電聲比(EPR)來衡量重現(xiàn)方法的穩(wěn)定性，EPR定義為揚聲器驅(qū)動信號時域有效均值與傳聲器聲壓時域有效均值的比值，其定義式如下：
(7)
2時域去卷積濾波網(wǎng)絡(luò)的求解
2.1病態(tài)問題的處理
在頻域內(nèi)，放音系統(tǒng)的特性可用傳遞函數(shù)來描述。記式(1)中g(shù)(t)和p(t)的傅里葉變換結(jié)果分別為G(f)=[G1(f),G2(f),…,GL(f)]T和P(f)=[P1(f),P2(f),…,PL(f)]T(其中f為頻率)。它們之間的關(guān)系用放音系統(tǒng)的傳遞函數(shù)H(f)可描述為：
P(f)=H(f)G(f)
(8)
式中，H(f)為M×L維復(fù)數(shù)矩陣，是無限沖激響應(yīng)函數(shù)矩陣在頻域的表達(dá)形式，可由線性掃頻法[12]直接測量得到。實際情況下H(f)往往呈病態(tài)性[9-10,13-14]，本文采用Tikhonov正則化方法[13]和奇異值截斷法[9]來求解去卷積網(wǎng)絡(luò)頻域特性。
對于上述兩種方法來說，門限值的選取是關(guān)鍵。無論M和L的大小關(guān)系如何，H(f)存在以下奇異值分解形式：
(9)
式中，U為M×M維酉矩陣，V為L×L維酉矩陣。Σ=diag[σ1,σ2,…,σD]為M×L維對角矩陣。當(dāng)M
λ=εσ1
(10)
式中，0<ε<1。對于Tikhonov方法來說，此種選取方式給對角矩陣Σ的對角元素帶來的增量為最大奇異值的ε倍；對奇異值截斷法來講，以最大奇異值為標(biāo)準(zhǔn)截斷了貢獻(xiàn)率小于其ε倍的成分。這種門限值的選取方法在一定程度上量化了正則化的影響。
因此根據(jù)上述思路，由奇異值截斷法求出的頻域特性可表示為：
(11)
式中，C為正整數(shù)。
基于Tikhonov方法求得的頻域特性可表示為：
(12)
2.2時域去卷積濾波器設(shè)計
去卷積網(wǎng)絡(luò)的頻域特性R(f)可寫成如下形式：
(13)
式(13)中的Rlm(f)為圖2中rlm(t)的頻域特性。Rlm(f)理論上應(yīng)定義在(-∞,+∞)的連續(xù)函數(shù)，其模關(guān)于x軸偶對稱，其相位關(guān)于x軸奇對稱。由連續(xù)傅里葉逆變換，可求得理想的時域去卷積濾波器rlm(t)：
(14)
式中，i2=-1。rlm(t)也定義在(-∞,+∞)上，而實際上不可能求解連續(xù)的Rlm(f)和rlm(t)，需要先進(jìn)行離散化。首先令rlm(t)的離散形式為rlm(τj)，采樣頻率fs，離散時間間隔Δt=1/fs。令實測傳遞函數(shù)的最低頻率和最高頻率分別為fL和fM(fs>2fM)，再令Rlm(f)的頻率離散間隔為Δf。Δf決定了濾波器rlm(τj)的時間長度(2T=1/Δf)。將Rlm(f)離散為Rlm(fn)后，根據(jù)2.1所述方法計算出Rlm(fn)在區(qū)間(fL,,fM)上的值后，利用快速傅里葉變換(FFT)可實現(xiàn)式(14)的離散快速計算。首先根據(jù)Rlm(fn)在區(qū)間(fL,fM)上的值在頻段[0,fs)內(nèi)構(gòu)造出奈奎斯特頻譜，記為B(fn)：
(15)
式中(·)*表示共軛，B(fn)的離散點數(shù)設(shè)置為偶數(shù)記為2N，離散頻率序數(shù)依次為0,1,…,2N-1。再將B(fn)轉(zhuǎn)換到時域：
b(τj)=Δf{FFT[B*(fn)]}*
(16)
b(τj)的離散時間序數(shù)為j=0,1,…,2N-1?？紤]到FFT的周期性，對b(τj)進(jìn)行時序調(diào)整后獲得rlm(τj)：
rlm(τj)=
(17)
式中:2T=1/Δf，為FFT變換中的時間周期，也是rlm(τj)的時間長度。此處rlm(τj)的時間范圍應(yīng)為-T≤τj≤T，離散時間序數(shù)為j=-N,-N+1,…,N-1。式(15)～式(17)為濾波器的快速計算方法。為了抑制rlm(τj)發(fā)生時域混疊現(xiàn)象需恰當(dāng)選擇濾波器的時間長度，即需要確定足夠小的Δf。
3快速濾波計算
當(dāng)時域去卷積濾波器的長度較長時，直接利用式(5)完成如圖2所示的卷積濾波運算將耗費大量時間。本文的實驗中M=18，L=21，fs=44 100 Hz，利用長度2 s的去卷積網(wǎng)絡(luò)對長度18 s的聲壓信號進(jìn)行濾波處理需要耗費大約20 h。為了提高運算效率，本文利用快速傅里葉變換來實現(xiàn)卷積的快速計算。式(5)中卷積運算的部分為ym(tk)，即令
(18)
式(18)為線卷積計算公式。設(shè)dm(tk)的長度為K(K>2N)，利用FFT實現(xiàn)ym(tk)的快速計算，如下式：
ym(tk)=iFFT[FFT[dm(tk),K]°FFT[rlm(τj),K]]
j=-N,-N+1,…,N-1;k=1,2,…,K
(19)
式中，iFFT[·]表示快速傅里葉逆變換。FFT[·,K]表示對數(shù)據(jù)進(jìn)行K點的快速傅里葉變換，若數(shù)據(jù)長度不足K點，則將數(shù)據(jù)后端添零至K點?！啊恪睘镠adamard乘積。式(19)中iFFT的結(jié)果點數(shù)為K點，其后端第2N到第K個數(shù)值(共K-2N+1個數(shù)值)為有效的運算結(jié)果，前2N-1個數(shù)值無意義。
4車內(nèi)加速噪聲的室內(nèi)重現(xiàn)
4.1實際噪聲的采集
為了驗證本文所述方法的有效性，對某汽車車內(nèi)加速噪聲進(jìn)行重現(xiàn)。將傳聲器陣列置于車內(nèi)測量加速時的噪聲，如圖3所示。傳聲器布置成環(huán)形包絡(luò)面，數(shù)量為M=18，平均間距Δ≈0.1 m，因此該陣列可有效的采集1 700 Hz以下的聲場。采集時采樣頻率為fs=44 100 Hz，第5個傳聲器(m=5)采集的聲信號的時頻圖如圖4所示，可以看出階次線的頻率隨時間遞增，噪聲的主要能量集中在低頻部分。
圖3 車內(nèi)加速噪聲的采集Fig.3 Interior noise measurement of an acceleration vehicle
圖4 車內(nèi)加速噪聲時頻圖Fig.4 Time-frequency map of the interior acceleration noise
4.2房間的布置
在重現(xiàn)實驗室內(nèi)，聲場重現(xiàn)的目標(biāo)區(qū)域處于房間中心(見圖5)，房間的墻壁內(nèi)部放置吸聲材料以改善房間的聲學(xué)特性。實驗室內(nèi)揚聲器的布置如圖6所示，所有揚聲器都指向目標(biāo)區(qū)域。為了獲得更寬頻帶的聲音共使用3種揚聲器，揚聲器的總數(shù)為L=21。
圖5 聲場重現(xiàn)實驗室圖Fig.5 Laboratory of sound field reproduction
圖6 實驗室內(nèi)揚聲器的布置Fig.6 Arrangement of loudspeakers in the laboratory
4.3實際時域去卷積濾波器設(shè)計
在求解時域去卷積網(wǎng)絡(luò)時，先測量房間的傳遞函數(shù)H(f)。當(dāng)l=17和m=9時測量的傳遞函數(shù)如圖7所示。由圖可見，由于房間的反射傳遞函數(shù)存在較大波動。
圖7 實測傳遞函數(shù)Fig.7 Measured transfer function
圖8顯示了傳遞函數(shù)的條件數(shù)，在噪聲能量較集中的100 Hz到300 Hz范圍內(nèi)，條件數(shù)已經(jīng)到達(dá)了相當(dāng)高的量級，分別利用Tikhonov正則化方法和奇異值截斷法處理H(f)的病態(tài)性，門限值選擇為0.01σ1。
圖8 傳遞函數(shù)的條件數(shù)Fig.8 Condition number of the measured transfer function
Δf的取值分別為1 Hz和0.5 Hz時，利用奇異值截斷法處理病態(tài)性后獲得的濾波器r10,10(t)如圖9所示，濾波器的長度分別為1 s和2 s。從局部放大圖中可以看出，當(dāng)Δf為1 Hz時設(shè)計出的濾波器邊緣還存在較大幅值，且該值與Δf為0.5 Hz時的不相同，這是由于此時濾波器出現(xiàn)了時域混疊；當(dāng)Δf為0.5 Hz時，濾波器邊緣的幅值為0，有效抑制了時域混疊。圖10為以上兩種情況下濾波器實際的頻域特性與目標(biāo)特性R10,10(f)的局部對比圖，從圖中可以看出，當(dāng)發(fā)生時域混疊時濾波器的特性發(fā)生畸變，而混疊較小的濾波器特性與目標(biāo)特性基本一致。也就是說，室內(nèi)聲場重現(xiàn)中由于反射的影響需要使用長時長的濾波器以抑制時域混疊。本例中Δf取0.5 Hz，此時對于兩種處理方法而言皆可有效抑制時域混疊。
圖9 利用不同Δf設(shè)計的濾波器Fig.9 Different filters designed with different Δf
圖10 濾波器的特性對比Fig.10 Property comparison of different filters
4.4濾波效率對比
式(18)和式(19)兩種方法進(jìn)行濾波時計算時間對比如表1所示(計算機(jī)CPU：i5 3230，2.6 GHz；計算軟件：Matlab)。原始聲場信號的總長度為18 s，用某一濾波器處理某通道聲壓信號時，直接濾波的耗時是快速濾波的2 130倍。若直接利用式(18)進(jìn)行濾波計算，進(jìn)行去卷積處理的總耗時約為20 h，快速濾波則只需要34 s。
表1 濾波耗時對比
4.5重現(xiàn)實驗結(jié)果
利用Tikhonov正則化方法和奇異值截斷法重現(xiàn)的聲場中第5個傳聲器的時頻圖分別如圖11和圖12所示。由于所采用的揚聲器不能夠發(fā)出低于40 Hz的低頻聲，因此不驗證低于40 Hz的聲場重現(xiàn)。對比圖4可知，兩種方法都很好的重現(xiàn)了原始聲場中的頻率成分。
圖11 Tikhonov方法重現(xiàn)聲場時頻圖Fig.11 Time-frequency map of the reproduction sound acquired by Tikhonov method
對原始聲場進(jìn)行40 Hz的高通濾波，然后計算聲場重現(xiàn)的平均聲壓級誤差，結(jié)果如圖13所示。從圖中可以看出，奇異值截斷方法與Tikhonov方法結(jié)果的精度相當(dāng)，大部分時間平均聲壓級誤差都在0.5 dB以下。在3.5 s左右兩種方法都有較大誤差，這是由于該時刻原始聲場劇烈變化所致(見圖4)。在整個時段，Tikhonov方法的重現(xiàn)結(jié)果較平穩(wěn)，TSVD方法的重現(xiàn)結(jié)果波動稍大。
圖12 TSVD方法重現(xiàn)聲場時頻圖Fig.12 Time-frequency map of the reproduction sound acquired by TSVD method
圖13 平均聲壓級誤差Fig.13 Averaged sound pressure level error
圖14 平均電聲比Fig.14 The ratio between averaged loudspeaker driving signal(RMS) and averaged reproduced sound pressure(RMS)
聲場重現(xiàn)的平均電聲比如圖14所示。從圖中可看出，雖然兩種方法的重現(xiàn)精度相當(dāng)，但整個時段奇異值截斷法的平均電聲比都為Tikhonov的兩倍，這意味著利用奇異值截斷法求出的揚聲器驅(qū)動信號幅值是Tikhonov方法的兩倍，由此說明針對本文所采集的樣本Tikhonov方法較奇異值截斷方法具有更好的重現(xiàn)穩(wěn)定性。
5結(jié)論
論文分析了聲場重現(xiàn)的時域模型，提出了利用長時長的時域濾波器對原始聲場的時域信號進(jìn)行濾波處理的方法，并采用快速卷積算法解決了耗時問題，進(jìn)而實現(xiàn)準(zhǔn)確高效地聲場重現(xiàn)。在這個方法的研究過程中可得出以下結(jié)論：
(1) 提出了時域去卷積濾波器的長度選擇方法及快速計算方法，有效地抑制了時域混疊現(xiàn)象。
(2) 提出的快速濾波算法提高了卷積濾波的運算效率，為去卷積方法的實際應(yīng)用提供了有利條件。
參 考 文 獻(xiàn)
[ 1 ] Gerzon M A. Periphony: with-height sound reproduction[J]. Journal of the Audio Engineering Society, 1972,21(1):2-10.
[ 2 ] Gerzon M A. Practical periphony: the reproduction of full-sphere sound[C]//Audio Engineering Society Convention 65. Audio Engineering Society, 1980.
[ 3 ] Poletti M A. Three-dimensional surround sound systems based on spherical harmonics[J]. Joutnal of the Audio Engineering Society, 2005,53(11):1004-1025.
[ 4 ] Daniel J, Moreau S. Further study of sound field coding with higher order ambisonics[C]//Audio Engineering Society Convention 116. Audio Engineering Society, 2004.
[ 5 ] Berkhout A J. A holographic approach to acoustic control[J]. Journal of the Audio Engineering Society, 1988,36(12):977-995.
[ 6 ] Berkhout A J, De Vries D, Vogel P. Acoustic control by wave field synthesis[J]. The Journal of the Acoustical Society of America, 1993,93(5):2764-2778.
[ 7 ] Blauert J. Spatial hearing: the psychophysics of human sound localization[M]. Cambridge, MA:MIT Press, 1997.
[ 8 ] Trevino J, Okamoto T, Iwaya Y, et al. High order Ambisonic decoding method for irregular loudspeaker arrays[C]//Proceedings of 20th International Congress on Acoustics,2010.
[ 9 ] Kolundzija M, Faller C, Vetterli M. Reproducing sound fields using MIMO acoustic channel inversion[J]. Journal of the Audio Engineering Society, 2011,59(10):721-734.
[10] Gauthier P A, Camier C, Gauthier O, et al. Aircraft sound environment reproduction: sound field reproduction inside a cabin mock-up using microphone and actuator arrays[C]//Proceedings of Meetings on Acoustics. Acoustical Society of America, 2013, 19(1): 055008.
[11] Nelson P A. Active control of acoustic fields and the reproduction of sound[J]. Journal of Sound and Vibration, 1994,177(4):447-477.
[12] Müller S, Massarani P. Transfer-function measurement with sweeps[J]. Journal of the Audio Engineering Society, 2001,49(6):443-471.
[13] Betlehem T, Withers C. Sound field reproduction with energy constraint on loudspeaker weights[J]. Audio, Speech, and Language Processing, IEEE Transactions on, 2012,20(8):2388-2392.
[14] Teal P D, Betlehem T, Poletti M A. An algorithm for power constrained holographic reproduction of sound[C]//Acoustics Speech and Signal Processing (ICASSP), 2010 IEEE International Conference on. Dallas,TX, 2010: 101-104.
Time domain deconvolution for reproducing sound field in an ordinary room
PENGBo,ZHENGSi-fa,LIAOXiang-ning,LIANXiao-min
(State Key Laboratory of Automotive Safety and Energy, Tsinghua University, Beijing 100084, China)
Abstract:In order to reproduce a desired sound field in an ordinary room with a reflection boundary, deconvolution of the original sound field signals is neccessary to get loudspeaker driving signals. Here, a time-domain model used for sound field reproduction in a room was established. In the model, a long-time deconvolution network was used to acquire a high accuracy and suppress the time-domain aliasing effect at the same time. Aiming at the massive computations of a convolution filtering process, the fast filtering method was applied to improve the computational efficiency. The reproduction test of the interior noise in an accelerating vehicle indicated that the proposed time-domain model is able to complete the deconvolution of sound pressure signals precisely and efficiently, and then realize the quantificational reproduction of the desired sound field in an ordinary room.
Key words:sound field reproduction; unfree field; deconvolution; filter design; fast convolution
中圖分類號:TB 535.3;TN912.27
文獻(xiàn)標(biāo)志碼：A
DOI:10.13465/j.cnki.jvs.2016.04.019
通信作者鄭四發(fā) 男，教授,博士生導(dǎo)師，1970年10月生
收稿日期：2014-09-16修改稿收到日期：2015-03-10
基金項目：國家自然科學(xué)資金資助項目(51275262);清華大學(xué)自主科研項目(2011Z02139)
第一作者 彭博 男，博士，1989年11月生
