高頻地波雷達射頻干擾慢時域抑制方法

高玉斌，岳顯昌，周 慶，吳雄斌，張 蘭

（武漢大學電子信息學院，湖北武漢 430072）

0 引 言

高頻地波雷達（High Frequency Surface Wave Radar,HFSWR）工作在短波段（3～30 MHz），極易接收到來自各種短波電臺、通信設備的射頻干擾信號（Radio Frequency Interference,RFI），導致雷達數據質量下降，嚴重影響高頻地波雷達目標探測以及海態信息的提取。

避免RFI 的一種直接方法就是利用實時頻譜監測結果，找到一個不被RFI污染的干凈頻段作為雷達的工作頻段［1-2］。但是在實際應用中，幾乎很難找到一段有足夠帶寬（＞30 kHz）且無明顯干擾的頻段。RFI 抑制方法主要有基于壓縮感知的方法［3］、經驗模態分解方法［4-5］、分數階傅里葉變換［6-8］等。近年來，正交子空間投影技術被應用于RFI的抑制并取得了較好的效果。正交子空間投影的關鍵在于如何準確估計出相互正交的信號子空間和干擾子空間。文獻［9］討論了RFI在“五域六圖”的特性，子空間類算法利用RFI的強距離相關性和方向特性［10］對其進行抑制，可以在慢時域［11］、多普勒域［12］、空域［13-14］等單一域處理，也可以聯合多個域的特征抑制RFI［15-16］。慢時域即雷達重復周期（Period），而快時域則指一個雷達重復周期內接收數據的采樣點，快時間維做一次快速傅里葉變換（Fast Fourier Transform,FFT）解距離得到距離域，慢時間維做一次FFT 解速度得到多普勒域。文獻［15］提出的基于高階奇異值分解（Higher-Order Singular Value Decomposition,HOSVD）的多域聯合子空間投影RFI 抑制算法，聯合RFI 在多普勒域的強距離相關性和空域的方向特性，能夠得到較單一域處理更好的干擾抑制效果。

現有的子空間類方法存在以下三個問題，難以應用到雷達數據的實時批處理中。首先，這類算法需要得到雷達數據的距離-周期（Range-Period,RP）譜或距離-多普勒（Range-Doppler,RD）譜，然后根據譜圖特征劃定干擾范圍，再進行干擾抑制。應用到雷達數據實時批處理中時，一般通過滑窗的方法對每段數據進行無差別抗干擾處理，這樣會大大降低數據處理效率，而且對于無污染的數據，強行估計干擾子空間作投影，也會影響回波數據質量。其次，這類算法估計干擾子空間沒有一個判決標準，通常約定最大的一個或者兩個特征值對應的特征矢量構成干擾子空間，數據批處理中往往會遇到不同的干擾情況，固定特征值的選取不利于干擾的準確消除，需要有一個選取標準，在數據批處理時能夠根據不同數據特征自適應地估計干擾子空間。最后，這類算法都只關注左奇異矩陣列向量信息，忽略了右奇異矩陣行向量所包含的信息，限制了干擾子空間投影的方式，如文獻［15］提出的HOSVD 方法在多普勒域消除RFI，由于RFI 在不同多普勒單元沒有相關性，導致HOSVD 后的第一展開模式矩陣不能估計干擾子空間，利用右奇異矩陣信息，三種展開模式矩陣都可以估計干擾子空間，有利于提高干擾子空間估計準確性。本文主要針對以上三個問題，提出解決方案，形成可以在實測數據批量實時處理中使用的RFI抑制算法。

1 RFI檢測

批處理雷達數據時，增加干擾檢測的預處理流程，可以大大提高數據處理效率。常見的干擾檢測算法有基于功率大小的門限檢測［17］、恒虛警檢測器［18］以及一些圖像處理算法［19-20］，以上方法都需要先得到雷達數據RD 譜，之后通過圖像特征作進一步分析，不適合用在雷達數據實時批處理中。將雷達數據在慢時域分段后，利用頻譜監測數據對每一段數據干擾情況進行判定，給出干擾有無的判定標志位，便于后續干擾抑制工作的進行，簡化處理流程。

在2021年東山、龍海雙站雙頻組網探測實驗中，增加了頻譜監測模塊，采集了實時環境噪聲數據，實時頻譜監測數據在雷達工作頻段的頻譜特征能很好地反映雷達數據是否受RFI的影響，以該實驗項目為背景介紹基于頻譜監測數據的RFI 慢時域檢測方法。頻譜監測實現方式為在48 MS/s的采樣率下，在接收脈沖末尾連續采集2 048 點環境噪聲數據，采樣時長約為43 μs。對該2 048點數據進行打包上傳，在上位機實現FFT 得到實時的環境噪聲頻譜。頻譜監測每隔2 s 實現一次，每次只實現單個天線的高低頻環境的同時監測。每隔兩秒切換下一個監測通道，實現完整的4根天線環境監測周期為8 s。

2 048 個采樣點的頻率分辨精度完全不夠，為了提高精度，需要將4 根天線的環境監測采樣數據按照時間先后順序進行拼接，在雷達相干積累時間5 min 內有150 場采樣數據，共有150×2 048=307 200個采樣點，以25 600個采樣點為一段進行拼接，共分為12段，此時頻率分辨率為

掃頻周期250 ms，5 min回波數據包含1 200個掃頻周期，將100 個掃頻周期劃為一段，即每段時長為25 s，在慢時域上將數據分為12 段，每段數據按照時間先后依次對應到12 段頻譜監測數據，通過對每段頻譜監測數據的分析判斷該段時間內是否存在RFI，對有RFI 影響的數據段進行定位，之后在慢時域上利用HOSVD 方法進行干擾的消除。圖1展示了雷達高頻范圍（11.5～13.5 MHz）與低頻范圍（7～9 MHz）的頻譜監測譜圖，將距離這兩個頻段都較遠的頻段（20～24 MHz）作為背景噪聲頻段，定義為底噪窗，檢測窗的頻率范圍設置取決于雷達工作頻率、掃頻帶寬以及掃頻模式，經過混頻、低通濾波后，只有載頻在雷達發射信號頻帶范圍內的干擾信號方能影響回波數據。一般存在RFI時，檢測窗內會有明顯尖峰，通過檢測窗峰值功率與底噪窗平均功率計算檢測信噪比，大于設定閾值時即認為該時間段存在RFI。

圖1 頻譜監測譜圖

2 慢時域分段RFI抑制

雷達原始采樣數據做第一次FFT 變換，解距離得到通道-距離-慢時間維數據，接著利用HOSVD 方法抑制RFI。干擾抑制前需要先將雷達數據在慢時域上分段，來降低HOSVD 算法的運算量，提高干擾抑制效率。與第1 節中頻譜監測數據分段檢測RFI相對應，將一個相干積累時間內雷達數據在慢時域分為12 段，即每段數據25 s。HOSVD通過矩陣展開的方式建立了張量和矩陣之間的聯系，定義了一種張量分解方法［21］，并已應用于圖像處理［22］、人臉識別［23］、雷達信號處理［24］等多個領域。

2.1 訓練張量構造

RFI分布于所有距離元，而海洋回波僅存在于近距離元，因此構造訓練張量時盡量選取不包含海洋回波的遠距離元或負距離元，有利于提高干擾子空間估計的準確性。假設在慢時間維上將每Q個掃頻周期分割為一段，每一段選取M個通道、P個遠距離元或者負距離元數據構造一個大小為P×Q×M的張量，稱之為訓練張量A，如圖2所示。

圖2 訓練張量構造示意圖

2.2 RFI子空間估計

將訓練張量A展開，得到三種展開模式矩陣如圖3所示，其中陰影部分代表RFI。分析三種模式展開矩陣的形式，利用其得到的左、右奇異矩陣包含的頻率信息估計RFI 子空間，其中左奇異矩陣列向量表征模式矩陣中各分量列與列之間的相關性，而右奇異矩陣行向量表征模式矩陣中各分量行與行之間的相關性，左、右奇異向量對應的奇異值矩陣中奇異值的大小表征了模式矩陣各分量列與列或行與行之間相關性的強弱。數的判定確定d1，d2，d3的取值，提高RFI 子空間估計的準確性，以免d選取偏小導致RFI 不能被完全消除或者d選取偏大將部分噪聲分量劃入RFI子空間，投影到其他距離元時影響有用信號分量。

對于第一展開模式矩陣A(1)∈CP×QM，每一列表示同一接收通道、同一掃頻周期所有距離元數據，每一行表示一個距離元內所有數據的拼接。由于RFI在距離維的強相關性，可以利用右奇異矩陣行向量來估計RFI 子空間，得到的投影矩陣為P1=，其中V1表示A(1)奇異值分解得到的右奇異矩陣，d1表示RFI分量對應右奇異矩陣的行數。

對于第二展開模式矩陣A(2)∈CQ×MP，每一列代表同一距離元、同一接收通道所有掃頻周期的數據，每一行表示一個掃頻周期內所有數據的拼接。只能通過A(2)對應的左奇異矩陣列向量來估計RFI 子空間，得到的投影矩陣為P2=，其中U2表示A(2)奇異值分解得到的左奇異矩陣，d2表示RFI分量對應左奇異矩陣的列數。

對于第三展開模式矩陣A(3)∈CM×PQ，每一列表示同一掃頻周期、同一距離元所有接收通道的數據，每一行表示一個接收通道內所有數據。左奇異矩陣列向量和右奇異矩陣行向量都可以用來估計RFI 子空間，雖然左、右奇異矩陣的自由度是相同的，都取決于奇異值矩陣的秩，但是左奇異矩陣的維度取決于接收通道的數目M，列向量點數過少，FFT 之后頻率分辨率過低，無法反映真實的頻率信息，相反右奇異矩陣的維度取決于訓練集所選距離元與每段數據選取的掃頻周期數的乘積，即P×Q，行向量FFT 之后可以準確反映相應分量的多普勒頻率信息，于是選擇右奇異矩陣行向量估計RFI 子空間，得到的投影矩陣為P3=，其中V3表示A(3)奇異值分解得到的右奇異矩陣，d3表示RFI分量對應右奇異矩陣的行數。

定義各展開模式矩陣對應左奇異矩陣列向量做FFT 得到左奇異頻域，對應右奇異矩陣行向量做FFT 得到右奇異頻域，分析RFI和噪聲在不同模式矩陣左奇異頻域和右奇異頻域的表現形式，利用該頻率信息劃定干擾子空間和噪聲子空間。通過對譜峰頻率、譜峰功率以及功率均值這三個參

2.3 處理張量子空間投影

估計出RFI子空間后，將每一段處理數據劃分成與訓練張量大小相等的若干處理張量D，將干擾子空間在處理張量上投影，得到處理張量中的RFI分量Dr：

這里的模n乘積與傳統的HOSVD 中的模n乘積有所不同，由于P1和P3是利用右奇異矩陣構造的投影矩陣，張量與其進行模n乘積時應使得其對應模式展開矩陣左乘以投影矩陣，以P1為例，有如下關系：

而對于P2而言，應使得對應模式展開矩陣右乘以投影矩陣。

將處理張量中的RFI 分量減去，就得到消除RFI后其余信號分量Ds：

2.4 算法步驟總結

基于HOSVD 的慢時域分段RFI抑制方法的具體步驟如下：

步驟1 分段：取第一次FFT 得到的數據，在慢時間維分段；

步驟2 構造訓練張量：利用不包含海洋回波的遠距離元或負距離元數據構造訓練張量；

步驟3 估計RFI 子空間：分別利用各展開模式矩陣奇異值分解后左、右奇異矩陣包含的頻率信息估計RFI子空間；

步驟4 處理張量子空間投影：將得到的干擾子空間投影到各處理張量D，得到只包含RFI的干擾張量Dr；

步驟5 消除RFI：利用處理張量D減去干擾張量Dr，得到抑制干擾后的數據張量Ds；

步驟6 重復步驟2～6，直到每一段存在RFI的數據都處理完成。

RFI抑制算法流程圖如圖4所示。

圖4 算法流程圖

3 實驗結果分析

本節首先通過抑制仿真RFI檢驗算法有效性，隨后利用2021年8月福建省東山、龍海高頻地波雷達站實測數據驗證算法實用性，最后統計分析該算法對一天雙站雙頻數據的抑制效率與抑制效果。

3.1 仿真實驗

選取一場不存在明顯干擾的雷達接收數據，這里以2017年1月5日00:20赤湖高頻地波雷達站接收數據為例，雷達參數如表1所示。

表1 仿真實驗雷達系統參數

實測數據距離-周期（Range-Period,RP）譜和距離-多普勒（Range-Doppler,RD）譜分別如圖5（a）、（b）所示，在該實測數據301 到400 掃頻周期內添加兩個射頻干擾，載波頻率分別為13 158 023.52 Hz、13 167 512.36 Hz，到達角相同，均為130°。繪制出添加仿真RFI 后的RP 譜和RD 譜分別如圖5（c）、（d）所示，由于載波頻率的不同，兩個射頻干擾分別位于RD 譜-0.48 Hz 與0.35 Hz 處，其中位于0.35 Hz處的RFI遮蓋了正一階Bragg峰區域。

圖5 原始數據與添加仿真RFI譜圖

利用本文方法抑制RFI，首先將每100 個掃頻周期分為一段，由于仿真RFI 僅存在于301 至400掃頻周期內，所以只需要處理該段數據來驗證算法有效性。該段數據為第301 到400 掃頻周期內所有距離元和接收通道的數據，選取遠距離元數據構造訓練張量，這里497 對應0 距離元，距離元是反的，選457 到466，共10 個距離元的數據，將訓練張量得到的三種展開模式矩陣依次進行SVD，歸一化奇異值分布如圖6所示，設置檢測閾值為0.1，認為歸一化奇異值大于0.1 則為大奇異值，從圖6中可以看到三種展開模式矩陣分別有2，4，3個大奇異值，依次將第一展開模式矩陣對應右奇異矩陣前兩行、第二展開模式矩陣對應左奇異矩陣前四列、第三展開模式矩陣對應右奇異矩陣前三行做FFT，如圖7所示，可以看到第一展開模式矩陣對應右奇異矩陣第一行、第二展開模式矩陣對應左奇異矩陣第一列、第三展開模式矩陣對應右奇異矩陣第一行FFT 譜圖中均有兩個明顯譜峰，譜峰頻率分別為-0.48 Hz 和0.35 Hz，與仿真RFI多普勒頻率相同。

圖6 歸一化奇異值分布圖

圖7 三種展開模式矩陣各行、列FFT譜圖

表2列出了各行各列FFT 譜圖功率均值、譜峰功率，從表中可以讀到各模式矩陣對應奇異矩陣第一行或第一列功率均值依次為-14.794 6，-12.629 9 和-16.002 4 dB，明顯低于0 dB，通過對大量數據的驗證分析，得到如下結論：若該行向量或列向量對應噪聲子空間，做FFT 之后，頻譜功率均值在0 dB 上下浮動，浮動范圍一般不超過5 dB。因此可以用譜峰功率、功率均值兩個指標來判斷該行向量或列向量對應噪聲子空間還是干擾子空間，而且可以通過譜峰對應頻率確定RFI 在RD 譜上的多普勒頻率位置。綜上所述，第一展開模式矩陣對應右奇異矩陣第一行、第二展開模式矩陣對應左奇異矩陣第一列以及第三展開模式矩陣對應右奇異矩陣第一行對應到仿真的兩個RFI 分量子空間，其余行、列向量對應噪聲子空間，由此確定第一、第二、第三展開模式矩陣估計的RFI 子空間投影矩陣的參數分別為d1=1，d2=1，d3=1。

表2 功率均值、譜峰功率統計表

估計出RFI子空間后，投影到各處理張量并消除RFI，畫出RP 譜和RD 譜分別如圖8（a）、（b）所示，從RP 譜中可以看到，在301 到400 掃頻周期添加的RFI 被完全消除，而且在RD 譜上也完全看不到RFI，正一階Bragg峰凸顯出來。

圖8 干擾抑制后RP譜與RD譜

3.2 實測實驗

利用2021年8月東山、龍海雙站雙頻組網探測的實測數據驗證本文RFI抑制算法的實用性，雷達系統參數如表3所示。

表3 實測實驗雷達系統參數

選取8月11日東山高頻07:40 的數據，該時間段數據被RFI 污染。將該時段頻譜監測數據拼接后均分為12 段，對每段數據做FFT 得到不同分段環境噪聲頻譜圖。根據表3可知,雷達工作頻率為12.500 MHz，掃頻周期為30 kHz，掃頻模式為上掃，因此檢測窗范圍設為12.500 MHz 至12.530 MHz，底噪窗范圍為20～24 MHz。求出檢測窗峰值功率與底噪窗功率均值做差，設定檢測閾值為20 dB。得到12 段RFI 標志位，其中第1～3、8～12 段數據標志位為1，即1 至300、701 至1 200 掃頻周期內的數據都被RFI 污染。畫出第2、第4 段環境噪聲頻譜圖如圖9所示，紅色虛線框表示檢測窗范圍，二者進行對比，可以看到第2段噪聲頻譜圖在檢測窗內有明顯的尖峰，而第4段噪聲頻譜圖在檢測窗內并無明顯尖峰。

圖9 環境噪聲頻譜圖

實測數據RP 譜和RD 譜分別如圖10（a）、（b）所示，可以看到，RP 譜中被RFI污染的掃頻周期分段與頻譜監測數據得到的檢測結果一致，在RD 譜上，RFI 呈現為一帶狀豎條紋，多普勒頻率大致為0.52 Hz，有一定程度的展寬，覆蓋了龍海發東山收得到海洋回波譜的負一階峰區域。本次實測雷達數據-80 至0 距離元是0 至80 距離元的搬移，即負距離元不再只包含RFI，而且由于是雙站組網，一場數據中包含了東山、龍海雙站的回波，圖10（b）中0 至15 距離元的Bragg 峰是東山站自發自收的回波，40 至55 距離元是龍海發東山收得到的海洋回波，而且由于高低頻發射都各有兩路，因此通過加多普勒偏置來區分兩路信號回波，導致兩路信號產生的海洋回波在RD 譜對稱分布。海洋回波占據了多數距離元，為保證RFI 抑制效果，需要自適應地選取不包含海洋回波的距離元進行訓練。設定訓練范圍是10 個距離元，首先通過參數配置文件讀取雙站雙頻的距離偏置，根據距離偏置確定雙站海洋回波所占據的距離元大致范圍，在該范圍外取10個距離元進行訓練。

圖10 實測數據RP譜與RD譜

以第二段數據，即101 至200 掃頻周期的數據為例進行分析，訓練距離元選取61 至70，得到訓練張量各展開模式矩陣歸一化奇異值分布如圖11所示，可以得到三種展開模式矩陣對應大奇異值個數分別為3，3，2，得到各行、列FFT 譜圖如圖12所示，功率均值、譜峰功率如表4所示，分析得出第一展開模式矩陣對應右奇異矩陣第一行、第二展開模式矩陣對應左奇異矩陣第一列以及第三展開模式矩陣對應右奇異矩陣第一行對應RFI子空間，從第三展開模式矩陣前兩行FFT 譜圖中可以看到有兩個明顯譜峰，其中一個譜峰頻率在0.52 Hz 左右，對應RD 譜中RFI 分量，而另一個譜峰頻率在0 Hz 左右，對應RD 譜中的零頻干擾（Zero Frequency Interference,ZFI）分量，可以將零頻干擾一并消除。確定第一、第二、第三展開模式矩陣估計的RFI 子空間投影矩陣的參數分別為d1=1，d2=1，d3=2，將干擾子空間投影到待處理張量，完成該段數據RFI 的抑制。

圖11 歸一化奇異值分布

圖12 三種展開模式矩陣各行、列FFT譜圖

表4 功率均值、譜峰功率統計表

根據RFI標志位依次處理各段數據，得到干擾抑制后的RP 譜和RD 譜分別如圖13（a）、（b）所示。從圖中可以看到，RP 譜中每個掃頻周期的RFI 都被抑制掉了，只有少量信噪比較低的殘余量，從RD 譜中可以看到，在第二次FFT 之后，這些殘余量基本降為底噪水平，RFI 被完全消除，被掩蓋的一階Bragg 峰凸顯出來。圖14（a）給出了RFI抑制前后第800 掃頻周期的距離譜對比圖，可以看到，RFI 被抑制，海洋回波并未受到影響，且信噪比抬升10 dB左右。圖14（b）為RFI抑制前后第42距離元多普勒譜對比圖，可以看到多普勒頻率位于0.52 Hz左右的射頻干擾以及0 Hz處的零頻干擾都被完全消除，海洋回波完整保留，Bragg峰的信噪比得到提升。

圖13 RFI抑制后RP譜與RD譜

圖14 RFI抑制前后距離譜與多普勒譜

3.3 統計分析

利用頻譜監測數據得到的RFI 標志位畫出2021年8月11日東山、龍海雙站雙頻回波數據被RFI 污染情況時間分布，如圖15所示，其中東山高頻RFI 存在時間占比12.93%，東山低頻RFI存在時間占比22.94%，龍海高頻RFI 存在時間占比53.68%，龍海低頻RFI 存在時間占比0.07%。總體來看，龍海低頻數據質量較好，其他數據受RFI影響較大。利用本文算法對2021年8月11日一整天雙站雙頻的數據進行無差別抑制，耗時3 494.531 373 s，加入RFI 檢測算法后，耗時1 670.540 791 s。測試所用電腦GPU型號：Intel（R）Core（TM）i5-9400 CPU@2.90 GHz。在數據批處理中，加入RFI 檢測算法可以有效降低干擾消除時間，提高數據批處理效率。

圖15 雙頻雙站RFI污染情況分布圖

通過對干擾抑制后RP 譜的查看，得到干擾抑制后剩余RFI時間分布如圖16所示，RFI存在時間占比依次為0.21%，2.66%，0.21%，0%。可以看到，干擾并未完全消除，而且在原本頻譜監測數據未檢測到RFI 的時間段也存在少量RFI，東山低頻數據仍有較長一段時間被RFI影響，分析原因主要有以下兩個：其一，由于RFI 檢測信噪比閾值的設置略大導致漏警，使得部分時間段RFI 并未有效檢測；其二，RFI并非總是由一系列單頻信號構成，它們可能經過各種調制，導致距離相關性減弱，子空間投影的抑制方法未能奏效。不過總體來說，大部分RFI被檢測并抑制，證明本文算法是有效的。

圖16 干擾抑制后雙頻雙站RFI污染情況分布圖

4 結束語

本文針對高頻地波雷達RFI抑制，提出了一種慢時域分段檢測與抑制方法，該方法解決了傳統子空間類方法無RFI檢測、無干擾子空間嚴格判決標準、未關注右奇異矩陣行向量信息的問題。在預處理階段，利用頻譜監測數據對慢時域分段數據進行RFI 檢測。在后處理階段，利用HOSVD 方法對雷達接收數據進行第一次FFT后得到的通道-距離-慢時間維數據進行處理，首先將待處理數據在慢時間維上分段，對于存在RFI 的數據段，自適應地選擇合適距離元數據構造訓練張量，通過分析訓練張量三種展開模式矩陣的形式，結合各展開模式矩陣奇異值分解后得到的左奇異矩陣列向量或右奇異矩陣行向量所包含的頻率信息，給出了干擾子空間估計的參數確定方法，準確估計出干擾子空間，進而得到信號子空間，在慢時域上完成RFI的抑制。實驗結果表明，該方法可以有效檢測并抑制RFI，提高了數據批處理的效率。