基于壓縮感知的時域抗射頻干擾方法

李宇環(huán), 岳顯昌, 張 蘭

(武漢大學(xué)電子信息學(xué)院，武漢 430072)

高頻地波雷達(dá)(high frequency surface wave radar, HFSWR)利用垂直極化的短波沿海洋表面繞射傳播衰減很小的特點(diǎn)，能實(shí)現(xiàn)對海洋表面動力學(xué)要素和海面運(yùn)動目標(biāo)的超視距探測[1-2]。在其工作的短波頻段(3～30 MHz)，用于各種業(yè)務(wù)的短波通信、廣播電臺等頻率占用率非常高，在自由空間存在著大量的短波無線電信號[3]。在合適的電磁波傳播環(huán)境下，來自于外部無線電設(shè)備的位于HFSWR系統(tǒng)工作帶寬內(nèi)的無線電波進(jìn)入雷達(dá)接收機(jī)，形成射頻干擾(radio frequency interference, RFI)，會極大地限制高頻地波雷達(dá)系統(tǒng)的探測性能。因此，抗射頻干擾技術(shù)的研究一直是高頻地波雷達(dá)應(yīng)用中的一個備受關(guān)注的課題。

傳統(tǒng)的射頻干擾抑制方法，例如自適應(yīng)波束形成[4]，極化濾波[5]和旁瓣對消[6]等方法，在窄波束HFSWR中取得了很好的應(yīng)用效果，然而由于角度分辨率的限制，該類方法在寬波束HFSWR中應(yīng)用效果無法令人滿意。中外學(xué)者通過分析射頻干擾在信號處理各個階段的表現(xiàn)形式和特點(diǎn)，對寬波束HFSWR中射頻干擾的抑制展開深入研究并提出了相應(yīng)算法。鑒于射頻干擾的強(qiáng)距離相關(guān)性和方向性，正交子空間投影技術(shù)[7-9]被應(yīng)用于射頻干擾的抑制并取得了良好的效果。“時域剔除類”算法[10-12]基于射頻干擾在原始時域上的短時性特征，對干擾數(shù)據(jù)段依次進(jìn)行置零和信號恢復(fù)處理。

壓縮感知(compressive sensing， CS)理論[13-14]是由Donoho、Candes及華裔科學(xué)家Tao等提出的一種新的信息獲取指導(dǎo)理論。該理論指出：對于稀疏或可壓縮(在某個表征域近似稀疏)的信號，采用遠(yuǎn)低于奈奎斯特標(biāo)準(zhǔn)的采樣率進(jìn)行采樣，就可以實(shí)現(xiàn)對信號的精確重構(gòu)。直至目前，該理論在圖像處理、醫(yī)療成像、模式識別、雷達(dá)信號處理等領(lǐng)域得到了廣泛應(yīng)用[15-17]。

在HFSWR回波數(shù)據(jù)中，目標(biāo)信號在距離域和多普勒域上是稀疏的，這為壓縮感知在射頻干擾抑制上的應(yīng)用提供了機(jī)會。文獻(xiàn)[18]分析了射頻干擾在時域上的短時性特征和艦船目標(biāo)在距離多普勒(range-Doppler, RD)譜上的稀疏性，提出了基于壓縮感知的射頻干擾抑制方法?；趬嚎s感知理論，該方法利用少量的無干擾時域數(shù)據(jù)即可完成對艦船目標(biāo)信號的距離信息和速度信息的提取，從而避免射頻干擾的影響。對于海態(tài)探測高頻地波雷達(dá)，海洋回波一般也只存在于距離譜的少部分距離元，然而由于其在距離域和多普勒域上分布的連續(xù)性，該方法在HFSWR海態(tài)探測中難以保證參數(shù)反演的精度和準(zhǔn)確度。

基于射頻干擾的時域短時性特征，提出一種新的壓縮感知射頻干擾抑制方法。該方法在原始時域數(shù)據(jù)中對射頻干擾位置進(jìn)行閾值檢測，通過對其他無干擾數(shù)據(jù)的隨機(jī)欠采樣，基于壓縮感知理論恢復(fù)被干擾的海洋回波時間采樣，以達(dá)到射頻干擾抑制的目的。為了驗(yàn)證算法的有效性和對海洋回波的影響，基于武漢大學(xué)OSMAR071高頻地波雷達(dá)實(shí)測數(shù)據(jù)，分別給出仿真和實(shí)測射頻干擾的處理結(jié)果。

1 射頻干擾的短時性特征

在高頻地波雷達(dá)系統(tǒng)中，短波、廣播通信等窄帶信號是射頻干擾的主要來源，其中影響雷達(dá)性能的主要是采用幅度(AM)調(diào)制的射頻干擾，且載波信號的頻率和能量一般遠(yuǎn)高于調(diào)制信號，因此總是忽略射頻干擾的調(diào)制部分而將其簡化為單頻信號進(jìn)行分析[19]。

以采用線性調(diào)頻中斷連續(xù)波(frequency modulate interrupted continuous wave, FMICW)體制的高頻地波雷達(dá)為例，單個掃頻周期內(nèi)雷達(dá)發(fā)射信號ST(t)可以表示為線性調(diào)頻信號：

(1)

式(1)中：f0是雷達(dá)系統(tǒng)的工作頻率；k是調(diào)頻斜率；φT是初相位；aT是信號幅度；Ts是掃頻周期。雷達(dá)接收機(jī)的本振信號與發(fā)射信號相同。

在第p個掃頻周期，一個單頻干擾信號可以寫成：

SI(p,t)=aIcos[2πfIt+φI(pTs)], 0≤t≤Ts

(2)

式(2)中：aI為射頻干擾信號的幅度；fI是干擾信號的載頻。

φI(pTs)=φI0+2πfI(p-1)Ts

(3)

式(3)中：φI0和φI(pTs)分別表示干擾在第1個和第p個掃頻周期的初相位。

圖1 RFI信號在接收機(jī)中的時頻示意圖Fig.1 The time-frequency diagram of RFI signal in the receiver

圖2 雷達(dá)原始數(shù)據(jù)Fig.2 Radar raw data

圖1為射頻干擾在接收機(jī)中的時頻示意圖。單頻射頻干擾信號SI進(jìn)入雷達(dá)接收機(jī)后，與本振信號ST進(jìn)行混頻，再經(jīng)過截止頻率為fL的低通濾波處理后，在雷達(dá)接收數(shù)據(jù)中第p個掃頻周期的表現(xiàn)形式為

φI(pTs)-φT)],t1≤t≤t2

(4)

式(4)中：fΔ=fI-f0。如圖1所示，由于fL一般遠(yuǎn)小于發(fā)射信號的掃頻帶寬B，低通濾波后干擾只存在于t1和t2時刻之間，由圖1中幾何關(guān)系可以求得：

(5)

也就是說，在雷達(dá)原始時域數(shù)據(jù)中，干擾只存在于每個掃頻周期的少數(shù)采樣點(diǎn)，當(dāng)干擾載頻在相干積累時間內(nèi)保持穩(wěn)定，則干擾數(shù)據(jù)位置也保持穩(wěn)定，且干擾瞬時功率一般遠(yuǎn)高于海洋回波，即射頻干擾在時域上的短時性特征，如圖2所示，海洋回波平穩(wěn)存在于雷達(dá)接收數(shù)據(jù)的每個采樣點(diǎn)，而射頻干擾只存在于陰影部分所示區(qū)域。

2 基于壓縮感知的射頻干擾抑制

2.1 壓縮感知理論簡介

壓縮感知理論主要包括3個部分：信號的稀疏表示，觀測矩陣的構(gòu)造和稀疏重構(gòu)算法。

對于一個自然信號的N點(diǎn)采樣序列x=[x1,x2, …,xN]T，找到一組合適的基向量{ψi∈CN×1, (i=1, 2, …,N)}對其進(jìn)行線性表示：

(6)

式(6)中，Ψ∈CN×N為稀疏基矩陣；s∈CN×1為系數(shù)向量；x和s是同一個信號在不同表征域下的等價體現(xiàn)。若向量s中非零元素或者大系數(shù)的個數(shù)K滿足K?N，則信號x是稀疏的，向量s是信號x的稀疏表示，K為信號x的稀疏度。

尋找一個合適的觀測矩陣Φ∈CM×N(K

y=Φx=ΦΨs

(7)

式(7)中，y∈CM×l是觀測信號，令D=ΦΨ，則得到：

y=Dx

(8)

式(8)為壓縮感知方程，其中D∈CM×N為傳感矩陣?？芍?8)中方程個數(shù)遠(yuǎn)小于未知數(shù)個數(shù)，不存在唯一確定解。但由于信號是K稀疏的，當(dāng)傳感矩陣D滿足有限等距性質(zhì)(restricted isometry property, RIP)[20]時，通過稀疏重構(gòu)算法[15]可唯一確定稀疏解s，根據(jù)式(6)則可精確重構(gòu)信號x。

2.2 射頻干擾抑制

使用的重構(gòu)算法是簡化后的正交匹配追蹤(orthogonal matching pursuit, OMP)算法。OMP算法的本質(zhì)思想是以貪婪迭代的方式選擇傳感矩陣中與當(dāng)前觀測信號相關(guān)性最大的列，通過最小二乘法估計(jì)觀測信號中已選擇的列集合的相關(guān)分量并減去，反復(fù)迭代，直到達(dá)到迭代終止條件。該算法主要包含兩個步驟：選擇最大相關(guān)列、最小二乘估計(jì)。根據(jù)高頻地波雷達(dá)信號處理特點(diǎn)，海洋回波在距離譜上的位置，即傳感矩陣中海洋回波對應(yīng)的列集合是可以先驗(yàn)確定的。

本文提出的基于壓縮感知的射頻干擾抑制方法根據(jù)雷達(dá)系統(tǒng)參數(shù)得到海洋回波在距離譜上的位置，直接利用最小二乘估計(jì)即可得到海洋回波保持的系數(shù)向量估計(jì)，減小了迭代帶來的計(jì)算量和相關(guān)性選擇帶來的誤差，增加了重構(gòu)準(zhǔn)確度。

基于壓縮感知的射頻干擾抑制方法具體步驟如下。

步驟1對每個掃頻周期的回波數(shù)據(jù)進(jìn)行求和，得到時間采樣和序列z∈CN×l：

(9)

式(9)中：xp∈CN×l表示第p個掃頻周期的接收數(shù)據(jù)。

步驟2設(shè)置閾值：T=d×med，其中d為經(jīng)驗(yàn)因子，一般取1.1～1.35，med為和序列z的中位數(shù)。

步驟3判定和序列z中幅值超過閾值T的數(shù)據(jù)存在射頻干擾，記其對應(yīng)的索引集合為Id。

步驟4選擇傅里葉變換基作為稀疏基矩陣Ψ，確定海洋回波對應(yīng)的列索引集合為Λ，集合中的元素總數(shù)就是稀疏度K，初始化p=1。

步驟5選擇回波數(shù)據(jù)xp，構(gòu)造觀測矩陣Φ為隨機(jī)欠采樣矩陣：

(10)

(11)

式(11)中，DΛ表示傳感矩陣D的第Λ列。

(12)

步驟8p=p+1，重復(fù)步驟5～步驟7，直到所有的掃頻周期處理結(jié)束。

其中，步驟1～步驟3完成射頻干擾定位，步驟4～步驟8完成射頻干擾抑制。

3 仿真和實(shí)測分析

利用布設(shè)在福建省海洋預(yù)報(bào)臺東山地波雷達(dá)站的OSMAR071陣列式高頻地波雷達(dá)實(shí)測數(shù)據(jù)，分別應(yīng)用本文算法對仿真的射頻干擾和實(shí)測的射頻干擾進(jìn)行抗干擾處理。

3.1 仿真射頻干擾處理

以2017年4月4日12:11東山地波雷達(dá)接收站接收數(shù)據(jù)為例，其原始時域數(shù)據(jù)和對應(yīng)的RD譜如圖3所示，該場數(shù)據(jù)無明顯的射頻干擾存在。相關(guān)雷達(dá)系統(tǒng)參數(shù)為：工作頻率f0=13.15 MHz，掃頻周期Ts=0.5 s，掃頻斜率k=60 kHz/s，相干積累時間為5 min，每個掃頻周期采樣點(diǎn)數(shù)N=994。

在此數(shù)據(jù)上，根據(jù)式(2)給出的單頻干擾模型，人為加入兩個仿真的射頻干擾信號，相關(guān)仿真參數(shù)分別為：兩個干擾載頻分別為fI1=13.167 5 MHz，fI2=13.158 0 MHz，低通濾波器截止頻率為fL=150 Hz，持續(xù)時間為200個掃頻周期，持續(xù)時間內(nèi)，每個掃頻周期受到干擾的時間采樣單元有212個，且設(shè)置兩個仿真射頻干擾信號的DOA為40°。添加射頻干擾后的原始時域數(shù)據(jù)和對應(yīng)RD譜如圖4所示。由圖4(b)可以看出，在RD譜上，正一階Bragg峰對應(yīng)的海洋回波已經(jīng)被仿真的射頻干擾完全覆蓋，對比圖3(b)，噪聲基底的色階也有顯著提高。

圖3 OSMAR071高頻雷達(dá)實(shí)測數(shù)據(jù)Fig.3 OSMAR071 HF radar measured data

圖4 實(shí)測數(shù)據(jù)添加干擾后Fig.4 The measured data with injected RFI

采用本文提出的基于壓縮感知的時域抗干擾方法，對圖4(a)中干擾存在的200個掃頻周期進(jìn)行處理，結(jié)果如圖5所示。時域數(shù)據(jù)中的短時性射頻干擾成分被完全去除，RD譜上射頻干擾對應(yīng)的豎條狀成分也不再顯現(xiàn)，且正一階Bragg峰從雜波背景凸顯出來?？芍獔D5和圖3越接近，表明抑制效果越好。為了更清晰地展現(xiàn)本方法的抑制效果，將圖5(b)分別和圖4(b)、圖3(b)的第6個距離單元的多普勒譜進(jìn)行對比，如圖6所示。圖6(a)對比了第6個距離單元的干擾抑制前后的多普勒譜，可以看出經(jīng)過抗干擾處理后，仿真的RFI被抑制約25 dB，同時射頻干擾旁瓣所帶來的噪聲也得到明顯改善；圖6(b)對比第6個距離單元干擾抑制后和干擾添加前的多普勒譜，可以看到正、負(fù)Bragg峰區(qū)域的海洋回波吻合程度非常高，表明海洋回波信號的強(qiáng)度未被衰減。

圖5 實(shí)測數(shù)據(jù)仿真干擾抑制后Fig.5 The measured data after the artificial RFI suppressed

圖6 第6個距離元的多普勒譜Fig.6 Doppler spectrums at the 6th range bin

為了驗(yàn)證該抗射頻干擾算法對回波譜相位(方位信息)的影響，選取第6個距離單元正Bragg峰區(qū)域的譜點(diǎn)，利用MUSIC算法[21]對其進(jìn)行DOA估計(jì)。MUSIC算法對距離譜進(jìn)行距離滑窗得到多快拍數(shù)據(jù)，設(shè)置窗長為距離譜總長度的一半，得到快拍數(shù)為40。

圖7為滑窗后得到的多普勒譜，其中正Bragg峰紅色標(biāo)記的區(qū)域?yàn)檫x擇的11個多普勒譜點(diǎn)。對原數(shù)據(jù)、添加仿真干擾后、干擾抑制后多普勒譜的這些譜點(diǎn)進(jìn)行DOA估計(jì)，結(jié)果如圖8所示。對比原數(shù)據(jù)和干擾抑制前(添加仿真射頻干擾后)的DOA結(jié)果可以看出，添加的射頻干擾完全掩蓋了原有海洋回波的方位信息；對比原數(shù)據(jù)和干擾抑制后的DOA估計(jì)結(jié)果可知，該抗干擾算法可以很大程度地保留海洋回波的DOA信息。但由于第1和第11個譜點(diǎn)信噪比過低，干擾抑制后其DOA信息受到較大影響。

圖7 第6個距離元滑窗后的多普勒譜Fig.7 Doppler spectrums at the 6th range bin after sliding window

圖8 選取的11個譜點(diǎn)的DOA估計(jì)結(jié)果Fig.8 DOA estimation results of the selected 11 spectrum points

3.2 實(shí)測數(shù)據(jù)處理

對2017年3月31日01：03船收東山的實(shí)測數(shù)據(jù)進(jìn)行抗干擾處理，其原始時域數(shù)據(jù)和RD譜如圖9所示?？梢钥闯?，在時域數(shù)據(jù)中存在有明顯的短時性射頻干擾，導(dǎo)致在RD譜上出現(xiàn)多個豎條紋狀的射頻干擾，正一階Bragg峰區(qū)域也被完全淹沒，對海洋回波的提取造成極大困難。通過閾值檢測，每個掃頻周期存在46個時間采樣單元被污染，采用本文提出的基于壓縮感知的時域抗干擾方法對其進(jìn)行抑制，干擾抑制之后的原始時域數(shù)據(jù)和RD譜如圖10所示。可以看出，在時域數(shù)據(jù)中，存在于每個掃頻周期的短時性射頻干擾成分已經(jīng)被去除，而在每個掃頻周期的第160個時間采樣單元附近仍存在較弱的短時性成分未被消除，但是在RD譜上，明顯的射頻干擾已經(jīng)被完全抑制，這是因?yàn)樵摱虝r性成分功率較低，對RD譜的影響很小。而且干擾抑制后，正Bragg峰區(qū)域的海洋回波也從RD譜中凸顯出來，背景噪聲也得到一定改善。圖11給出了第4個距離單元抑制RFI前后的多普勒譜對比，可以看出，射頻干擾被抑制了大約10 dB，受干擾影響較小的負(fù)Bragg峰區(qū)域海洋回波在干擾抑制前后十分吻合，表明海洋回波得到完整的保留。

圖9 2017年3月31日01：03分船收東山的實(shí)測數(shù)據(jù)Fig.9 The measured data of ship from Dongshan at 01:03 am on Mar 31, 2017

圖10 實(shí)測數(shù)據(jù)干擾抑制后Fig.10 The measured data after RFI suppressed

圖11 第4個距離單元干擾抑制前后的多普勒譜Fig.11 The Doppler spectrum at the 4th range bin before and after RFI suppression

4 結(jié)論

通過分析射頻干擾在時域數(shù)據(jù)上的短時性特征，利用海洋回波信號在傅里葉變換基上的稀疏性，提出了基于壓縮感知的時域射頻干擾抑制方法。該方法對雷達(dá)接收數(shù)據(jù)中含有干擾的掃頻周期分別進(jìn)行處理，通過對部分未受干擾的采樣數(shù)據(jù)進(jìn)行隨機(jī)欠采樣，恢復(fù)被干擾的海洋回波時間采樣數(shù)據(jù)，從而在時域上完成射頻干擾抑制。結(jié)合雷達(dá)信號處理特點(diǎn)對OMP算法進(jìn)行簡化，極大地節(jié)省了算法的時間成本并提高了稀疏重構(gòu)的準(zhǔn)確性。通過實(shí)測數(shù)據(jù)和仿真分析，得到以下結(jié)論。

(1)本文提出的基于壓縮感知的時域抗干擾算法能有效抑制射頻干擾。

(2)當(dāng)射頻干擾與海洋回波在RD譜上重疊時，本方法在有效抑制射頻干擾的同時能保持海洋回波不受影響。

(3)重構(gòu)和替換過程會帶來誤差，少數(shù)極低信噪比的譜點(diǎn)容易受到抗干擾算法的影響，但在實(shí)際海態(tài)參數(shù)反演中，過低信噪比的譜點(diǎn)一般不進(jìn)行反演或反演結(jié)果可信度低，因此不會影響該算法的實(shí)際應(yīng)用效果。