高頻地波雷達(dá)海雜波抑制的時(shí)頻－矩陣聯(lián)合法*

王祎鳴，張 杰，紀(jì)永剛，于長(zhǎng)軍，毛興鵬

（1.國(guó)家海洋局第一海洋研究所，山東 青島 266061；2.哈爾濱工業(yè)大學(xué)，黑龍江 哈爾濱 150001；3.海洋遙測(cè)工程技術(shù)研究中心（國(guó)家海洋局與航天科技集團(tuán)共建），山東 青島 266061）

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

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

本文將利用時(shí)頻分析與矩陣分解所得信息，提出一種新的海雜波抑制方法。首先，針對(duì)高頻地波雷達(dá)的特點(diǎn)，分析海雜波信號(hào)的頻率、幅度以及時(shí)頻特性；然后，通過海雜波信號(hào)統(tǒng)計(jì)模型與高頻雷達(dá)回波狀態(tài)方程表明海雜波及船只信號(hào)的頻率變化可在此基礎(chǔ)上進(jìn)行估計(jì)，結(jié)合雜波特性提出時(shí)頻－矩陣聯(lián)合的海雜波抑制方法；最后，采用實(shí)測(cè)數(shù)據(jù)驗(yàn)證該方法的有效性。

1 高頻海雜波特性分析

1.1 海雜波頻率

高頻地波雷達(dá)對(duì)海探測(cè)期間，沿海面?zhèn)鞑サ碾姶挪ㄅc其波長(zhǎng)尺度相近的海浪作用發(fā)生Bragg散射，在回波頻譜中表現(xiàn)為一階Bragg峰及二階連續(xù)譜。諧振產(chǎn)生的一階峰幅度通常遠(yuǎn)強(qiáng)于連續(xù)譜，成為船只檢測(cè)應(yīng)用中海雜波的主要成分。一階峰［10］由2個(gè)對(duì)稱于零頻的尖峰構(gòu)成，多普勒頻移的理論值為：

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

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

1.2 海雜波幅度

海洋是一個(gè)不斷變化的動(dòng)態(tài)系統(tǒng)，雷達(dá)觀測(cè)到的回波能量隨波浪的運(yùn)動(dòng)產(chǎn)生變化［13］。海雜波譜主要與雷達(dá)工作頻率和海態(tài)（風(fēng)速、風(fēng)向）有關(guān)，根據(jù)Barrick［14］的一階和二階雷達(dá)散射截面方程進(jìn)行數(shù)值模擬，仿真計(jì)算出雷達(dá)工作頻率為8.9MHz，4級(jí)風(fēng)速，風(fēng)向?yàn)?0°，45°以及20°時(shí)各多普勒頻率的海浪散射系數(shù)。圖1中的多普勒頻率以海雜波峰值點(diǎn)為準(zhǔn)進(jìn)行了歸一化，風(fēng)向?yàn)楹Ｃ骘L(fēng)與產(chǎn)生諧振海浪的夾角?？梢?，單由風(fēng)向改變，海雜波幅度就可發(fā)生較大的波動(dòng)。

圖1 不同風(fēng)向條件下的海浪散射系數(shù)Fig.1 Wave scattering coefficients under different wind direction

1.3 海雜波的時(shí)頻表征

為獲得有利于船只檢測(cè)的信噪比，高頻地波雷達(dá)相干處理的積累時(shí)間通常達(dá)數(shù)分鐘。傳統(tǒng)的二次傅立葉處理過程，相當(dāng)于在此期間對(duì)回波進(jìn)行平均處理；而海面處于連續(xù)的運(yùn)動(dòng)狀態(tài)，海水波動(dòng)頻率也隨之緩慢改變［15］，長(zhǎng)時(shí)間的積累不能準(zhǔn)確反映海雜波幅度與頻率的變化情況。

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

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

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

（2）海雜波能量分布是隨頻率變化的，其能量在幾十秒的時(shí)間內(nèi)呈現(xiàn)較強(qiáng)連續(xù)性且伴隨有強(qiáng)弱交替，這是由海面回波的相干性以及海面風(fēng)向與雷達(dá)波束夾角的改變所致。

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

圖2 高頻雷達(dá)海雜波時(shí)頻分布Fig.2 Time-frequency distribution of sea clutter

2 海雜波信號(hào)模型及抑制方法

2.1 海雜波信號(hào)模型

海浪可被看作是由各個(gè)波長(zhǎng)的正弦波疊加形成，高頻雷達(dá)海浪回波主成分服從高斯分布，其所對(duì)應(yīng)的回波多普勒頻譜幅度符合瑞利分布。經(jīng)數(shù)字化采樣的一階海雜波通常可用復(fù)正弦信號(hào)模擬［17］，即

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

船只等目標(biāo)回波，同樣可以用復(fù)正弦信號(hào)模型化；高頻地波雷達(dá)接收的海面回波信號(hào)可表示為雷達(dá)分辨單元內(nèi)海浪與船只回波的疊加，即

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

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

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

據(jù)方程（4）和（5），通過最小二乘法可得到狀態(tài)反饋矩陣［18］，進(jìn)而得到雷達(dá)回波信號(hào)分量的瞬時(shí)頻率w1、w2及w3。由海雜波信號(hào)統(tǒng)計(jì)模型與雷達(dá)回波狀態(tài)方程表明，海雜波及船只信號(hào)的頻率變化可在此基礎(chǔ)上進(jìn)行估計(jì)。

2.2 海雜波抑制方法

海雜波抑制是一個(gè)將雜波分量從復(fù)合回波中剔除的過程。利用海雜波在時(shí)頻域分析中表現(xiàn)出的特征，結(jié)合信號(hào)各分量的瞬時(shí)頻率估計(jì)算法，從矩陣構(gòu)造與時(shí)頻域聯(lián)合的角度進(jìn)行海雜波的抑制。

首先采用矩陣分解的方法，如SVD（Sigular Value Decomposition），通過對(duì)特定奇異值置零，移除相對(duì)應(yīng)的特征向量。將高頻雷達(dá)回波信號(hào)分段構(gòu)造成長(zhǎng)方陣，對(duì)N點(diǎn)時(shí)域采樣數(shù)據(jù)y（n），以Hankel矩陣的形式表示，即：

其中：L≥3r；r為數(shù)據(jù)所含的慢時(shí)變信號(hào)分量個(gè)數(shù)，即海雜波與船只回波的數(shù)量。

式中：d為L(zhǎng)的中值；F（d＋k）為狀態(tài)反饋矩陣，寫作

該狀態(tài)反饋矩陣可由方程（10）用最小二乘法獲解，信號(hào)分量的瞬時(shí)頻率可根據(jù)矩陣特征值的相位角估計(jì)獲得。

短時(shí)傅立葉變換（STFT）作為一種經(jīng)典的時(shí)頻分析方法，可通過窗函數(shù)的選擇，獲取不同的時(shí)間分辨率和頻率分辨率。數(shù)字信號(hào)x（n）的STFT表達(dá)式為

式中：k＝0，1，…，M－1，M為傅立葉點(diǎn)數(shù)；N為窗函數(shù)g的步長(zhǎng)。

海雜波信號(hào)波形變化較為平緩，即低頻信息較多，分析窗口可適當(dāng)延長(zhǎng)，獲取較長(zhǎng)的時(shí)間采樣，提高頻率分辨率。將變換值求平方，得到采樣數(shù)據(jù)的功率譜。

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

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

通過時(shí)域序列數(shù)據(jù)構(gòu)造矩陣時(shí)，也可先進(jìn)行時(shí)頻分析，只截取海雜波能量較強(qiáng)時(shí)段的雷達(dá)回波數(shù)據(jù)進(jìn)行處理，以降低運(yùn)算量。本方法將矩陣SVD和時(shí)頻分析抑制算法有機(jī)結(jié)合，克服了前者雜波奇異值不易確定的問題，也避免了后者的逆變換與補(bǔ)償重構(gòu)問題。

3 實(shí)測(cè)數(shù)據(jù)分析與處理

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

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

圖4（b）為矩陣SVD估計(jì)的4個(gè)奇異值所對(duì)應(yīng)的瞬時(shí)頻率，頻率按相應(yīng)的奇異值依次表示為w1，w2，w3及w4。其中，w1和w3在±0.3Hz附近小幅波動(dòng)，在時(shí)頻分析識(shí)別的一階海雜波頻譜范圍之內(nèi)，頻率w1和w3對(duì)應(yīng)的奇異值即為海雜波奇異值。

圖4 回波信號(hào)時(shí)頻分析Fig.4 Time-frequency analysis of received signals

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

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

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

圖6 抑制前后檢測(cè)結(jié)果Fig.6 Detection results before/after suppression

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

4 結(jié)語(yǔ)

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

［1］紀(jì)永剛，張杰，王祎鳴，等.雙頻地波雷達(dá)船只目標(biāo)點(diǎn)跡關(guān)聯(lián)與融合處理 ［J］.系統(tǒng)工程與電子技術(shù)，2014，36（2）：266-271.

［2］凡俊梅，焦培南，肖景明.海洋雜波對(duì)高頻雷達(dá)檢測(cè)海面上低速目標(biāo)的影響 ［J］.電波科學(xué)學(xué)報(bào)，1997，12（2）：205-210.

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

［4］楊煉，孫合敏，潘新龍.基于擴(kuò)展Prony算法的海雜波循環(huán)對(duì)消法［J］.現(xiàn)代雷達(dá)，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］周文瑜，焦培南.超視距雷達(dá)技術(shù)［M］.北京：電子工業(yè)出版社，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］董英凝，張寧，許榮慶.高頻地波雷達(dá)工作環(huán)境對(duì)系統(tǒng)性能影響的分析 ［J］.電波科學(xué)學(xué)報(bào)，2007，22（2）：325-330.

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

［17］仇永斌，張寧，張樹春.雙基地高頻地波雷達(dá)海雜波抑制［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，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檢測(cè)器性能分析 ［J］.空軍雷達(dá)學(xué)院學(xué)報(bào)，2009，23（5）：356-358.