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

1 引 言

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

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

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

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

2 單通道多分量激光微多普勒信號(hào)模型

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

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

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

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

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

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

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

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

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

對(duì)于包含多個(gè)微多普勒分量的回波信號(hào),若各分量幅度相似,則一次SVR可估計(jì)出各分量的微動(dòng)頻率;若各分量信號(hào)幅度相差較大,則一次SVR只得到主周期分量的微動(dòng)頻率,繼續(xù)估計(jì)該分量剩余其他參數(shù)后,從回波信號(hào)中將其重構(gòu)信號(hào)去除,用同樣方法可估計(jì)其他分量的微動(dòng)頻率.通過仿真得到各分量信號(hào)幅度比與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分別表示奇異值比譜中周期較長(zhǎng)的分量和周期較短的分量對(duì)應(yīng)的峰值大小(即σ1/σ2的值),橫坐標(biāo)代表SVR1和SVR2對(duì)應(yīng)分量的幅值比.當(dāng)信號(hào)幅值比超出圖中范圍時(shí),SVR中只有一個(gè)明顯的峰值.圖1的關(guān)系可為下一節(jié)中設(shè)計(jì)適合激光微多普勒信號(hào)的似然函數(shù)提供參數(shù)選擇的依據(jù).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2)t時(shí)刻,根據(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個(gè)樣本矢量burn-in,對(duì)剩下的樣本求均值得到參數(shù)的估計(jì)為

和

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

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

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

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

式中

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

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

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

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

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

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

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

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

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

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

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

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

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

圖5 分量1的微動(dòng)參數(shù)估計(jì)結(jié)果 (a)奇異值比譜;(b)參數(shù)D v1和ρ01的概率密度分布;(c)振動(dòng)幅度;(d)振動(dòng)初始相位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的微動(dòng)參數(shù)估計(jì)結(jié)果 (a)奇異值比譜;(b)參數(shù)D v1和ρ01的概率密度分布;(c)振動(dòng)幅度;(d)振動(dòng)初始相位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é)果和相對(duì)誤差如表1所列.表中RE=(?ψ?ψ)×100%/ψ,表示相對(duì)估計(jì)誤差;γ=|E[sk(n)?sk(n)]|·|E[sk(n)]2·E[?sk(n)]2|?1/2表示波形相似度,通過對(duì)比重構(gòu)信號(hào)?sk(n)與實(shí)際信號(hào)在波形上的相似程度來(lái)反映參數(shù)估計(jì)的準(zhǔn)確度,γ越接近1則波形相似度越高,表明參數(shù)估計(jì)越接近真實(shí)值.從表1可以看出,本文方法對(duì)目標(biāo)微動(dòng)參數(shù)估計(jì)的相對(duì)誤差都在10?5量級(jí),具有較高的精度.

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

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

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

圖7 微動(dòng)參數(shù)估計(jì)均方誤差與克拉美羅界對(duì)比 (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)證與對(duì)比

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

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

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

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

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

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

5 結(jié) 論

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

[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]