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

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

1 引 言

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

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

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

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

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

根據電磁散射理論,電大尺寸目標在光學區的散射具有局部效應,目標總的電磁散射可等效為一些局部散射點的疊加,這些局部的散射源就叫作等效散射中心.對于激光探測的目標微多普勒效應回波信號,其波長遠小于目標尺寸,采用電磁散射在光學區的等效散射中心理論是合理的.

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

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

3 時頻交疊激光微多普勒信號的參數分離估計

3.1 基于SVR譜的微動頻率估計

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

如果直接用完整的信號構建矩陣,矩陣列數與信號周期之間的誤差會隨著行數增加不斷累加,這常常會導致SVR方法失效.所以,這里提出固定列數的周期掃描的改進方法來實現周期的估計.構建矩陣可寫為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

因為P=P2=PH,所以(12)式可化簡為

一般直接用(14)式的似然函數形式構建參數ψ的概率分布函數,再用MCMC方法產生服從分布的樣本實現對參數的最大似然估計,這對于傳統微波探測或是平穩信號處理是可行的.但是對于激光微多普勒信號,由于探測波長極短,信號由雷達探測中的弱調制變為深度調制,似然函數的非線性程度急劇增加,概率分布變為復雜的多峰形狀.這將導致馬爾可夫鏈在產生平穩分布時會錯誤收斂至局部最大值,最終估計錯誤.微多普勒信號參數估計概率分布與探測波長之間的關系如圖2.

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

圖2 不同波長下傳統似然函數分布 (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,并以此作為先驗信息計算Dv1和ρ01參數域上的似然函數分布.當估計值?Dv1和?ρ01與信號分量1實際的參數值接近時,sl1(n)的頻譜能量集中在零頻附近很窄的帶寬里,能量峰值極高;而此時,由于分量2與分量1的微動參數不同,在參數域中sl2(n)的頻譜能量只能分布在一個較寬的帶寬內;slnoi(n)則均勻分布在整個頻譜中.當估計值完全與實際值相同時sl1(n)為常數,頻譜為沖擊響應形狀.所以,可以把似然函數頻譜中占總能量的比例為η時的能帶寬度作為評價參數估計值是否接近真值的標準,構建新的似然函數,

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

A1/A2的值可根據圖1中的關系進行確定.新似然函數各分量頻譜關系和參數示意如圖3所示.

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

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

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

3.2.2 MCMC算法實現

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

當候選狀態使目標分布概率密度增大時,以概率1更新為這個狀態,否則,以概率π(x?i)/π(xti)進行更新.所以,t+1時刻狀態的實際轉移概率為

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

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

2)t時刻,根據條件建議分布的轉移概率得到參數 ρ0k的候選狀態 ρ?0k;

3)計算接受概率

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

5)根據條件建議分布的轉移概率得到參數Dvk的候選狀態D?vk;

6)計算接受概率

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

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

和

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

3.3 信號幅度和初相的估計

根據3.1節和3.2節估計的第k分量的微動參數重構信號

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

式中

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

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

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

3.4 激光M D信號參數估計的克拉美羅界

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

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

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

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

4 仿真分析和實驗驗證

4.1 仿真結果與分析

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

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

從總的回波信號中減去分量1,按照相同的步驟對剩余分量進行參數估計.剩余分量的SVR如圖6(a)所示.圖中峰值在第391點處,可求出對應分量2的微動頻率為300 Hz.圖6(b)為分量2對應微動參數的概率密度分布,剩余信號的概率分布仍然是平滑單峰形狀,可以保證分量2參數的正確收斂.圖6(c)和圖6(d)為分量2微動參數的MCMC估計對比結果,整體結果與圖5一致,利用傳統的似然函數得到的馬氏鏈被困在初始值附近,而改進的似然函數確保了正確的收斂.此外,與圖5對比還可發現分量2估計的馬爾可夫鏈收斂到真值的速度較分量1更快.這是因為分量1參數被精確地估計,從原始信號中移除重構的分量1后,分量2成為剩余信號的主要成分,沒有其他信號分量的影響,所以收斂更快.

圖5 分量1的微動參數估計結果 (a)奇異值比譜;(b)參數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的微動參數估計結果 (a)奇異值比譜;(b)參數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ξ.

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

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

從圖7中可以看出,由于兩個分量中參數值的差異導致各信噪比下的CRB也存在差異.但隨著信噪比的增加,分量1和分量2的微動參數估計MSE都逐漸接近各自的克拉美羅界,說明本文基于MCMC和重新定義的似然函數來實現微動參數的最大似然估計具有最優的估計性能.此外,仔細對比可發現,各信噪比下分量2中參數的估計誤差都略大于分量1的,這是因為對分量2的估計依賴于分量1的結果,會存在一定誤差傳遞的干擾.

表1 參數估計結果Table 1.Parameter estimation result.

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

4.2 實驗驗證與對比

利用實驗數據估計目標微動參數來驗證算法有效性,實驗設置如圖8.

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

圖8 實驗結構 (a)系統框圖;(b)實物圖Fig.8.Experiment system:(a)Structure diagram;(b)experiment setup.

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

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

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

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

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

圖11 實驗數據分量1的微動參數估計結果 (a)奇異值比譜;(b)參數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 實驗數據分量2的微動參數估計結果 (a)奇異值比譜;(b)參數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 重構信號波形相似度Table 2.Waveform similarity comparison.

從表2可以看出,利用最大似然估計得到的微動參數重構波形,其相似度達到0.9以上,遠高于基于傳統時頻分布的逆Radon變換重構波形的相似度,驗證了參數估計的準確性.說明本文所提方法更有利于實現對微動參數的精確估計,這為基于微動特征的目標分類和精細識別奠定了基礎.

5 結 論

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

[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博士學位論文 (西安:西安電子科技大學)]

[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)[王童,童創明,李西敏2015物理學報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物理學報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博士學位論文 (西安:西安電子科技大學)]

[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)[霍凱,游鵬,姜衛東2010電子與信息學報32 355]

[14]Deng D H,Zhang Q,Luo Y 2013 Acta Electronica Sinica 41 2339(in Chinese)[鄧冬虎,張群,羅迎 2013電子學報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物理學報 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統計信號處理(合肥:中國科學技術大學出版社)第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武漢理工大學學報 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物理學報 63 130701]

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