高頻地波雷達海雜波抑制的時頻－矩陣聯合法*

王祎鳴，張 杰，紀永剛，于長軍，毛興鵬

（1.國家海洋局第一海洋研究所，山東 青島 266061；2.哈爾濱工業大學，黑龍江 哈爾濱 150001；3.海洋遙測工程技術研究中心（國家海洋局與航天科技集團共建），山東 青島 266061）

高頻地波雷達具有探測距離遠、范圍廣、持續性強等特點，已成為海洋遙感、遙測的一種重要手段。其利用高頻電磁波沿海面繞射傳播的原理，能檢測到海平面視距以外的船只與低空飛機等目標。但是，作為傳播媒介的海洋中始終存在與雷達波長尺度相近的海浪，在諧振作用下疊加產生很強的回波。同時，由于海水運動和雷達體制等原因導致了這種回波的能量具有一定程度的拓展。對于CFAR（Constant False Alarm Rate）等基于閾值門限的海面目標檢測，相當于在更大的范圍內抬高了檢測基底，導致雜波被誤檢為船只目標，而部分船只無法被測出，造成了虛警與漏警現象。相比于雷電、電臺及電離層雜波等干擾，這種影響的持續性更強，空間跨度更大。作為一種嚴重制約高頻雷達海面目標探測能力的主要干擾［1－2］，海雜波的有效抑制成為高頻雷達領域的重要問題。

海雜波抑制需要考慮兩個關鍵方面：一是保持目標特性，避免有效信號在處理域內的損失，提高信雜噪比；二是雜波抑制手段應為后續的信號與數據處理保留更多信息提取的自由度。高頻超視距雷達海雜波抑制有對消法［3－4］、線性預測法［5］、圖像處理法［6］以及奇異值分解法［7－8］等。其中迭代對消法需要對表征海雜波的復正弦信號的幅度、頻率和初相進行精確估計，否則將致使雜波剩余和擴散，性能也會急劇惡化；線性預測法中的參數估計復雜度高，適應性不強；圖像處理法在去除海雜波的同時，回波局部特征易受破壞；而基于奇異值分解的方法可以從雷達回波信號中有選擇性的去除雜波分量，在避免上述問題的同時，能最大程度上提高信雜噪比。文獻［9］對此類方法進行了比較，指出其還具有對海雜波相關性要求低，不需高維求逆等大運算量計算等優點，但存在的問題是難以準確估計雜波奇異值，導致處理效果的下降。

本文將利用時頻分析與矩陣分解所得信息，提出一種新的海雜波抑制方法。首先，針對高頻地波雷達的特點，分析海雜波信號的頻率、幅度以及時頻特性；然后，通過海雜波信號統計模型與高頻雷達回波狀態方程表明海雜波及船只信號的頻率變化可在此基礎上進行估計，結合雜波特性提出時頻－矩陣聯合的海雜波抑制方法；最后，采用實測數據驗證該方法的有效性。

1 高頻海雜波特性分析

1.1 海雜波頻率

高頻地波雷達對海探測期間，沿海面傳播的電磁波與其波長尺度相近的海浪作用發生Bragg散射，在回波頻譜中表現為一階Bragg峰及二階連續譜。諧振產生的一階峰幅度通常遠強于連續譜，成為船只檢測應用中海雜波的主要成分。一階峰［10］由2個對稱于零頻的尖峰構成，多普勒頻移的理論值為：

式中：fB為布拉格頻率；f0為雷達工作頻率。

當海流存在時，海浪疊加到運動的海面，引起多普勒譜中的2個布拉格尖峰朝同一方向等量偏移［11－12］。由海雜波的這種頻譜性質可知正負布拉格頻率理論位置雖有所偏離，但是頻率間隔相對固定，即2fB，這是海雜波的顯著頻率特性。

1.2 海雜波幅度

海洋是一個不斷變化的動態系統，雷達觀測到的回波能量隨波浪的運動產生變化［13］。海雜波譜主要與雷達工作頻率和海態（風速、風向）有關，根據Barrick［14］的一階和二階雷達散射截面方程進行數值模擬，仿真計算出雷達工作頻率為8.9MHz，4級風速，風向為90°，45°以及20°時各多普勒頻率的海浪散射系數。圖1中的多普勒頻率以海雜波峰值點為準進行了歸一化，風向為海面風與產生諧振海浪的夾角?？梢姡瑔斡娠L向改變，海雜波幅度就可發生較大的波動。

圖1 不同風向條件下的海浪散射系數Fig.1 Wave scattering coefficients under different wind direction

1.3 海雜波的時頻表征

為獲得有利于船只檢測的信噪比，高頻地波雷達相干處理的積累時間通常達數分鐘。傳統的二次傅立葉處理過程，相當于在此期間對回波進行平均處理；而海面處于連續的運動狀態，海水波動頻率也隨之緩慢改變［15］，長時間的積累不能準確反映海雜波幅度與頻率的變化情況。

因此，需采用時頻分析的技術手段，描述信號頻譜含量相對于時間的變化。通過建立時變頻譜的分布，方便在時間和頻率上同時表示信號的局部化信息及能量分布［16］。圖2為實測高頻地波雷達海雜波信號的時頻分布。圖中2條橫向的亮線分別對應于正負一階Bragg峰。

由圖2可見，使用時頻分析的方法觀測時變海面對高頻電磁波的作用，反映出海面回波的更多有效信息：

（1）實測海雜波多普勒頻率的位置較理論公式（1）計算得到的fB＝±0.304 3Hz，產生同向偏移，但峰值頻率間距保持相對固定。據高頻海洋遙感反演原理，表明海面存在徑向流且遠離雷達接收站。

（2）海雜波能量分布是隨頻率變化的，其能量在幾十秒的時間內呈現較強連續性且伴隨有強弱交替，這是由海面回波的相干性以及海面風向與雷達波束夾角的改變所致。

（3）海雜波頻譜能量主要集中于10－2Hz量級的頻率范圍內。從趨勢看，攜帶主要能量的海雜波頻率變化并不劇烈（時間窗口為30s）。

圖2 高頻雷達海雜波時頻分布Fig.2 Time-frequency distribution of sea clutter

2 海雜波信號模型及抑制方法

2.1 海雜波信號模型

海浪可被看作是由各個波長的正弦波疊加形成，高頻雷達海浪回波主成分服從高斯分布，其所對應的回波多普勒頻譜幅度符合瑞利分布。經數字化采樣的一階海雜波通?？捎脧驼倚盘柲M［17］，即

其中：n＝0，1，…，N－1為時間采樣點數；Ts為采樣間隔；fp（n），fn（n）為正負多普勒頻率；A（n），A′（n）為滿足高斯分布的振幅。

船只等目標回波，同樣可以用復正弦信號模型化；高頻地波雷達接收的海面回波信號可表示為雷達分辨單元內海浪與船只回波的疊加，即

其中：B（n），ft（n）分別為船只目標的幅度和頻率。

因此，回波信號中將含有多個有效頻率分量。通常情況下，海雜波及船只的頻率在一定的時間內變化緩慢，可將接收信號建模為振蕩器的輸出。狀態空間模型［5］為

其中：x（n）為狀態向量；F（n）為狀態反饋矩陣；h為輸出向量。

據方程（4）和（5），通過最小二乘法可得到狀態反饋矩陣［18］，進而得到雷達回波信號分量的瞬時頻率w1、w2及w3。由海雜波信號統計模型與雷達回波狀態方程表明，海雜波及船只信號的頻率變化可在此基礎上進行估計。

2.2 海雜波抑制方法

海雜波抑制是一個將雜波分量從復合回波中剔除的過程。利用海雜波在時頻域分析中表現出的特征，結合信號各分量的瞬時頻率估計算法，從矩陣構造與時頻域聯合的角度進行海雜波的抑制。

首先采用矩陣分解的方法，如SVD（Sigular Value Decomposition），通過對特定奇異值置零，移除相對應的特征向量。將高頻雷達回波信號分段構造成長方陣，對N點時域采樣數據y（n），以Hankel矩陣的形式表示，即：

其中：L≥3r；r為數據所含的慢時變信號分量個數，即海雜波與船只回波的數量。

式中：d為L的中值；F（d＋k）為狀態反饋矩陣，寫作

該狀態反饋矩陣可由方程（10）用最小二乘法獲解，信號分量的瞬時頻率可根據矩陣特征值的相位角估計獲得。

短時傅立葉變換（STFT）作為一種經典的時頻分析方法，可通過窗函數的選擇，獲取不同的時間分辨率和頻率分辨率。數字信號x（n）的STFT表達式為

式中：k＝0，1，…，M－1，M為傅立葉點數；N為窗函數g的步長。

海雜波信號波形變化較為平緩，即低頻信息較多，分析窗口可適當延長，獲取較長的時間采樣，提高頻率分辨率。將變換值求平方，得到采樣數據的功率譜。

利用海雜波特性，尤其是在時頻域所表現出的特征，可以進一步確定回波信號中海雜波分量隨時間及頻率的變化。將降秩SVD方法所得到的主分量奇異值與估計的瞬時頻率相匹配，實現海雜波奇異值的判別。反推海雜波奇異值并置零，剔除其中的海雜波分量，得到新的Hankel矩陣。依據原始數據矩陣的構造方式，還原出信號的時間序列。圖3為提出的海雜波抑制方法流程。

圖3 海雜波抑制流程Fig.3 Flowchart of sea clutter suppression

通過時域序列數據構造矩陣時，也可先進行時頻分析，只截取海雜波能量較強時段的雷達回波數據進行處理，以降低運算量。本方法將矩陣SVD和時頻分析抑制算法有機結合，克服了前者雜波奇異值不易確定的問題，也避免了后者的逆變換與補償重構問題。

3 實測數據分析與處理

為驗證本文算法的有效性，利用在某高頻地波雷達站采集的實測數據開展了海雜波抑制實驗。實驗系統的通道數據經插值處理后距離單元分辨率是1.5km，時域采樣點間隔為0.14s。截取130s的雷達回波進行STFT變換，滑窗采用漢明窗，按距離單元分別對數據中所包含的信息進行分析處理。

圖4（a）給出第21距離單元的STFT處理結果，圖中以亮度越亮表示能量越強。由圖4（a）可見，回波能量主要集中于4個頻率變化區間內，即：1區間0.26～0.306 7Hz，2 區間 0.193 3～0.226 7Hz，3 區 間－0.093～－0.126 7Hz及4區間－0.28～－0.326 7 Hz。此時，雷達發射頻率為8.9MHz，據布拉格頻率公式計算一階峰位置為fB＝±0.304 3Hz，正負一階Bragg峰的頻差為0.608 6Hz，頻譜展寬約0.04Hz。另外，參照與之相間隔距離單元的情況，在頻寬較小的頻率范圍2和3內無較強能量分布，滿足船只目標在頻譜中沿距離向的延展性不強的特征。由此可判定譜中頻率范圍1和4為一階海雜波，范圍2和3為船只目標。

圖4（b）為矩陣SVD估計的4個奇異值所對應的瞬時頻率，頻率按相應的奇異值依次表示為w1，w2，w3及w4。其中，w1和w3在±0.3Hz附近小幅波動，在時頻分析識別的一階海雜波頻譜范圍之內，頻率w1和w3對應的奇異值即為海雜波奇異值。

圖4 回波信號時頻分析Fig.4 Time-frequency analysis of received signals

圖5給出海雜波抑制前后的時域波形及速度譜。圖5（a）為只有海雜波存在時的時域信號，圖5（b）為海雜波和2個船只目標同時存在的時域信號，圖5（c）為船只目標所在距離單元的頻譜。從圖5（a）和圖5（b）中可以看出，雜波抑制后海雜波主成分得到了顯著的抑制。由圖5（c）可見，相應距離單元中的目標信噪比提高了約20dB。此外，除海雜波區域外，目標的頻譜特征沒有受到任何影響。

圖5 海雜波抑制前后對比Fig.5 Results before/after suppression

為更加全面地衡量高頻地波雷達作用距離內海雜波的抑制效果，將所有時域回波進行海雜波抑制，利用抑制前后的數據分別形成距離－多普勒（R-D）譜。采用CFAR檢測算法［19］對此二維圖像進行船只目標檢測驗證。首先，分別在距離向、頻率向進行單元平均CFAR目標檢測，綜合虛警概率設為10－6。然后綜合距離向、頻率向檢測結果，保留在距離域和頻域均超過檢測門限的目標。最后，經峰值檢測輸出檢測結果。抑制前檢測到的目標用圓圈標出，抑制后檢測到的用小正方形標出（見圖6）。

圖6 抑制前后檢測結果Fig.6 Detection results before/after suppression

由圖6（a）可見，海雜波抑制前的R-D譜中2條海雜波峰的幅度非常高，強于附近的點狀船只回波散射強度。經檢測，海雜波區出現很多虛假目標（由橢圓框框出），而其周圍的點狀船只目標沒有被檢測到。從圖6（b）給出的抑制后的R-D譜可以看出，海雜波抑制后的船只回波得到凸顯，由海雜波引起的虛假目標被剔除，同時檢測到之前漏檢的點狀船只目標（由矩形框框出），降低了海雜波對檢測背景的影響，提高了船只目標的檢測性能。

4 結語

本文分析了高頻雷達海雜波的頻率、幅度以及時頻特性，提出了基于時頻分析和矩陣分解聯合的海雜波抑制方法。該方法綜合利用多個處理域，豐富了海雜波抑制方法的信息量，提高了海雜波抑制效果。剔除海雜波后信雜比大大提高，而原始數據的其它時域及頻域特征都得到了保留。經雷達實測數據驗證，本文提出的方法有效提高了CFAR船只目標檢測的能力，具有很好的應用前景。

［1］紀永剛，張杰，王祎鳴，等.雙頻地波雷達船只目標點跡關聯與融合處理 ［J］.系統工程與電子技術，2014，36（2）：266-271.

［2］凡俊梅，焦培南，肖景明.海洋雜波對高頻雷達檢測海面上低速目標的影響 ［J］.電波科學學報，1997，12（2）：205-210.

［3］錢文振，紀永剛，王祎鳴，等.一種改進的地波雷達鄰近距離單元格一階海雜波對消方法 ［J］.海洋科學進展，2013，01：138-144.

［4］楊煉，孫合敏，潘新龍.基于擴展Prony算法的海雜波循環對消法［J］.現代雷達，2011，33（6）：53-57.

［5］Khan R.Ocean-Clutter model for high frequency radar［J］.IEEE Journal of Oceanic Engineering，1991，16（2）：181-188.

［6］WANG Jian.High Resolution Methods for Small Target Detection and Estimation in High Frequency Radar［D］.Canada：University of Victoria，2004.

［7］Poon W Y，Khan R，Le-ngoc S.A singular value decomposition based method for suppressing ocean clutter in high frequency radar［J］.IEEE Trans on Signal Processing，1993，41（3）：1421-1424.

［8］Khan R，Power D，Walsh J.Ocean clutter suppression for an HF ground wave radar［C］.IEEE conf on Electrical and Computer Engineering，1997，2：512-515.

［9］LI Yang，JI Zhenyuan，XIE Junhao，et al.Detection of HF firstorder sea clutter and its splitting peaks with image feature：results in strong current shear environment［C］.Advanced Concepts for Intelligent Vision Systems，Berlin：Springer，2012：503-514.

［10］周文瑜，焦培南.超視距雷達技術［M］.北京：電子工業出版社，2008.

［11］Barrick D E，Headrick J M，Bogle R W，et al.Sea backscatter at HF：interpretation and utilization of the echo［J］.Proc of the IEEE，1974，62（6）：673-680.

［12］Panagopoulos S，Soraghan J.Small target detection in sea clutter［J］.IEEE Trans on Geoscience and Remote Sensing，2004，42（7）：1355-1361.

［13］Haykin S，Bakker R，Curriew B.Uncovering nonlinear dynamics-The case study of sea clutter［J］.Proc IEEE，2002，90（5）：860-881.

［14］Lipa B J，Barrick D E.Extraction of sea state from HF radar sea echo：mathematical theory and modeling ［J］.Radio Science，1986，21（1）：81-100.

［15］董英凝，張寧，許榮慶.高頻地波雷達工作環境對系統性能影響的分析 ［J］.電波科學學報，2007，22（2）：325-330.

［16］熊新農，萬顯榮，柯亨玉，等.基于時頻分析的高頻地波雷達電離層雜波抑制 ［J］.系統工程與電子技術，2008，30（8）：20-25.

［17］仇永斌，張寧，張樹春.雙基地高頻地波雷達海雜波抑制［J］.哈爾濱工業大學學報，2012，44（1）：71-77.

［18］Dimonte C，Arun K.Tracking the frequencies of superimposed time-varying harmonics ［C］.Albuquerque：Proc ICASSP 90，1990，5：2536-2542.

［19］王騰，徐向東，范叢旺，等.兩種檢波方式CA-CFAR檢測器性能分析 ［J］.空軍雷達學院學報，2009，23（5）：356-358.