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

（武漢大學(xué)電子信息學(xué)院，湖北武漢 430072）

0 引 言

高頻地波雷達(dá)（High Frequency Surface Wave Radar,HFSWR）工作在短波段（3～30 MHz），極易接收到來(lái)自各種短波電臺(tái)、通信設(shè)備的射頻干擾信號(hào)（Radio Frequency Interference,RFI），導(dǎo)致雷達(dá)數(shù)據(jù)質(zhì)量下降，嚴(yán)重影響高頻地波雷達(dá)目標(biāo)探測(cè)以及海態(tài)信息的提取。

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

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

1 RFI檢測(cè)

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

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

2 048 個(gè)采樣點(diǎn)的頻率分辨精度完全不夠，為了提高精度，需要將4 根天線(xiàn)的環(huán)境監(jiān)測(cè)采樣數(shù)據(jù)按照時(shí)間先后順序進(jìn)行拼接，在雷達(dá)相干積累時(shí)間5 min 內(nèi)有150 場(chǎng)采樣數(shù)據(jù)，共有150×2 048=307 200個(gè)采樣點(diǎn)，以25 600個(gè)采樣點(diǎn)為一段進(jìn)行拼接，共分為12段，此時(shí)頻率分辨率為

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

圖1 頻譜監(jiān)測(cè)譜圖

2 慢時(shí)域分段RFI抑制

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

2.1 訓(xùn)練張量構(gòu)造

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

圖2 訓(xùn)練張量構(gòu)造示意圖

2.2 RFI子空間估計(jì)

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

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

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

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

定義各展開(kāi)模式矩陣對(duì)應(yīng)左奇異矩陣列向量做FFT 得到左奇異頻域，對(duì)應(yīng)右奇異矩陣行向量做FFT 得到右奇異頻域，分析RFI和噪聲在不同模式矩陣左奇異頻域和右奇異頻域的表現(xiàn)形式，利用該頻率信息劃定干擾子空間和噪聲子空間。通過(guò)對(duì)譜峰頻率、譜峰功率以及功率均值這三個(gè)參

2.3 處理張量子空間投影

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

這里的模n乘積與傳統(tǒng)的HOSVD 中的模n乘積有所不同，由于P1和P3是利用右奇異矩陣構(gòu)造的投影矩陣，張量與其進(jìn)行模n乘積時(shí)應(yīng)使得其對(duì)應(yīng)模式展開(kāi)矩陣左乘以投影矩陣，以P1為例，有如下關(guān)系：

而對(duì)于P2而言，應(yīng)使得對(duì)應(yīng)模式展開(kāi)矩陣右乘以投影矩陣。

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

2.4 算法步驟總結(jié)

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

步驟1 分段：取第一次FFT 得到的數(shù)據(jù)，在慢時(shí)間維分段；

步驟2 構(gòu)造訓(xùn)練張量：利用不包含海洋回波的遠(yuǎn)距離元或負(fù)距離元數(shù)據(jù)構(gòu)造訓(xùn)練張量；

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

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

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

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

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

圖4 算法流程圖

3 實(shí)驗(yàn)結(jié)果分析

本節(jié)首先通過(guò)抑制仿真RFI檢驗(yàn)算法有效性，隨后利用2021年8月福建省東山、龍海高頻地波雷達(dá)站實(shí)測(cè)數(shù)據(jù)驗(yàn)證算法實(shí)用性，最后統(tǒng)計(jì)分析該算法對(duì)一天雙站雙頻數(shù)據(jù)的抑制效率與抑制效果。

3.1 仿真實(shí)驗(yàn)

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

表1 仿真實(shí)驗(yàn)雷達(dá)系統(tǒng)參數(shù)

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

圖5 原始數(shù)據(jù)與添加仿真RFI譜圖

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

圖6 歸一化奇異值分布圖

圖7 三種展開(kāi)模式矩陣各行、列FFT譜圖

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

表2 功率均值、譜峰功率統(tǒng)計(jì)表

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

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

3.2 實(shí)測(cè)實(shí)驗(yàn)

利用2021年8月東山、龍海雙站雙頻組網(wǎng)探測(cè)的實(shí)測(cè)數(shù)據(jù)驗(yàn)證本文RFI抑制算法的實(shí)用性，雷達(dá)系統(tǒng)參數(shù)如表3所示。

表3 實(shí)測(cè)實(shí)驗(yàn)雷達(dá)系統(tǒng)參數(shù)

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

圖9 環(huán)境噪聲頻譜圖

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

圖10 實(shí)測(cè)數(shù)據(jù)RP譜與RD譜

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

圖11 歸一化奇異值分布

圖12 三種展開(kāi)模式矩陣各行、列FFT譜圖

表4 功率均值、譜峰功率統(tǒng)計(jì)表

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

圖13 RFI抑制后RP譜與RD譜

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

3.3 統(tǒng)計(jì)分析

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

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

通過(guò)對(duì)干擾抑制后RP 譜的查看，得到干擾抑制后剩余RFI時(shí)間分布如圖16所示，RFI存在時(shí)間占比依次為0.21%，2.66%，0.21%，0%?？梢钥吹?，干擾并未完全消除，而且在原本頻譜監(jiān)測(cè)數(shù)據(jù)未檢測(cè)到RFI 的時(shí)間段也存在少量RFI，東山低頻數(shù)據(jù)仍有較長(zhǎng)一段時(shí)間被RFI影響，分析原因主要有以下兩個(gè)：其一，由于RFI 檢測(cè)信噪比閾值的設(shè)置略大導(dǎo)致漏警，使得部分時(shí)間段RFI 并未有效檢測(cè)；其二，RFI并非總是由一系列單頻信號(hào)構(gòu)成，它們可能經(jīng)過(guò)各種調(diào)制，導(dǎo)致距離相關(guān)性減弱，子空間投影的抑制方法未能奏效。不過(guò)總體來(lái)說(shuō)，大部分RFI被檢測(cè)并抑制，證明本文算法是有效的。

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

4 結(jié)束語(yǔ)

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