基于時間-調(diào)頻斜率分布的多普勒調(diào)頻率估計

汪金真 宿紹瑩 陳曾平

1 引言

在ISAR成像中，尤其在高分辨成像中，正確估計出目標(biāo)回波線性調(diào)頻（Linear Frequency Modulated, LFM）信號的多普勒調(diào)頻斜率后對成像結(jié)果進行補償，可以有效提高成像質(zhì)量。最大似然法[1]基于2維搜索具有最優(yōu)的性能，但是運算量大。為了克服這個問題。文獻[2]提出基于相位域的調(diào)頻斜率估計算法，簡化了運算。文獻[3]采用相位線性擬合的方法大大提高了計算效率。文獻[4]提出了短時傅里葉變換調(diào)頻直線線性擬合的方法，較精確地估計出LFM信號的調(diào)頻斜率；但是它們都只能用于單分量LFM信號的場合。在多分量LFM信號應(yīng)用場合，Wigner-Ville分布[5]對各分量LFM信號具有最佳的時頻聚集性，但存在交叉項的影響，通過Hough變換[6,7]可以一定程度上抑制交叉項。分數(shù)階傅里葉變換（Fractional Fourier Transform, FrFT）[8]和離散chirp傅里葉變換[9]都是比較有效的方法，但它們都需要 2維搜索。離散多項式變換[10,11]計算量小，但要求較高的信噪比。文獻[12,13]提出時間-調(diào)頻斜率分布（Time Chirp Distribution, TCD）的定義并作了改進，文獻[14,15]將TCD應(yīng)用在多分量線性調(diào)頻信號的檢測與參數(shù)估計當(dāng)中并提出了多分量線性調(diào)頻信號交叉項的抑制方法。文獻[16]改進了三次相位函數(shù)計算方法，較快實現(xiàn)了多分量LFM信號參數(shù)估計，但它們都是直接針對線性調(diào)頻信號進行研究的，對目標(biāo)回波線性調(diào)頻信號則研究甚少。

本文建立了單散射點和多散射點目標(biāo)回波模型并分析其瞬時相位特點，提出了基于 TCD的多普勒調(diào)頻斜率估計算法。針對多散射點目標(biāo)回波的TCD的自項和交叉項均在多普勒調(diào)頻斜率處取得極值的特點，采用垂直調(diào)頻斜率軸投影積分有效抑制非多普勒調(diào)頻斜率交叉項和噪聲，增強多普勒調(diào)頻斜率項，通過二次搜索投影積分的最大值位置估計多普勒調(diào)頻斜率，減少了計算耗時。仿真表明，該方法具有很好的單散射點和多散射點目標(biāo)回波多普勒調(diào)頻斜率估計性能。

2 目標(biāo)回波模型

在高斯白噪聲環(huán)境下目標(biāo)回波模型可表示為

式中，t是時間；P是散射點數(shù)； { sm（t） }1是與目標(biāo) P個散射點相對應(yīng)的回波線性調(diào)頻信號；w（t）是與信號無關(guān)的均值為0、方差為 δ2的高斯白噪聲。

距離雷達為Rm的散射點回波LFM信號經(jīng)距離壓縮后的信號為

式中，A為幅度；fi為起始頻率；Tp為脈寬；k為調(diào)頻率；c為電磁波傳播速度；散射系數(shù)為σm; λ為波長。對于單散射點回波，記瞬時相位為 θm（t），

從式（3）可以看出，瞬時相位 θm（t）是時間t的二次函數(shù)，對 θm（t）進行多項式擬合可以估計出多普勒調(diào)頻斜率k。對于多散射點目標(biāo)回波 ∑ sm（t ），雖然各分量相位的二次項系數(shù)均為 k的π倍，但回波相位是各分量的合成相位，不等于各分量相位的線性疊加，因此，必須采用新方法估計多普勒調(diào)頻斜率k。

3 時間-調(diào)頻斜率（TCD）分布

3.1 TCD分布及自項分析

類似WVD定義，信號x（t）的TCD分布定義[12]為

考慮距離雷達 R、散射系數(shù)為σ的單散射點的回波 LFM 信號 x（t），將 x（t）代入式（4），可得

從式（5）可以看出，對于任意時刻t，當(dāng)u=k, T C （ t, u）均具有極大值，并且以極值點為中心沿著調(diào)頻斜率軸急劇下降，而加性高斯白噪聲不具有這種特性。

圖 1（a）是單散射點目標(biāo)模型，圖 1（b）是該單散射點目標(biāo)回波LFM信號的TCD，采樣點數(shù)2501，采樣頻率為 10 MHz，調(diào)頻斜率為 - 1.6 × 10-4MHz/μs，信噪比為1 dB。從圖1（b）可以看出，單散射點回波LFM信號的TCD在時間-調(diào)頻斜率平面投影為一條垂直于調(diào)頻斜率軸的直線，任意時刻的極大值位置均落在該直線上。

3.2 交叉項分析

考慮多散射點目標(biāo)上任意兩個散射點的回波LFM 信號 sm1（t）, sm2（t）, sm1（t）與 sm2（t）之間的 TCD的交叉項為

將式（2）代入式（6），計算可得

圖1 單散射點目標(biāo)模型及回波LFM信號的TCD

由式（7）可以看出，任意時刻 t，當(dāng) u=k 時，sm1（t）與sm2（t）之間的TCD的交叉項絕對值TCm1,m2（t, u）將取得極大值。式（7）取絕對值可得

記多散射點目標(biāo)回波LFM信號的TCD交叉項為 T Ccross（t, u），有

將式（7）代入式（9）中，可得

由式（10）可知，在任意時刻 t，多散射點目標(biāo)回波LFM信號的 TCD交叉項TCcross（t, u）也將在u=k處取得極大值。

3.3 多散射點目標(biāo)回波LFM信號的TCD

綜合3.1節(jié)和3.2節(jié)的分析，可以得到多散射點目標(biāo)回波 LFM 信號 x（t）的 TCD，記為 TCx（t,u） ，如式（11）所示。式中，噪聲w（t）的TCD 自項以及噪聲 w（t）與各散射點回波 LFM 信號sm（t）的交叉項一起記為（,）t uξ。從式（11）可以看出，對于無噪的多散射點目標(biāo)回波LFM信號，其TCD分布將沒有與噪聲有關(guān)的最后一項，任意時刻t, T Cx（t, u）將在u = k 處取得峰值；對于含噪情況，噪聲項 ξ（t, u）出現(xiàn)，會使得x（t）的TCD噪底增強，峰值相對不再尖銳，此外，還有可能出現(xiàn)一些虛假偽峰。

圖 2（a）是多散射點飛機目標(biāo)模型，圖 2（b）是該多散射點目標(biāo)回波的TCD，采樣頻率為10 MHz，采樣點數(shù) 2501，調(diào)頻斜率為 - 1 .6× 1 0-4MHz/μs，信噪比為5 dB。從圖2（b）可以看出，在多散射點目標(biāo)模型下，交叉項TCcross（t, u）在LFM信號調(diào)頻斜率值處取得極大值，該極大值位置與自項極大值位置大致重疊在一起，TCD在時間-調(diào)頻斜率平面投影仍然為一條近似垂直于調(diào)頻斜率軸的直線，受噪聲和交叉項的影響，該直線不是很緊湊，極值點不是嚴格地分布在該直線上，而是分布在該直線具有一定鄰域的條帶內(nèi)，該直線上各點的分布強度顯然不如圖 1（b）中單散射點情況的分布強度均勻，但從統(tǒng)計學(xué)上，該直線條帶的中心線為 k =-1.6 × 10-4MHz/μs，任意時刻的極值點對應(yīng)的調(diào)頻斜率都將趨近 1 .6×10-4MHz/μs。

4 多普勒調(diào)頻斜率估計

4.1 參數(shù)估計

記 x（t）的 TCD 分布為離散形式為 TCx（n,u），定義垂直調(diào)頻斜率軸投影積分累積量（）uΛ

由于多散射點目標(biāo)回波LFM信號的TCD的交叉項和自項均在u=K處取得極值，且以極值點為中心沿著調(diào)頻斜率軸急劇下降，通過式（12）可以很好累積多普勒調(diào)頻斜率處的自項和交叉項，即增強多普勒調(diào)頻斜率項，同時抑制非多普勒調(diào)頻斜率處的交叉項和噪聲。 Λ（u）相當(dāng)于TCx（n, u ） 對整個時間軸進行積分，可看作 T Cx（n, u ） 在零度角上的離散Radon變換[17]，由于只在一個角度上進行Radon變換，因而運算量小。

圖2 多散射點目標(biāo)模型及回波LFM信號的TCD

圖3 給出了單散射點和多散射點目標(biāo)回波TCD在垂直調(diào)頻斜率軸上進行投影積分的情況，單散射點和多散射點目標(biāo)回波TCD分別如圖1（b），圖2（b）所示。從圖3中可以看出，經(jīng)過垂直調(diào)頻斜率軸投影積分累積后，多普勒調(diào)頻斜率項得到了有效增強，非多普勒調(diào)頻斜率項得到了有效抑制，因此，可通過搜索（）uΛ最大值位置估計出多普勒調(diào)頻斜率。為了較少計算耗時，本文采用一次粗搜索方法獲取積分累積量（）uΛ最大值粗略位置基礎(chǔ)上進行二次精搜索的方法來精確估計多普勒調(diào)頻斜率，具體方法為：

（1）根據(jù)式（11）計算 x（t）的 TCD 分布 TCx（t, u），其離散形式記為 TCx（n,u）；

（2）根據(jù)式（12）計算 TCx（n,u）的垂直調(diào)頻斜率軸投影積分，得到（）uΛ；

（3）設(shè)置 u的初級搜索步長 Nstep1，在u的離散取值范圍 - 1 / N + Nstep1, 1/N - Nstep1內(nèi)搜索 Λ（u）的最大值以及相鄰的次最大值、次次最大值的位置M1, M2, M3；

（4）取M2, M3中較小值賦給 Mmin，較大值賦給Mmax，設(shè)置 u的精細步長 Nstep2, u =Mmin, Mmin+Nstep2, Mmin+ 2 Nstep2, … ,Mmax，重新計算 TCx（n,u），重復(fù)步驟（2），在u的離散取值范圍[Mmin, Mmax]內(nèi)搜索Λ（u）的最大值的精細位置M︿；

（5）根據(jù)數(shù)字多普勒調(diào)頻斜率公式 K = k Ts2，得到多散射點目標(biāo)回波 LFM 信號多普勒調(diào)頻斜率精估計值?k 為

4.2 計算量分析

表1，表2分別分析了N為奇數(shù)和偶數(shù)情況下計算TCD的所需的復(fù)數(shù)乘法和復(fù)數(shù)加法量，Λ（u）計算所需的復(fù)數(shù)加法量以及 Λ （ u ）最大值搜索總次數(shù)，最大值搜索次數(shù)包括第1次搜索次數(shù)和第2次搜索次數(shù)，均與Nstep1, Nstep2有關(guān)。

表1 TCD計算量分析（N為奇數(shù)）

表2 DCD計算量分析（N為偶數(shù)）

其中復(fù)數(shù)乘法最耗時，記復(fù)數(shù)乘法次數(shù)為N*，在N一定的情況下，N為奇數(shù)時，

圖3 單散射點和多散射點目標(biāo)回波TCD在垂直調(diào)頻斜率軸上的投影積分

N為偶數(shù)時，

當(dāng) 且 僅 當(dāng) 1 /（N Nstep1） = （ Mmax- Mmin） /Nstep2， 即Nstep2=NNstep1（ Mmax- Mmin） 時，式（21）和式（22）取得最小值，此時 N*min= 2 （N + 1 ）（N - 1）/（N Nstep1）或 N*min= 2 （N + 2 ）/Nstep1。取 Nstep1=1/（100N ） ，則Nstep2= （ Mmax- Mmin）/100, N*min= 2 00（N + 1）?（N-1）或N*min= 2 00（N + 2 ）N。

5 實驗結(jié)果與分析

為了驗證本文算法的有效性，對不同模型的數(shù)據(jù)進行仿真實驗。單散射點目標(biāo)模型和多散射點目標(biāo)模型分別如圖 1（a），圖 2（a）所示，回波 LFM 信號采樣點數(shù) N=2501，采樣頻率 fs=10 MHz，調(diào)頻斜率均為 k =-1 .6× 1 0-4MHz/μs 。

圖4給出了單散射點情況下的多普勒調(diào)頻斜率估計性能曲線，并與二次相位擬合法（Quadratic Phase Fitting, QPF）、平滑二次相位擬合法（Smoothed Quadratic Phase Fitting, SQPF）, FrFT法進行了比較。實驗中TCD法中 Nstep1=1/（100N ）,Nstep2= （ Mmax- Mmin）/100，平滑二次相位擬合法的平滑窗選用點數(shù)為11的中值濾波器。從圖中可以看出，QPF在較高的信噪比下具有較好的性能，當(dāng)SNR≥11 dB時，QPF的性能優(yōu)于FrFT；但當(dāng)SNR減小時，二次相位關(guān)系變差，QPF的性能急劇惡化，明顯次于FrFT。SQPF在一定信噪比范圍內(nèi)能夠改善二次相位關(guān)系，SNR≥9 dB時，SQPF的性能優(yōu)于FrFT，當(dāng)SNR降低時，對相位的平滑也無濟于改善二次相位關(guān)系，SQPF失效。FrFT具有較穩(wěn)定的性能，RMSE與SNR近似成線性 dB關(guān)系。當(dāng)SNR ≥ -8 dB 時，本文 TCD法的性能明顯優(yōu)于QPF, SQPF, FrFT，而且當(dāng)SNR ≥ -5 dB ，均方根誤差保持在一個很低的穩(wěn)定水平。

圖4 單散射點目標(biāo)回波LFM信號多普勒調(diào)頻斜率估計性能

圖5 給出了多散射點情況下的多普勒調(diào)頻斜率估計性能曲線，并與 FrFT法進行了比較。從圖中可以看出，SNR在-7～0 dB范圍內(nèi)時，本文TCD法的性能優(yōu)于FrFT法的性能；當(dāng)SNR≥1 dB時，本文TCD法的性能趨于穩(wěn)定。在-7～15 dB范圍內(nèi)，本文算法較 FrFT法具有更好的抗噪性能，這主要是因為多散射點目標(biāo)回波LFM信號的TCD分布在調(diào)頻斜率處進行時間上的累積可以有效抑制噪聲和非多普勒調(diào)頻斜率處的交叉項，而多散射點目標(biāo)回波LFM信號的FrFT峰值受噪聲影響很大。

6 結(jié)束語

本文針對單散射點和多散射點目標(biāo)回波相位自身的特點，提出了基于 TCD的多普勒調(diào)頻斜率估計算法，采用垂直調(diào)頻斜率軸投影積分可以有效抑制 TCD的非多普勒調(diào)頻斜率交叉項及噪聲、增強多普勒調(diào)頻斜率項，計算量分析表明，通過改變數(shù)字調(diào)頻斜率的步長，可以調(diào)節(jié)總計算量，使得復(fù)數(shù)乘法次數(shù)在一定條件下達到最小。仿真實驗驗證了該算法的有效性，單散射點情況下，SNR 3 dB≥-時，該算法的性能優(yōu)于FrFT, QPF, SQPF；多散射點情況下，該算法的整體抗噪性能優(yōu)于FrFT。雖然TCD計算量較大，然而，隨著計算機的發(fā)展，其處理能力提升，該問題非常容易得到解決，因此，本文算法具有較好的工程指導(dǎo)意義。

圖5 多散射點目標(biāo)回波LFM信號多普勒調(diào)頻斜率估計性能

[1] Tao R, Feng Y, and Wang Y. Theoretical and experimental study of a signal feature extraction algorithm for measuring propeller acceleration in a port surveillance system[J]. IET Radar, Sonar & Navigation, 2011, 5（2）: 172-181.

[2] 金勝, 王峰, 鄧振淼, 等. 一種 LFM 信號相位域快速高精度參數(shù)估計算法[J]. 系統(tǒng)工程與電子技術(shù), 2011, 33（2）: 264-267.Jin Sheng, Wang Feng, Deng Zhen-miao, et al.. Fast and accurate estimator on parameters of chirp signals in phase domain[J]. Systems Engineering and Electronics, 2011, 33（2）:264-267.

[3] 牛萌, 李釗, 朱曉光, 等. LFM信號參數(shù)估計的最大似然改進算法[J]. 無線電通信技術(shù), 2009, 35（3）: 59-61.Niu Meng, Li Zhao, Zhu Xiao-guang, et al.. Improved maximum likelihood estimation algorithm to estimate parameters of LFM signal[J]. Radio Communications Technology, 2009, 35（3）: 59-61.

[4] Wang J Z, Su S Y, and Chen Z P. Parameter estimation of chirp signal under low SNR[J]. Science China Information Sciences, 2015, 58: 020307（13）.

[5] Chen G H, Ma S W, Qin T H, et al.. Geometryinformation-aided efficient motion parameter estimation for moving-target imaging and location[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12（1）: 155-159.

[6] Zeng D, Zeng X, and Tang B. Automatic modulation classification of radar signals using the Rihaczek distribution and Hough transform[J]. IET Radar, Sonar & Navigation,2012, 6（5）: 322-331.

[7] Jia S Y, Wang G H, and Tan S C. Radial acceleration estimation of multiple high maneuvering targets[J]. Journal of Systems Engineering and Electronics, 2014, 25（2）: 183-193.[8] Miah K and Potter D. Geophysical signal parameterization and filtering using the fractional Fourier transform[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2014, 7（3）: 845-852.

[9] Wu L, Wei X Z, Yang D G, et al.. ISAR Imaging of targets with complex motion based on discrete chirp Fourier transform for cubic chirps[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50（10）: 4201-4212.

[10] 高春霞, 張?zhí)祢U, 金翔, 等. 基于離散多項式變換的寬帶線性調(diào)頻信號波達方向估計[J]. 計算機應(yīng)用, 2011, 31（10）:2872-2875.Gao Chun-xia, Zhang Tian-qi, Jin Xiang, et al.. Arrival direction estimation of wideband linear frequency modulation signal based on DPT[J]. Journal of Computer Application,2011, 31（10）: 2872-2875.

[11] 龐存鎖. 基于離散多項式相位變換和分數(shù)階傅里葉變換的加速目標(biāo)檢測算法[J]. 電子學(xué)報, 2012, 40（1）: 184-188.Pang Cun-suo. An accelerating target detection algorithm based on DPT and fractional Fourier transform[J]. Acta Electronica Sinica, 2012, 40（1）: 184-188.

[12] Peter O. A fast algorithm for estimating the parameters of a quadratic FM signal[J]. IEEE Transactions on Signal Processing, 2004, 52（2）: 385-393.

[13] Maree F and Peter O. Extending the performance of the cubic phase function algorithm[J]. IEEE Transactions on Signal Processing, 2007, 55（10）: 4767-4774.

[14] 杜雨洺, 楊建宇. 基于時間-調(diào)頻斜率分布的多線性調(diào)頻信號檢測與參數(shù)估計[J]. 電子與信息學(xué)報, 2007, 29（3）: 631-634.Du Yu-ming and Yang Jian-yu. Superimposed LFM signals detection and parameter estimation based on time-frequency rate distribution[J]. Journal of Electronics & Information Technology, 2007, 29（3）: 631-634.

[15] 杜雨洺, 楊建宇. LFM 信號時間調(diào)-頻斜率分布交叉項抑制方法[J]. 系統(tǒng)工程與電子技術(shù), 2010, 32（2）: 229-231.Du Yu-ming and Yang Jian-yu. Elimination of cross terms of time-frequency rate of LFM signals[J]. Systems Engineering and Electronics, 2010, 32（2）: 229-231.

[16] 李利, 司錫才, 張雯雯, 等. 改進的多分量 LFM 信號參數(shù)估計算法及快速實現(xiàn)[J]. 系統(tǒng)工程與電子技術(shù), 2009, 31（11）:2560-2562.Li Li, Si Xi-cai, Zhang Wen-wen, et al.. Improved estimation algorithm of multi-component LFM signal parameters and its fast implementation[J]. Systems Engineering and Electronics,2009, 31（11）: 2560-2562.

[17] 于媛. 基于Radon變換的LFM信號檢測與參數(shù)估計[J]. 現(xiàn)代防御技術(shù), 2014, 41（1）: 136-141.Yu Yuan. Detection and parameter estimation of linear frequency modulated signals based on Radon transform[J].Modern Defense Technology, 2013, 41（1）: 136-141.