郝國成， 張必超， 鍋 娟， 張雅冰， 石光耀， 王盼盼， 張 薇

(1. 中國地質(zhì)大學(xué)(武漢) 機(jī)械與電子信息學(xué)院，武漢 430074； 2. 杜克大學(xué) 數(shù)學(xué)系，美國Durham 27708； 3. 中國地質(zhì)大學(xué)(武漢) 復(fù)雜系統(tǒng)先進(jìn)控制與智能自動化湖北重點(diǎn)實(shí)驗(yàn)室， 武漢 430074)

時頻分析是處理非平穩(wěn)信號的重要工具之一，能夠?qū)⒁痪S時序信號變換到時間頻率坐標(biāo)軸并提供信號頻率隨時間變化的聯(lián)合分布關(guān)系[1].目前，時頻分析方法已經(jīng)廣泛應(yīng)用于地震信號[2-4]、旋轉(zhuǎn)機(jī)械故障信號[5]、金屬破裂信號[6]等.海雜波信號也具有非平穩(wěn)隨機(jī)的特征，如何提高海面小目標(biāo)(如浮冰、小船、蛙人以及飛機(jī)殘骸)的檢測能力是一項(xiàng)長期且艱巨的任務(wù)[7].傳統(tǒng)時頻分析方法包括短時Fourier變換(STFT)[8]、小波變換(WT)[9]及Wigner-Ville分布(WVD)[10]等.其中，STFT和WT屬于線性時頻分析方法，根據(jù)海森堡不確定準(zhǔn)則，線性時頻方法不能同時實(shí)現(xiàn)時間和頻率的高分辨率.二次型時頻分析方法WVD在處理單分量信號時能夠得到高分辨率的時頻分布，但是在處理多分量信號會出現(xiàn)交叉項(xiàng)干擾，這限制了WVD在現(xiàn)實(shí)生活中的應(yīng)用.

近些年，針對傳統(tǒng)時頻分析方法的不足，許多優(yōu)秀的時頻分析方法被提出以解決上述的問題.根據(jù)選擇的窗函數(shù)如果和待分析信號的內(nèi)部特征相近，就會得到聚集度較高的時頻表示這種思想.參數(shù)化時頻分析方法朝著識別信號的內(nèi)部特征和構(gòu)造窗函數(shù)的方向發(fā)展.線調(diào)頻小波變換(CT)[11]通過引入額外的調(diào)頻率(CR)參數(shù)，在處理線性調(diào)制信號時能夠得到聚集度較高的時頻分布.但是，當(dāng)CR與信號的瞬時頻率差別較大時，通過該方法獲得的時頻表示存在時頻能量模糊現(xiàn)象，不利于信號時頻特征的提取[12].后處理時頻分析方法通過對時頻譜的進(jìn)一步處理，使得傳統(tǒng)時頻分析方法的能量譜聚集度得到極大的改善.再賦值方法(RM)[13]通過在時間軸和頻率軸方向重排時頻點(diǎn)，實(shí)現(xiàn)了時頻分布的高分辨率.然而，RM是完全依賴時頻譜的，導(dǎo)致該方法失去了信號重構(gòu)能力.文獻(xiàn)[14]首先提出同步壓縮理論并應(yīng)用到連續(xù)小波變換(CWT).文獻(xiàn)[15]將同步壓縮理論應(yīng)用到STFT中，基于STFT的同步壓縮變換(SST)極大地改善了STFT時頻模糊現(xiàn)象并且保留了信號的重構(gòu)能力.SST在處理強(qiáng)時變信號時，其時頻分布結(jié)果會出現(xiàn)分辨率降低的現(xiàn)象.為了克服這一問題，涌現(xiàn)了許多新的時頻分析方法，比如2階SST[16-17]、3階SST[18-19]、高階SST[20]等.隨著階數(shù)的提高，這些方法的計(jì)算量也在不斷地增加[21].文獻(xiàn)[21]提出一種新的時頻分析方法，局部最大值同步變換(LMSST).該算法通過在頻率方向檢測時頻譜的局部最大值，構(gòu)造一個全新的頻率再賦值規(guī)則.通過這種再賦值操作，時頻分布的分辨率與SST和RM相比進(jìn)一步得到提高.針對CT存在的問題，本文將LMSST的頻率重排規(guī)則推廣到CT中，提出局部最大值同步壓縮線調(diào)頻小波變換(LMSCT)算法，解決了調(diào)頻率參數(shù)選取對CT時頻表示分辨率的影響.同時，仿真實(shí)驗(yàn)表明，LMSCT比CT具有更好的抗噪性能.最后，將LMSCT應(yīng)用到海雜波信號中，為海雜波背景下的小目標(biāo)檢測提供一個可參考的技術(shù)手段.

1 算法介紹

1.1 同步壓縮理論

一個具有調(diào)幅和調(diào)頻規(guī)律的多分量信號f(t)可以表示為

(1)

式中：Ak(t)和φk(t)分別為第k個分量的瞬時幅度和瞬時相位；t為時間變量.

對于一個待分析信號f滿足f∈L2(R)，L為勒貝格積分，對應(yīng)的窗函數(shù)g為實(shí)函數(shù)且滿足g∈L2(R)，那么此信號的STFT可以定義為

(2)

式中：ω為信號的瞬時頻率；u為時間變量.

SST通過獲得STFT的瞬時頻率信息，然后進(jìn)行頻率點(diǎn)壓縮來提高時頻譜的聚集度.要想獲得STFT的瞬時頻率，首先求取式(2)對時間t的1階導(dǎo)數(shù)，則有：

iωS(t,ω)-Sg′(t,ω)

(3)

STFT的瞬時頻率信息ω0(t,ω)可由如下公式求得：

(4)

(5)

式中：δ(·)為狄拉克函數(shù).通過這一后處理操作，可以得到比原STFT聚集度更高的時頻分布.

1.2 LMSST算法

作為SST算法的改進(jìn)方法，LMSST算法定義了一個新的頻率再賦值規(guī)則，規(guī)則如下式所示：

ωm(t,ω)=

(6)

式中：Δ為離散頻率間隔；|S(t,ω)|為STFT的時頻譜.假如任意兩個分量相隔足夠的距離，由于窗函數(shù)的Fourier變換在0處達(dá)到最大值，式(6)可以簡化為

ωm(t,ω)=

(7)

為了獲得理想時頻分布，所有的模糊時頻系數(shù)都應(yīng)該沿頻率方向分配到時頻軌跡.因此，LMSST可以實(shí)現(xiàn)更高的時頻能量聚集度，可以表示為

(8)

1.3 LMSCT算法

CT同樣是一個非常有效的時頻分析方法，通過引入CR參數(shù)β，其本身就成為了STFT的推廣形式，公式定義為[22]

(9)

如果參數(shù)β能夠很好地匹配待分析信號的調(diào)頻率，就能夠獲得能量聚集度較高的時頻分布結(jié)果.但是對于非線性調(diào)頻信號而言，此時參數(shù)β不能夠很好地匹配變化的調(diào)頻率，因此CT算法在處理這類信號時就會受到限制.

根據(jù)LMSST方法的思想，通過對CT的結(jié)果進(jìn)一步處理，在頻率方向重新分配時頻系數(shù)，再賦值規(guī)則可以表示為

ωL(t,ω)=

(10)

因此，提出一個新的時頻分析方法并命名為局部最大值同步壓縮線調(diào)頻小波變換，算法的定義如以下式所示：

TLMSCT(t,η)=

(11)

傳統(tǒng)高階SST隨著階數(shù)的增加，計(jì)算量也在相應(yīng)增加.原始SST需要執(zhí)行1次STFT步驟，當(dāng)階數(shù)增加到4時，執(zhí)行STFT的次數(shù)就變?yōu)?1次.而LMSCT方法只需要執(zhí)行1次CT，計(jì)算量也相對較小.

2 數(shù)值實(shí)驗(yàn)

為了驗(yàn)證本文提出算法的有效性，需要進(jìn)行2組仿真信號的實(shí)驗(yàn) ，將本文提出的算法與前文引用的STFT、CT和SST方法進(jìn)行對比.首先將文獻(xiàn)[23]中的仿真信號進(jìn)行驗(yàn)證，信號1的函數(shù)表達(dá)式為

f1(t)=sin[2π(10t+5t2/4+t3/9-t4/160)]

(12)

信號1的瞬時頻率是時變的并且是單分量的，實(shí)驗(yàn)中該信號的抽樣頻率為100 Hz.根據(jù)式(12)，在圖1(a)中繪制了信號的真實(shí)頻率分布，圖1(b)展示的是該信號STFT結(jié)果，STFT存在嚴(yán)重的能量擴(kuò)散，不適合分析信號頻率隨時間變化的關(guān)系.在圖2中分別使用CT、SST和LMSCT對信號1做進(jìn)一步處理.

圖1 仿真信號1的時頻分布Fig.1 Time-frequency distribution of simulation Signal 1

圖2(a)的時頻分布圖由CT得到，且β=4π.同STFT的結(jié)果相比，CT在一定程度上提高了該信號時頻分布的聚集度.但是從圖2(a)中的兩處標(biāo)記可以看出，在瞬時頻率較小標(biāo)記處，CR的值與信號的瞬時頻率相近，其時頻分布的聚集度相對較好，在瞬時頻率較大標(biāo)記處，CR的值與信號瞬時頻率的差別較大，時頻能量擴(kuò)散幅度較大.這兩處的對比驗(yàn)證了前文提到的CR值的選取對信號時頻分布的影響.圖2(b) 和2(c)分別由SST與LMSCT得到.與SST的結(jié)果相比，LMSCT的聚集度更高.為了更加直觀地比較，分別截取了SST和LMSCT的局部效果圖進(jìn)行展示，截取部分為圖2(b)和2(c)中的標(biāo)記，局部效果如圖3所示.圖3(a)和3(b)分別為SST與LMSCT的局部展示.通過局部結(jié)果的展示可知，LMSCT的聚集度更高.

圖2 不同算法時頻分布Fig.2 Time frequency distributions of different algorithms

圖3 時頻分布局部效果Fig.3 Local effects of time-frequency distribution

為了系統(tǒng)地比較幾種常見的時頻分析方法聚集度的高低，在這里引入Rényi熵的概念[24].文獻(xiàn)[25]提出將3階Rényi熵用于提供時頻分布的信息度量，文獻(xiàn)[24]詳細(xì)研究了Rényi熵的性質(zhì)和潛在的應(yīng)用價值.Rényi熵是評價時頻分布能量聚集度的指標(biāo)之一.Rényi熵越小代表時頻聚集度越高，相反Rényi熵越大，代表時頻聚集度越低.α階Rényi熵的計(jì)算公式如下式所示：

(13)

式中：T(t,ω)為信號時頻分布.一般情況下，α默認(rèn)為3.其中，信號1在不同算法下的Rényi熵值如表1所示.

表1 不同算法處理信號1得到的Rényi熵

接下來選取一個二分量信號進(jìn)行仿真驗(yàn)證，信號2的函數(shù)表達(dá)式為

f2(t)=sin{2π[40t+sin(1.5t)]}+

sin{2π[17t+6sin(1.5t)]}

(14)

仿真實(shí)驗(yàn)中，設(shè)置信號抽樣頻率為100 Hz，抽樣時間為4 s.信號實(shí)際頻率分布如圖4(a)所示，圖4(b)為WVD的時頻分布結(jié)果.從圖4中可以觀察到，對于多分量信號，由于交叉項(xiàng)干擾的存在，不能有效提取頻率隨時間變化的信息，所以WVD不適合處理這類型的信號.將主流的幾種時頻分析方法與本文提出的方法進(jìn)行比較，對比結(jié)果如圖5所示.由圖5可知，由LMSCT算法得到的時頻分布的能量聚集度比STFT、CT和SST更高.將圖5的時頻分布轉(zhuǎn)換到三維中，如圖6所示，其中：A為幅值；E為能量譜.由圖6可以清晰地觀察到，STFT和CT的時頻分布存在能量擴(kuò)散的情況；SST與LMSCT能量聚集度較高，然而相較于SST的結(jié)果，LMSCT的幅度比較一致，結(jié)果更加準(zhǔn)確.

圖4 仿真信號2時頻結(jié)果Fig.4 Time-frequency distribution of simulation Signal 2

圖5 信號2時頻分布Fig.5 Time-frequency distribution of Signal 2

圖6 時頻分布的三維展示Fig.6 Three-dimensional representation of time-frequency distribution

信號2在不同算法下的Rényi熵值如表2所示.為了更直觀地比較，在圖7中繪制了Rényi熵值的散點(diǎn)圖，其中：R為Rényi熵.

表2 不同算法處理信號2得到的Rényi熵

圖7 不同算法的Rényi熵值Fig.7 Rényi entropies of different algorithms

通過表2和圖7列出的幾種算法Rényi熵值大小的比較，可以確定本文提出的算法具有更高的時頻聚集度.通常情況下，現(xiàn)實(shí)生活中的真實(shí)信號都會含有噪聲，將會影響信號時頻分布的可讀性.接下來對信號2添加高斯白噪聲，然后再進(jìn)行時頻分析算法處理，評估不同算法的抗噪聲性能.

在信號2中加入信噪比為7 dB的加性高斯白噪聲，不同算法的處理結(jié)果如圖8所示.從圖8(a)和8(b)可以看出，受噪聲影響，STFT與CT的時頻分布不是連續(xù)的，嚴(yán)重影響信號時頻變化特征的提取.圖8(c)和8(d)分別為SST與LMSCT的結(jié)果，同STFT、CT結(jié)果相比，時頻譜能量聚集度都有較大的改善.

圖8 不同算法抗噪效果Fig.8 Anti-noise effects of different algorithms

為了進(jìn)一步比較SST與LMSCT的抗噪性能，分別截取SST和LMSCT的部分時頻表示進(jìn)行對比，截取部分為圖8(c)和8(d)標(biāo)記部分，局部結(jié)果如圖9所示.圖9(a)和9(b)分別為SST與LMSCT的局部效果.通過對比可以看出，在相同的頻率點(diǎn)范圍，LMSCT的頻率變化曲線比SST更加平滑，而且沒有出現(xiàn)能量消失的情況，受噪聲影響較小.綜上所述，LMSCT算法的時頻聚集度更高，且抗噪性能更好.

圖9 抗噪性能局部圖Fig.9 Local diagram of anti-noise performance

3 IPIX雷達(dá)數(shù)據(jù)分析

IPIX雷達(dá)由加拿大McMaster大學(xué)Adaptive Systems實(shí)驗(yàn)室設(shè)計(jì)，Haykin教授帶領(lǐng)其團(tuán)隊(duì)分別于1993年和1998年利用IPIX采集并公開了大量高分辨海雜波數(shù)據(jù)，該數(shù)據(jù)已經(jīng)成為測試?yán)走_(dá)檢測算法的重要基準(zhǔn)數(shù)據(jù)，在海雜波特性研究方面也做出了重要的貢獻(xiàn)[26].IPIX雷達(dá)是一個全相參的X波段雷達(dá)，具有I通道和Q通道兩路收發(fā)信號.

實(shí)驗(yàn)中，雷達(dá)架設(shè)在加拿大東海岸一個高出海平面25 ～ 30 m的固定位置，雷達(dá)朝大西洋海面照射，待檢測目標(biāo)是被鋁絲包裹直徑1 m的漂浮圓球.雷達(dá)工作頻率為9.3 GHz，波束寬度為0.9°，距離分辨率為30 m.雷達(dá)工作在駐留模式，連續(xù)接收來自某一確定方向的海面回波，脈沖重復(fù)頻率為 1 kHz，駐留時間約為131 s，每組數(shù)據(jù)包含14個距離單元的回波信號.由于雷達(dá)以低掠射角照射目標(biāo)，目標(biāo)物體隨海面起伏和擺動導(dǎo)致目標(biāo)能量擴(kuò)散，并且在進(jìn)行數(shù)據(jù)采集時采取了距離過采樣，所以目標(biāo)所在單元周圍的臨近單元會受到目標(biāo)能量的影響，記為受影響單元.每個數(shù)據(jù)文件都由純海雜波、目標(biāo)所在單元回波和受影響目標(biāo)單元回波組成.

測量該數(shù)據(jù)時的雷達(dá)環(huán)境參數(shù)及數(shù)據(jù)組成如表3所示.回波信號幅度隨距離和時間的變化關(guān)系如圖10所示，其中：U為距離單元.每個距離單元由 131 072 個時域數(shù)據(jù)點(diǎn)組成，目標(biāo)位于第8個單元，距離為2.66 km，由于目標(biāo)的起伏和漂移，7、9、10這3個單元成為受影響單元.為了進(jìn)一步分析，本文選取水平發(fā)射水平接收(HH)極化下的I通道數(shù)據(jù)進(jìn)行處理與分析.

表3 IPIX雷達(dá)54#文件環(huán)境參數(shù)及數(shù)據(jù)組成

圖10 海雜波時間-距離-能量圖Fig.10 Didgram of sea clutter time-distance-energy

單元1不是單元8的鄰近單元，在數(shù)據(jù)采集過程中不會受到目標(biāo)起伏和漂移的影響.首先處理該單元的數(shù)據(jù)，分析純凈海雜波信號的時頻分布結(jié)果.單元1純凈海雜波使用3種算法的處理結(jié)果如圖11所示，其中：Sp為抽樣點(diǎn).從圖11(a)和11(b)中可以看出，時頻分布結(jié)果包含了大量的背景噪聲，并且從圖11(c)的標(biāo)記部分展示的結(jié)果可知，LMSCT算法較好地抑制了噪聲的影響，可以為后續(xù)的含噪聲信號的處理提供參考依據(jù).接下來將對目標(biāo)單元及其鄰近單元的數(shù)據(jù)進(jìn)行處理，對目標(biāo)單元和其他單元數(shù)據(jù)的時頻分布結(jié)果進(jìn)行分析，提取有效的時變特征信息.

圖11 純凈海雜波(單元1)Fig.11 Pure sea clutter (Unit 1)

圖12、13、15和16為單元8的鄰近單元處理結(jié)果，分別為單元6、單元7、單元9和單元10.圖14為含有目標(biāo)的單元8處理結(jié)果.由圖14可知，對含目標(biāo)的單元8進(jìn)行處理后，與單元1的處理結(jié)果相比，可以觀察到在0頻附近出現(xiàn)清晰的頻率曲線，并且通過對比圖14(a)～14(c)的時頻結(jié)果發(fā)現(xiàn)，LMSCT的頻率曲線更加精細(xì)，更有利于時頻變化特征的提取.從圖13、15和16的處理結(jié)果中可以看出，由于距離目標(biāo)單元較近，所以也會出現(xiàn)相似的頻率曲線.由圖12可知，同樣距離目標(biāo)單元較近的單元6卻沒有受到影響.這個現(xiàn)象說明，此影響并不是均勻向兩邊擴(kuò)散.若不能提前確定單元8是含目標(biāo)的數(shù)據(jù)，很難從單元7、8、9、10中確定出目標(biāo)單元.為了能夠進(jìn)一步確定目標(biāo)所在單元的位置，采用前文提到的Rényi熵.熵代表信號的無序性，也就是信息量越大，不確定性越小，對應(yīng)熵值越小.時頻分布Rényi熵越小的單元，其含有的信息量越大，即出現(xiàn)目標(biāo)的可能性越大.單元1、6、7、8、9和10分別使用LMSCT算法得到的時頻分布后，計(jì)算取得的Rényi熵值如圖17所示.從每個單元Rényi熵的散點(diǎn)圖可以得出，在單元8的Rényi熵值最小，此單元含有的信息量最大，因此可以判定單元8就是含目標(biāo)的單元.

圖12 純凈海雜波(單元6)Fig.12 Pure sea clutter (Unit 6)

圖13 受影響海雜波(單元7)Fig.13 Affected sea clutter (Unit 7)

圖14 含目標(biāo)海雜波(單元8)Fig.14 Target sea clutter (Unit 8)

圖15 受影響海雜波(單元9)Fig.15 Affected sea clutter (Unit 9)

圖16 受影響海雜波(單元10)Fig.16 Affected sea clutter (Unit 10)

圖17 由LMSCT得到的Rényi熵Fig.17 Rényi entropies obtained by LMSCT

4 結(jié)語

為提高信號時頻分布的聚集度及抑制噪聲干擾的性能，本文提出一種新的高質(zhì)量的時頻分析算法，即LMSCT.該算法主要思想是將LMSST與CT算法結(jié)合，首先對待分析信號進(jìn)行CT處理，采用局部最大值同步壓縮規(guī)則對CT結(jié)果進(jìn)行頻率點(diǎn)的重新分配，進(jìn)而可得到目標(biāo)信號較為精確的時頻分布圖.數(shù)值實(shí)驗(yàn)結(jié)果表明，LMSCT與其他時頻分析方法相比，該算法時頻分布的能量聚集度更高且抗噪聲性能更好.實(shí)際IPIX雷達(dá)信號分析應(yīng)用中，與目前主流的時頻分析算法相比，LMSCT有較小的Rényi熵與瞬時頻率聚集度更高的優(yōu)點(diǎn)以及大幅抑制噪聲干擾，通過比較目標(biāo)信號的Rényi熵，可以有效地確定目標(biāo)所在的距離單元.