高質量LMSCT時頻分析算法及其在雷達信號目標檢測中的應用

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

(1. 中國地質大學(武漢) 機械與電子信息學院，武漢 430074； 2. 杜克大學 數學系，美國Durham 27708； 3. 中國地質大學(武漢) 復雜系統先進控制與智能自動化湖北重點實驗室， 武漢 430074)

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

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

1 算法介紹

1.1 同步壓縮理論

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

(1)

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

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

(2)

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

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

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

(3)

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

(4)

(5)

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

1.2 LMSST算法

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

ωm(t,ω)=

(6)

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

ωm(t,ω)=

(7)

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

(8)

1.3 LMSCT算法

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

(9)

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

根據LMSST方法的思想，通過對CT的結果進一步處理，在頻率方向重新分配時頻系數，再賦值規則可以表示為

ωL(t,ω)=

(10)

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

TLMSCT(t,η)=

(11)

傳統高階SST隨著階數的增加，計算量也在相應增加.原始SST需要執行1次STFT步驟，當階數增加到4時，執行STFT的次數就變為11次.而LMSCT方法只需要執行1次CT，計算量也相對較小.

2 數值實驗

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

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

(12)

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

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

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

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

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

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

(13)

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

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

接下來選取一個二分量信號進行仿真驗證，信號2的函數表達式為

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

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

(14)

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

圖4 仿真信號2時頻結果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熵值的散點圖，其中：R為Rényi熵.

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

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

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

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

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

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

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

3 IPIX雷達數據分析

IPIX雷達由加拿大McMaster大學Adaptive Systems實驗室設計，Haykin教授帶領其團隊分別于1993年和1998年利用IPIX采集并公開了大量高分辨海雜波數據，該數據已經成為測試雷達檢測算法的重要基準數據，在海雜波特性研究方面也做出了重要的貢獻[26].IPIX雷達是一個全相參的X波段雷達，具有I通道和Q通道兩路收發信號.

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

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

表3 IPIX雷達54#文件環境參數及數據組成

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

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

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

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

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

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

圖14 含目標海雜波(單元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 結語

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