王龍崗 岳顯昌 吳雄斌 易先洲
(武漢大學電子信息學院,武漢 430072)
短波波段垂直極化的電磁波沿海洋表面?zhèn)鞑r具有衰減小且能夠繞射傳播的特點,因此高頻地波雷達(high frequency surface wave radar, HFSWR)可克服地球曲率的影響,實現(xiàn)對海洋表面和低空目標的超視距、全天候探測,具有探測范圍廣、成本低、實時性好等特點,成為探測海面和低空目標的重要手段[1].
HFSWR進行目標探測時,除了目標回波信號外,回波當中還存在大量的干擾信號,其中由海浪散射產生的回波信號稱為海雜波. 回波中一階海雜波的強度往往比目標信號強得多[2],而且艦船等低速目標的多普勒頻率與海雜波的Bragg峰相近,從而被淹沒在海雜波中,很難被檢測到[3]. 因此對海雜波進行抑制,是HFSWR實現(xiàn)對一階海雜波譜中目標檢測的關鍵.
目前,HFSWR檢測艦船目標的方法主要有海雜波對消法和現(xiàn)代譜估計方法[4]. 其中海雜波對消方法包括基于模型的對消方法、循環(huán)對消方法和預測對消方法. 基于模型的對消方法是通過對海雜波進行建模,然后將利用模型得到的海雜波作為背景,進行對消處理[5-6]. 這種方法要求建立的模型與真實的海雜波特性較好地吻合,否則無法達到理想的效果. 循環(huán)對消方法主要有Root[7]、郭欣等[8]提出的基于快速傅里葉變換(fast Fourier transform, FFT)的幅度、相位分析的對消方法和基于擴展Prony算法的循環(huán)對消算法[9]. 預測對消法[10]是借助神經網絡通過數(shù)據(jù)訓練來提取海雜波的參數(shù),進行海雜波對消. 循環(huán)對消法和預測對消法是將海雜波的主要分量近似看作正弦信號,估計對應的頻率、幅度和相位等參數(shù),對消強的海雜波及其諧波分量. 現(xiàn)代譜估計方法主要有小波變換法和基于子空間估計的海雜波抑制算法.在小波變換方法中,基函數(shù)的選擇對抑制效果有較大的影響,需要選擇合適的基函數(shù)才能得到較好效果,并且這種方法會破壞回波信號的相位信息[11-12].子空間估計法是根據(jù)雜波在子空間的聚集特性來實現(xiàn)海雜波的抑制,如基于特征值分解(eigen value decomposition, EVD)的方法[13]、基于奇異值分解(singular value decomposition, SVD)的方法[14-15]、改進特征值分解(modified eigen value decomposition, MEVD)的方法[16]. 這些基于子空間估計的海雜波抑制方法,只能在目標與海雜波多普勒頻率差別較大的情況下實現(xiàn)海雜波的抑制. 當兩者頻率接近或者重疊時,目標會被誤消,造成漏警. 而高階奇異值分解(high order SVD, HOSVD)方法,從方位維、距離維和慢時間維三維聯(lián)合來估計海雜波子空間,但是由于一階回波的復雜性,不容易實現(xiàn)海雜波子空間與噪聲子空間的準確劃分,對海雜波的抑制效果不理想. 此外還有斜投影(oblique projection, OP)方法,在已知海雜波方位的情況下,通過空間OP來抑制海雜波[17],然而在實際探測目標時,海雜波的方位無法預先得知.
針對上述問題,本文提出了基于SVD的空域海雜波抑制方法(以下簡稱空域SVD算法),充分利用天線陣列的信息,在海雜波方位未知的情況下,通過SVD來估計空域的海雜波子空間和噪聲子空間,然后利用子空間之間的正交性,通過投影實現(xiàn)在空間域對海雜波的抑制. 空域SVD算法是在空域對海雜波進行抑制的,一般來說距離-多普勒(range-Doppler,RD)譜中單個譜點信號源只存在有限的幾個方位,在對陣列信號的協(xié)方差矩陣進行SVD時,得到的奇異值數(shù)目較少,比較容易準確劃分海雜波子空間與噪聲子空間,海雜波的抑制效果更好. 文中首先研究了HFSWR海雜波產生的機理,然后闡述了基于SVD的空域海雜波抑制算法的原理,最后利用仿真與實測數(shù)據(jù)對算法性能進行驗證,并將結果與現(xiàn)有的空域海雜波抑制方法進行了對比,驗證了本文算法的有效性.
海浪的運動狀態(tài)是一種復雜和混沌的隨機過程[18],由不同頻率、波高和傳播方向的海浪疊加而成,這些海浪近似呈現(xiàn)正弦波動,如圖1所示[4].干涉. 設海浪的波長為L,雷達發(fā)射電磁波的波長為
圖1 海雜波形成機理Fig. 1 Formation mechanism of sea clutter
HFSWR發(fā)射的高頻無線電波與海浪發(fā)生相互作用,使電磁波發(fā)生散射,散射的電磁波之間會發(fā)生λ,當電磁波束照射到圖1所示的海浪上時,入射波到達海浪各散射點的路程差為Lcosα,當滿足式(1)時,海浪后向散射產生的回波信號間的波程差為λ,回波信號同相加強,產生Bragg諧振:
式中,α為入射波的擦地角. 當α趨于0時,即當無線電波波長等于海浪波長的兩倍時,發(fā)生Bragg諧振,在回波多普勒譜中產生兩個尖峰,稱為一階Bragg峰.
根據(jù)深水區(qū)動力學原理,重力波的波長與流速之間滿足:
海浪相對雷達站的徑向速度為
Bragg峰對應的多普勒頻率為
將式(1)、(2)和(3)代入式(4)可得
式中:g表 示重力加速度;fB表示Bragg峰的多普勒頻率,單位是Hz;f0表示雷達工作的中心頻率,單位是MHz;±是指海浪朝向和背離雷達. 通常雷達的擦地角α很小,式(5)可近似為
式中,fB表示Bragg峰的多普勒頻率的理論值. 實際上由于海流的存在,實測的多普勒頻率相比理論值會存在海流引起的多普勒頻率fc左右的偏移,這是造成海雜波一階譜展寬的重要原因. 因此,在對海雜波進行建模仿真的過程中,必須考慮海流的影響[19].
在一個積累周期內,對接收陣的每個陣元接收的回波信號進行傅里葉變換,得到對應的RD譜[20]. 不同距離、不同速度的散射源的回波信號將被分離開. 對于一階譜中的目標,無法從距離維和多普勒維將其與海雜波區(qū)分開,考慮在空域對他們進行分離,并抑制海雜波.
假設空間中有K個窄帶遠場信號s1(t),s2(t),···,sK(t)入射到如圖2所示線性均勻陣列,設陣元的個數(shù)為N(N>K),陣元間距為d,信號源的入射方向分別為?i(i=1,2,···,K),其中前P個方位為海雜波的方位.N個陣元在t時刻的輸出數(shù)據(jù)矢量為
式中:x(t)= [x1(t),x2(t),···,xN(t)]T;A表 示陣列流型,A=[a(?1),a(?2),···,a(?K)],其中,,;s表示信號矢量,s=[s1,s2,···,sK]T,si表示第i個 信號源的復振幅;n表示均值為0、方差為σ2的高斯白噪聲,n=[n1,n2,···,nN]T;C表示海雜波在均勻線性陣列N個陣元的輸出數(shù)據(jù)矢量,C=A1s1,A1=[a(?1),a(?2),···,a(?P)]為海雜波的陣列流型,s1=[s1,s2,···,sP]T為 海雜波的復振幅;T表示目標在均勻線性陣列的N個陣元的輸出數(shù)據(jù)矢量,T=A2s2,A2=[a(?P+1),(?P+2),···,a(?K)]為目標的陣列流型,s2=[sP+1,sP+2,···,sK]T為目標的復振幅.
圖2 均勻線性陣列Fig. 2 Uniform linear array
出數(shù)據(jù)矢量可記為x,假設各信號源方位不同,其協(xié)方差矩陣為
對于RD譜中的一個譜點,陣列接收該信號的輸
式中:A1和A2是相互正交的滿秩矩陣;IN表示N維的單位矩陣.
為了估計空域的海雜波子空間和噪聲子空間,對陣列接收信號的協(xié)方差矩陣C進行SVD運算:
式中:U代表左奇異矩陣;S代表由奇異值構成的對角陣;V代表右奇異矩陣. 由于海雜波的強度遠大于噪聲強度,因此根據(jù)奇異值的大小進行子空間的劃分,大奇異值對應的矢量構成的子空間為海雜波子空間U1,小奇異值對應的矢量構成的子空間為噪聲子空間U2,由于海雜波與目標在不同方位,U1和U2是正交的. 海雜波子空間的投影矩陣為U1U1H,噪聲子空間的投影矩陣為將接收的回波信號x向U1空間進行投影,得到投影分量為
進而得到與海雜波子空間正交的子空間內的投影向量為
由于噪聲子空間與海雜波子空間正交,根據(jù)矢量運算法則,經過子空間投影后,海雜波得到了抑制.
在實際運用中,由于目標回波信號與海雜波信號混雜在一起,無法從接收的回波數(shù)據(jù)中準確估計出雜波子空間. 因此在估計子空間時,利用海雜波在相近的距離元之間具有較強的相關性[17],可以選擇感興趣譜點前后鄰近的距離元對應譜點作為參考單元,來估計雜波子空間.
基于SVD的空域海雜波抑制的基本步驟如下.
1)對雷達接收機的每個天線單元接收的數(shù)據(jù)進行兩次FFT得到RD譜,并選擇感興趣譜點作為待處理單元.
2)選擇待處理單元前后鄰近的譜點作為參考單元,對參考單元的數(shù)據(jù)進行短時傅里葉變換,并計算參考單元陣列輸出數(shù)據(jù)矢量xf的協(xié)方差矩陣C.
3)對參考單元的協(xié)方差矩陣C進行SVD,并根據(jù)奇異值的大小,將空間化分為雜波子空間U1和目標與噪聲子空間U2.
4)利用式(10)和(11),將待處理單元陣列輸出矢量x向U1空間進行投影,并從陣列輸出矢量x中減去其在U1空間的投影分量x′,從而達到抑制海雜波的目的.
算法流程圖如圖3所示.
圖3 算法流程圖Fig. 3 Flowchart of the algorithm
利用仿真實驗,將空域SVD算法對海雜波的抑制性能與現(xiàn)有的空域海雜波抑制算法進行對比. 在仿真海雜波時,設置載頻為13.15 MHz,脈沖重復周期為0.5 s,積累周期數(shù)目為600,Bragg峰的多普勒頻率為fB=±0.37 Hz,加入的噪聲為高斯白噪聲,設置海流的速度為0.5 m/s,正的一階譜區(qū)的多普勒頻率為[0.32,0.41]. 接收天線為八元天線陣,近似線性. 加入三個速度分別為4 m/s、4.5 m/s和?4 m/s的目標,對應的多普勒頻率為0.35 Hz、0.39 Hz和?0.35 Hz,目標距離雷達站均為50 km,相對于雷達站的方位分別是北偏東65°、北偏東35°、北偏東10°,初始的信雜噪比(signal to clutter plus noise ratio, SCNR)為?8 dB、?15 dB和?10 dB,仿真得到的第一個天線單元的10號距離元的多普勒譜如圖4所示. 可以看出,由于海流的存在,海洋回波的一階譜發(fā)生了展寬,在[0.32, 0.41] Hz和[?0.41, ?0.32] Hz兩個頻率區(qū)域內,回波信號強度較大,對應的正是展寬的一階譜區(qū). 由于仿真的目標回波強度低于海雜波,目標被海雜波淹沒,箭頭所指處僅是目標所對應的多普勒頻率. 要想檢測被淹沒在一階譜區(qū)的目標,需要對海雜波進行抑制.
圖4 海雜波抑制前的多普勒譜(仿真)Fig. 4 Doppler spectrum without sea clutter suppression (simulation)
圖5給出了利用空域SVD算法、HOSVD算法和OP算法三種方法對海雜波進行抑制后的多普勒譜. 利用空域SVD算法處理數(shù)據(jù)時,由于實測實驗中目標回波信號的旁瓣會泄漏到相鄰距離元,因此在實際估計海雜波子空間時參考距離單元與感興趣的距離單元中間間隔一個距離單元. 對一階譜內的每一個譜點的陣列向量x減去其在海雜波子空間的投影,從而實現(xiàn)在空間維對海雜波的抑制. HOSVD算法則是利用一階譜內多普勒維、距離維和方位維的三維數(shù)據(jù)構造三階張量,通過HOSVD來抑制海雜波. OP算法是在已知海雜波方位情況下計算其陣列流型,并對一階譜內的每一個譜點的陣列向量通過OP來抑制海雜波.
圖5 抑制海雜波后的多普勒譜(仿真)Fig. 5 Doppler spectrum with sea clutter suppression(simulation)
從圖5可以看出,在仿真實驗中,三種算法均能抑制海雜波,使目標凸顯出來,提高了SCNR,有利于一階譜內目標的檢測. 在海雜波方位已知的情況下,OP算法與空域SVD算法在海雜波抑制方面性能幾乎是相同的,而且輸出的SCNR比HOSVD算法至少提高了5 dB,使目標更加突出. 主要是由于空域海雜波抑制算法進行SVD的維數(shù)較低,RD譜上一個譜點對應的方位維信源個數(shù)較少,容易準確估計出海雜波子空間,對海雜波的抑制效果更好. OP算法預先已知海雜波方位,可以計算出其陣列流型和海雜波子空間,因此這兩種算法通過投影消除海雜波均能有比較好的效果,而且對目標的影響較小. 而HOSVD算法對空間維、距離維和多普勒維進行了聯(lián)合處理,SVD算法得到的奇異值個數(shù)較多,難以利用奇異值的大小進行子空間的準確劃分,存在海雜波子空間向目標子空間的泄露問題,在消除海雜波的同時,也會造成目標能量被削弱.
為了驗證算法的有效性,做了100次蒙特卡洛實驗,圖6為利用三種算法抑制海雜波之后的輸出SCNR. 可以看出,空域SVD算法和OP算法輸出的SCNR平均高于HOSVD算法約7 dB,而且具有較好的穩(wěn)定性,說明在抑制海雜波方面,空域海雜波抑制算法具有較好的性能.
圖6 抑制海雜波之后的輸出SCNR(仿真)Fig. 6 Output SCNR with sea clutter suppression (simulation)
通過基于實測數(shù)據(jù)的仿真結果來對比分析HOSVD算法、OP算法和空域SVD算法對海雜波的抑制性能. 實驗采用的是2016年12月26日東山站單站發(fā)射、單站接收的實驗數(shù)據(jù),雷達站接收天線共有八個單元,且近似成線性排列,雷達的工作頻率為13.15 MHz,脈沖重復周期為0.5 s,Bragg峰的多普勒頻率為fB=±0.37Hz,脈沖積累周期數(shù)為600,積累時間為5 min,距離分辨率為5 km. 實驗過程中,在實測數(shù)據(jù)中加入了三個仿真目標,目標均在距離雷達站50 km處,相對于雷達站的方位分別是北偏東30°、北偏東60°、北偏東45°,速度分別為正北5 m/s、正西5 m/s、正北6 m/s,對應的多普勒頻率為?0.38 Hz、0.38 Hz、?0.36 Hz,初始的SCNR均為?5 dB,目標被完全淹沒,加入仿真目標后的10號距離元的多普勒譜如圖7所示. 把海流流速為± 0.5 m/s對應的多普勒區(qū)域當作是一階譜區(qū),對應的多普勒頻率區(qū)間分別為[0.32,0.41] Hz和[?0.41,?0.32] Hz.
圖7 加入仿真目標后的多普勒譜(實測)Fig. 7 Doppler spectrum with simulation target (measured)
分別利用HOSVD算法、OP算法和空域SVD算法對數(shù)據(jù)進行處理. 在運用空域SVD算法抑制海雜波時,參考單元與感興趣距離單元之間間隔一個距離元,經過SVD后,將前5個大奇異值對應的子空間設為海雜波對應的子空間;HOSVD算法中三階張量的構造依然是利用前20個距離元、一階譜區(qū)的多普勒元與所有陣列單元的數(shù)據(jù);對于實測數(shù)據(jù),無法預先知道海雜波的方位,通過MUSIC算法為OP算法提供海雜波的估計方位,進而得到海雜波的子空間. 目標所在的10號距離單元利用三種算法對海雜波進行抑制后的多普勒譜如圖8所示.
圖8 抑制海雜波之后的多普勒譜(實測)Fig. 8 Doppler spectrum with sea clutter suppression(measured)
從圖8可以看出,HOSVD算法、OP算法和空域SVD算法都能抑制較強的海雜波分量,在目標處形成一個尖峰,然而空域SVD算法抑制海雜波后輸出的SCNR比HOSVD算法至少高4 dB,OP算法輸出的SCNR介于兩者之間,因此空域SVD算法在圖8中形成的尖峰比較明顯. 產生上述現(xiàn)象的原因是實測數(shù)據(jù)的一階譜區(qū)的海洋回波分量比較復雜,根據(jù)經驗劃分的海雜波子空間與噪聲子空間不準確,海雜波子空間向噪聲子空間的擴散造成海雜波不能被有效消除,旁瓣較高,子空間劃分不準確,對HOSVD算法海雜波的抑制性能影響較大. 而空域SVD算法進行SVD時,奇異值數(shù)目較少,比較容易進行子空間劃分,且子空間投影算法利用兩個相距較近的距離元的空域信息,海雜波的相關性比較強,通過投影法海雜波消除得更充分,對目標影響較小. OP算法要求預先知道的海雜波子空間信息在實測中無法得到,通過MUSIC算法估計海雜波的方位,進而獲取子空間,使得誤差較大,因此實測中對于海雜波的抑制效果OP算法比空域SVD算法穩(wěn)定性差. 在實測中三種算法的性能都比仿真結果差,可能是實測中海雜波子空間和噪聲子空間不完全正交,海雜波不能完全被消除造成的.
為了驗證算法的有效性,利用三種算法進行了基于24場2個小時實測數(shù)據(jù)的仿真實驗,圖9為利用三種算法抑制海雜波之后的輸出SCNR.
圖9 抑制海雜波之后的輸出SCNR(實測)Fig. 9 Output SCNR with sea clutter suppression (measured)
從圖9可以看出:空域SVD算法的輸出SCNR比HOSVD高約5 dB,OP算法的輸出SCNR介于兩者之間;但OP算法輸出SCNR波動比較劇烈,因為通過MUSIC算法估計海雜波方位,進而獲取海雜波對應的陣列流型和子空間會存在較大誤差,進而影響OP算法的性能. 這三種算法的輸出SCNR隨著時間變化都會波動,在抑制海雜波時都有一定效果,但是不能完全抑制海雜波,目標周圍的雜波殘余分量會變化,最終造成輸出SCNR的起伏. 圖10給出了這2個小時內海流方位與目標方位之間的角度差.可以看出,雖然目標的方位不變,但是海流的方向會發(fā)生變化,海流與目標之間的方位差會隨時間變化.由于空域SVD算法是在空域對海雜波進行抑制,當目標與海雜波的來波方向相同或相近時,目標存在被誤消的可能,也會造成輸出SCNR的起伏.
圖10 海流和目標方位之間的差值(實測)Fig. 10 Difference of direction between current and target(measured)
本文在SVD的基礎上提出了空域SVD算法,用于一階譜區(qū)海雜波的抑制,進而實現(xiàn)一階譜中目標的檢測. 對算法進行了理論分析和推導,并利用仿真與實測數(shù)據(jù)驗證了算法的有效性. 理論仿真和基于實測數(shù)據(jù)的海雜波抑制結果均表明,與HOSVD算法、OP算法相比,空域SVD算法在未能預先知道海雜波的子空間的情況下,能夠比較準確地估計海雜波子空間,產生比較好的海雜波抑制效果;該算法利用海雜波在距離上的相關性,通過在空間域子空間投影抑制海雜波,不受目標多普勒頻率與Bragg頻率差值大小的影響;在空間域進行SVD,奇異值數(shù)目較少,容易劃分雜波子空間與噪聲子空間,對海雜波有比較好的抑制效果. 但是該算法的性能受到海洋回波在距離上相關性的影響,應用有一定的局限,而且只能將海雜波分量消弱,并不能完全消去海雜波;當目標與海雜波處于相同或相近方位時,存在目標被誤消的情況.