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

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

0 引言

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

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

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

1 MBD時(shí)頻分布

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

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

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

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

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

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

2.1 時(shí)頻圖像二值化

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

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

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

時(shí)頻圖像進(jìn)行二值化后，信號(hào)分量被進(jìn)一步展寬，本文進(jìn)一步采用數(shù)學(xué)形態(tài)學(xué)方法消除時(shí)頻面上的噪聲，細(xì)化信號(hào)分量，最終估計(jì)出信號(hào)的IF。形態(tài)學(xué)處理是應(yīng)用具有一定形態(tài)的結(jié)構(gòu)元素對(duì)集合進(jìn)行腐蝕和膨脹的操作，膨脹使時(shí)頻圖連通域擴(kuò)張，腐蝕使時(shí)頻圖連通域收縮。開運(yùn)算先腐蝕再膨脹，用于濾除圖像中區(qū)域小于結(jié)構(gòu)元素的時(shí)頻獨(dú)立點(diǎn)或明顯區(qū)別于信號(hào)分量的斑點(diǎn)，而保留相應(yīng)時(shí)頻聚集面積大于結(jié)構(gòu)元素的時(shí)頻點(diǎn)，從而使信號(hào)在時(shí)頻分布平面對(duì)應(yīng)的自分量的輪廓變得光滑，消除時(shí)頻分布平面上少量噪聲對(duì)應(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ū)域的中心線，提取骨架的目的是減少圖像成分，只留下時(shí)頻圖像最基本信息，要求最大限度地細(xì)化原圖形，并且要求原圖像中屬于同一連通域的像素不出現(xiàn)斷裂。本文通過對(duì)時(shí)頻圖像骨骼化，找出時(shí)頻能量脊線，由于時(shí)頻能量沿著IF 曲線方向聚集，因而時(shí)頻能量脊線和IF曲線方向是一致的。圖像A 的骨骼化表示如下：

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

2.3 IF估計(jì)

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

圖2 時(shí)頻圖像處理流程

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

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

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

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

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

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

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

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

4 結(jié)束語

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

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