基于稀疏協(xié)方差矩陣迭代的單快拍氣流速度估計算法

虞 飛陶建武陳 誠錢立林

①(海軍航空工程學(xué)院控制工程系 煙臺 264001)

②(空軍航空大學(xué)飛行器控制系 長春 130022)

③(94638部隊 南昌 330201)

基于稀疏協(xié)方差矩陣迭代的單快拍氣流速度估計算法

虞 飛*①②陶建武②陳 誠③錢立林②

①(海軍航空工程學(xué)院控制工程系 煙臺 264001)

②(空軍航空大學(xué)飛行器控制系 長春 130022)

③(94638部隊 南昌 330201)

該文研究基于聲傳感器陣列的單快拍氣流速度估計問題。首先，根據(jù)聲波在亞音速和超音速氣流中的傳播特性，針對特定的測量裝置，建立了聲傳感器線性陣列的輸出模型。在此基礎(chǔ)上，提出一種稀疏協(xié)方差矩陣迭代的單快拍(Sparse Covariance Matrix Iteration with a Single Snapshot, SCMISS)氣流速度估計算法，與其他稀疏估計方法相比，該文提出的SCMISS算法無需正則化參數(shù)選擇，計算量更低，具有更強(qiáng)的實時性，且只需單快拍采樣數(shù)據(jù)就可對亞音速和超音速氣流速度進(jìn)行統(tǒng)一估計。最后，為了評價所提算法的估計性能，推導(dǎo)了氣流速度估計的克拉美-羅界(Cramér-Rao Bound, CRB)表達(dá)式。仿真實驗驗證了該算法的有效性。

陣列信號處理；氣流速度估計；單快拍；稀疏參數(shù)估計；聲傳感器陣列；克拉美-羅界

1 引言

當(dāng)前，對高速、高機(jī)動性飛行器的設(shè)計研究引起了世界各主要軍事強(qiáng)國的極大關(guān)注，而傳統(tǒng)的以空速管為代表的空速測量方法無法滿足這類先進(jìn)飛行器的設(shè)計要求。現(xiàn)有的空速測量方法大都可以測量飛行器的亞音速空速，但是在超音速飛行狀態(tài)下，空氣流動特性發(fā)生質(zhì)的變化而產(chǎn)生激波，使得目前的測量方法要么難以適用，如空速管，要么需要通過參數(shù)校正才能適用，如嵌入式大氣數(shù)據(jù)傳感(Flush Air Data Sensing, FADS)系統(tǒng)[1]。因此，有必要建立統(tǒng)一的亞音速和超音速空速測量模型，使模型在各個飛行狀態(tài)都能輸出正確的空速值，同時避免因亞音速和超音速飛行狀態(tài)的氣流差異帶來的高成本的參數(shù)校正，但是，如何建立統(tǒng)一的跨音速測量模型是當(dāng)前大氣數(shù)據(jù)測量研究的一個難點。

聲矢量傳感器作為一種新型傳感器，可以同步測量聲場中同一點處的聲波質(zhì)點振速和聲壓[2]。文獻(xiàn)[3]采用聲矢量傳感器構(gòu)成了一種新型的氣流速度測量方法，它可以利用聲矢量傳感器的測量值計算氣流速度，再根據(jù)測量裝置所在的機(jī)體表面形狀及飛行姿態(tài)進(jìn)行解算，可由氣流速度成功導(dǎo)出飛行器的空速。它的缺點在于僅采用單個聲矢量傳感器來測量氣流速度，使得系統(tǒng)在抗干擾能力和測量精度上難以滿足實際要求。文獻(xiàn)[4]則建立了聲矢量傳感器陣列輸出數(shù)據(jù)與氣流速度之間的數(shù)學(xué)模型，并研究了一種數(shù)據(jù)缺失情形下的氣流速度估計算法，此算法可以有效提高氣流速度的估計精度。然而，前面提出的氣流速度測量模型和算法都只適用于亞音速氣流速度的估計，難以滿足超音速飛行器的使用要求。針對超音速氣流速度測量問題，文獻(xiàn)[5]提出了一種聲波在超音速氣流中傳播的分析方法——等效聲源法，并由此近似求得馬赫錐內(nèi)一點的質(zhì)點振速表達(dá)式，建立了基于聲矢量傳感器的信號輸出模型，但是該模型只適用于超音速氣流速度測量。為此，本文基于聲傳感器，提出了一種能同時適用于亞音速和超音速氣流速度測量的新方法。該方法綜合考慮了亞音速和超音速兩種狀態(tài)下氣流流動的差異，建立了統(tǒng)一的跨音速測量模型，從原理上，實現(xiàn)了從亞音速到高超音速整個范圍內(nèi)氣流速度的統(tǒng)一測量。

前面提到的氣流速度估計算法都是基于多快拍采樣數(shù)據(jù)的方法，而在實際應(yīng)用中，由于飛行器作復(fù)雜機(jī)動飛行動作，以及氣流環(huán)境本身復(fù)雜多變等因素的影響，使得待測氣流速度可能是實時變化的。這時，只有少數(shù)快拍甚至單快拍數(shù)據(jù)可以利用。雖然文獻(xiàn)[6]提出的基于魯棒H∞濾波的氣流速度估計算法只需單快拍采樣數(shù)據(jù)，但該算法無法用于超音速氣流速度測量。本文在建立統(tǒng)一的跨音速測量模型的基礎(chǔ)上，提出了一種基于稀疏協(xié)方差矩陣迭代的單快拍氣流速度估計算法。此算法只需單快拍采樣數(shù)據(jù)就可以對亞音速和超音速氣流進(jìn)行統(tǒng)一測量，這樣不僅節(jié)省了數(shù)據(jù)存儲空間，而且能實現(xiàn)對快變氣流速度的跟蹤測量。

2 陣列接收模型

氣流速度測量裝置如圖1所示，管路內(nèi)徑為D，一般可要求D在1～2 cm范圍內(nèi)，整個管路測量裝置可以嵌入式安裝在飛機(jī)機(jī)頭或機(jī)翼前緣。在管路內(nèi)剖面建立xOy直角坐標(biāo)系。在坐標(biāo)原點O處安裝一個聲波發(fā)射換能器，在管路內(nèi)壁正上方平行于x軸方向布置L個聲波接收換能器，分別編號1,2,…, L，設(shè)第l(l=1,2,…,L )個接收換能器相對于發(fā)射換能器的方位角為θl，則θl∈(-π/2,π/2)(從y軸測量)。

如圖1所示，設(shè)聲波波陣面由發(fā)射換能器到達(dá)第l個接收換能器的傳播距離為Rl，則Rl=D/ cosθl。假設(shè)來流從x軸正向吹來。當(dāng)氣流速度v=0時，聲波波陣面由發(fā)射換能器發(fā)射到達(dá)第l個接收換能器的傳播時間為

圖1 測量裝置剖面圖

其中，c為聲波傳播速度。當(dāng)0vc＜＜時，對應(yīng)亞音速情形，聲波波陣面?zhèn)鞑r間變?yōu)?/p>

當(dāng)vc≥時，對應(yīng)超音速情形。這時，波陣面一方面擴(kuò)展，一方面順流而下，所有波陣面被擠壓而聚集在一個被稱為馬赫錐的錐面內(nèi)。馬赫錐的半頂角與氣流速度之間有如下關(guān)系：

式中α稱為馬赫角。只有當(dāng)接收換能器位于y軸的左半平面且方位角滿足π/2+θl＜α?xí)r，換能器才能接收到聲波信號，且聲波波陣面由發(fā)射換能器傳播到該接收換能器的時間為

由π/2+θl＜α得cosθl＜c/v，則有v＜c/ cosθl。這個不等式給出了方位角為θl的接收換能器可測量到的超音速氣流速度的上界。當(dāng)待測氣流速度v≥c/cosθl時，該接收換能器將接收不到發(fā)射的聲波。

以下在建立接收換能器基陣的陣列接收模型時，假設(shè)在速度為v的氣流作用下，L個接收換能器都能接收到發(fā)射換能器發(fā)射的聲波信號。考慮發(fā)射換能器發(fā)射一個連續(xù)單頻聲波信號()st= Aej(ωt+φ(t))，則第l(l=1,2,…,L)個接收換能器的接收信號可以表示為

式中，A為發(fā)射換能器發(fā)射的聲波信號的幅值，ω=2πf , f為聲波信號的頻率，φ(t)為在區(qū)間[0,2π]之間均勻分布的隨機(jī)相位。σs(t)為零均值，單位方差復(fù)高斯隨機(jī)過程，“”表示取復(fù)數(shù)模值。nl(t)為第l個接收換能器上的零均值，方差為σl的加性復(fù)高斯白噪聲。假設(shè)信號與噪聲互不相關(guān)，以x1(t)為參照，將式(5)變換得

式中，s1(t)表示第1個接收陣元接收的聲波信號。Dl1=τl-τ1為聲波到達(dá)接收換能器1和l的延時差值，其中包含了待估計的氣流速度信息v,R1/Rl為相應(yīng)接收信號的幅度比值。將各個接收換能器的接收信號聯(lián)立起來，可得到接收換能器基陣的陣列輸出模型。

式中，x(t)=[x(t),x(t),…,x(t )]T,a(v)=[1,(R/

12L1為接收換能器基陣的L×1維陣列流形矢量，為加性噪聲矢量。

3 稀疏協(xié)方差矩陣迭代的單快拍氣流速度估計算法

3.1 協(xié)方差矩陣的稀疏表示

將待處理的氣流速度范圍V進(jìn)行N次均勻網(wǎng)格劃分{v1,v2,…,vN}，假設(shè)網(wǎng)格劃分密度足夠大，以至于真實的氣流速度v始終可以落在其中某一網(wǎng)格上，則可得如下稀疏表示模型

式中，Φ=[a(v1),a(v2),…,a(vN)]為L×N維稀疏基矩陣，a(vn)表示第n個網(wǎng)格對應(yīng)的陣列流形矢量。α(t)為N×1維稀疏信號矢量。應(yīng)注意的是，由于在亞音速氣流作用下y軸左、右兩端的接收換能器都可以接收到發(fā)射的聲波信號，而在超音速情形下，只有y軸左半平面的接收換能器才可以接收到聲波，因此，兩種情形下氣流速度范圍的網(wǎng)格劃分應(yīng)分開考慮。在亞音速時，只需對整個亞音速域[0,)c進(jìn)行均勻網(wǎng)格劃分；在超音速時，只需對超音速域[c,c/cosθl),(θl＜0)進(jìn)行均勻網(wǎng)格劃分，而無需對整個超音速域進(jìn)行劃分，這里l表示在待測超音速氣流作用下可用的接收陣元個數(shù)。

一旦求出了α(t)的最稀疏解?α(t)，則α?()的t歸一化稀疏譜的譜峰對齊的速度網(wǎng)格就是當(dāng)前氣流速度的估計值v?。目前，求解()tα的最稀疏解使用最廣泛的方法是1-范數(shù)最小化稀疏恢復(fù)方法。該方法將1-范數(shù)最小化問題轉(zhuǎn)化為二階錐規(guī)劃問題的一般形式，并在二階錐規(guī)劃框架下求解。但在實際測量氣流速度時，速度網(wǎng)格劃分的數(shù)量級將達(dá)到310，這將明顯降低二階錐規(guī)劃求解速度。同時，在很多情況下對正則化參數(shù)β的準(zhǔn)確選擇是很困難的。

下面提出的稀疏協(xié)方差矩陣迭代的單快拍氣流速度估計算法將極大地解決上述這些問題。設(shè)n(t)內(nèi)各個分量互不相關(guān)，則E[n(t)nH(t)]=diag[σ1,…,σL]。假設(shè)稀疏信號矢量α(t)與n(t)不相關(guān)，且α(t)內(nèi)各個分量也互不相關(guān)，則有陣列接收數(shù)據(jù)的稀疏協(xié)方差矩陣

式中，pn=E[αn(t)(t )]表示α(t)中第n個分量的信號功率，IL為L×L維單位陣。

于是

本文提出的稀疏協(xié)方差矩陣迭代的氣流速度估計算法就是利用單快拍陣列接收數(shù)據(jù)，通過迭代算法來估計功率矩陣P的對角元素，中最大元素的下標(biāo)即為真實氣流速度v在網(wǎng)格中的位置下標(biāo)，由此來實現(xiàn)氣流速度的估計。

3.2 單快拍氣流速度估計算法

對于單快拍陣列接收數(shù)據(jù)()tx，有如下協(xié)方差矩陣擬合準(zhǔn)則[7]

其中，tr(·)表示求矩陣的跡，且由式(12)可得

同時考慮到協(xié)方差矩陣擬合準(zhǔn)則式(13)中目標(biāo)函數(shù)f的自變量是R,?R可作為常量處理，則目標(biāo)函數(shù)f的最小化可以等效為式(16)所示的目標(biāo)函數(shù)g的最小化問題：

其證明過程類似于文獻(xiàn)[8]，本文不再贅述。為了節(jié)省計算量，這里參照文獻(xiàn)[8]轉(zhuǎn)而求解一個與優(yōu)化問題式(17)相關(guān)的優(yōu)化問題來間接求優(yōu)化問題式(17)的最優(yōu)解，則歸納出SCMISS算法的更新公式：

本文算法與文獻(xiàn)[7]中所述方法的主要區(qū)別有兩點：(1)在應(yīng)用場合方面，本文算法可以應(yīng)用到部分測量數(shù)據(jù)缺失的場合；而文獻(xiàn)[7]算法僅適用于具有完整測量數(shù)據(jù)的場合。對于不同的待估計氣流速度，有效的接收陣元個數(shù)會有所變化。因此，有效的測量數(shù)據(jù)是隨著氣流速度動態(tài)變化的。如果將文獻(xiàn)[7]中的算法直接用于本文的氣流速度估計，會因有些測量數(shù)據(jù)的缺失而導(dǎo)致估計性能嚴(yán)重下降甚至算法失效。而本文算法解決了在測量數(shù)據(jù)缺失的場合下參數(shù)估計問題。(2)在網(wǎng)格劃分方面，本文采用了動態(tài)網(wǎng)格劃分方法，而文獻(xiàn)[7]采用的是靜態(tài)網(wǎng)格劃分方法。本文算法在氣流速度范圍進(jìn)行網(wǎng)格劃分時，充分考慮了亞音速和超音速階段的氣流流動差異，采用了動態(tài)網(wǎng)格劃分方法。這樣既保證了算法的估計精度，又有效降低算法的計算復(fù)雜度。

4 算法性能分析

4.1 CRB推導(dǎo)

由給出的陣列輸出模型式(7)可知，觀測信號x(t )符合復(fù)高斯隨機(jī)分布，記為x(t )～CN(μ,Rx)，其中Rx和μ分別表示觀測矢量x(t)的L×L維協(xié)方差矩陣和L×1維均值矢量，且μ=0, Rx=(A2/定義未知的實值參數(shù)矢量，其中，則關(guān)于參數(shù)矢量η的Fisher信息矩陣的第(i,j)個元素

式中，ηi表示矢量η的第i個元素。在實際應(yīng)用中，一般只關(guān)注氣流速度v，其他參數(shù)為多余參數(shù)，且氣流速度v與其他參數(shù)不是互耦的，則關(guān)于氣流速度v的Fisher信息函數(shù)為可表示為[9]

故相對于氣流速度v的CRB為

4.2 計算量分析與比較

下面分析SCMISS算法的計算量，并與l1-范數(shù)最小化方法作比較。l1-范數(shù)最小化方法一般是在二階錐規(guī)劃的框架下采用內(nèi)點法來求最優(yōu)解，其計算復(fù)雜度為O(N3)[10]。設(shè)MDN表示復(fù)數(shù)乘和除的次數(shù)，由于SCMISS算法計算量主要集中于功率迭代更新計算中，因此通過分析迭代更新的計算量可以得到該算法的計算復(fù)雜度。在每次迭代的第1步中：

計算R(i)共需要η1=2L2(N+L)次MDN；計算R-1(i)需要2L3/3次MDN，計算需要L2次MDN，計算需要2L2次MDN，因此計算共需要次 MDN；計算N+L共需要η3=3(N+L)次MDN。因此，SCMISS算法每次迭代的第1步所需的MDN次數(shù)近似為

每次迭代的2步所需的MDN次數(shù)為

設(shè)算法達(dá)到要求的估計精度所需迭代次數(shù)為cn，則SCMISS算法中迭代更新計算所需的MDN總次數(shù)大致為

由于接收換能器基陣陣元數(shù)L數(shù)量級一般為10，而速度網(wǎng)格劃分?jǐn)?shù)N的數(shù)量級高達(dá)310，同時SCMISS算法一般經(jīng)過不超過15次迭代即可得到很高精度的全局最優(yōu)解，因此，SCMISS算法的計算復(fù)雜度明顯低于1-范數(shù)最小化方法，從而具有更好的實時性。

5 仿真實驗

考慮到管路測量裝置是擬用于機(jī)載嵌入式大氣數(shù)據(jù)傳感系統(tǒng)，設(shè)管路內(nèi)徑D=1.5 cm ，靜止聲場中的聲速c=340 m/s，發(fā)射換能器發(fā)射的聲波信號頻率取6.8 kHz。SCMISS算法和L1-Min算法速度網(wǎng)格劃分間隔均取1 m/s。在如下各組實驗中，L1-Min算法正則化參數(shù)的選擇依據(jù)是：選擇正則化參數(shù)β，使得以一個較小的概率成立，從而使殘留量與噪聲范數(shù)的期望盡量匹配，具體參考文獻(xiàn)[10]。

實驗1 亞音速情形下的算法統(tǒng)計性能實驗。在均勻穩(wěn)定氣流中發(fā)射換能器發(fā)射一單頻聲波信號入射到由10個接收換能器構(gòu)成的基陣，各陣元方位角分別為θ=[-88°,-80°,-70°,-60°,-50°,50°,60°, 70°,80°,88°]T。亞音速氣流速度v=240 m/s，對單快拍陣列觀測數(shù)據(jù)分別采用SCMISS算法、L1-Min算法和魯棒H∞濾波氣流速度估計算法[6]進(jìn)行100次Monte Carlo仿真實驗，得到氣流速度估計v?的均方根誤差(Root-Mean-Square Error, RMSE)隨信噪比的變化曲線，如圖2所示。由圖2可以看出，本文提出的SCMISS算法的估計精度明顯高于另外兩種算法，且當(dāng)信噪比高于2.5 dB時，算法對氣流速度估計誤差的方差接近CRB，說明此時估計誤差接近于最小值。

實驗2 超音速情形下的算法統(tǒng)計性能。取7個接收換能器構(gòu)成的基陣，且各陣元的方位角為θ=[-88°,-87°,-86°,-85°,-82°,-80°,-78°]T。設(shè)超音速氣流速度v=500 m/s，分別采用SCMISS算法和L1-Min算法進(jìn)行50次Monte Carlo仿真實驗，得到氣流速度估計v?的RMSE隨信噪比的變化曲線，如圖3所示。由圖3可以看出，采用這兩種算法估計的氣流速度v?的RMSE隨著信噪比的增大，與CRB越來越接近，說明算法對氣流速度的估計精度越來越高。其中，本文提出的SCMISS算法的估計精度要高于L1-Min算法。

實驗3 容錯性實驗。假設(shè)接收換能器基陣的第2個陣元失效，則x2(t)=0。亞音速氣流速度v= 240 m/s, SNR=7.5 dB，圖4(a)給出了氣流速度估計的歸一化稀疏功率譜圖。圖4(b)為無陣元失效時的稀疏功率譜圖。圖4(c)給出了在速度為v= 1435 m/s的超音速氣流作用下，氣流速度估計的稀疏功率譜圖。從圖4可以發(fā)現(xiàn)，當(dāng)接收換能器基陣的某個陣元失效時，本文算法仍可以對亞音速氣流速度進(jìn)行有效估計，但無法正確估計出超音速氣流速度。

6 結(jié)論

本文建立了穩(wěn)定氣流作用下聲傳感器陣列的輸出模型，提出了一種基于稀疏協(xié)方差矩陣迭代的單快拍氣流速度估計算法(SCMISS)，分析了SCMISS算法的計算量，并與L1-Min算法的計算量進(jìn)行了比較，推導(dǎo)了氣流速度估計的CRB表達(dá)式。理論分析和仿真實驗表明：本文提出的SCMISS算法可同時適用于亞音速和超音速氣流速度的高精度估計；當(dāng)某個接收陣元失效時，算法仍可以對亞音速氣流速度進(jìn)行有效估計；SCMISS算法比L1-Min算法具有更低的計算量。

圖2 氣流速度估計的RMSE隨信噪比變化曲線

圖3 氣流速度估計的RMSE隨信噪比變化曲線

圖4 氣流速度估計的稀疏功率譜圖

[1] Whitmore S A, Cobleigh B R, and Haering E A. Design and calibration of the X-33 flush airdata sensing (FADS) system [R]. NASA/TM-1998-206540, 1998.

[2] Wu Y I, Wong K T, and Lau S K. The acoustic vectorsensor’s near-field array-manifold[J]. IEEE Transactions on Signal Processing, 2010, 58(7): 3946-3951.

[3] 王立新, 陶建武. 一種新型的大氣數(shù)據(jù)測量方法[J]. 計量學(xué)報, 2011, 32(1): 31-35.

Wang Li-xin and Tao Jian-wu. A new measuring airdata algorithm[J]. Acta Metrologica Sinica, 2011, 32(1): 31-35.

[4] 陳誠, 陶建武, 曾賓. 數(shù)據(jù)缺失情形下的基于聲矢量傳感器陣列的空氣流動速度估計算法[J]. 電子學(xué)報, 2014, 42(3): 491-497.

Chen Cheng, Tao Jian-wu, and Zeng Bin. Estimation of airspeed in case of missing data[J]. Acta Electronica Sinica, 2014, 42(3): 491-497.

[5] Zeng Bin, Tao Jian-wu, Yu Fei, et al.. A new estimation method of airspeed based on acoustic vector sensor[C]. 2013 3th Intermational Conference on Signal Processing, Communications and Computing, Kunming, 2013: 307-310.

[6] 陳誠, 陶建武. 基于聲矢量傳感器陣列的魯棒H∞空氣速度估計算法[J]. 航空學(xué)報, 2013, 34(2): 361-370.

Chen Cheng and Tao Jian-wu. Robust H∞estimation of airspeed based on acoustic vector sensor array[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(2): 361-370.

[7] Stoica P, Babu P, and Li J. New method of sparse parameter estimation in separable models and its use for spectral analysis of irregularly sampled data[J]. IEEE Transactions on Signal Processing, 2011, 59(1): 35-47.

[8] Stoica P, Babu P, and Li J. SPICE: a sparse covariance-based estimation method for array processing[J]. IEEE Transactions on Signal Processing, 2011, 59(2): 629-638.

[9] Korso M N El, Boyer R, Renaux A, et al.. Conditional and unconditional Cramer-Rao bounds for near-field source localization[J]. IEEE Transactions on Signal Processing, 2010, 58(5): 2901-2907.

[10] Malioutov D. A sparse signal reconstruction perspective for source localization with sensor arrays[D]. [Master dissertation], Massachusetts Institute of Technology, 2003.

虞 飛： 男，1987年生，博士生，研究方向為傳感器與智能測量系統(tǒng)、陣列信號處理.

陶建武： 男，1959年生，博士，教授，博士生導(dǎo)師，從事陣列信號處理及應(yīng)用、非平穩(wěn)信號處理等方面的研究工作.

Single Snapshot Airspeed Estimation Based on Sparse Covariance Matrix Iteration

Yu Fei①②Tao Jian-wu②Chen Cheng③Qian Li-lin②①(Department of Control Engineering, Naval Aeronautical and Astronautical University, Yantai 264001, China)
②(Department of Aircraft Control Engineering, Aviation University of Air Force, Changchun 130022, China)
③(94638 Department, Nanchang 330201, China)

The issue of single snapshot airspeed estimation is researched based on acoustic sensor array. According to the propagation property of acoustic waves in subsonic and supersonic air current, the output model of acoustic sensor array is constructed for a given measuring equipment. Then an airspeed estimation algorithm based on Sparse Covariance Matrix Iteration with a Single Snapshot (SCMISS) is presented. SCMISS has several unique features not shared by other sparse estimation methods: it does not require the user to make any difficult selection of regularization parameters, and it has lower computational complexity and better real-time. What is more, the proposed algorithm can be applied to both subsonic and supersonic circumstances with single snapshot measurement. Finally, a compact expression for the Cramér-Rao Bound (CRB) on the estimation error of airspeed is derived to evaluate the performance of the proposed algorithm. Simulations are implemented to show the effectiveness of SCMISS.

Array signal processing; Airspeed estimation; Single snapshot; Sparse parameter estimation; Acoustic sensor array; Cramér-Rao Bound (CRB)

TN911.7; TN06

A

1009-5896(2015)03-0574-06

10.11999/JEIT140668

2014-05-21收到，2014-09-24改回

國家自然科學(xué)基金(61172126)和吉林省自然科學(xué)基金(20140101073 JC)資助課題

*通信作者：虞飛 yufei19871128@163.com