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

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

0 引言

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

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

圖1 本文IF估計(jì)方法流程圖

1 MBD時頻分布

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

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

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

2 形態(tài)學(xué)的時頻圖像處理

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

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

2.1 時頻圖像二值化

時頻圖像二值化實(shí)際上是對圖像進(jìn)行閾值處理，將圖像上的灰度值置為0或1，將256個亮度等級的灰度圖像通過適當(dāng)?shù)拈撝颠x取轉(zhuǎn)化為仍然可以反映圖像整體和局部特征二值圖像，同時也減少后期圖像處理的計(jì)算量和存儲空間，圖像的二值化處理可以描述如下：

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

2.2 時頻圖像形態(tài)學(xué)處理

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

式中，B1和B2分別為腐蝕和膨脹的結(jié)構(gòu)元素，Θ 表示腐蝕運(yùn)算，⊕表示膨脹運(yùn)算。本文中B1選擇半徑為5的圓盤型結(jié)構(gòu)元素，B2擇半徑為3的菱型結(jié)構(gòu)元素。

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

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

2.3 IF估計(jì)

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

圖2 時頻圖像處理流程

3 仿真實(shí)驗(yàn)

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

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

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

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

表1 瞬時頻率估計(jì)性能統(tǒng)計(jì) （dB）

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

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

4 結(jié)束語

本文提出了一種將時頻分析和圖像處理相結(jié)合的雷達(dá)信號IF 估計(jì)方法，將時頻分布轉(zhuǎn)化為二值圖像后，采用圖像處理中的形態(tài)學(xué)方法估計(jì)出信號的IF，實(shí)驗(yàn)結(jié)果表明該方法在較低信噪比條件下能夠獲得質(zhì)量較好的瞬時頻率曲線，適用于多分量調(diào)頻信號的IF估計(jì)。隨著信噪比的降低，許多信號處理中的IF估計(jì)方法效果會變得很差，而采用圖像處理方法能有效降低噪聲對估計(jì)性能的影響。IF包含了豐富的信號調(diào)制信息，因此本文的研究內(nèi)容對于雷達(dá)信號參數(shù)估計(jì)和調(diào)制方式識別具有較為重要的參考價值?！?/p>

［1］劉潮，李政杰，童寧寧.改進(jìn)的相位建模法在瞬時頻率估計(jì)中的應(yīng)用［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］尚海燕，水鵬朗，張守宏，等.基于時頻形態(tài)學(xué)濾波的能量積累 檢 測［J］.電 子 與 信 息 學(xué) 報(bào)，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.