張宏欣,周穗華,馮士民
(海軍工程大學(xué)兵器工程系,湖北武漢 430033)
高斯粒子流濾波器
張宏欣,周穗華,馮士民
(海軍工程大學(xué)兵器工程系,湖北武漢 430033)
粒子流濾波器以粒子流速度場(chǎng)描述隨機(jī)樣本從先驗(yàn)分布到后驗(yàn)分布的演化,實(shí)現(xiàn)對(duì)系統(tǒng)狀態(tài)的貝葉斯估計(jì).針對(duì)其一般解計(jì)算復(fù)雜、難于濾波求解的問(wèn)題,導(dǎo)出一種高斯假設(shè)條件下的粒子流濾波器.在線性高斯條件下推導(dǎo)了速度場(chǎng)的解析解;證明了當(dāng)演化步長(zhǎng)趨近于0時(shí),該解析解與Kalman-Bucy濾波器的解具有一致的形式;基于該解導(dǎo)出了非線性高斯系統(tǒng)速度場(chǎng)的表達(dá)式,并進(jìn)一步利用Unscented變換近似求解.通過(guò)若干仿真算例表明,高斯粒子流濾波器放寬了系統(tǒng)噪聲為高斯型的限制,其精度優(yōu)于經(jīng)典非線性高斯濾波器,計(jì)算復(fù)雜度低于一般粒子濾波器,且具有良好的穩(wěn)定性.
非線性濾波;貝葉斯估計(jì);粒子流濾波器;速度場(chǎng);Unscented變換
貝葉斯非線性濾波基于非線性動(dòng)態(tài)系統(tǒng)模型,通過(guò)帶噪聲的觀測(cè)采樣序列近似計(jì)算出目標(biāo)狀態(tài)的后驗(yàn)分布統(tǒng)計(jì)量,來(lái)實(shí)現(xiàn)對(duì)系統(tǒng)狀態(tài)的估計(jì)[1],近幾十年來(lái)已經(jīng)在目標(biāo)跟蹤、導(dǎo)航等工程領(lǐng)域得到了廣泛應(yīng)用.根據(jù)近似方法的不同,可分為解析近似方法和仿真近似方法.解析法基于卡爾曼濾波器(Kalman Filter,KF)框架,假設(shè)后驗(yàn)分布為高斯型,近似計(jì)算非線性變換后的狀態(tài)均值和誤差矩陣,利用卡爾曼濾波框架進(jìn)行求解,其中,擴(kuò)展卡爾曼濾波(Extended Kalman Filter,EKF)基于一階偽線性化,對(duì)于高階模型精度不高;以無(wú)跡卡爾曼濾波(Unscented Kalman Filter,UKF)為代表的σ點(diǎn)卡爾曼濾波器[2,3]采用一組確定的樣本點(diǎn)更精確的捕捉均值和方差信息,對(duì)高斯分布可達(dá)到三階精度;但對(duì)于某些高階模型(如指數(shù)函數(shù),三角函數(shù)等),在高維情況下仍可能引入較大誤差,而捕捉更高階矩信息則需要更龐大的Sigma點(diǎn)集和克服數(shù)值問(wèn)題[4].
仿真近似方法采用大量隨機(jī)樣本經(jīng)過(guò)特定“變換”來(lái)近似后驗(yàn)分布.粒子濾波[5](Particle Filter,PF)從建議分布中抽取滿足后驗(yàn)分布樣本(重采樣),但并未改變樣本本身,因此建議分布樣本需足量覆蓋后驗(yàn)分布樣本,否則將導(dǎo)致濾波性能不穩(wěn)定甚至發(fā)散[6],該問(wèn)題在狀態(tài)維數(shù)較高時(shí)尤為明顯,通常需要大量隨機(jī)樣本(隨問(wèn)題維數(shù)呈指數(shù)增長(zhǎng)[7]).針對(duì)此問(wèn)題,反饋型的粒子濾波算法通過(guò)對(duì)每個(gè)先驗(yàn)分布樣本疊加反饋量來(lái)構(gòu)造建議分布樣本,普遍做法是對(duì)每個(gè)樣本進(jìn)行非線性卡爾曼濾波,從而衍生出EKF-PF及UKF-PF[8]濾波器,這類算法的時(shí)間復(fù)雜度較高,且在重要性權(quán)值計(jì)算中很容易出現(xiàn)數(shù)值不穩(wěn)定情況.
粒子流濾波器[9~12](Particle Flow Filter,PFF)是一種無(wú)需構(gòu)造建議分布和重采樣過(guò)程的反饋式粒子濾波器,通過(guò)構(gòu)造同倫(Homotopy)函數(shù),將先驗(yàn)分布粒子到后驗(yàn)分布粒子的變化過(guò)程描述為粒子的“流動(dòng)”,結(jié)合描述概率分布隨時(shí)間變化的Fokker-Planck方程(FPE),得到粒子流速度場(chǎng)(以下簡(jiǎn)稱速度場(chǎng))的偏微分方程,求出速度場(chǎng)即可利用數(shù)值方法得到后驗(yàn)分布樣本;其優(yōu)勢(shì)在于無(wú)需構(gòu)造建議分布和重采樣過(guò)程,且適于求解高維問(wèn)題[6].文獻(xiàn)[9]中通過(guò)直接求解泊松方程給出了粒子流速度場(chǎng)理論表達(dá)式的Monte-Carlo (MC)積分解,其計(jì)算復(fù)雜度過(guò)高,難于實(shí)現(xiàn).
與文獻(xiàn)[9]不同,本文從粒子流濾波器的基本理論出發(fā),推導(dǎo)出一種高斯型粒子流濾波器(Gaussian Particle Flow Filter,GPFF).假設(shè)系統(tǒng)噪聲為高斯分布的基礎(chǔ)上,根據(jù)速度場(chǎng)所滿足的偏微分方程,在線性高斯條件下推導(dǎo)了速度場(chǎng)的解析解,證明了該解與連續(xù)時(shí)間條件下的Kalman-Bucy濾波解具有一致的形式;并將該解推廣到非線性系統(tǒng),推導(dǎo)了速度場(chǎng)的期望形式解,結(jié)合無(wú)跡變換(Unscented Transform,UT)進(jìn)行近似求解.由于GPFF采用粒子進(jìn)行非線性傳播,放寬了對(duì)系統(tǒng)噪聲的高斯型限制[13],相比EKF和UKF而言具有較高的精度;保留了PFF無(wú)需構(gòu)造建議分布及重采樣的優(yōu)點(diǎn),相對(duì)PF來(lái)講降低了所需的粒子數(shù)量,減小了計(jì)算量,且穩(wěn)定性得到提高.通過(guò)對(duì)若干問(wèn)題模型進(jìn)行仿真驗(yàn)證了GPFF的有效性.
考慮如下濾波問(wèn)題模型:
若將狀態(tài)向量x視為R Rd空間標(biāo)準(zhǔn)基下的某一點(diǎn)的位置坐標(biāo),則f(x,λ)可以看做λ時(shí)刻該點(diǎn)上的速度,若f對(duì)所有粒子都有定義,則它表示粒子從先驗(yàn)分布到后驗(yàn)分布“流動(dòng)”的速度場(chǎng).因此,若能求得f,即可通過(guò)數(shù)值積分方法計(jì)算后驗(yàn)粒子.
假設(shè)變化過(guò)程不存在噪聲,根據(jù)連續(xù)時(shí)間隨機(jī)濾波理論,q(x,λ)滿足零擴(kuò)散項(xiàng)(diffusion term)的Fokker-Planck方程:
其中▽·為散度算子.對(duì)式(4)取自然對(duì)數(shù)得,
對(duì)式(7)求關(guān)于λ的偏導(dǎo)數(shù),并結(jié)合式(6)可得:
式(9)是難處理的,這是因?yàn)椋海╝)獲得合適的β值是困難的[9];(b)q(x,λ)求值需要采用分布擬合算法,計(jì)算代價(jià)過(guò)高;(c)格林公式僅適用于d>2情況;(d)已經(jīng)證明[14],這種形式的解在數(shù)值上是不穩(wěn)定的.因此需要尋求其他方法.事實(shí)上,文獻(xiàn)[9]中僅給出了使用式(9)計(jì)算速度場(chǎng)的結(jié)果,并未給出基于該式的濾波結(jié)果.
3.1解析解推導(dǎo)
對(duì)于線性系統(tǒng),假設(shè)先驗(yàn)概率密度 g(x)滿足高斯分布;似然函數(shù)h(x)為高斯型,即
分別為λ條件下的后驗(yàn)均值及協(xié)方差矩陣.在線性高斯條件下,設(shè)粒子流速度場(chǎng)f(x,λ)具有如下形式
其中A(λ)∈R Rd×d,b(λ)∈R Rd.記f=f(x,λ),注意到式(8)右端的形式,則可將其重寫(xiě)為
將式(11)、(12)和(14)代入式(8)可得:
可得
根據(jù)式(18)右端的形式,可設(shè)
將式(19)和(20)代入式(15),可以得到如下關(guān)系
為減小計(jì)算量,利用如下關(guān)系
則A(λ)可僅由先驗(yàn)協(xié)方差矩陣表示如下
將式(25)和式(26)代入式(20),則b(λ)可僅由先驗(yàn)信息表示為
式(13)、(24)和(27)給出了粒子流速度場(chǎng)在線性高斯條件下的解析解.
3.2速度場(chǎng)與Kalman-Bucy濾波器的一致性
由考慮如下的連續(xù)時(shí)間線性隨機(jī)系統(tǒng)
其中xt為連續(xù)時(shí)間狀態(tài)向量,F(xiàn)t為漂移項(xiàng);Ht為觀測(cè)矩陣,wt,vt分別為滿足零均值的高斯過(guò)程的狀態(tài)噪聲和觀測(cè)噪聲,且
對(duì)于連續(xù)時(shí)間線性高斯系統(tǒng),Kalman-Bucy濾波器[16]以線性微分方程形式給出了最小均方無(wú)偏估計(jì)解及其協(xié)方差矩陣Pt,即
利用式(26)求Pλ關(guān)于λ的導(dǎo)數(shù),可得
ε為一任意小的正數(shù).將式(37)代入式(36),化簡(jiǎn)得
利用式(23),得到Pλ關(guān)于λ的導(dǎo)數(shù)為
對(duì)式(22)取關(guān)于x的期望可得,
綜上可知線性高斯條件下λ變化過(guò)程中,速度場(chǎng)與零漂移項(xiàng)和零過(guò)程噪聲系統(tǒng)的Kalman-Bucy濾波器解是一致的.
設(shè)非線性條件下速度場(chǎng)的解具有如下形式:
將h(x)在狀態(tài)預(yù)測(cè)值^xk|k-1附近做一階泰勒展開(kāi),即
將式(41)和式(42)代入式(14),利用2.1節(jié)的推導(dǎo)過(guò)程可得
為將式(43)表示為統(tǒng)計(jì)量形式,利用式(23),并注意到
其中
式(41)和(46)給出了速度場(chǎng)在非線性高斯條件下的近似解.
采用無(wú)跡變換(UT)近似求解式(44)及式(45),UT通過(guò)一組確定的采樣點(diǎn)來(lái)近似計(jì)算隨機(jī)變量經(jīng)過(guò)非線性變換后的均值和方差,首先生成Sigma點(diǎn)以及相應(yīng)用于計(jì)算均值的權(quán)值和用于計(jì)算方差的權(quán)值即:
關(guān)于狀態(tài)的統(tǒng)計(jì)量直接通過(guò)粒子實(shí)現(xiàn),則先驗(yàn)統(tǒng)計(jì)量為
其中x(i)k|k-1為第i個(gè)先驗(yàn)粒子,N為粒子個(gè)數(shù).利用導(dǎo)出的速度場(chǎng)計(jì)算公式,對(duì)狀態(tài)在λ∈(0,1]內(nèi)進(jìn)行數(shù)值積分即可求出后驗(yàn)粒子,相應(yīng)的后驗(yàn)統(tǒng)計(jì)量可根據(jù)式(51)和(52)類似地得到.
為驗(yàn)證 GPFF的有效性,對(duì)三個(gè)典型算例,即連續(xù)Wiener過(guò)程加速度(CWPA)模型、一維非穩(wěn)定增長(zhǎng)(UNG)模型以及只測(cè)角跟蹤(BOT)模型進(jìn)行仿真,CWPA模型為線性模型,用于在高斯噪聲和非高斯噪聲條件下將GPFF與KF進(jìn)行比較;UNG模型用于在精度上將GPFF與EKF、UKF及PF等經(jīng)典非線性濾波器進(jìn)行比較.BOT模型考察多維和非理想初值條件下各濾波器性能.GPFF的粒子數(shù)量均為100個(gè),步長(zhǎng)Δλ=0.05.
5.1平面CWPA模型
平面CPWA模型為目標(biāo)加速度受隨機(jī)擾動(dòng)的連續(xù)時(shí)間運(yùn)動(dòng)模型,以隨機(jī)微分方程描述如下:
其中,
k=1,2,…,T,T為觀測(cè)點(diǎn)數(shù),
噪聲協(xié)方差矩陣
Δt為離散化時(shí)間間隔.
以平面位置坐標(biāo)作為觀測(cè)量,觀測(cè)方程為
觀測(cè)矩陣H=(I2×202×4).設(shè)Δt=0.1,T=80.初始狀態(tài)為x0=06×1,濾波器初始分布為N(x0,P0|0),初始誤差矩陣P0|0=diag(0.1,0.1,0.1,0.1,0.5,0.5).分別在高斯噪聲和混合高斯噪聲條件下進(jìn)行仿真,即
為直觀比較,采用由下式定義的維度無(wú)關(guān)均方誤差(MSE)
圖1(a)~(b)分別給出了在不同類型觀測(cè)噪聲下,KF與GPFF的100次MC計(jì)算MSE對(duì)比.由圖可知,在高斯噪聲條件下,GPFF與KF誤差近似相同,這與2.2中證明的結(jié)論是一致的;而在非高斯噪聲條件下,KF誤差增幅較大,且變得不穩(wěn)定,而GPFF仍具有較小的誤差,這表明GPFF保持了粒子型算法對(duì)系統(tǒng)噪聲的適應(yīng)性.
圖2給出了非高斯噪聲條件下GPFF與KF的跟蹤結(jié)果對(duì)比,可知在非高斯觀測(cè)噪聲下,KF跟蹤軌跡與真實(shí)軌跡產(chǎn)生了較大偏差,而GPFF仍可準(zhǔn)確地跟蹤目標(biāo).
5.2UNG模型
該模型是非線性濾波中的基準(zhǔn)問(wèn)題[5],其狀態(tài)空間模型為
圖3給出了各個(gè)算法 MSE的100次 MC仿真計(jì)算結(jié)果,其中PF粒子數(shù)N=100.由于UNG是一個(gè)不含采樣時(shí)間的完全離散模型,且包含高階非線性項(xiàng),導(dǎo)致EKF的一階線性化誤差較大,穩(wěn)定性較差;GPFF在MSE和穩(wěn)定性上均顯著地優(yōu)于EKF以及同樣采用UT變換的UKF.
表1進(jìn)一步給出了各個(gè)算法100次MC計(jì)算的MSE均值和方差以及單次平均計(jì)算時(shí)間.其中,對(duì)于PF 和GPFF分別給出了 N=100,N=500和N=1000時(shí)的計(jì)算結(jié)果(如括號(hào)內(nèi)所示),執(zhí)行時(shí)間在MATLAB R2013(8.1.0.604)版本和2.5GHz雙核計(jì)算機(jī)條件下測(cè)得.由表可知,作為高斯型濾波器,GPFF精度優(yōu)于EKF和UKF;而作為粒子型算法,GPFF的計(jì)算量隨粒子數(shù)量的增長(zhǎng)變化是較慢的;PF采用重采樣方法,其計(jì)算時(shí)間隨粒子數(shù)量增長(zhǎng)較快,因此GPFF的計(jì)算效率高于PF.綜合圖表結(jié)果看,由于PF對(duì)樣本分布不基于任何假設(shè),對(duì)于本例的情況(一維模型和理想初始條件),其精度最高.而在算例3中將會(huì)看到,對(duì)于多維問(wèn)題,PF的性能將會(huì)大幅下降.
表1100 次MC計(jì)算結(jié)果:MSE均值,方差,計(jì)算時(shí)間
5.3BOT模型
BOT模型[15]是只以方位角測(cè)量值作為觀測(cè)量的跟蹤模型,其狀態(tài)模型為離散Wiener速度模型,
觀測(cè)噪聲vk~N(02×1,R),噪聲協(xié)方差矩陣R= 0.052I2×2.觀測(cè)點(diǎn)坐標(biāo),設(shè)Δt=0.1s,時(shí)刻k=1,2,…,T,T= 200.狀態(tài)初值x0=(0,36,0.2,-1)T,濾波器初始分布為N(^x0,P0|0),初始估計(jì)值^x0=(-6,17,0,0)T,初始誤差矩陣P0|0=(25,25,1,1)T.PF粒子個(gè)數(shù)N=500.
采用各算法進(jìn)行100次Monte-Carlo計(jì)算,對(duì)于k時(shí)刻的某一個(gè)狀態(tài)分量εk,按下式計(jì)算各個(gè)狀態(tài)分量的MSE
其中Nm為Monte-Carlo仿真次數(shù),^εk為k時(shí)刻狀態(tài)估計(jì)值.由上式可知,MSEkε是100次Monte-Carlo計(jì)算的平方誤差在特定時(shí)刻上的均值,可衡量算法隨觀測(cè)值信息增加的收斂性能.其計(jì)算結(jié)果如圖4(a)~(d)所示,由圖可知,EKF、UKF和GPFF均能夠收斂到較小的誤差范圍,而PF則出現(xiàn)了顯著偏差,且隨時(shí)間變化的趨勢(shì)是起伏的,可見(jiàn)在多維模型和不理想初始值條件下,經(jīng)典PF的性能較差;進(jìn)一步地可以看出,相較EKF和 UKF,GPFF能夠更快的收斂,且誤差較小.
圖5(a)~(d)給出了各個(gè)算法MSE的100次MC計(jì)算結(jié)果,可以看到,PF對(duì)于目標(biāo)縱軸坐標(biāo)的估計(jì)產(chǎn)生了很大誤差,而對(duì)于目標(biāo)橫軸速度估計(jì)的MSE則較小,且PF穩(wěn)定性較差,這種現(xiàn)象是由于y的估計(jì)初值與真實(shí)初值偏差最大,而.x的估計(jì)初值與真實(shí)初始偏差很小造成的,這從圖4中也可得出.
圖4和圖5表明一定粒子數(shù)量條件下,對(duì)于多維估計(jì)問(wèn)題,PF的性能對(duì)于初始條件是十分敏感的,若初始條件不夠理想的,其平均性能甚至?xí)佑贓KF算法.而GPFF對(duì)于每個(gè)分量估計(jì)的MSE均明顯地低于其他算法,具有較高的精度,且在所有仿真中都保持了穩(wěn)定的性能.
表2中給出了各個(gè)算法條件下,100次MC計(jì)算所得的狀態(tài)分量估計(jì)MSE均值,方差和執(zhí)行時(shí)間,其中對(duì)于PF分別給出了粒子數(shù)時(shí)的結(jié)果.可知各算法中GPFF 的MSE均值為最小,且穩(wěn)定性好;UKF精度高于 EKF;PF算法在粒子數(shù)量為4000時(shí)給出了較好的估計(jì)結(jié)果,但此時(shí)其精度仍低于GPFF,表明在不理想初值和/或多維模型條件下,PF若無(wú)足夠粒子來(lái)代表多維樣本空間,或建議分布本身就與后驗(yàn)分布相去甚遠(yuǎn)(不理想的初值),則重采樣方法此時(shí)是低效的.
表2 MSE均值、標(biāo)準(zhǔn)差和執(zhí)行時(shí)間
針對(duì)粒子流濾波器中速度場(chǎng)難于計(jì)算的問(wèn)題,基于粒子流濾波器基本理論,導(dǎo)出一種高斯粒子流濾波器.在線性高斯條件下推導(dǎo)了速度場(chǎng)的解析解;證明了當(dāng)粒子流演化步長(zhǎng)趨近于0時(shí)速度場(chǎng)解析解與Kalman-Bucy濾波器解的一致性;并將算法推廣到非線性系統(tǒng),利用Unscented變換近似求解.仿真算例表明,算法放寬了系統(tǒng)噪聲為高斯型的限制,其精度優(yōu)于EKF 和UKF等經(jīng)典非線性高斯濾波器,計(jì)算復(fù)雜度低于一般粒子濾波器,且性能穩(wěn)定,具有一定的應(yīng)用價(jià)值.
[1]James C.Bayesian Signals Processing[M].Hoboken:John Wiley&Sons Press,2009.3-10.
[2]Simon J,Jeffery U,Hugh D.A new method for the nonlinear transformation of means and covariance in filters and estimators[J].IEEE Transactions on Auto Control,2000,45(3):477-481.
[3]Arasaratnam I,Haykin S.Cubature Kalman filters[J]. IEEE Transactions on Auto Control,2009,56(6):1254 -1269.
[4]Rudolph M.Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Models[D].Oregon:Oregon Health&Science University,2004.251-256.
[5]Gordon N,Salmond D,Smith A.Novel approach to nonlinear/non-Gaussian Bayesian state estimation[J].IEE Radar and Signal Processing,1993,140(2):107-113.
[6]Fred D.Nonlinear filters:Beyond Kalman filter[J].IEEE Aerospace&Electronics Systems Magazine,2005,43(8):57-69.
[7]Synder C,Bengtsson T,Bickel P,et al.Obstacles to highdimensional particle filtering[J].Monthly Weather Review,2008,136(1):4629-4640.
[8]Rudolph M.The Unscented Particle Filter[R].Cambridge:Cambridge University,2000.
[9]Daum F.Coulomb's law particle flow for nonlinear filter[A].Proceedings of SPIE on Signal Processing and Sensor Fusion[C].USA:SPIE Press,2011.3351-3362.
[10]Daum F.Numerical experiments with Coulomb’s law for particle flow in nonlinear filters[A].Proceedings of SPIE Conference[C].San Diego:SPIE Press,2011.2137 -2146.
[11]Fred D,Jim H.Particle flow with non-zero diffusion for nonlinear filters[A].Proceedings of SPIE on Signal Processing and Sensor Fusion[C].USA:SPIE Press,2013. 1745-1755.
[12]Tao D,Mark C.Implementation of the daum-huang exactflow particle filter[A].IEEE Statistical Signal Processing Workshop[C].Ann Arbor:IEEE Press,2012.257-261.
[13]Kotecha J,Djuric P.Gaussian particle filtering[J].IEEE Transactions on Signal Processing,2003,51(10):2592 -2601.
[14]Tao Y,Mehta P G,Meyn S P.Feedback particle filter [J].IEEE Transactions on Automatic Control,2013,58 (10):2037-2043.
[15]Bar-Shalom Y,Xiaorong Li.Estimation with Applications to Tracking and Navigation[M].New York:John Wiley &Sons Press,2005.50-54.
[16]Simon D.Optimal State Estimation[M].New York:John Wiley&Sons Press,2006.22-27.
張宏欣 男,1987年12月出生,陜西漢中人.2010年畢業(yè)于西安理工大學(xué),現(xiàn)為海軍工程大學(xué)博士生,從事統(tǒng)計(jì)信號(hào)處理及目標(biāo)跟蹤相關(guān)研究.
E-mail:mylifeforthebattle@hotmail.com
周穗華 男,1962年10月出生,廣東五華人,1984年畢業(yè)于海軍工程學(xué)院,1990年在海軍工程學(xué)院獲得博士學(xué)位.現(xiàn)為海軍工程大學(xué)教授,從事軍用目標(biāo)特性信息處理及武器系統(tǒng)總體設(shè)計(jì)方面研究.
Gaussian Particle Flow Filter
ZHANG Hong-xin,ZHOU Sui-hua,F(xiàn)ENG Shi-min
(Department of Weapon Engineering,Naval University of Engineering,Wuhan,Hubei 430033,China)
Particle flow filter formulated the dynamics from prior samples with posterior samples with particle flow velocity field to perform Bayesian estimation of system state.To address difficulties of particle velocity field computation in present particle flow filter,a novel particle flow filter based on Gaussian assumption was proposed.The analytical solution of velocity field under linear Gaussian condition was derived.The consistency of this analytical solution with Kalman-Bucy filter for continuous system,when discrete dynamic step goes to zero,was proved.The solution was finally extended to obtain the nonlinear Gaussian velocity field expression which can be approximated by using unscented transformation.Several simulations revealed the effectiveness over classic nonlinear Gaussian assumption on accuracy and particle filter on efficiency and stability.
nonlinear filtering;Bayesian estimation;particle flow filter;velocity field;Unscented transformation
TP202
A
0372-2112(2016)04-0795-09
電子學(xué)報(bào)URL:http://www.ejournal.org.cn 10.3969/j.issn.0372-2112.2016.04.007
2014-10-28;
2015-03-23;責(zé)任編輯:孫瑤