運(yùn)動(dòng)目標(biāo)激光微多普勒效應(yīng)平動(dòng)補(bǔ)償和微動(dòng)參數(shù)估計(jì)?

郭力仁 胡以華 董驍 李敏樂

(國防科技大學(xué)電子對(duì)抗學(xué)院,脈沖功率激光技術(shù)國家重點(diǎn)實(shí)驗(yàn)室,合肥 230037)(2017年12月28日收到;2018年4月26日收到修改稿)

1 引 言

運(yùn)動(dòng)目標(biāo)都會(huì)存在由于發(fā)動(dòng)機(jī)振動(dòng)、旋翼轉(zhuǎn)動(dòng)等引起的相對(duì)于整體平動(dòng)而言的微運(yùn)動(dòng).這些微動(dòng)會(huì)對(duì)探測信號(hào)造成多普勒效應(yīng)外的附加頻率調(diào)制,產(chǎn)生微多普勒效應(yīng)[1].不同目標(biāo)的微動(dòng)特征具有惟一性,開辟了探測和識(shí)別的新途徑.現(xiàn)有的微多普勒研究主要用于目標(biāo)分類,但通過對(duì)目標(biāo)運(yùn)動(dòng)和微動(dòng)參數(shù)的精確估計(jì),再結(jié)合充足的目標(biāo)先驗(yàn)信息,還可開創(chuàng)性地實(shí)現(xiàn)對(duì)同類目標(biāo)的精細(xì)識(shí)別.實(shí)際的遙感目標(biāo)振動(dòng)幅度在微米量級(jí),如汽車、飛機(jī)發(fā)動(dòng)機(jī)運(yùn)轉(zhuǎn)引起的表面振動(dòng)等.微波探測對(duì)此類微多普勒效應(yīng)并不敏感,但采用波長更短的激光探測手段,可以得到更為顯著的微多普勒頻移[2],具有更高的靈敏度和分辨率,有利于提高參數(shù)估計(jì)及目標(biāo)識(shí)別的準(zhǔn)確性.

目標(biāo)運(yùn)動(dòng)產(chǎn)生的多普勒效應(yīng)一般可用多項(xiàng)式相位信號(hào)(PPS)進(jìn)行建模,而目標(biāo)微動(dòng)廣泛采用正弦運(yùn)動(dòng)模擬,其回波是典型的正弦調(diào)頻(SFM)信號(hào)[3,4].兩者在回波中形成PPS-SFM混合信號(hào),這在目標(biāo)探測的工程應(yīng)用中普遍存在[5].但是,目前多數(shù)對(duì)目標(biāo)微多普勒特征的研究只考慮回波為SFM模型[6],忽略了目標(biāo)整體運(yùn)動(dòng)的影響.對(duì)于混合信號(hào),傳統(tǒng)的適用于純SFM信號(hào)或PPS的估計(jì)算法都不再適用.而現(xiàn)有的針對(duì)混合信號(hào)參數(shù)估計(jì)的研究又存在估計(jì)誤差大、計(jì)算過程復(fù)雜的不足.一般而言,直接估計(jì)平動(dòng)微動(dòng)混合信號(hào)參數(shù)難以實(shí)現(xiàn),實(shí)際操作中往往先將兩項(xiàng)分離再分別進(jìn)行估計(jì).分離方法可大致分為兩類[7],即基于時(shí)頻分布類和基于信號(hào)模型變換類.前者對(duì)信號(hào)進(jìn)行時(shí)頻分析,提取信號(hào)的瞬時(shí)頻率特征[8,9],利用多項(xiàng)式擬合補(bǔ)償平動(dòng)分量,再對(duì)剩余微動(dòng)分量進(jìn)行參數(shù)估計(jì)[10,11].該類方法有利于直觀地觀測信號(hào)調(diào)頻特征及各步處理效果,但對(duì)信噪比要求較高,且受限于時(shí)頻分布分辨率、峰值提取準(zhǔn)確度、多項(xiàng)式擬合精度等多重影響,存在多層誤差傳遞,估計(jì)精度和時(shí)效性較差,甚至?xí)霈F(xiàn)運(yùn)動(dòng)階數(shù)的錯(cuò)誤估計(jì).后者則根據(jù)兩類信號(hào)的不同特性,利用高階模糊函數(shù)[6]和延時(shí)共軛相乘法[12,13]逐階消去信號(hào)中多項(xiàng)式相位,只保留SFM再進(jìn)行參數(shù)估計(jì),但都需要目標(biāo)運(yùn)動(dòng)階數(shù)這一先驗(yàn)信息,且運(yùn)算過程復(fù)雜.此外還有通過粒子濾波(PF)[14]、廣義周期[15]法從信號(hào)層面對(duì)混合信號(hào)進(jìn)行估計(jì),但只考慮了固定載頻的情況,相當(dāng)于勻速運(yùn)動(dòng)產(chǎn)生的一階PPS,模型比較簡單,不適合目標(biāo)實(shí)際運(yùn)動(dòng)情況.

在激光探測的微多普勒信號(hào)參數(shù)估計(jì)中,由于探測波長短、采樣率高,對(duì)SFM部分參數(shù)估計(jì)時(shí)將產(chǎn)生較微波探測非線性程度更高的代價(jià)函數(shù)[16],密集的局部峰值極易使傳統(tǒng)的迭代類方法、統(tǒng)計(jì)類方法不收斂或錯(cuò)誤收斂[17],同時(shí)計(jì)算量也會(huì)隨參數(shù)數(shù)量的增加而呈指數(shù)上升.現(xiàn)有的SFM變換[18]和循環(huán)自相關(guān)[19]等參數(shù)估計(jì)方法都只針對(duì)兩參數(shù)的簡化模型,不適合處理實(shí)際的多參數(shù)信號(hào).對(duì)于以上的混合信號(hào)分離和微動(dòng)參數(shù)估計(jì)問題,本文提出基于分?jǐn)?shù)階傅里葉變換(FrFT)的平動(dòng)補(bǔ)償方法,通過設(shè)計(jì)FrFT參數(shù)域的帶寬搜索方法,可以從混合信號(hào)中精確估計(jì)平動(dòng)參數(shù),實(shí)現(xiàn)平動(dòng)和微動(dòng)的分離.對(duì)補(bǔ)償后的信號(hào),通過改進(jìn)的靜態(tài)參數(shù)PF方法估計(jì)其中的微動(dòng)參數(shù),所提方法不受噪聲類型限制,可同時(shí)估計(jì)多個(gè)參數(shù),避免誤差傳遞的影響,而且估計(jì)參數(shù)數(shù)量的增加并不會(huì)過多地增加算法的計(jì)算量.利用累積殘差定義新的權(quán)重函數(shù),并用馬爾可夫-蒙特卡羅(MCMC)保持粒子多樣性,保證了靜態(tài)參數(shù)模型下PF的快速有效收斂.通過對(duì)算法估計(jì)性能進(jìn)行仿真和實(shí)驗(yàn)分析,本文方法可實(shí)現(xiàn)對(duì)平動(dòng)和微動(dòng)混合信號(hào)參數(shù)的準(zhǔn)確估計(jì).

2 基于FrFT的平動(dòng)參數(shù)估計(jì)

2.1 運(yùn)動(dòng)目標(biāo)微多普勒信號(hào)模型

利用點(diǎn)散射模型對(duì)微多普勒效應(yīng)進(jìn)行建模[1],經(jīng)過后續(xù)化簡和整合,各類微動(dòng)的點(diǎn)散射模型都具有SFM信號(hào)的基本形式.若考慮目標(biāo)微動(dòng)特征平穩(wěn),信號(hào)處理中一般只需要一到兩個(gè)周期的信號(hào)長度就可以實(shí)現(xiàn)對(duì)微動(dòng)參數(shù)的估計(jì).實(shí)際目標(biāo)的微動(dòng)

頻率多在百赫茲量級(jí)甚至更高,所以要求處理的信號(hào)時(shí)長很短.對(duì)于平動(dòng),人們常采用PPS進(jìn)行模擬,在較短的時(shí)間內(nèi),平動(dòng)相對(duì)于微動(dòng)是緩變量,用勻加速運(yùn)動(dòng)的模型可完全滿足對(duì)平動(dòng)的模擬,這時(shí)信號(hào)平動(dòng)引起的多普勒效應(yīng)表現(xiàn)為線性調(diào)頻(LFM)信號(hào)的形式.所以運(yùn)動(dòng)目標(biāo)微多普勒效應(yīng)回波基帶信號(hào)可表示為

式中,ν表示目標(biāo)運(yùn)動(dòng)初速,a表示加速度,λ表示激光波長,Dv表示目標(biāo)微動(dòng)幅度,fv表示目標(biāo)微動(dòng)頻率,ρ0為微動(dòng)初始相位,f0和μ0分別對(duì)應(yīng)目標(biāo)的速度和加速度項(xiàng),fmD表示微多普勒頻率項(xiàng),為簡便起見,模型中考慮散射點(diǎn)振動(dòng)方位角和俯仰角及目標(biāo)相對(duì)于雷達(dá)的方位角和俯仰角為.θ0為信號(hào)初始相位,noi(t)為信號(hào)噪聲.

目前對(duì)微多普勒效應(yīng)的研究大多只針對(duì)微動(dòng)引起的SFM模型展開,而不考慮平動(dòng)對(duì)微動(dòng)特征的影響.事實(shí)上,未補(bǔ)償平動(dòng)或只用測速數(shù)據(jù)粗略補(bǔ)償?shù)男盘?hào)都是一個(gè)LFM-SFM的混合形式,平動(dòng)分量的存在會(huì)破壞微多普勒特征的周期性,導(dǎo)致現(xiàn)有微動(dòng)參數(shù)估計(jì)方法的失效.而且平動(dòng)還會(huì)使頻譜展寬,對(duì)采樣率和探測器帶寬都提出了更高的要求,對(duì)于同樣時(shí)長的信號(hào),需要處理的數(shù)據(jù)點(diǎn)數(shù)會(huì)極大地增加,導(dǎo)致運(yùn)算量劇增,所以有必要先對(duì)平動(dòng)參數(shù)進(jìn)行精確估計(jì)和補(bǔ)償.

2.2 平動(dòng)參數(shù)估計(jì)原理

FrFT可看作是在一組正交的LFM信號(hào)基上展開信號(hào)[20].所以,一個(gè)確定的LFM信號(hào)可以在FrFT的某一階次上得到最好的能量聚集效果[21].定義信號(hào)s(t)的p階FrFT為

式中階次p可以為任意實(shí)數(shù),角度α=pπ/2;Kα(t,u)為FrFT的變換核,可表示為

圖1 LFM-SFM信號(hào)FrFT示意圖 (a)能量聚集原理;(b)(α,u)域分布Fig.1.FrFT on LFM-SFM signal:(a)Energy accumulation principle;(b)FrFT distribution in parameter f i eld.

本文研究的運(yùn)動(dòng)目標(biāo)微多普勒信號(hào)表現(xiàn)為LFM-SFM混合的信號(hào)形式,經(jīng)過FrFT的結(jié)果如圖1所示.

從圖1可以看出,對(duì)于混合信號(hào)進(jìn)行FrFT在對(duì)應(yīng)的α下不是向單純的LFM信號(hào)一樣聚集到一點(diǎn),而是在u域有最窄的帶寬,從圖1(b)可以看到繼續(xù)采用傳統(tǒng)的最大值搜索法來確定LFM項(xiàng)參數(shù)顯然已經(jīng)不再合適.這里提出雙向閾值搜索法來精確計(jì)算各u域的帶寬.用M×N維矩陣F表示FrFT參數(shù)域分布,其中M為對(duì)α離散化掃描的次數(shù),與步進(jìn)長度成反比;N為信號(hào)長度.

式中,Fαm表示FrFT在第m行的分布向量,對(duì)每行的Fαm取均值作為該行的閾值th(αm).帶寬搜索原理如圖2.

圖2 雙向閾值搜索法u域確定帶寬Fig.2. Bandwidth determination by bidirectional threshold searching method.

圖2中首先對(duì)Fαm行進(jìn)行正向搜索,將第一次出現(xiàn)大于閾值的u作為帶寬的左邊界,有然后再對(duì)進(jìn)行反向搜索,確定帶寬的右邊界其中表示對(duì)括號(hào)中的向量進(jìn)行左右反轉(zhuǎn)操作. 此時(shí)第m行的u域帶寬等于所以通過搜索矩陣F各行最窄的帶寬,再代入(5)式就可以確定平動(dòng)參數(shù),注意此時(shí)這里應(yīng)注意閾值的選擇應(yīng)與信號(hào)信噪比關(guān)聯(lián),當(dāng)信噪比較低(<10 dB)時(shí),為避免噪聲對(duì)帶寬邊界搜索的干擾,閾值可取為均值的兩倍.

得到平動(dòng)參數(shù)后重構(gòu)相應(yīng)的多普勒分量sD(t),并將其從混合信號(hào)中去除,實(shí)現(xiàn)平動(dòng)補(bǔ)償,此時(shí)剩余信號(hào)中只包含微多普勒效應(yīng).

3 基于靜態(tài)參數(shù)PF的微動(dòng)參數(shù)估計(jì)

3.1 靜態(tài)參數(shù)PF模型

基于PF的參數(shù)估計(jì)方法不受信號(hào)模型的限制,可實(shí)現(xiàn)多參數(shù)同時(shí)估計(jì),避免誤差傳遞,保證每個(gè)參數(shù)的估計(jì)精度.SFM參數(shù)估計(jì)在PF中屬于靜態(tài)參數(shù)估計(jì)問題,以(7)式中的補(bǔ)償結(jié)果smD(t)作為系統(tǒng)的觀測方程,微動(dòng)參數(shù)作為狀態(tài)方程,有

式中Yt為t時(shí)刻觀測值,Xt為t時(shí)刻粒子狀態(tài),t時(shí)刻第i個(gè)粒子可寫為對(duì)于靜止參數(shù)模型,粒子多樣性不會(huì)增加,只能在初始化形成的粒子組中尋找最優(yōu)解,無法保證正確收斂,所以需要加入抖動(dòng)vt,來維持粒子的多樣性.此時(shí),vt相當(dāng)于粒子狀態(tài)的更新量,當(dāng)vt過大時(shí),粒子更新幅度也較大,不利于參數(shù)的精細(xì)搜索,容易出現(xiàn)無效的迭代計(jì)算;當(dāng)vt太小時(shí),粒子每次的更新幅度也較小,需要大量迭代才能收斂,效率較低.所以vt決定了整個(gè)算法的效率.本文采用MCMC[22]算法更新粒子狀態(tài),算法的建議分布方差就等效于模型中添加的抖動(dòng)vt,效果相當(dāng)于在隨機(jī)更新的基礎(chǔ)上加入了方向的選擇,使粒子始終向真值方向更新,有效減少了收斂所需的迭代次數(shù).

3.2 多維參數(shù)權(quán)值計(jì)算

根據(jù)序貫重要性采樣方法[23],重要性密度函數(shù)取此時(shí)各粒子權(quán)重可化簡為

3.3 粒子狀態(tài)更新

對(duì)于多參數(shù)情況,同時(shí)產(chǎn)生所有參數(shù)的候選狀態(tài)時(shí),由于各參數(shù)粒子間隨機(jī)匹配,難以得到高權(quán)值的粒子組合,這會(huì)導(dǎo)致極低的接收概率,粒子可能長時(shí)間得不到更新,不僅對(duì)改善粒子多樣性沒有任何幫助,還會(huì)降低運(yùn)算效率.對(duì)此,本文提出自適應(yīng)方差Gibbs算法來具體實(shí)現(xiàn),雖然增加了單次抽樣的計(jì)算量,但接收概率更高,整體上是加快了收斂,具有更高的效率.

自適應(yīng)方差是指根據(jù)實(shí)時(shí)的估計(jì)誤差來調(diào)整算法建議分布的方差,在初始時(shí),粒子往往與真值差距較大,這時(shí)選用較大的方差來提高候選狀態(tài)在整個(gè)支撐域上的搜索范圍,以便于快速更新到真值附近,相當(dāng)于粗搜索;當(dāng)幾次迭代搜索之后,殘差變小,則減小方差,使其在真值附近精細(xì)搜索,提高樣本接受概率和參數(shù)估計(jì)精度.自適應(yīng)方差法可表示為

3.4 靜態(tài)參數(shù)PF參數(shù)估計(jì)流程

PF前首先要解決多維參數(shù)初始化的問題.對(duì)于未知參數(shù)值,一般在取值范圍內(nèi)隨機(jī)產(chǎn)生初始值,這需要大量的粒子數(shù)目來提高產(chǎn)生優(yōu)質(zhì)初始粒子的概率,但過多的粒子會(huì)導(dǎo)致后續(xù)繁重的計(jì)算量.對(duì)此提出變粒子數(shù)目的PF方法:首先設(shè)置較大的初始粒子數(shù)目,計(jì)算各粒子組的權(quán)重,只保留權(quán)重較高的少數(shù)粒子作為之后迭代的初值.該方法可提高初始值的質(zhì)量,并避免在低權(quán)重的粒子上進(jìn)行無謂的計(jì)算,優(yōu)化了算法效率.應(yīng)當(dāng)說明的是,算法只在初始化后進(jìn)行一次減少粒子數(shù)的操作.

PF的終止條件可通過設(shè)置累積誤差閾值Rstop來進(jìn)行判斷.本文利用靜態(tài)參數(shù)PF進(jìn)行目標(biāo)微動(dòng)參數(shù)估計(jì)的流程圖如圖3.

具體的算法步驟可以總結(jié)為:

1)確定參數(shù)范圍,設(shè)定粒子數(shù)為Np,初始化粒子迭代次數(shù)t=0;

5)根據(jù)(12)式確定建議分布方差,采用MCMC算法更新粒子狀態(tài),保持多樣性;

圖3 靜態(tài)參數(shù)PF算法流程圖Fig.3.Flow chart of static parameter PF algorithm.

4 仿真分析與實(shí)驗(yàn)驗(yàn)證

4.1 仿真信號(hào)參數(shù)估計(jì)

4.1.1 平動(dòng)補(bǔ)償

設(shè)目標(biāo)為4缸發(fā)動(dòng)機(jī)的汽車,發(fā)動(dòng)機(jī)轉(zhuǎn)速為3000 min,每缸點(diǎn)火一次產(chǎn)生一次振動(dòng),根據(jù)發(fā)動(dòng)機(jī)工作原理,對(duì)應(yīng)的振動(dòng)頻率為fv=3000×2/60=100 Hz,目標(biāo)表面振動(dòng)幅度Dv為5μm,振動(dòng)對(duì)信號(hào)調(diào)制的初始相位ρ0為π/6 rad,激光波長為1.55μm,信噪比為20 dB.設(shè)目標(biāo)平動(dòng)引起的多普勒參數(shù)為f0=104Hz(認(rèn)為是初步的平動(dòng)補(bǔ)償后的殘余分量),μ0=4×106Hz,對(duì)應(yīng)的加速度為3 m/s2,仿真信號(hào)的時(shí)頻分布如圖4所示.

圖4中LFM-SFM信號(hào)的瞬時(shí)頻率特征為兩者的融合,兩者相互影響破壞了兩種信號(hào)原有的時(shí)頻特征,傳統(tǒng)的針對(duì)單純的SFM或LFM信號(hào)特征的提取方法已不再使用.根據(jù)2.2節(jié)所提方法對(duì)平動(dòng)參數(shù)進(jìn)行估計(jì),得到FrFT參數(shù)域分布及帶寬搜索結(jié)果如圖5所示.

圖5(a)中顏色較亮的區(qū)域?qū)?yīng)能量聚集區(qū)域,兩側(cè)亮暗交界線即對(duì)應(yīng)各旋轉(zhuǎn)角下的u域帶寬.當(dāng)α與信號(hào)中LFM對(duì)應(yīng)的瞬時(shí)頻率曲線垂直時(shí),信號(hào)在u域有最好的能量聚集,此時(shí)帶寬最窄.圖5(b)為雙向閾值搜索法估計(jì)出的帶寬隨α的變化趨勢,與圖5(a)中的亮暗交界線包絡(luò)的變化趨勢一致,證明了帶寬估計(jì)方法的正確性.確定帶寬最小時(shí)的α值,代入(5)式計(jì)算LFM參數(shù),并和目前典型的基于時(shí)頻特征多項(xiàng)式擬合的平動(dòng)參數(shù)估計(jì)方法進(jìn)行對(duì)比.

圖4 LFM-SFM混合信號(hào)時(shí)頻分布Fig.4.Time-Frequency distribution of LFM-SFM signal.

圖5 FrFT對(duì)混合信號(hào)處理結(jié)果 (a)FrFT在(α,u)域投影;(b)u域帶寬隨α的變化Fig.5.FrFT results of LFM-SFM signal:(a)Projection on(α,u)domain;(b)bandwidth changes with α.

表1 平動(dòng)參數(shù)估計(jì)結(jié)果對(duì)比Table 1.Comparison of translational parameter estimation results.

表1中處理的信號(hào)時(shí)長為兩個(gè)微動(dòng)周期,可以看出擬合法得到的參數(shù)估計(jì)精度受微動(dòng)時(shí)頻特征初始相位影響極大,這可通過圖1(a)進(jìn)行解釋,不同的初相導(dǎo)致瞬時(shí)頻率曲線的起始位置和中心頻率之間距離不同,而這個(gè)距離會(huì)直接影響擬合的結(jié)果.當(dāng)初相為π/2時(shí),起始點(diǎn)恰和中心頻率重合,不存在擬合誤差,所以此時(shí)平動(dòng)參數(shù)估計(jì)精度最高,而在其他相位下,估計(jì)結(jié)果會(huì)存在不同程度的誤差,且遠(yuǎn)高于本文參數(shù)化方法的估計(jì)結(jié)果.相比之下,本文采用了基于投影原理的參數(shù)估計(jì)方法,不受信號(hào)初相的影響,具有穩(wěn)定的估計(jì)精度.

對(duì)目標(biāo)平動(dòng)的補(bǔ)償效果會(huì)直接影響到后續(xù)微動(dòng)參數(shù)的估計(jì)精度.這里用平動(dòng)補(bǔ)償后信號(hào)與實(shí)際SFM信號(hào)的相位波形相似度作為標(biāo)準(zhǔn)對(duì)LFM參數(shù)的估計(jì)精度提出要求.定義波形相似度為圖6為真實(shí)SFM相位和經(jīng)過平動(dòng)補(bǔ)償后的信號(hào)相位在γ=0.9999時(shí)的波形對(duì)比.圖7為波形相似度隨LFM參數(shù)估計(jì)相對(duì)誤差的變化情況.

從圖6可以看出,γ=0.9999時(shí)兩波形基本重合,此時(shí)可以認(rèn)為后續(xù)的微動(dòng)參數(shù)估計(jì)不受影響.圖7中γ隨參數(shù)相對(duì)誤差的增加逐漸降低,通過對(duì)γ的限制可以確定可容忍的LFM參數(shù)估計(jì)精度.從圖7可以看到,在時(shí),要求REf0<10%就可保證γ>0.9999,而本文方法求出的REf0=0.45%,完全滿足要求,而擬合法只有在微動(dòng)初相為π/2時(shí)才能滿足要求.當(dāng)然,通過增加處理信號(hào)的長度可以減小擬合法的誤差,結(jié)果列于表2.

表2中,T=1/fv表示一個(gè)微動(dòng)周期的信號(hào)長度,Tope表示算法運(yùn)行時(shí)間.可以看出,隨著信號(hào)長度的增加,估計(jì)誤差不斷降低,這是因?yàn)樾盘?hào)越長,端點(diǎn)的誤差對(duì)擬合的影響越小.但過長的信號(hào)也增加了時(shí)頻分析算法的運(yùn)算時(shí)間.而文中提出的FrFT方法只需兩個(gè)周期長的信號(hào)就可實(shí)現(xiàn)更低誤差的參數(shù)估計(jì),而且運(yùn)算時(shí)間遠(yuǎn)低于相同條件下的時(shí)頻分析擬合方法.

圖6 γ=0.9999時(shí)的相位波形對(duì)比Fig.6.Waveform comparison at γ=0.9999.

圖7 γ隨LFM參數(shù)相對(duì)估計(jì)誤差的變化Fig.7.Relation between γ and relative error.

表2 不同信號(hào)長度估計(jì)結(jié)果對(duì)比Table 2.Comparison of estimation results under dif f erent signal lengths.

4.1.2 微動(dòng)參數(shù)估計(jì)

令初始時(shí)粒子數(shù)Np為400,計(jì)算各粒子組的權(quán)重,只保留權(quán)重較高的400個(gè)粒子作為下一次迭代的初值.殘差累積長度為200個(gè)采樣點(diǎn),將觀測估計(jì)值的相對(duì)誤差連續(xù)穩(wěn)定且小于0.001作為迭代結(jié)束的條件.

首先分析粒子初始化對(duì)算法結(jié)果的影響,對(duì)于未知參數(shù),根據(jù)經(jīng)驗(yàn)設(shè)置初始化范圍:[50,150]Hz,Dv∈ [2,10]μm,ρ0∈ [0,π/2]rad.不同初始粒子數(shù)下的仿真結(jié)果如圖8.

從圖8(a)可以看出,對(duì)于多維參數(shù)而言,較少的粒子數(shù)不足以產(chǎn)生優(yōu)質(zhì)的初始參數(shù)組合,所有粒子的累積殘差都處在一個(gè)較大的水平(>350).當(dāng)粒子數(shù)增多時(shí),在相同的參數(shù)取值范圍下,多維參數(shù)組合的多樣性增加,產(chǎn)生接近真實(shí)值的概率也隨之增加,如圖8(b),出現(xiàn)了殘差小于200的初始粒子.根據(jù)3.4節(jié)的變粒子數(shù)初始化方法,只取權(quán)值高的部分粒子進(jìn)行PF,并與圖8(a)代表的固定粒子數(shù)方法進(jìn)行對(duì)比,用估計(jì)殘差反映估計(jì)參數(shù)和真實(shí)值的近似程度,得到圖8(c).圖8(b)代表的變粒子數(shù)方法只在第一次迭代時(shí)處理較多的粒子,計(jì)算量較大,之后迭代的計(jì)算量與傳統(tǒng)方法相同,兩算法整體計(jì)算量幾乎一致,但從對(duì)比結(jié)果可見前者的估計(jì)殘差明顯低于固定粒子數(shù)方法,且不到1次迭代就實(shí)現(xiàn)了收斂,而固定粒子數(shù)方法則由于初始值質(zhì)量較差,PF過程中殘差一直處于較高水平,不能確保正確的收斂,仿真結(jié)果體現(xiàn)了本文初始方法的有效性和必要性.

圖8 初始化值影響 (a)Np=400;(b)Np=4000;(c)不同初始值的估計(jì)結(jié)果對(duì)比Fig.8.The inf l uence of the initial value:(a)Np=400;(b)Np=4000;(c)etimation results of dif f erent initial values.

考慮到算法的收斂效率與建議分布σ2的取值密切相關(guān),在參數(shù)估計(jì)前,應(yīng)先確定方差σ2.根據(jù)3.3節(jié)對(duì)方差的分析,本文對(duì)方差和算法效率之間的關(guān)系進(jìn)行了仿真,并與常用Metropolis-Hastings(MH)算法的性能進(jìn)行了對(duì)比.最高迭代次數(shù)設(shè)為100,對(duì)各σ2值進(jìn)行100次蒙特卡羅仿真取平均.由于不同參數(shù)對(duì)應(yīng)方差不同,這里主要對(duì)進(jìn)行自適應(yīng)部分的方差進(jìn)行仿真,它對(duì)各參數(shù)的作用相同,可以概括地反映方差與迭代次數(shù)的關(guān)系,結(jié)果如圖9所示.

圖9 方差與算法迭代次數(shù)的關(guān)系Fig.9.Relationship between variance and iteration times.

仿真中用算法穩(wěn)定所需要的迭代次數(shù)反映算法效率,從圖9可以看出,相比于文獻(xiàn)[14]中采用的MH算法,在同樣條件下本文方法實(shí)現(xiàn)收斂的迭代次數(shù)更少,說明本文方法效率更高.MCMC建議分布方差存在最優(yōu)選擇獲得最佳的算法效率,對(duì)于本文設(shè)置的參數(shù),最優(yōu)方差在10?12.當(dāng)方差逐漸大于最優(yōu)值時(shí),算法仍可收斂,但迭代次數(shù)增加,因?yàn)榇蟮姆讲顚?dǎo)致粒子候選狀態(tài)的接受概率降低,但至少在范圍上可以覆蓋真值;當(dāng)方差逐漸小于最優(yōu)值時(shí),算法甚至在100次的迭代中無法實(shí)現(xiàn)穩(wěn)定,這是因?yàn)闋顟B(tài)更新幅度太小,粒子需要大量迭代才能接近真實(shí)值.

從圖10可以看出,本文所提算法采用了變粒子數(shù)、Gibbs多維更新算法和自適應(yīng)方差后,只用了不到10次迭代便使估計(jì)值收斂至真值,而只采用傳統(tǒng)MH算法去實(shí)現(xiàn)PF,需要40次左右的迭代才能實(shí)現(xiàn)收斂.如果未使用最佳σ2進(jìn)行估計(jì),算法效率將更低,需要更多次迭代才能實(shí)現(xiàn)收斂[14].下面對(duì)比兩種方法單次迭代運(yùn)算的時(shí)間復(fù)雜度.本文方法為其中τ為通過累積殘差計(jì)算粒子權(quán)重時(shí)設(shè)定的累積時(shí)間長度(具體分析見3.1節(jié));p為估計(jì)參數(shù)個(gè)數(shù);第二項(xiàng)表示重采樣帶來的時(shí)間復(fù)雜度;第三項(xiàng)表示狀態(tài)更新的時(shí)間復(fù)雜度.經(jīng)典的MH算法一次迭代的時(shí)間復(fù)雜度為兩者差別主要體現(xiàn)在多維粒子狀態(tài)的更新運(yùn)算上.記錄算法每次迭代的運(yùn)行時(shí)間,其中本文方法為0.0353 s,MH方法為0.0211 s,綜合考慮算法收斂所需的迭代次數(shù),前者耗時(shí)0.353 s(1次迭代),后者0.844 s(4次迭代),本文方法整體上用時(shí)更少,效率更高.從計(jì)算復(fù)雜度的角度來看,即使待估計(jì)參數(shù)的數(shù)量增加,也并不會(huì)過多地增加本文算法的計(jì)算量.

圖10 參數(shù)估計(jì)結(jié)果隨迭代次數(shù)的變化 (a)振動(dòng)幅度;(b)振動(dòng)頻率;(c)振動(dòng)初相;(d)相對(duì)誤差Fig.10.Estimation results changing with the iteration times:(a)Vibration amplitude;(b)vibration frequency;(c)vibration initial phase;(d)relative estimated residual.

4.2 實(shí)驗(yàn)數(shù)據(jù)參數(shù)估計(jì)

對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理估計(jì)目標(biāo)微動(dòng)參數(shù)來驗(yàn)證算法有效性.實(shí)驗(yàn)設(shè)置如圖11所示.

實(shí)驗(yàn)裝置為全光纖結(jié)構(gòu),采用波長1550 nm連續(xù)波激光器,輸出功率40 mW,線寬小于0.1 kHz.激光器輸出光通過保偏光纖分束器,分為90%和10%兩路.其中光強(qiáng)強(qiáng)的一路作為信號(hào)光,經(jīng)過口徑為20 mm的光纖準(zhǔn)直器照射到目標(biāo)上;光強(qiáng)弱的一路則經(jīng)過可調(diào)衰減器,作為本振光.信號(hào)接收采用口徑為80 mm的透射式望遠(yuǎn)鏡,將接收信號(hào)光和本振光通過2×2保偏光纖耦合器后輸入帶寬為80 MHz的InGaAs平衡探測器.輸出的中頻電信號(hào)由12 bit的A/D采集卡采集,采樣率在10—500 MHz可調(diào).目標(biāo)為振膜揚(yáng)聲器模擬微動(dòng)目標(biāo),揚(yáng)聲器驅(qū)動(dòng)信號(hào)的振動(dòng)頻率設(shè)為256 Hz,振幅約60μm,振動(dòng)初始相位由截取信號(hào)的時(shí)刻決定.利用本文算法對(duì)回波信號(hào)微動(dòng)參數(shù)進(jìn)行估計(jì),結(jié)果如圖12所示.

圖11 實(shí)驗(yàn)系統(tǒng)結(jié)構(gòu)圖Fig.11.Structure diagram of experiment system.

圖12 實(shí)驗(yàn)數(shù)據(jù)參數(shù)估計(jì)結(jié)果 (a)振動(dòng)幅度;(b)振動(dòng)頻率;(c)振動(dòng)初相Fig.12.Parameter estimation for experiment data:(a)Vibration amplitude;(b)vibration frequency;(c)vibration initial phase.

圖13 時(shí)頻分布IRT方法參數(shù)估計(jì)結(jié)果 (a)回波信號(hào)時(shí)頻分布;(b)IRT變換結(jié)果Fig.13.Parameter estimation results of IRT:(a)Time-frequency distribution;(b)IRT results.

圖14 重構(gòu)信號(hào)與實(shí)驗(yàn)信號(hào)對(duì)比 (a)IRT方法;(b)本文方法Fig.14.Comparison of reconstructed signal and experiment signal:(a)IRT;(b)proposed method.

從圖12(a)—(c)可以看出,在經(jīng)過約30次迭代后,本文算法對(duì)實(shí)驗(yàn)數(shù)據(jù)微動(dòng)參數(shù)的估計(jì)開始收斂,實(shí)現(xiàn)了對(duì)目標(biāo)微動(dòng)參數(shù)的估計(jì).與傳統(tǒng)的基于時(shí)頻分布逆Radon變換(IRT)的方法[24]進(jìn)行對(duì)比,驗(yàn)證本文的參數(shù)估計(jì)效果.

從圖13的IRT結(jié)果可以看出,基于時(shí)頻分布的參數(shù)估計(jì)只能把微動(dòng)參數(shù)的估計(jì)定位到一個(gè)大致區(qū)域,并不精確到具體值,這難以保證估計(jì)精度.而且,IRT定位的參數(shù)區(qū)域大小嚴(yán)重依賴于時(shí)頻分布的分辨率,會(huì)存在嚴(yán)重的誤差傳遞影響.圖14為利用估計(jì)參數(shù)重構(gòu)的信號(hào)和實(shí)驗(yàn)測得信號(hào)的對(duì)比,從圖14(a)可以看到,IRT估計(jì)的參數(shù)重構(gòu)信號(hào)與實(shí)驗(yàn)觀測信號(hào)差別較大,而圖14(b)中本文方法得到的兩波形基本重合.用兩者波形相似度γ和平均絕對(duì)誤差來對(duì)比參數(shù)估計(jì)的性能,結(jié)果列于表3.

從表3中可看出,參數(shù)化的PF方法在兩個(gè)評(píng)價(jià)標(biāo)準(zhǔn)上都優(yōu)于基于非參數(shù)化時(shí)頻分析的IRT,體現(xiàn)了所提方法在參數(shù)估計(jì)精度上的優(yōu)勢(STFT表示短時(shí)傅里葉變換).

表3 估計(jì)性能對(duì)比Table 3.Estimation performance comparison.

實(shí)驗(yàn)結(jié)果驗(yàn)證了算法對(duì)微動(dòng)參數(shù)估計(jì)的正確性.實(shí)際估計(jì)的微動(dòng)頻率在257 Hz附近,與設(shè)置值256 Hz有差別,這是由于驅(qū)動(dòng)源自身的誤差和揚(yáng)聲器自身對(duì)信號(hào)存在的響應(yīng)誤差,導(dǎo)致了實(shí)際振動(dòng)頻率和設(shè)定頻率有微小差距.

5 結(jié) 論

本文以基于微多普勒特征參數(shù)的目標(biāo)識(shí)別這一新應(yīng)用為研究背景,對(duì)實(shí)際運(yùn)動(dòng)目標(biāo)探測中遇到的平動(dòng)和微動(dòng)混合信號(hào)模型提出了一種精確的分離和參數(shù)估計(jì)方法.通過定義平動(dòng)補(bǔ)償信號(hào)與理想微多普勒信號(hào)波形相似度參數(shù),給出了對(duì)平動(dòng)參數(shù)估計(jì)精度的要求.利用FrFT理論結(jié)合參數(shù)域最窄帶寬搜索,可以實(shí)現(xiàn)滿足精度要求的平動(dòng)參數(shù)估計(jì),完成對(duì)LFM-SFM混合信號(hào)的精確平動(dòng)補(bǔ)償.以累積殘差為基礎(chǔ)定義權(quán)重計(jì)算函數(shù),再結(jié)合MCMC采樣方法對(duì)補(bǔ)償后的信號(hào)進(jìn)行靜態(tài)參數(shù)PF,實(shí)現(xiàn)了多維參數(shù)的同時(shí)估計(jì),避免了誤差傳遞影響,而且算法收斂效率較傳統(tǒng)MH方法有明顯提高.仿真和實(shí)驗(yàn)結(jié)果都證明了算法的有效性,體現(xiàn)了該算法對(duì)微動(dòng)參數(shù)的精確估計(jì)能力,這有利于對(duì)微弱振動(dòng)幅度目標(biāo)的探測以及對(duì)微動(dòng)參數(shù)接近的目標(biāo)微動(dòng)特征進(jìn)行分辨.實(shí)際中往往存在多目標(biāo)或多散射點(diǎn)產(chǎn)生的多分量混合信號(hào),下一步將在此基礎(chǔ)上針對(duì)多分量混合信號(hào)的參數(shù)估計(jì)開展研究.