基于時頻圖像形態學的雷達信號IF估計

白 航，李傳兵，李富軍，徐 文

（中國人民解放軍61906部隊，河北 廊坊065001）

0 引言

雷達信號是一種非平穩信號，傳統的Fourier頻率不能對其進行有效描述，瞬時頻率（IF）反映了信號頻率隨時間的變化規律，是表征非平穩信號的重要參數。因此IF估計是分析和研究雷達信號的一項基本任務。時頻分析是目前IF 估計研究中較為普遍和有效的方法［1-7］，根據信號的時頻聚集特性，信號的能量在時頻面上將沿著IF的方向產生聚集，從而可以采用時頻分布峰值檢測法對信號進行IF 估計［3-4］。在信噪比較高的情況下，時頻分布峰值法能準確估計出的信號IF，但隨著信噪比的降低，時頻分布不能準確地描述信號的能量分布，因而對于IF的估計性能也大大降低。文獻［5］提出了自適應窗長WVD（Wigner-Ville Distribution）峰值檢測法，在計算WVD 時選擇一個最優窗長使得IF 估計的均方根誤差最小。為了避免初始估計中出現大的偏差，一般采用較窄的窗長，然而隨著信噪比的降低，較窄的數據窗長會使得時頻分布的峰值不在信號自項區域的概率增大，從而導致IF估計出現較大的偏差。文獻［6］綜合利用Gabor變換WVD 的特點，提出一種自適應時頻分布，然后采用時頻分布一階矩方法估計信號IF，該方法有效降低了運算量，但在低信噪比條件下信號IF估計效果也不是很理想。

將圖像處理技術和信號處理相結合，為調頻信號IF估計提供了新的視角。針對較低信噪比條件下以及多分量信號的IF估計，本文首先對信號進行時頻變換，并將其轉化為灰度圖；然后檢測信號的起止頻率，剪切出信號時頻分布的有效區域［8］；最后采用形態學圖像處理方法估計出信號的IF。將信號轉化為時頻圖像后，可以進一步利用圖像處理中的降噪算法，增強信號的前景像素，降低噪聲對IF估計的影響。具體流程如圖1所示。

圖1 本文IF估計方法流程圖

1 MBD時頻分布

雷達信號是一種典型的非平穩信號，傳統的時域和頻域分析方法只能獲取有限的信號信息，時頻分析反映了信號能量隨時間和頻率的分布，能在時頻域更為精確地描述信號。本文應用時頻分析方法估計信號的IF曲線，為了提升估計效果，需要時頻分布具有較高的時頻聚集性，使得信號時頻能量聚集在IF曲線附近，同時又能盡可能消除交叉項的影響。Cohen類時頻分布通過核函數對WVD 進行平滑，在抑制交叉項和保持高時頻分辨率之間做一個折中，其定義為：

Hussainn和Boashash 于2002年提出一種改進的B分布（MBD）的時頻分布方法［7］，其核函數為：

式中，kβ＝Γ（2β）／（22β-1Γ2（β）），Γ 為gamma函數，即：MBD 能滿足時頻分布的大多數特性要求，其核函數滿足二維低通特性。從時頻分布的時頻聚集性、交叉項抑制能力、時頻分辨率和噪聲抑制能力等綜合指標來看，相對其他的二次時頻分布，MBD 分布性能最優。

2 形態學的時頻圖像處理

信號的時頻分布可以看作是一幅二維圖像，因而可以采用圖像處理方法對時頻圖像做進一步處理。首先將時頻圖像轉化為灰度圖，圖像中像素點的不同灰度值對應時頻點的能量值。從時頻圖像中可以看出，信號的時頻分布區域聚集在IF曲線周圍，而噪聲的時頻點散布在整個時頻面上。時頻面上信號的自項可以看作是圖像中的“對象”，而噪聲和交叉項則構成了圖像的“背景”。本文從圖像處理角度對信號的時頻表示結果進行處理，實際上就是在灰度圖像中去除背景而保留對象的過程。

各個信號的時頻圖像灰度值的動態范圍是不一樣的，為了減少數據間的不平衡性，首先對時頻圖像灰度值進行歸一化，然后使用自適應維納濾波器去除時頻圖像的噪聲點，對圖像進行增強。從時頻圖中可以看出，并非所有區域都分布有信號，對此可以檢測分析信號的起止頻率，將沒有信號分布的圖像區域剪切掉［8］，減小冗余信息，更有利于下一步對時頻圖像的分析。接著本文對時頻圖像依次進行二值化、形態學處理和標注連接分量，最終得到信號的IF估計。

2.1 時頻圖像二值化

時頻圖像二值化實際上是對圖像進行閾值處理，將圖像上的灰度值置為0或1，將256個亮度等級的灰度圖像通過適當的閾值選取轉化為仍然可以反映圖像整體和局部特征二值圖像，同時也減少后期圖像處理的計算量和存儲空間，圖像的二值化處理可以描述如下：

選擇合理的閾值Thr是時頻圖像二值化的關鍵，對時頻圖像二值化時，應盡量保留信號在時頻圖中對應的像素點，并盡可能去除噪聲。本文閾值選取參照文獻［9］中方法。

2.2 時頻圖像形態學處理

時頻圖像進行二值化后，信號分量被進一步展寬，本文進一步采用數學形態學方法消除時頻面上的噪聲，細化信號分量，最終估計出信號的IF。形態學處理是應用具有一定形態的結構元素對集合進行腐蝕和膨脹的操作，膨脹使時頻圖連通域擴張，腐蝕使時頻圖連通域收縮。開運算先腐蝕再膨脹，用于濾除圖像中區域小于結構元素的時頻獨立點或明顯區別于信號分量的斑點，而保留相應時頻聚集面積大于結構元素的時頻點，從而使信號在時頻分布平面對應的自分量的輪廓變得光滑，消除時頻分布平面上少量噪聲對應的細的突出物，經形態學開操作處理后的二值圖像可表示為：

式中，B1和B2分別為腐蝕和膨脹的結構元素，Θ 表示腐蝕運算，⊕表示膨脹運算。本文中B1選擇半徑為5的圓盤型結構元素，B2擇半徑為3的菱型結構元素。

形態學骨骼化可以把二值圖像區域縮成單像素的線條，以逼近區域的中心線，提取骨架的目的是減少圖像成分，只留下時頻圖像最基本信息，要求最大限度地細化原圖形，并且要求原圖像中屬于同一連通域的像素不出現斷裂。本文通過對時頻圖像骨骼化，找出時頻能量脊線，由于時頻能量沿著IF 曲線方向聚集，因而時頻能量脊線和IF曲線方向是一致的。圖像A 的骨骼化表示如下：

式中，Sk（A）為骨骼子集，（AΘkB）表示對A 連續腐蝕k 次，。表示開運算。時頻圖像骨骼化后會出現許多毛刺，對此可以采用去毛刺算法［10］，平滑所得到的時頻脊線。

2.3 IF估計

二值時頻圖像是由以前景像素為基本單位組成的若干個連接分量構成的，因此找出信號項對應的時頻圖像上的連接分量，即可確定信號的IF。而對于像素點元素比較少的連接分量可以認為是噪聲，可以通過統計各連接分量像素點的個數，剔除噪聲分量。對此文中采用標注連接分量方法［11］得到信號的連接分量，統計連接分量前景像素點的行索引和列索引，即可估計出信號的IF。當信號中存在多個分量時，采用連接分量標記算法同樣可以區分出各個信號分量。圖2中以LFM 信號（SNR＝-3dB）為例，說明了本文中時頻圖像的處理流程，由左至右分別是信號的時頻圖、經剪切后時頻灰度圖、二值化時頻圖、開運算后的時頻圖、骨骼化后時頻圖以及去毛刺后的時頻圖。

圖2 時頻圖像處理流程

3 仿真實驗

本節通過Matlab仿真對本文IF估計方法性能進行驗證。實驗中分別生成LFM（線性調頻）、SFM（正弦頻率調制）、BFSK（二進制頻移鍵控）和EQFM（偶二次調頻）四種信號。其中LFM 信號載頻設為25MHz，SFM 載 頻 設 為15MHz，BFSK 上 邊 頻 為10MHz、下邊頻為20MHz，EQFM 載頻設為10MHz，采樣頻率均為200MHz，脈沖寬度為11μs，為了簡化分析和計算，信號幅度設為1，仿真時噪聲采用高斯白噪聲。時頻分布窗函數采用改進的B 分布函數，核函數參數β取為0.05，時頻窗長設為161。

首先采用所提出方法分別對LFM、BFSK 和EQFM 信號的IF進行估計，信噪比為-3dB。圖3（a）為LFM 信號的IF估計曲線，可以看出IF估計曲線比較平滑，較為準確地描述了信號真實的IF 變化規律；圖3（b）為FSK 信號的IF 估計曲線，同樣得到了該信號的有效估計，表明本文IF估計方法適用于頻率突變信號；圖3（c）為EQFM 信號的估計曲線，由于該信號的時頻能量主要聚集在IF曲線波谷處，信號兩端能量分布較少，因而在信號兩端的IF 估計會出現偏差，但總體上來看其IF估計值是較為準確的。

圖3 本文方法的IF估計曲線（SNR＝-3dB）

進一步將所提出方法同WVD 峰值檢測法［4］以及時頻分布一階矩法［6］的IF估計效果進行比較。信噪比變化范圍為-9～15dB，分別對三種IF 估計方法做500次Monte-Carlo實驗。CRLB為正弦頻率調制信號IF估計的克拉美羅界。表1對SFM 信號IF估計的MSE 隨信噪比變化情況進行了統計。總的來說，時頻分布一階矩法的估計性能最差，在無噪聲的信號環境下，采用時頻分布一階矩可得到無偏估計，但當信號中有噪聲干擾時，該方法的估計性能急劇下降；當SNR≥6dB時，WVD 峰值檢測法和本文方法的MSE 都接近克拉美羅界，而當SNR≤3dB 時，本文方法明顯優于WVD 峰值檢測法，表明本文方法在較低信噪比的IF估計性能要優于其它兩種方法。從統計特性上來看，本文方法的估計性能有一定程度的提升。

表1 瞬時頻率估計性能統計 （dB）

圖5 多分量信號的IF估計方法比較（SNR＝-3dB）

圖4為SNR＝-3dB時單分量信號的IF 估計效果圖，由圖中可以看出，時頻分布一階矩法得到的IF估計和真實IF曲線明顯差別最大；WVD 峰值檢測法得到的IF估計曲線也不是很理想，受噪聲的影響，許多信號時頻分布的峰值點遠離IF曲線，產生了許多突變點；本文方法得到的IF曲線與真實的IF比較接近，表明采用圖像處理方法能有效降低噪聲對信號IF 估計的影響。圖5為SNR＝-3dB 時多分量信號的IF估計效果圖，由圖中可以看出，傳統的時頻分布一階矩和WVD 峰值法將會失效，而本文方法仍能得到有效的IF估計，主要是因為采用圖像處理中的標注連接分量方法可以有效區分出時頻面上的各個信號分量，所以本文方法也適用于多分量信號的IF 估計。實驗結果驗證了本文方法的有效性。

4 結束語

本文提出了一種將時頻分析和圖像處理相結合的雷達信號IF 估計方法，將時頻分布轉化為二值圖像后，采用圖像處理中的形態學方法估計出信號的IF，實驗結果表明該方法在較低信噪比條件下能夠獲得質量較好的瞬時頻率曲線，適用于多分量調頻信號的IF估計。隨著信噪比的降低，許多信號處理中的IF估計方法效果會變得很差，而采用圖像處理方法能有效降低噪聲對估計性能的影響。IF包含了豐富的信號調制信息，因此本文的研究內容對于雷達信號參數估計和調制方式識別具有較為重要的參考價值。■

［1］劉潮，李政杰，童寧寧.改進的相位建模法在瞬時頻率估計中的應用［J］.航天電子對抗，2011，27（3）：23-25.

［2］Orovic′I，Stankovic′S，Thayaparan T，et al.Multiwindow S-method for instantaneous frequency estimation and its application in radar signal analysis［J］.IET Signal Processing，2010，4（4）：363-370.

［3］Khan NA，Taj IA，Jaffri M.Instantaneous frequency estimation using fractional Fourier transform and Wigner distribution［C］∥Bangalore：ICSAP’10 International Conference on Signal Acquisition and Processing，2010：319-321.

［4］Chen Guanghua，Ma Shiwei，Liu Ming，et al.Wigner-Ville distribution and cross Wigner-Ville distribution of noisy signals［J］.Journal of Systems Engineering and Electronics，2008，19（5）：1053-1057.

［5］Lerga J，Sucic V.Nonlinear IF estimation based on the pseudo WVD adapted using the improved sliding pairwise ICI rule［J］.IEEE Signal Processing Letters，2009，16（11）：953-956.

［6］Baraniuk RG，Coates M，Steeghs P.Hybrid linear／quadratic time–frequency attributes［J］.IEEE Trans.Signal Processing，2001，49（4）：760-766.

［7］Hussain Z，Boashash B.Adaptive instantaneous frequency estimation of multi-component FM signals using quadratic time frequency distributions［J］.IEEE Trans.Signal Processing，2002，50（8）：1866-1876.

［8］Zilberman ER，Pace PE.Autonomous time-frequency morphological feature extraction algorithm for LPI radar modulation classification［C］∥Atlanta，GA：Proc.of the IEEE International Conference on Image Processing，2006：2321-2324.

［9］尚海燕，水鵬朗，張守宏，等.基于時頻形態學濾波的能量積累 檢 測［J］.電 子 與 信 息 學 報，2007，29（6）：1416-1420.

［10］Gonzalez RC，Woods RE，Eddins SL.Digital image processing using MATLAB［M］.2nd ed.Gatemark Publishing，Inc.，2009.

［11］Lagrange M，Marchand S.Estimating the instantaneous frequency of sinusoidal components using phase-based methods［J］.Journal of the Audio Engineering Society，2007，1（55）：385-399.