窄帶雷達頻譜混疊微動信號的平動補償

王 超, 文樹梁, 葉春茂

(北京無線電測量研究所, 北京 100854)

?

窄帶雷達頻譜混疊微動信號的平動補償

王 超, 文樹梁, 葉春茂

(北京無線電測量研究所, 北京 100854)

針對存在頻譜混疊的窄帶雷達微動信號平動補償問題,從高階和低階頻譜混疊兩種情況出發,給出了一種平動補償方法。在對兩種頻譜混疊情況分析的基礎上,針對高階混疊將混疊的瞬時頻率近似為線性調頻鋸齒波信號,通過對頻率以及調頻率的估計,實現信號混疊程度的快速降階。對于低階混疊,首先將時頻圖在頻率維周期擴展,再基于形態學圖像處理方法提取出無模糊的瞬時頻率曲線,最后對其進行多項式的最小二乘擬合完成平動補償。仿真實驗驗證了所提方法的有效性以及對不同的信號形式和噪聲的魯棒性。

微動; 平動補償; 頻譜混疊

0 引 言

微動目標產生的微多普勒特征可用于對復雜目標的分類、識別和成像等領域[1-3]。對于微動目標的分析,人們通常將視線集中在特征提取方面,并假設目標沒有平動或平動已經補償[4-6]。在處理實際問題時,目標平動帶來的影響往往不能忽略,因此平動補償問題難以避免。國內外在該方面已有一部分成果,但仍顯不夠,下面進行簡要總結。

文獻[7-8]針對寬帶雷達的距離-慢時間像提出了基于經驗模式分解(empirical mode decomposition, EMD)的平動補償方法。具體做法是將距離-慢時間像中多個“傾斜的”正弦曲線分解到不同的本征模態函數(intrinsic mode function, IMF)中,通過去除分解余量實現了有效的平動補償。但對于窄帶雷達時頻(time-frequency, TF)曲線的平動補償,這類方法可能因頻譜混疊問題而失效。文獻[9-10]將微動回波信號視為附加了正弦調頻的多項式相位信號(sinusoidal frequency modulation-polynomial phase signal, SFM-PPS),利用高階模糊函數(high-order ambiguity function,HAF)實現了多項式相位項的補償。這類參數化方法對時域信號直接進行處理,不受頻譜混疊的影響且具有一定的抗噪性能,但僅適合處理單分量信號,在多分量時將產生交叉項干擾。另一類常用的平動補償方法是以頻譜最小熵或最小頻譜寬度為準則,對加速度、加加速度等平動參數直接進行搜索得到參數的估計值[11-12]。這類方法不考慮頻率的瞬時狀態,而將重點放在整體的聚集程度上,因此對不同的信號形式以及噪聲具有更強的魯棒性。但這類方法的補償精度通常不高,且計算復雜度與補償精度對搜索間隔與范圍很敏感,在沒有先驗知識的情況下合理設置這些參數將變得很困難。

本文針對存在頻譜混疊的窄帶雷達微動信號平動補償問題,首先依據混疊的程度分為高階和低階兩種情況。對于高階混疊,通過將混疊的瞬時頻率近似為線性調頻(linear frequency modulation, LFM)鋸齒波信號,提出了一種快速降階方法;對于低階混疊,在頻率維周期擴展的時頻圖上基于形態學圖像處理給出了一種無模糊瞬時頻率提取及平動補償方法。所提方法對目標的運動形式以及噪聲有較好的魯棒性,且參數設置簡單。

1 微動回波信號的兩種頻譜混疊情況

微動目標相對雷達的運動可分為平動和轉動,則不同時刻目標相對雷達的距離可表達為

r(tm)=rT(tm)+rR(tm)

(1)

式中,tm為慢時間;rT(tm)為距離的平動分量;rR(tm)為距離的轉動分量。

為簡化問題而又不失一般性,假設目標的平動由速度、加速度和加加速度組成,則距離的平動分量可表達為

(2)

式中,r0為目標與雷達的初始距離;v為平動速度;a為平動加速度;g為平動加加速度。

對于單散射點的旋轉運動,產生的距離分量可表示為

rR(tm)=Lksin(2πfotm+φk)

(3)

式中,Lk為散射點旋轉幅度;fo為旋轉頻率;φk為旋轉初相角。

對應的無模糊瞬時多普勒頻率可表達為

(4)

式中,λ為雷達波長;fT(tm)為瞬時平動頻率;fR(tm)為瞬時轉動頻率;Ak=4πf0Lk/λ為轉動頻率的幅度。

對于窄帶微動信號,假設沒有發生越距離單元徙動,則在目標所在的距離單元上利用時頻分析方法(未作特殊說明則默認采用短時傅里葉變換)可獲得混疊的二維時頻像X(tm,f)。定義時頻峰值曲線為

(5)

相應的峰值差分序列為

(6)

式中,fr為脈沖重復頻率(pulse repetition frequency,PRF)。

圖1表示出了兩種典型的頻譜混疊情況。

(1) 高階頻譜混疊:回波瞬時頻率變化劇烈且呈現明顯的周期性跳變,如圖1(a)～圖1(c)。

(2) 低階頻譜混疊:回波的瞬時頻率整體相對變化緩慢,觀察時間內只會跨越較少幾個脈沖重復頻率,如圖1(d)～圖1(f)。

圖1 頻譜混疊示意圖Fig.1 Spectrum aliasing scheme

從圖1(a)和圖1(b)中可發現,高階頻譜混疊的回波時頻峰值曲線基本由平動表征。由于以PRF為間隔被周期性的截斷,而在峰值差分序列中產生周期性的峰值,如圖1(c)。而低階頻譜混疊的回波時頻峰值曲線由平動和轉動共同表征,如圖1(d)和圖1(e)。此時峰值差分序列不再顯示出明顯的周期性,如圖1(f)。但相比于高階頻譜混疊,低階的混疊現象僅出現在相鄰的少數幾個PRF內。

2 頻譜混疊微動信號的平動補償

對于頻譜混疊的窄帶雷達微動信號,首先進行無模糊瞬時頻率的提取,再基于無模糊瞬時頻率趨勢項的最小二乘擬合實現平動補償。當微動回波信號出現高階頻譜混疊時,需先進行預處理將信號降為低階甚至無頻譜混疊的情況。

對于無越距離單元徙動的窄帶多分量微動信號,假設多個轉動分量的幅度Ak都相對PRF較小,且微動信號的不同分量在時頻面上是相交的。

2.1 高階頻譜混疊的快速降階方法

在高階頻譜混疊情況下,模糊的時頻峰值曲線主要由共同的平動加速度a和加加速度g表征,而多個轉動分量可視為附加在平動項上的“偽隨機噪聲”ε(tm)。此時,無模糊的時頻峰值曲線可表達為

(7)

受到信號回波PRF的限制,實際的時頻峰值曲線可近似為鋸齒波,具有如下形式：

(8)

式中,f0為鋸齒波信號fsw(tm)基頻分量的頻率。

由式(8)可知,鋸齒波的各次諧波分量能量依次遞減。這里我們僅考慮最強分量,即基頻分量。在時頻像中,鋸齒波基頻頻率f0相當于峰值曲線fpk(tm)變化一個PRF所需要的時間,即

(9)

式中,f′(·)為求導運算。

基于式(9),當加加速度為g=0時

(10)

因此,fsw(tm)的最強分量為一單頻信號,其頻率與目標平動加速度a相對應。

而當加加速度g≠0時,有

(11)

此時,基頻信號的頻率具有線性調頻的形式。我們熟悉的LFM信號是正弦包絡的,而這里fsw(tm)相當于鋸齒波包絡的LFM信號,在時頻域其最強分量對應一條具有一定斜率的斜線。

綜上所述,當加加速度g=0時,平動加速度的估計相當于頻域的譜峰估計問題,可直接利用FFT求峰值或其他更高精度的估計方法;當加加速度g≠0時,平動加速度和加加速度的估計相當于LFM信號初始頻率與調頻率的估計問題,可利用基于時頻面的Radon-Wigner變換或基于乘積型高階模糊函數(product high-order ambiguity function,PHAF)的參數化估計方法?；诩铀俣?、加加速度的估值進行平動的初次補償即可實現降階。由于平動速度僅會使瞬時頻率產生一個恒定的偏移量,這里采用一個簡單的圓移位實現補償[13]。

2.2 時頻面擴展

當信號為低階混疊時,為恢復瞬時頻率的連續性先將混疊的時頻圖進行正負方向的周期重復,得到擴展的時頻圖Xex(tm,f)。這里基于峰值差分序列來確定重復周期數

num+=Num{tm|Δfpk(tm)<-α·fr}

(12)

num-=Num{tm|Δfpk(tm)>α·fr}

(13)

式中,num+為正頻率上重復的周期數;num-為負頻率上重復的周期數;Num{.}為集合中的元素個數;α為接近但小于1的數,這里為了減小峰值頻率在多分量間跳變的影響而引入。實際情況中,還可對差分序列中的峰值最小間隔Δtmin或最小正負向擴展周期數nmin進行設定。

2.3 無模糊瞬時頻率提取

首先,將擴展的時頻圖進行二值化,可得

(14)

式中,Thr為二值化門限。

變換后的時頻圖Xex2(tm,f)可視為一幅二值化圖像,利用數字圖像處理的知識可實現噪聲環境下連續瞬時頻率曲線的提取。下面給出將用到幾種特殊的形態學圖像處理方法的表達式,具體說明參見文獻[14],這里不再贅述。

(1) 開操作

用結構元素B對集合A進行開操作可表達為

A°B=(A⊙B)⊕B

(15)

式中,A⊕B為結構元素B對集合A的膨脹操作;A⊙B為相應的腐蝕操作。

(2) 連通區域提取

令Y是包含于集合A中的連通分量,并假設Y中一點p是已知的,則連通分量可用下面的表達式迭代的得到:

Xk=(Xk-1⊕B)∩A,k=1,2,3…

(16)

式中,X0=p;B為一個適當大小的結構元素。當Xk=Xk-1時,迭代終止;Y=Xk即為提取出的連通分量。

當存在噪聲時,在二值化的擴展時頻圖中將出現個別離散的小區域,并且信號瞬時頻率曲線也將出現“毛刺”。開操作可使對象的輪廓變得光滑,消除細長的突出物。為了剔除“毛刺”并避免與鄰近的噪聲區域錯誤連接,這里先對圖像進行一次開操作,再由混疊時頻圖中幅度最強點出發進行連通區域提取來實現無模糊瞬時頻率的提取。其處理過程可表達為

Xda2=(Xk-1⊕B2)∩(Xex2°B1)

(17)

式中,Xda2為無混疊的二值化時頻圖;B1為開操作中用到的結構元素;B2為連通分量提取時用到的結構元素。

將二值化時頻圖視為一個模板,可將擴展時頻圖中一條完整而無模糊的時頻曲線提取出來

Xda=Xda2·Xex

(18)

2.4 基于無模糊瞬時頻率的平動補償

對于無模糊瞬時頻率的平動補償,可將Xda按行求重心作為多分量信號的平均瞬時頻率fme(tm),并對其采用基于最小二乘(least square, LS)的二次多項式擬合來實現補償。即

fme=D0·x+ε

(19)

式中

fme=(fme(t1),…,fme(tn))T

x=(v,a,g)T,ε=(ε(t1),…,ε(tn))T

需要指出的是,若得到的無混疊瞬時頻率可近似為單個正弦信號與多項式信號的和信號,采用文獻[15]中的方法能在提高補償精度的同時實現正弦信號周期的估計。

3 實驗仿真與分析

3.1 高階頻譜混疊的快速降階

假設目標圍繞中心點做二維轉動,以中心點為圓心,轉動平面為XOY面建立目標坐標系,目標中包含的3個分量在目標坐標系中的坐標分別為(0,1,0),(-1.2,0,0),(1.5,0,0)。平動加速度a=75.8 m/s2,加加速度g=5.45 m/s3,這里不考慮速度的影響,因此令平動速度v=0 m/s。同時,目標進行二維轉動,對應的轉動頻率為frot=0.5 Hz。雷達的仿真參數有載頻fc=3 GHz,脈沖重復頻率fr=800 Hz,雷達視線對應的單位方向矢量為nradian=(0.5,-0.866,0),仿真時間T=5 s,信噪比(signal to noise ratio, SNR)SNR=10 dB。

此時,信號的時頻圖、時頻峰值曲線和峰值差分序列如圖1(a)～圖1(c)所示。將時頻峰值曲線做Wigner-Ville變換,如圖2所示。觀察圖2可知,在時頻面中時頻峰值曲線可近似表征為一組平行的“斜線”,且強度隨頻率的增加依次遞減,與之前鋸齒波包絡LFM信號的時頻特征相符。

圖2 高階頻譜混疊時頻峰值曲線時頻圖Fig.2 TF distribution of TF peak curve with high-order spectrum aliasing

圖3給出了基于本文的快速降階方法得到的平動初步補償結果。從圖中可知,通過降階方法去除了信號由平動造成的頻譜混疊,且處理時間僅需要1.5 s。這里平動加速度、加加速度的估計方法采用PHAF。

圖3 高階頻譜混疊的快速降階Fig.3 Fast degrade of high-order spectrum aliasing

3.2 低階頻譜混疊的平動補償

令平動加速度a=13.4 m/s2,加加速度g=1.24 m/s3,則仿真時間內信號的頻譜混疊將由高階變為低階,信號的時頻圖、時頻峰值曲線和峰值差分序列如圖1(d)～圖1(f)所示。由峰值差分序列可得到負頻率擴展周期數為3,正頻率擴展周期數為1,其中α=0.8,Δtmin=T/20,nmin=1。擴展后的時頻圖如圖4所示。將擴展時頻圖二值化,并利用形態學圖像處理提取出的無模糊信號時頻曲線如圖5所示。其中,二值化門限取0.1倍時頻圖最大值,若時頻圖換算成了dB值則為最大值-10 dB。采用的結構元素為B1=ones(3,3),B2=ones(11,11),ones(n1,n2)為n1行、n2列值全為1的矩陣。無模糊的信號時頻曲線的重心線如圖6所示,它能對平動產生的趨勢項進行有效近似。

圖4 擴展時頻圖Fig.4 Expanded TF distribution

圖5 無模糊時頻曲線提取Fig.5 TF curve without frequency ambiguity

圖6 時頻曲線重心線Fig.6 Center line of TF curve

圖7 低階頻譜混疊的補償結果Fig.7 Compensate result of low-order spectrum aliasing

3.3 平動補償方法的魯棒性分析

本文提出的頻譜混疊下的平動補償方法對不同的信號形式以及噪聲具有較好的魯棒性。下面,從SNR和大轉動分量兩方面對該方法的魯棒性進行分析。

3.3.1 信噪比

在之前高/低階頻譜混疊的仿真場景下,將SNR從10 dB降低到3 dB和-3 dB,得到的平動參數估計值示于表1中。從表中數據可知,本方法對噪聲具有較好的魯棒性。

高階頻譜混疊時噪聲的影響等價于時頻峰值曲線的鋸齒波近似中噪聲項ε(tm)的增加。由于PHAF對噪聲有較好的魯棒性,因而在SNR較低時仍能實現有效的降階。對于低階頻譜混疊時,圖8顯示出了3 dB SNR下開操作對二值化圖像的去噪效果。可發現,開運算不僅能濾除小的噪聲斑,還能增加剩余噪聲彼此或與信號間的離散度,與相互連通的信號分量呈現兩種不同的狀態?；诖?可實現無模糊時頻曲線的重心線提取,如圖9所示。通過與圖6對比可發現,隨著SNR的降低提取出的重心線產生了“毛刺”。

表1 不同SNR下平動參數估計結果

圖8 開操作去噪Fig.8 Denoising with open operation

圖9 無模糊時頻曲線的重心線Fig.9 Center line of TF curve without frequency ambiguity

3.3.2 特殊信號形式—單個大轉動分量

轉動產生的頻率變化受到分量在目標坐標系中的位置、雷達波長以及轉速的影響。這里將單個信號分量的坐標設為(0,10,0),其他參數同第3.2節,則轉動頻率的幅度Ak將超過1倍PRF,如圖10所示。采用低價頻譜混疊平動補償方法,由重心近似的瞬時頻率曲線如圖11所示。

圖10 單個大轉動分量的時頻圖Fig.10 TF distribution of single component with large frequency variation

圖11 提取的重心線Fig.11 Extracted center line

圖12 最小二乘均方誤差與周期的關系曲線Fig.12 Relation curve between LS mean square error and rotate period

值得注意的是,大轉動分量情況下基于最小熵的平動補償將不再適用,其補償結果如圖13所示。其中,加速度搜索范圍為-100～100m/s2,均勻搜索800點,加加速度搜索范圍為-20～20m/s3,均勻搜索400點。

圖13 基于最小熵的大轉動分量平動補償結果Fig.13 Translational compensation result based on minimum entropy

4 結 論

本文針對窄帶雷達頻譜混疊微動信號的平動補償問題,從高階和低階頻譜混疊兩種情況出發,給出了一種無模糊瞬時頻率提取及平動補償方法。提出的方法對噪聲以及一些特殊的信號形式具有較好的魯棒性?；谔崛〕龅臒o模糊瞬時頻率,本方法也可以與已有的一些針對無模糊情況的平動補償方法聯合使用,如之前提到的基于EMD的平動補償方法。

值得指出的是,為了簡化問題,本文采用了點目標的仿真模型。它能較好的表征彈道類等具有簡單散射特性的目標。而針對電磁散射特性更為復雜的目標,由于時頻面重心線通常能夠表征平動引起的整體趨勢變化,因此保證了算法的有效性。

[1]ChenVC,LiF,HoSS,etal.Micro-Dopplereffectinradar:phenomenon,modelandsimulationstudy[J]. IEEE Trans.on Aerospace and Electrical System, 2006, 42(1): 2-21.

[2]LuoY,ZhangQ,QiuCW,etal.Micro-DopplereffectanalysisandfeatureextractioninISARimagingwithstepped-frequencychirpsignals[J]. IEEE Trans.on Geoscience and Remote Sen-sing, 2010, 48(4): 2087-2098.

[3]MaL,LiuJ,WangT,etal.Micro-Dopplercharacteristicsofsliding-typescatteringcenteronrotationallysymmetrictarget[J]. Science in China Scientia Sinica Information Sciences,2011,41(5):605-616.

[4]SunHX,LeiZ.Nutationfeatureextractionofballisticmissilewarhead[J]. Electronics Letters, 2011, 47(13):770-772.

[5]LeiP,WangJ,SunJ.Analysisofradarmicro-Dopplersignaturesfromrigidtargetsinspacebasedoninertialparameters[J]. IET Radar Sonar Navigation, 2011, 5(3): 93-102.

[6]StankvicL,DjurovicI,ThayaparanT.Separationoftargetrigidbodyandmicro-DopplereffectsinISARimaging[J]. IEEE Trans.on Aerospace and Electrical System, 2006, 42(4):1496-1506.

[7]LuoY,BaiYQ,ZhangQ,etal.Translationalmotioncompensationandmicro-Dopplerfeatureextractionofballistictargets[J]. Journal of Electronics & Information Technology, 2012, 34(3):602-608.

[8]ZhaoYQ,WangGZ,BaiYQ,etal.AnewalgorithmoftranslationalcompensationforspacespinningtargetbasedonEMD[C]∥Proc.of IEEE International Conference on Signal Processing, Communications and Computing, 2013: 1-5.

[9]YangYC,TongNN,FengCQ,etal.Translationcompensationandmicro-Dopplerextractionoftheechofromballistictargetsinmidcourse[J]. Scientia Sinica(Informationis), 2013, 43(9):1172-1182.

[10]LiuJ.Radarsignalparameterestimationandphysicalfeatureextractionofmicro-motiontargets[D].Changsha:NationalUniversityofDefenseTechnology, 2010.(劉進.微動目標雷達信號參數估計與物理特征提取[D].長沙:國防科學與技術大學,2010.)

[11]YuWY,WangJF.Phaseadjustmentforextractionofmicro-motioninformationofballistictargets[C]∥Proc.of the IEEE International Congress on Image and Signal Processing, 2012: 1837-1840.

[12]GaoHW,XieLG,WenSL,etal.Researchontheinfluenceofaccelerationonmicro-Doppleranditscompensation[J]. Journal of Astronautics, 2009, 30(2): 705-711. (高紅衛, 謝良貴,文樹梁,等.加速度對微多普勒的影響及其補償研究[J].宇航學報,2009,30(2):705-711.)

[13]GaoHW,XieLG,WenSL,etal.Researchontheinfluenceofvelocityonmicro-Doppleranditscompensation[J]. Aerospace Electronic Warfare, 2008, 24(4):46-50. (高紅衛,謝良貴,文樹梁,等.速度對微多普勒的影響及其補償研究[J].航天電子對抗, 2008, 24(4):46-50.)

[14]ConzalezRC,WoodsRE. Digital image processing[M].YuanQT,trans. 2nded.Beijing:PublishingHouseofElectronicsIndustry,2006:423-437.(岡薩雷斯伍茲.數字圖像處理[M].阮秋琦,譯.2版. 北京:電子工業出版社,2006:423-437.)

[15]WuYH,YangJF,ChengJY,etal.Highaccuracyestimationalgorithmfor4parametersofsinewave[J].Journal of University of Science and Technology of China, 2006, 36(6):625-629.(吳義華,楊俊峰,程敬原,等.正弦信號四參數的高精度估計算法[J].中國科學技術大學學報,2006,36(6):625-629.)

Translational compensation of micro-motion with spectrum aliasing for narrow-band radar

WANG Chao, WEN Shu-liang, YE Chun-mao

(BeijingInstituteofRadioMeasurement,Beijing100854,China)

A translational compensation algorithm of micro-motion with spectrum aliasing for narrow-band radar is proposed. It is done in high-order and low-order spectrum aliasing cases. Based on the analysis of these two aliasing cases, the translational can be fast degraded by approximating the aliasing frequency as a linear frequency modulation saw-tooth wave in high-order spectrum aliasing cases. To deal with the low-order spectrum aliasing, the time-frequency image is expanded in frequency, and then the time-frequency curve without frequency ambiguity can be extracted based on morphological image processing. Finally, the least-square fitting of polynomial is utilized in the curve for translational compensation. Simulations are given to validate the effectiveness and robustness of the proposed algorithm.

micro-motion; translational compensation; spectrum aliasing

2015-12-24;

2016-08-08;網絡優先出版日期:2016-08-22。

國家自然科學基金(61271417)資助課題

TN 95

A

10.3969/j.issn.1001-506X.2016.12.08

王 超(1990-),男,博士研究生,主要研究方向為雷達目標識別技術。

E-mail:whhitwa@126.com

文樹梁(1971-),男,研究員,博士研究生導師,主要研究方向為雷達系統設計與信號處理。

E-mail:wenshul@sina.com

葉春茂(1981-),男,高級工程師,碩士研究生導師,主要研究方向為雷達ISAR成像與目標識別技術。

E-mail:danielgodman@163.com

網絡優先出版地址:http:∥www.cnki.net/kcms/detail/11.2422.TN.20160822.1005.004.html