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

郭力仁 胡以華 董驍 李敏樂

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

1 引 言

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

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

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

2 基于FrFT的平動參數(shù)估計

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

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

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

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

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

2.2 平動參數(shù)估計原理

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

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

圖1 LFM-SFM信號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)動目標(biāo)微多普勒信號表現(xiàn)為LFM-SFM混合的信號形式,經(jīng)過FrFT的結(jié)果如圖1所示.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

PF的終止條件可通過設(shè)置累積誤差閾值Rstop來進(jìn)行判斷.本文利用靜態(tài)參數(shù)PF進(jìn)行目標(biāo)微動參數(shù)估計的流程圖如圖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 仿真分析與實驗驗證

4.1 仿真信號參數(shù)估計

4.1.1 平動補(bǔ)償

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

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

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

圖4 LFM-SFM混合信號時頻分布Fig.4.Time-Frequency distribution of LFM-SFM signal.

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

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

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

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

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

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

圖6 γ=0.9999時的相位波形對比Fig.6.Waveform comparison at γ=0.9999.

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

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

4.1.2 微動參數(shù)估計

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

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

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

圖8 初始化值影響 (a)Np=400;(b)Np=4000;(c)不同初始值的估計結(jié)果對比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ù)估計前,應(yīng)先確定方差σ2.根據(jù)3.3節(jié)對方差的分析,本文對方差和算法效率之間的關(guān)系進(jìn)行了仿真,并與常用Metropolis-Hastings(MH)算法的性能進(jìn)行了對比.最高迭代次數(shù)設(shè)為100,對各σ2值進(jìn)行100次蒙特卡羅仿真取平均.由于不同參數(shù)對應(yīng)方差不同,這里主要對進(jìn)行自適應(yīng)部分的方差進(jìn)行仿真,它對各參數(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算法,在同樣條件下本文方法實現(xiàn)收斂的迭代次數(shù)更少,說明本文方法效率更高.MCMC建議分布方差存在最優(yōu)選擇獲得最佳的算法效率,對于本文設(shè)置的參數(shù),最優(yōu)方差在10?12.當(dāng)方差逐漸大于最優(yōu)值時,算法仍可收斂,但迭代次數(shù)增加,因為大的方差導(dǎo)致粒子候選狀態(tài)的接受概率降低,但至少在范圍上可以覆蓋真值;當(dāng)方差逐漸小于最優(yōu)值時,算法甚至在100次的迭代中無法實現(xiàn)穩(wěn)定,這是因為狀態(tài)更新幅度太小,粒子需要大量迭代才能接近真實值.

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

圖10 參數(shù)估計結(jié)果隨迭代次數(shù)的變化 (a)振動幅度;(b)振動頻率;(c)振動初相;(d)相對誤差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ù)據(jù)參數(shù)估計

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

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

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

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

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

圖14 重構(gòu)信號與實驗信號對比 (a)IRT方法;(b)本文方法Fig.14.Comparison of reconstructed signal and experiment signal:(a)IRT;(b)proposed method.

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

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

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

表3 估計性能對比Table 3.Estimation performance comparison.

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

5 結(jié) 論

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