基于最大似然的單通道交疊激光微多普勒信號參數(shù)分離估計(jì)?

郭力仁1)2) 胡以華1)2)? 王云鵬1)2) 徐世龍1)2)

1 引 言

在對運(yùn)動目標(biāo)進(jìn)行遙感探測時,目標(biāo)發(fā)動機(jī)振動、葉片轉(zhuǎn)動等微運(yùn)動的存在會對回波信號產(chǎn)生附加的頻率調(diào)制,產(chǎn)生所謂的“微多普勒效應(yīng)”[1].典型的目標(biāo)如坦克、汽車、飛機(jī)等都有自己獨(dú)特的微動形式和特定的微動參數(shù)范圍[2,3],通過準(zhǔn)確地提取微多普勒特征和精確估計(jì)微動參數(shù)可以實(shí)現(xiàn)對目標(biāo)的分類和精確識別.就發(fā)動機(jī)引起的振動而言,目標(biāo)的振動幅度往往在微米量級甚至更小.由于微多普勒調(diào)制效應(yīng)與波長成反比[4],一般雷達(dá)探測很難獲得明顯的微動特征,所以激光探測微多普勒效應(yīng)在準(zhǔn)確反演目標(biāo)特征和精確估計(jì)微動參數(shù)中就有著不可替代的優(yōu)勢.

微多普勒效應(yīng)的應(yīng)用關(guān)鍵在于目標(biāo)微動參數(shù)的精確估計(jì).在激光探測中,經(jīng)常存在多目標(biāo)或目標(biāo)等效為多個散射點(diǎn)的情況,而且目標(biāo)微動參數(shù)往往比較接近[5],在補(bǔ)償了目標(biāo)的主體運(yùn)動后,回波信號就成為單通道多分量(single channel multicomponent,SCMC)時頻域交疊的混合信號,目前還沒有有效的估計(jì)方法.對SCMC微多普勒信號的研究主要分為非參數(shù)化和參數(shù)化兩類[6].前者主要基于時頻分析方法,由于從時頻分布圖(timefrequency distribution,TFD)上可以直接對瞬時頻率變化規(guī)律進(jìn)行觀察,便于理解,使其成為過去十幾年間研究提取微多普勒特征的主要方向.但是多分量的時頻分布受交叉項(xiàng)和窗函數(shù)的影響嚴(yán)重[7],改進(jìn)的時頻分析方法[8,9]提高了時頻分辨率,但計(jì)算過程復(fù)雜,運(yùn)算量大.從包含多個分量的TFD中提取微動特征,傳統(tǒng)方法如峰值提取等[10]不再適用.文獻(xiàn)[11]根據(jù)微動瞬時頻率(instantaneous frequency,IF)曲線平滑無跳變的特點(diǎn),提出基于Viterbi的曲線跟蹤類方法,可以逐次提取時頻域上混疊的多條IF曲線,但該類方法對信噪比要求較高,對曲線斷點(diǎn)、交點(diǎn)處的處理性能差,且嚴(yán)重依賴前面時頻分析的效果.文獻(xiàn)[12]采用Radon變換從TFD中估計(jì)微動參數(shù),提高了算法的魯棒性,但是這類方法只能將參數(shù)定位到一個模糊的區(qū)域,估計(jì)精度無法保證,而且模型僅適用于兩參數(shù)估計(jì)的情況,參數(shù)增多計(jì)算量將呈指數(shù)增加,不適合對實(shí)際中的多參數(shù)情況進(jìn)行分離和估計(jì).

為提高目標(biāo)微動參數(shù)的估計(jì)精度,近幾年來參數(shù)化方法開始興起.文獻(xiàn)[13]提出循環(huán)平穩(wěn)方法，估計(jì)了微動頻率,但由于信號長度及循環(huán)周期性的限制,無法避免出現(xiàn)1/2或雙倍周期的估計(jì)誤差,不適合多分量信號.文獻(xiàn)[14]提出基于高階矩的方法估計(jì)微動參數(shù),雖然較圖像處理和正交匹配等方法耗時短,但依然不具備實(shí)時處理能力,且魯棒性較差.文獻(xiàn)[15]和文獻(xiàn)[16]則提出了正弦調(diào)頻匹配空間濾波法,通過依次估計(jì)各分量參數(shù),對信號進(jìn)行重建,實(shí)現(xiàn)分離,但其在參數(shù)估計(jì)中采用網(wǎng)格搜索,導(dǎo)致估計(jì)精度和計(jì)算量存在不可調(diào)和的矛盾.文獻(xiàn)[17]提出基于粒子濾波的SCMC參數(shù)估計(jì),但是對于多分量多參數(shù)的情況,對粒子數(shù)量要求量巨大,難以保證多維粒子的理想更新.文獻(xiàn)[18—20]分別引用盲源分離、總體經(jīng)驗(yàn)分解、局域均值分解等信號處理方法,這對于多通道或頻域不混疊的多分量微多普勒信號有一定效果,但是對SCMC這種極端欠定情況或時頻域存在交疊的信號則無法有效分離.Setlur等[21]基于最大似然估計(jì)的框架提出了迭代加權(quán)非線性最小二乘估計(jì)方法,得到了近似克拉美羅下界的估計(jì)精度,但算法步驟復(fù)雜,每次迭代計(jì)算量很大.最大似然方法已被證明是一種無偏漸進(jìn)最優(yōu)估計(jì)[22],在數(shù)據(jù)足夠多、信噪比較高時,MLE的性能可達(dá)到估計(jì)下界,遠(yuǎn)高于常用的非參數(shù)化方法.但是,目前的似然函數(shù)都是面向雷達(dá)信號建立,對于適用于微弱振幅探測的激光微多普勒信號,傳統(tǒng)方法構(gòu)建的似然函數(shù)非線性程度急劇增加[23],似然函數(shù)分布形狀由平滑的單峰變?yōu)槊芗亩喾錥24],使傳統(tǒng)迭代類或統(tǒng)計(jì)類的MLE在具體實(shí)現(xiàn)中無法得到有效的收斂.

為了從激光探測得到的單通道多分量且時頻交疊信號中精確估計(jì)目標(biāo)微動參數(shù),本文提出基于奇異值分解和最大似然框架的微動參數(shù)分離估計(jì)方法.通過精細(xì)化掃描信號周期點(diǎn)數(shù),改進(jìn)了信號矩陣的構(gòu)建方法,使基于奇異值比(singular value ratio,SVR)譜的微動頻率估計(jì)精度得到了數(shù)量級上的提高;分析對比了激光與微波探測微弱振幅微多普勒信號的特點(diǎn),基于頻域能量分布設(shè)計(jì)了新的似然函數(shù),在效果上相當(dāng)于對傳統(tǒng)似然函數(shù)分布進(jìn)行了平滑,而且由于噪聲在頻域上均勻分布,利用頻譜能量構(gòu)建似然函數(shù)有效提高了算法的抗噪性能;利用均值思想給出了最大似然估計(jì)的解析表達(dá),并用馬爾可夫-蒙特卡羅(Markov Chain Monte Carlo,MCMC)采樣方法解決了估計(jì)中的高維積分問題,實(shí)現(xiàn)了對參數(shù)進(jìn)行聯(lián)合估計(jì),提高了算法效率,同時避免了誤差傳遞問題.用該算法對仿真和實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理,得到了接近克拉美羅界的參數(shù)估計(jì)精度,驗(yàn)證了算法的可行性.通過與傳統(tǒng)基于時頻分布的逆Radon變換得到的估計(jì)結(jié)果對比,體現(xiàn)了本文方法在精確估計(jì)混合微動參數(shù)上的優(yōu)勢.

2 單通道多分量激光微多普勒信號模型

根據(jù)電磁散射理論,電大尺寸目標(biāo)在光學(xué)區(qū)的散射具有局部效應(yīng),目標(biāo)總的電磁散射可等效為一些局部散射點(diǎn)的疊加,這些局部的散射源就叫作等效散射中心.對于激光探測的目標(biāo)微多普勒效應(yīng)回波信號,其波長遠(yuǎn)小于目標(biāo)尺寸,采用電磁散射在光學(xué)區(qū)的等效散射中心理論是合理的.

Chen[4]基于散射中心模型對振動、轉(zhuǎn)動、錐動等典型微動的微多普勒效應(yīng)建立了回波數(shù)學(xué)模型.經(jīng)過后續(xù)化簡和整合,這些微動的散射中心模型都具備基本的正弦調(diào)頻形式.對于實(shí)際探測中常遇到的單通道多目標(biāo)或目標(biāo)具有多個等效散射中心的多分量情況,包含微多普勒效應(yīng)的基帶信號離散形式可寫為

式中 λ表示激光波長;f0k=fvk/fs,fs為采樣率;Dvk,fvk,ρ0k分別表示目標(biāo)第k個散射中心對應(yīng)的微動幅度、頻率、初始相位;Ak和θk表示信號強(qiáng)度和初相;K 表示回波中疊加的分量個數(shù);e(n)表示無噪聲干擾信號;w(n)為噪聲;α,β是目標(biāo)相對于雷達(dá)的方位角和俯仰角,對于目標(biāo)微動參數(shù)來說是慢變化量,在信號處理過程中認(rèn)為是常數(shù).另外考慮到激光探測發(fā)散角較小,認(rèn)為照射光斑內(nèi)各散射點(diǎn)的α,β相同.為簡便起見,在計(jì)算過程中可將散射點(diǎn)振動方位角和俯仰角,目標(biāo)方位角和俯仰角都設(shè)為0.

3 時頻交疊激光微多普勒信號的參數(shù)分離估計(jì)

3.1 基于SVR譜的微動頻率估計(jì)

根據(jù)奇異值分解理論,若矩陣S的各個行向量間的比值為一個常數(shù),則S的奇異值中只有第一項(xiàng)不為0.所以對于周期信號,若將其構(gòu)建為矩陣形式,則當(dāng)矩陣每行長度與信號周期長度接近時,奇異值比σ1/σ2就會是一個較大值.奇異值比譜方法就是通過不斷改變構(gòu)建矩陣的列數(shù),尋找最大的σ1/σ2來實(shí)現(xiàn)對信號周期的檢測[25].激光微多普勒信號具有正弦調(diào)頻信號的形式,調(diào)制信號的周期性使回波基帶信號也具備周期性,所以具備通過奇異值比譜法來估計(jì)調(diào)制周期的條件.

如果直接用完整的信號構(gòu)建矩陣,矩陣列數(shù)與信號周期之間的誤差會隨著行數(shù)增加不斷累加,這常常會導(dǎo)致SVR方法失效.所以,這里提出固定列數(shù)的周期掃描的改進(jìn)方法來實(shí)現(xiàn)周期的估計(jì).構(gòu)建矩陣可寫為

式中dm=round(m·d),d為假設(shè)的一個周期內(nèi)的點(diǎn)數(shù),根據(jù)實(shí)際情況在目標(biāo)微動周期的范圍[Tmin,Tmax]內(nèi)對d進(jìn)行掃描,構(gòu)建不同列數(shù)的矩陣;m=0,1,···,M,M=floor(N/d)?1,表示總信號長度N內(nèi)包含的周期數(shù);L為每個周期內(nèi)截取的固定長度,可取L=round(d)以充分利用信號信息.通過SVR估計(jì)的信號周期為

(3)式表示取奇異值比最大時對應(yīng)的d作為周期長度,?T代表估計(jì)值,Ts=1/fs表示采樣時間間隔.若d在掃描中取整數(shù),即掃描步長?d=1,則周期估計(jì)誤差范圍在[0—Ts/2],可見信號采樣率越高,則估計(jì)誤差越小.對于確定fs的信號,為進(jìn)一步減小估計(jì)誤差,可減小?d,經(jīng)過取整運(yùn)算后,SVR中可能存在最大值對應(yīng)連續(xù)多個位置的情況,則取其中間位置作為最終的點(diǎn)數(shù)估計(jì).對應(yīng)的微動頻率為

對于包含多個微多普勒分量的回波信號,若各分量幅度相似,則一次SVR可估計(jì)出各分量的微動頻率;若各分量信號幅度相差較大,則一次SVR只得到主周期分量的微動頻率,繼續(xù)估計(jì)該分量剩余其他參數(shù)后,從回波信號中將其重構(gòu)信號去除,用同樣方法可估計(jì)其他分量的微動頻率.通過仿真得到各分量信號幅度比與SVR峰值比的關(guān)系如圖1所示.

圖1 分量幅值比和SVR峰值比的關(guān)系Fig.1.Relationship between the signal amplitude ratio and the peak ratio of SVR spectrum.

圖1 中縱坐標(biāo)SVR1和SVR2分別表示奇異值比譜中周期較長的分量和周期較短的分量對應(yīng)的峰值大小(即σ1/σ2的值),橫坐標(biāo)代表SVR1和SVR2對應(yīng)分量的幅值比.當(dāng)信號幅值比超出圖中范圍時,SVR中只有一個明顯的峰值.圖1的關(guān)系可為下一節(jié)中設(shè)計(jì)適合激光微多普勒信號的似然函數(shù)提供參數(shù)選擇的依據(jù).

3.2 基于最大似然的微動幅度和初相估計(jì)

在3.1節(jié)估計(jì)了主周期分量微動頻率的基礎(chǔ)上,繼續(xù)利用MCMC方法搜索估計(jì)該分量的微動幅度和初相,用ψ=[Dvk,ρ0k]T表示第k個分量的微動參數(shù)矢量.MCMC是產(chǎn)生服從特定概率分布的樣本,所以首先對參數(shù)的概率密度函數(shù)進(jìn)行分析.

3.2.1 激光微多普勒信號的均值似然函數(shù)

認(rèn)為(1)式中的噪聲是服從N(0,σ2)的加性高斯白噪聲,那么信號矢量服從s～N(μ(ψ),C(ψ)),信號的概率密度函數(shù)可寫為[26]:

L(s;ψ)即為待估計(jì)參數(shù)的似然函數(shù),可對其求對數(shù)進(jìn)行化簡,

(8)式中,等號右邊前兩項(xiàng)為常數(shù)項(xiàng),求似然函數(shù)的最大值等同于求第三項(xiàng)的最小值,此時參數(shù)的似然估計(jì)值為:

令J(ψ)=(s?e)H(s?e)/2σ2,把信號e分解為振幅項(xiàng)和相位項(xiàng)兩個部分,(9)式可改寫為:

在求解(10)式時,Hk(ψ)和Ak的估計(jì)過程解耦分別進(jìn)行估計(jì)與直接求聯(lián)合MLE是等價的.利用加權(quán)最小二乘估計(jì)A[24],有

將(11)式代入(10)式得到

因?yàn)镻=P2=PH,所以(12)式可化簡為

一般直接用(14)式的似然函數(shù)形式構(gòu)建參數(shù)ψ的概率分布函數(shù),再用MCMC方法產(chǎn)生服從分布的樣本實(shí)現(xiàn)對參數(shù)的最大似然估計(jì),這對于傳統(tǒng)微波探測或是平穩(wěn)信號處理是可行的.但是對于激光微多普勒信號,由于探測波長極短,信號由雷達(dá)探測中的弱調(diào)制變?yōu)樯疃日{(diào)制,似然函數(shù)的非線性程度急劇增加,概率分布變?yōu)閺?fù)雜的多峰形狀.這將導(dǎo)致馬爾可夫鏈在產(chǎn)生平穩(wěn)分布時會錯誤收斂至局部最大值,最終估計(jì)錯誤.微多普勒信號參數(shù)估計(jì)概率分布與探測波長之間的關(guān)系如圖2.

圖2中的振動幅度設(shè)置為5μm,與實(shí)驗(yàn)中目標(biāo)實(shí)際振動幅度量級相同.圖2(a)為激光波段下的似然函數(shù),其分布為密集的多峰形狀.隨著波長的增加似然函數(shù)峰值分布由密變疏,如圖2(a)—圖2(d)所示.類似圖2(c)平滑的單峰的似然函數(shù)可以保證傳統(tǒng)的迭代搜索方法實(shí)現(xiàn)對參數(shù)的高精度估計(jì),但在圖2(a)所示的似然函數(shù)下,這些方法將失效.此外,隨著波長的進(jìn)一步增加,似然函數(shù)峰值消失,這意味著過長的波長失去了對微弱振動的探測和估計(jì)能力.所以,處理激光信號時需要對似然函數(shù)進(jìn)行改變.以兩分量為例對(14)式中的似然函數(shù)進(jìn)行分析:

圖2 不同波長下傳統(tǒng)似然函數(shù)分布 (a)λ=10?6 m;(b)λ=10?5 m;(c)λ=10?4 m;(d)λ=10?3 mFig.2.Traditional PDF distribution in diff erent wavelength:(a)λ =10?6 m;(b)λ =10?5 m;(c)λ =10?4 m;(d)λ =10?3 m.

首先通過SVR方法確定出一個分量的微動頻率?f01,并以此作為先驗(yàn)信息計(jì)算Dv1和ρ01參數(shù)域上的似然函數(shù)分布.當(dāng)估計(jì)值?Dv1和?ρ01與信號分量1實(shí)際的參數(shù)值接近時,sl1(n)的頻譜能量集中在零頻附近很窄的帶寬里,能量峰值極高;而此時,由于分量2與分量1的微動參數(shù)不同,在參數(shù)域中sl2(n)的頻譜能量只能分布在一個較寬的帶寬內(nèi);slnoi(n)則均勻分布在整個頻譜中.當(dāng)估計(jì)值完全與實(shí)際值相同時sl1(n)為常數(shù),頻譜為沖擊響應(yīng)形狀.所以,可以把似然函數(shù)頻譜中占總能量的比例為η時的能帶寬度作為評價參數(shù)估計(jì)值是否接近真值的標(biāo)準(zhǔn),構(gòu)建新的似然函數(shù),

式中F{·}表示傅里葉變換,分母為計(jì)算頻譜總能量,η為設(shè)定sl1(n)部分占總能量的比例,W表示E[sl1(n)]所占頻帶寬帶.由于信號時域和頻域能量守恒,η的表達(dá)式可寫為

A1/A2的值可根據(jù)圖1中的關(guān)系進(jìn)行確定.新似然函數(shù)各分量頻譜關(guān)系和參數(shù)示意如圖3所示.

圖3 改進(jìn)似然函數(shù) (a)參數(shù)示意圖;(b)頻譜能量累積分布函數(shù);(c)不同相對誤差(RE)下的對比Fig.3 The improved likelihood function:(a)The schematic diagram of parameters;(b)cumulative distribution on spectrum energy;(c)comparison between diff erent relative error(RE).

圖3 (b)為圖3(a)的累積分布,圖中參數(shù)估計(jì)值與真實(shí)值的相對誤差為5%,信噪比為5 d B.圖3(c)對比結(jié)果說明隨著相對誤差的減小,E[sl1(n)]對應(yīng)的寬度W也將逐漸減小,所以用W反映參數(shù)估計(jì)的精度是合理的.由于噪聲在頻域上均勻分布,利用頻譜能量構(gòu)建似然函數(shù)避免了個別時刻噪聲的突起對傳統(tǒng)基于時域構(gòu)建似然函數(shù)的影響,有效提高了算法的抗噪性能.這時,引入均值似然估計(jì)給出多參數(shù)MLE的解析表達(dá)形式:

式中p(Dvk,ρ0k)為歸一化壓縮似然函數(shù),其表達(dá)式為

由于(16)式對似然函數(shù)形式進(jìn)行了改進(jìn),得到了平滑單峰形狀的分布,不再需要專門設(shè)置壓縮指數(shù)γ來突出全局最大值以保證收斂,降低了工作量,這里γ取1.

3.2.2 MCMC算法實(shí)現(xiàn)

對于多維參數(shù)的MCMC抽樣,采用Gibbs方法來具體實(shí)現(xiàn).用π(x)表示目標(biāo)分布函數(shù),x=(x1,···,xm)是m維參數(shù)矢量,q(x)表示建議分布函數(shù). 在產(chǎn)生馬爾可夫鏈的過程中,t時刻第i個參數(shù)的狀態(tài)xti根據(jù)條件建議分布q(→ x?i|) 轉(zhuǎn)移到狀態(tài)x?i, 產(chǎn)生t+1時刻的候選樣本,=(,···,,,···,);樣本的接受概率A(→ x?i);產(chǎn)生服從0—1均勻分布的隨機(jī)數(shù)u與A(→|)比較,若u小于接受概率,則馬爾可夫鏈t+1時刻狀態(tài)更新為,否則=.接受概率由概率轉(zhuǎn)移函數(shù)的細(xì)致平衡等式推出[27]:

當(dāng)候選狀態(tài)使目標(biāo)分布概率密度增大時,以概率1更新為這個狀態(tài),否則,以概率π(x?i)/π(xti)進(jìn)行更新.所以,t+1時刻狀態(tài)的實(shí)際轉(zhuǎn)移概率為

由于T(xti,x?i)滿足細(xì)致平衡條件,所以在足夠多次的迭代后,馬爾可夫鏈最終將收斂于目標(biāo)分布π(x).將(20)式作為目標(biāo)概率分布,對激光微多普勒參數(shù)ψ進(jìn)行似然估計(jì)的MCMC算法具體實(shí)現(xiàn)步驟為:

1)初始化[D1vk,ρ10k];

2)t時刻,根據(jù)條件建議分布的轉(zhuǎn)移概率得到參數(shù) ρ0k的候選狀態(tài) ρ?0k;

3)計(jì)算接受概率

4)產(chǎn)生u ～ U(0,1), 若u

5)根據(jù)條件建議分布的轉(zhuǎn)移概率得到參數(shù)Dvk的候選狀態(tài)D?vk;

6)計(jì)算接受概率

7)產(chǎn)生u′～ U(0,1),若則

8)將步驟(2)—(7)重復(fù)M次,將達(dá)到收斂之前的τ?1個樣本矢量burn-in,對剩下的樣本求均值得到參數(shù)的估計(jì)為

和

算法建議分布q(x)的選擇直接決定了馬氏鏈的收斂效果.若q(x)分布過于分散,則馬氏鏈長期得不到更新,需要很多次迭代才能搜索到最優(yōu)值,降低算法效率;若分布過于集中雖然能使接受概率高,但每次狀態(tài)更新跨度太小,同樣需要大量的迭代才能實(shí)現(xiàn)馬氏鏈的收斂.Gelman等建議將接受概率控制在0.15—0.5比較合適,所以具體操作中可通過對接受概率的監(jiān)視來調(diào)制.

3.3 信號幅度和初相的估計(jì)

根據(jù)3.1節(jié)和3.2節(jié)估計(jì)的第k分量的微動參數(shù)重構(gòu)信號

則(13)式的似然函數(shù)形式可改寫為

式中

表示觀測混合信號中的第k個分量單位調(diào)制信號.當(dāng)估計(jì)值與實(shí)際參數(shù)相等時,重構(gòu)分量與實(shí)際分量抵消,等式最右邊第一項(xiàng)等于AkN exp(jθk);由于重構(gòu)分量與其他分量參數(shù)不同,所以第2項(xiàng)相當(dāng)于是一個寬帶正弦調(diào)頻信號求和,在一個調(diào)制周期內(nèi)其值正好為0,N的值一般根據(jù)SVR結(jié)果選取整周期長度;噪聲和重構(gòu)信號不相關(guān),經(jīng)過長時間累積后,第3項(xiàng)也近似為0.此時,可得到信號幅度和相位的估計(jì)為:

估計(jì)了信號幅度和相位后得到了分量k所有參數(shù),可重建k分量?Ak?sk(n)ej?θk,將其從初始混合信號中去除,用同樣方法繼續(xù)對剩余信號參數(shù)進(jìn)行估計(jì).完整的SCMC混合信號參數(shù)估計(jì)和分離流程如圖4所示.

圖4 SCMC激光微多普勒信號分離和參數(shù)估計(jì)流程Fig.4.The fl ow chart of the estimation and separation method for the SCMC signal.

3.4 激光M D信號參數(shù)估計(jì)的克拉美羅界

克拉美羅界(Cramer-Rao bounds,CRB)是所有無偏估計(jì)方法所能達(dá)到的估計(jì)方差的下限,其取值只與信號本身有關(guān),不受估計(jì)方法的影響,經(jīng)常被用作評價參數(shù)估計(jì)精度的標(biāo)準(zhǔn).對于矢量參數(shù)的CRB,由Fisher信息矩陣求逆后的對角線元素確定,即

式中var(?ψi)表示第i個參數(shù)估計(jì)的方差;ψ表示包含待估計(jì)微動參數(shù)的矢量;I(ψ)為Fisher信息矩陣,其定義為[28]

將微多普勒信號模型代入(28)式,可得到各參數(shù)的Fisher信息為:

(29)式—(31)式構(gòu)成微動參數(shù)的Fisher信息矩陣I(ψ),其中ω0=2πfv/fs,對I(ψ)求逆可得到參數(shù)的CRB.上式中的Fisher信息取值與波長成反比,這意味著CRB與波長成正比.所以,探測中的波長越短,參數(shù)估計(jì)的CRB越低,能夠達(dá)到的估計(jì)的精度也就越高,這證明了激光微多普勒探測較微波探測在精確估計(jì)上具有優(yōu)勢.

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

4.1 仿真結(jié)果與分析

以兩分量混合的單通道時頻域交疊激光微多普勒信號為例,對本文所提分離和參數(shù)估計(jì)方法進(jìn)行仿真驗(yàn)證.分量1的微動參數(shù):振動頻率、振動幅度和初始相位分別設(shè)為256 Hz,5×10?5m和π/3 rad,信號幅度和相位分別為A1=3,θ1=π/4 rad;分量2對應(yīng)的微動參數(shù)設(shè)置為300 Hz,2×10?5m和π/6 rad,信號幅度和相位設(shè)為A2=1,θ2=π/5 rad.設(shè)信號噪聲為零均值高斯白噪聲,信噪比已知為20 d B,激光波長λ=1550 nm,采樣率為fs=312 kHz.

首先根據(jù)3.1節(jié)介紹的方法計(jì)算信號的奇異值比譜,結(jié)果如圖5(a)所示.掃描周期點(diǎn)數(shù)從n=1001點(diǎn)開始,?d取0.1,圖中只有一個明顯的峰值在第2178點(diǎn),所以由SVR估計(jì)的微動頻率為fs·[1001+(2178?1)×?d]?1=256.0105 Hz,定義為fv1.將fv1作為已知量計(jì)算,利用改進(jìn)的似然函數(shù)計(jì)算Dv1和ρ01的概率密度分布,如圖5(b).與圖2(a)相比可以看出,對激光微多普勒信號的似然函數(shù)進(jìn)行重新設(shè)計(jì)后,可以把傳統(tǒng)密集多峰的分布形狀變?yōu)槠交瑔畏宓姆植?有利于參數(shù)的搜索和算法的收斂.此外還降低了對參數(shù)初始化精度的要求,只要設(shè)置范圍合理即可,不需要提前進(jìn)行精確計(jì)算.圖5(c)和圖5(d)給出了本文方法和傳統(tǒng)方法對激光微動參數(shù)進(jìn)行最大似然估計(jì)的收斂結(jié)果,圖中曲線標(biāo)注“LF”表示似然函數(shù).對比圖中四種方法可以看出利用本文設(shè)計(jì)的似然函數(shù),結(jié)合Gibbs算法具有最高的效率,在本文方法下參數(shù)Dv1和ρ01收斂到真值需要的迭代次數(shù)最少.而MH方法由于在多參數(shù)情況下的更新效率低,所以收斂速度較慢.相比之下,利用傳統(tǒng)方法計(jì)算激光微動信號的似然函數(shù),由于密集局部峰值的存在,不管利用哪種算法都不能正確估計(jì)出微動參數(shù),收斂結(jié)果停留在初始值附近的某個極大值處.把估計(jì)出的微動參數(shù)代入(23)式重構(gòu)分量1的單位調(diào)制信號,再根據(jù)(25)式和(26)式計(jì)算分量1的真實(shí)幅度和相位.

從總的回波信號中減去分量1,按照相同的步驟對剩余分量進(jìn)行參數(shù)估計(jì).剩余分量的SVR如圖6(a)所示.圖中峰值在第391點(diǎn)處,可求出對應(yīng)分量2的微動頻率為300 Hz.圖6(b)為分量2對應(yīng)微動參數(shù)的概率密度分布,剩余信號的概率分布仍然是平滑單峰形狀,可以保證分量2參數(shù)的正確收斂.圖6(c)和圖6(d)為分量2微動參數(shù)的MCMC估計(jì)對比結(jié)果,整體結(jié)果與圖5一致,利用傳統(tǒng)的似然函數(shù)得到的馬氏鏈被困在初始值附近,而改進(jìn)的似然函數(shù)確保了正確的收斂.此外,與圖5對比還可發(fā)現(xiàn)分量2估計(jì)的馬爾可夫鏈?zhǔn)諗康秸嬷档乃俣容^分量1更快.這是因?yàn)榉至?參數(shù)被精確地估計(jì),從原始信號中移除重構(gòu)的分量1后,分量2成為剩余信號的主要成分,沒有其他信號分量的影響,所以收斂更快.

圖5 分量1的微動參數(shù)估計(jì)結(jié)果 (a)奇異值比譜;(b)參數(shù)D v1和ρ01的概率密度分布;(c)振動幅度;(d)振動初始相位Fig.5.The estimation results of component 1 in the simulation data:(a)SVR;(b)the PDF of D v1 andρ01;(c)vibration amplitude;(d)vibration initial phase.

圖6 分量2的微動參數(shù)估計(jì)結(jié)果 (a)奇異值比譜;(b)參數(shù)D v1和ρ01的概率密度分布;(c)振動幅度;(d)振動初始相位Fig.6.The estimation results of component 2 in the simulation data:(a)SVR;(b)the pdf of D v1 andρ01;(c)vibration amplitude;(d)vibration initial phaseξ.

各分量參數(shù)估計(jì)的具體結(jié)果和相對誤差如表1所列.表中RE=(?ψ?ψ)×100%/ψ,表示相對估計(jì)誤差;γ=|E[sk(n)?sk(n)]|·|E[sk(n)]2·E[?sk(n)]2|?1/2表示波形相似度,通過對比重構(gòu)信號?sk(n)與實(shí)際信號在波形上的相似程度來反映參數(shù)估計(jì)的準(zhǔn)確度,γ越接近1則波形相似度越高,表明參數(shù)估計(jì)越接近真實(shí)值.從表1可以看出,本文方法對目標(biāo)微動參數(shù)估計(jì)的相對誤差都在10?5量級,具有較高的精度.

為定量分析本文方法的參數(shù)估計(jì)精度,計(jì)算不同信噪比下的估計(jì)均方誤差MSE=并與CRB對比.令信噪比變化范圍為?5—35 dB,每隔5 d B進(jìn)行100次蒙特卡羅實(shí)驗(yàn)計(jì)算MSE,對比結(jié)果如圖7.

從圖7中可以看出,由于兩個分量中參數(shù)值的差異導(dǎo)致各信噪比下的CRB也存在差異.但隨著信噪比的增加,分量1和分量2的微動參數(shù)估計(jì)MSE都逐漸接近各自的克拉美羅界,說明本文基于MCMC和重新定義的似然函數(shù)來實(shí)現(xiàn)微動參數(shù)的最大似然估計(jì)具有最優(yōu)的估計(jì)性能.此外,仔細(xì)對比可發(fā)現(xiàn),各信噪比下分量2中參數(shù)的估計(jì)誤差都略大于分量1的,這是因?yàn)閷Ψ至?的估計(jì)依賴于分量1的結(jié)果,會存在一定誤差傳遞的干擾.

表1 參數(shù)估計(jì)結(jié)果Table 1.Parameter estimation result.

圖7 微動參數(shù)估計(jì)均方誤差與克拉美羅界對比 (a)分量1參數(shù)估計(jì)性能;(b)分量2參數(shù)估計(jì)性能Fig.7.Comparison of MSE and CRB:(a)Parameters of component 1;(b)parameters of component 2.

4.2 實(shí)驗(yàn)驗(yàn)證與對比

利用實(shí)驗(yàn)數(shù)據(jù)估計(jì)目標(biāo)微動參數(shù)來驗(yàn)證算法有效性,實(shí)驗(yàn)設(shè)置如圖8.

實(shí)驗(yàn)結(jié)構(gòu)為全光纖相干激光探測結(jié)構(gòu),采用波長1550 nm連續(xù)波激光器,輸出功率40 mW,線寬小于0.1 kHz.激光通過90/10的保偏光纖分束器,9%一路作為信號光,經(jīng)過光纖擴(kuò)束系統(tǒng)后照射到兩個目標(biāo)上;另一路經(jīng)過可調(diào)衰減器,作為本振光.接收端采用口徑為80 mm的透射式望遠(yuǎn)鏡接收兩個目標(biāo)的散射光信號,回波光和本振光接入平衡探測器,然后由A/D采集卡采集.用振膜揚(yáng)聲器和電驅(qū)動音叉模擬微動目標(biāo),振動頻率分別設(shè)為300和256 Hz,振幅在微米量級,振動初始相位由截取信號的時刻決定.實(shí)驗(yàn)信號的時域波形與時頻分布如圖9.

圖8 實(shí)驗(yàn)結(jié)構(gòu) (a)系統(tǒng)框圖;(b)實(shí)物圖Fig.8.Experiment system:(a)Structure diagram;(b)experiment setup.

圖9 實(shí)驗(yàn)信號 (a)歸一化時域波形;(b)STFT時頻分布Fig.9.Experiment data:(a)Normalized time-domain waveforms;(b)STFT time-frequency distribution.

從圖9(a)中可以看出目標(biāo)振動微多普勒效應(yīng)對信號的調(diào)制效應(yīng)使兩個分量的混合信號在時域上出現(xiàn)了周期性,該周期與目標(biāo)振動周期一致,所以具備利用SVR譜來估計(jì)目標(biāo)振動頻率的條件.圖9(b)為混合信號的時頻分布,兩個分量的微多普勒特征在時頻域相互交疊,這一現(xiàn)象將導(dǎo)致傳統(tǒng)的信號分離和估計(jì)方法失效,體現(xiàn)了研究SCMC時頻交疊信號處理方法的必要性.首先利用傳統(tǒng)逆Radon變換對時頻圖進(jìn)行處理,估計(jì)各分量微動參數(shù),作為對比.逆Radon變換估計(jì)結(jié)果如圖10所示.

從圖10的逆Radon變換結(jié)果可以看出,基于時頻分布的參數(shù)估計(jì)只能把微動參數(shù)的估計(jì)定位到一個大致區(qū)域,并不精確到具體值,這難以保證估計(jì)精度.而且,逆Radon定位的參數(shù)區(qū)域大小嚴(yán)重依賴于時頻分布的分辨率,會存在嚴(yán)重的誤差傳遞影響.下面用本文方法對實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理,驗(yàn)證算法的有效性.參數(shù)估計(jì)結(jié)果如圖11和圖12所示.

從圖11(a)奇異值比譜中除了得到各分量的振動周期外,還可以得到最高峰和次高峰值比SVR1/SVR2=1.5098,根據(jù)圖1給出的關(guān)系可以得到對應(yīng)分量的幅值比約為2.1,代入(17)式計(jì)算得η=0.815,由此得到的概率分布是理想的平滑單峰形狀,如圖11(b)所示,可以保證算法快速準(zhǔn)確的收斂.利用MCMC方法實(shí)現(xiàn)對參數(shù)Dv1和ρ01的最大似然估計(jì),結(jié)果如圖11(c)和圖11(d),Markov鏈迅速得到了收斂,表明鎖定了信號真實(shí)的參數(shù),整個參數(shù)最大估計(jì)過程不到1 s,具備實(shí)時處理能力.從探測信號中去除分量1,繼續(xù)對剩余信號進(jìn)行估計(jì).分量2的周期在圖12(a)中有明顯的峰值,但SVR譜中除了分量2的主峰值還有分量1的殘留分量,這是由于實(shí)驗(yàn)中驅(qū)動源自身的不穩(wěn)定和目標(biāo)對驅(qū)動響應(yīng)不穩(wěn)定使實(shí)驗(yàn)?zāi)繕?biāo)對實(shí)際振動模擬不理想造成的.實(shí)驗(yàn)中的這些不確定性相當(dāng)于給信號增加了頻率噪聲,會影響參數(shù)估計(jì)的精度.但從圖12(b)的概率分布形狀以及圖12(c)和圖12(d)的收斂情況看,這并不會影響該方法對參數(shù)2中分量的估計(jì)效果.算法對實(shí)驗(yàn)數(shù)據(jù)處理得到的收斂效果與仿真分析一致,驗(yàn)證了算法的有效性.為定量分析對比參數(shù)估計(jì)的準(zhǔn)確度,這里用波形相似度γ來驗(yàn)證參數(shù)估計(jì)的準(zhǔn)確度,并與傳統(tǒng)非參數(shù)化方法進(jìn)行對比,結(jié)果如表2.

圖10 逆Radon變換估計(jì)結(jié)果 (a)分量1參數(shù);(b)分量2參數(shù)Fig.10.Results of iRadon:(a)Parameters of component 1;(b)parameters of component 2.

圖11 實(shí)驗(yàn)數(shù)據(jù)分量1的微動參數(shù)估計(jì)結(jié)果 (a)奇異值比譜;(b)參數(shù)D v1和ρ01的概率密度分布;(c)振動幅度;(d)振動初始相位Fig.11.The estimation results of component 1 in the experiment data:(a)SVR;(b)the PDF of D v1 andρ01;(c)vibration amplitude;(d)vibration initial phase.

圖12 實(shí)驗(yàn)數(shù)據(jù)分量2的微動參數(shù)估計(jì)結(jié)果 (a)奇異值比譜;(b)參數(shù)D v1和ρ01的概率密度分布;(c)振動幅度;(d)振動初始相位Fig.12.The estimation results of component 2 in the experiment data:(a)SVR;(b)the PDF of D v1 andρ01;(c)vibration amplitude;(d)vibration initial phase.

表2 重構(gòu)信號波形相似度Table 2.Waveform similarity comparison.

從表2可以看出,利用最大似然估計(jì)得到的微動參數(shù)重構(gòu)波形,其相似度達(dá)到0.9以上,遠(yuǎn)高于基于傳統(tǒng)時頻分布的逆Radon變換重構(gòu)波形的相似度,驗(yàn)證了參數(shù)估計(jì)的準(zhǔn)確性.說明本文所提方法更有利于實(shí)現(xiàn)對微動參數(shù)的精確估計(jì),這為基于微動特征的目標(biāo)分類和精細(xì)識別奠定了基礎(chǔ).

5 結(jié) 論

本文針對激光微多普勒探測中存在時頻域交疊的單通道多分量混合信號,提出了基于最大似然框架的參數(shù)估計(jì)方法和基于信號重構(gòu)的分離方法.激光對微弱振動微多普勒效應(yīng)的探測有不可替代的優(yōu)勢,但激光微多普勒效應(yīng)對信號的深度調(diào)制也使傳統(tǒng)似然函數(shù)性質(zhì)發(fā)生了本質(zhì)變化.對此文中給出了適合激光信號的似然函數(shù)計(jì)算方法,可以得到理想的概率密度分布形式,同時降低了初始化要求,提高了MLE的抗噪能力.通過改進(jìn)奇異值比譜方法精細(xì)掃描混合信號的周期性,得到目標(biāo)微動頻率信息.利用均值似然函數(shù)給出了MLE的解析表達(dá),并在求得微動頻率的基礎(chǔ)上利用MCMC方法實(shí)現(xiàn)對剩余微動參數(shù)的最大似然估計(jì),解決了復(fù)雜的高維積分的問題.利用估計(jì)的微動參數(shù)重構(gòu)單位調(diào)制信號,估計(jì)了信號的幅度和初始相位,通過對各個微動分量所包含參數(shù)的分離估計(jì),實(shí)現(xiàn)了對單通道混合信號中微動特征參量的精確提取.仿真和實(shí)驗(yàn)數(shù)據(jù)的處理結(jié)果驗(yàn)證了所提算法的有效性.該方法可實(shí)現(xiàn)對單通道混合微多普勒信號參數(shù)的精確估計(jì),為基于微動參數(shù)的目標(biāo)分類和識別以及對目標(biāo)振動精細(xì)成像提供了可能.

[1]Chen V C 2006 IEEE Trans.Aerosp.Electron.Syst.42 2

[2]Jiang Y 2014 Ph.D.Dissertation(Xi’an:Xidian University)(in Chinese)[姜悅 2014博士學(xué)位論文 (西安:西安電子科技大學(xué))]

[3]Yang J,Liu C,Wang Y 2015 IEEE Trans.Geosci.Remote Sens.53 920

[4]Chen V C 2011 The Micro-Doppler Eff ect in Radar(Fitchburg:Artech House)pp15–17

[5]Wang T,Tong C M,Li X M 2015 Acta Phys.Sin.64 058401(in Chinese)[王童,童創(chuàng)明,李西敏2015物理學(xué)報(bào)64 058401]

[6]Hong L,Dai F,Wang X 2016 IEEE Geosci.Remote Sens.Lett.13 1349

[7]Zhu H,Zhang S N,Zhao H C 2014 Acta Phys.Sin.63 058401(in Chinese)[朱航,張淑寧,趙惠昌2014物理學(xué)報(bào)63 058401]

[8]Simeunovic M,Popovic-Bugarin V,Djurovic I 2017 IEEE Trans.Aerosp.Electron.Syst.53 1273

[9]Tan R,Lim H S,Smits A B 2016 IEEE Region 10 Conference,TENCON,2016 p730

[10]Chen G F 2014 Ph.D.Dissertation(Xi’an:Xidian University)(in Chinese)[陳廣鋒 2014博士學(xué)位論文 (西安:西安電子科技大學(xué))]

[11]Zhao M M,Zhang Q,Luo Y 2017 IEEE Geosci.Remote Sens.Lett.14 174

[12]Yang Q,Deng B,Wang H 2014 EURASIP J.Wirel.Comm.1 61

[13]Huo K,You P,Jiang W D 2010 Jounal of Electronics&Information Technology 32 355(in Chinese)[霍凱,游鵬,姜衛(wèi)東2010電子與信息學(xué)報(bào)32 355]

[14]Deng D H,Zhang Q,Luo Y 2013 Acta Electronica Sinica 41 2339(in Chinese)[鄧冬虎,張群,羅迎 2013電子學(xué)報(bào)41 2339]

[15]Zhu H,Zhang S N,Zhao H C 2015 Digital Signal Process.40 224

[16]Sun Z G,Chen J,Cao X 2016 J.Syst.Engin.Electron.10 1973

[17]Zhang S N,Zhao H C,Xiong G,Guo C Y 2014 Acta Phys.Sin.63 158401(in Chinese)[張淑寧,趙惠昌,熊剛,郭長勇2014物理學(xué)報(bào) 63 158401]

[18]Sharafi nezhad S R,Alizadeh H,Eshghi M 2014 Elect.Eng.22nd Iranian Conference on IEEE Iran,2014 p1673

[19]Wang Y,Wu X,Li W 2016 Neurocomputing 171 48

[20]Yuan B,Chen Z,Xu S 2014 IEEE Trans.Geosci.Remote Sens.52 1285

[21]Setlur P,Fauzia A,Moeness A 2011 IET Signal Proc.5 194

[22]Ye Z F 2009 Statistical Signal Processing(Hefei:China University of Science and Technology Press)pp241–246(in Chinese)[葉中付 2009統(tǒng)計(jì)信號處理(合肥:中國科學(xué)技術(shù)大學(xué)出版社)第241–246頁]

[23]Guo L R,Hu Y H,Wang Y P 2016 Proceedings of the SPIE,Photonics Asia Beijng,China,October 11–14,2016 p21

[24]Hu Y,Guo L,Dong X 2016 Ubiquitous Positioning,Indoor Navigation and Location Based Services(UPINLBS)Fourth Int.Conf.IEEE Shanghai,China,November 2–4,2016 p264

[25]Hou Z F,Yang J,Zhang X 2011 Journal Wuhan University of Technology 1 142(in Chinese)[侯者非,楊杰,張雪2011武漢理工大學(xué)學(xué)報(bào) 1 142]

[26]Kay S M 2006 Fundamentals of Statistical Signal Processing:Estimation Theory(Prentice Hall PTR:Upper Saddle River)pp142–150

[27]Li j,Zhao Y J,Li D H 2014 Acta Phys.Sin.63 130701(in Chinese)[李晶,趙擁軍,李冬海 2014物理學(xué)報(bào) 63 130701]

[28]Guo L R,Hu Y H,Wang Y P 2017 Infrared and Laser Engineering 46 17(in Chinese)[郭力仁,胡以華,王云鵬2017紅外與激光工程 46 17]