基于時間-調頻斜率分布的多普勒調頻率估計

汪金真 宿紹瑩 陳曾平

1 引言

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

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

2 目標回波模型

在高斯白噪聲環境下目標回波模型可表示為

式中，t是時間；P是散射點數； { sm（t） }1是與目標 P個散射點相對應的回波線性調頻信號；w（t）是與信號無關的均值為0、方差為 δ2的高斯白噪聲。

距離雷達為Rm的散射點回波LFM信號經距離壓縮后的信號為

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

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

3 時間-調頻斜率（TCD）分布

3.1 TCD分布及自項分析

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

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

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

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

3.2 交叉項分析

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

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

圖1 單散射點目標模型及回波LFM信號的TCD

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

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

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

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

3.3 多散射點目標回波LFM信號的TCD

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

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

4 多普勒調頻斜率估計

4.1 參數估計

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

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

圖2 多散射點目標模型及回波LFM信號的TCD

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

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

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

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

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

（5）根據數字多普勒調頻斜率公式 K = k Ts2，得到多散射點目標回波 LFM 信號多普勒調頻斜率精估計值?k 為

4.2 計算量分析

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

表1 TCD計算量分析（N為奇數）

表2 DCD計算量分析（N為偶數）

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

圖3 單散射點和多散射點目標回波TCD在垂直調頻斜率軸上的投影積分

N為偶數時，

當 且 僅 當 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 實驗結果與分析

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

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

圖4 單散射點目標回波LFM信號多普勒調頻斜率估計性能

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

6 結束語

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

圖5 多散射點目標回波LFM信號多普勒調頻斜率估計性能

[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 信號相位域快速高精度參數估計算法[J]. 系統工程與電子技術, 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信號參數估計的最大似然改進算法[J]. 無線電通信技術, 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] 高春霞, 張天騏, 金翔, 等. 基于離散多項式變換的寬帶線性調頻信號波達方向估計[J]. 計算機應用, 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] 龐存鎖. 基于離散多項式相位變換和分數階傅里葉變換的加速目標檢測算法[J]. 電子學報, 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] 杜雨洺, 楊建宇. 基于時間-調頻斜率分布的多線性調頻信號檢測與參數估計[J]. 電子與信息學報, 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 信號時間調-頻斜率分布交叉項抑制方法[J]. 系統工程與電子技術, 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 信號參數估計算法及快速實現[J]. 系統工程與電子技術, 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信號檢測與參數估計[J]. 現代防御技術, 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.