劉 鵬, 王 盾, 彭 博
(1. 北京衛(wèi)星信息工程研究所, 北京 100095; 2. 天地一體化信息技術(shù)國家重點實驗室, 北京 100095)
衛(wèi)星導(dǎo)航系統(tǒng)具備全天候提供精確的導(dǎo)航、定位和授時信息的能力,在國民生活和國防建設(shè)中發(fā)揮著越來越重要的作用[1-3]。衛(wèi)星導(dǎo)航信號是一種擴頻信號,其在到達地面時極其微弱,僅靠擴頻增益本身無法保證接收機在強干擾環(huán)境下正常工作,必須借助其他的抗干擾措施[4-7]。自適應(yīng)陣列處理具有良好的抗干擾能力,在衛(wèi)星導(dǎo)航接收機的研制過程中得到廣泛應(yīng)用,常見的算法包含空時自適應(yīng)處理(space-time adaptive processing, STAP)[8-10]算法以及空頻自適應(yīng)處理[11-13]算法。
自適應(yīng)陣列處理的核心工作在于利用某種最優(yōu)準(zhǔn)則,求取陣列天線的濾波權(quán)值。濾波權(quán)值的求取方法包括塊自適應(yīng)處理算法和連續(xù)自適應(yīng)處理算法[14]。兩種算法分別利用當(dāng)前數(shù)據(jù)塊或數(shù)據(jù)快拍的二階統(tǒng)計特征,求取濾波權(quán)值,用于對當(dāng)前及后續(xù)若干數(shù)據(jù)塊或數(shù)據(jù)快拍進行自適應(yīng)濾波,達到保留有用信號、濾除干擾信號的目的。按照干擾的信號統(tǒng)計特征是否隨時間變化,干擾可以分為平穩(wěn)信號和非平穩(wěn)信號,其中非平穩(wěn)信號泛指具有時變能量譜的確定性信號或具有時變功率譜的隨機信號[15-16]。當(dāng)陣列接收到的干擾為平穩(wěn)信號時,相鄰數(shù)據(jù)的統(tǒng)計特征基本一致,自適應(yīng)處理算法可獲得良好的干擾抑制效果。然而,當(dāng)陣列接收到的干擾為非平穩(wěn)信號時,相鄰數(shù)據(jù)之間的統(tǒng)計特征可能存在較大差異,則此時自適應(yīng)處理算法的抗干擾性能便會顯著下降。
針對非平穩(wěn)信號,必須采用能夠表現(xiàn)信號局域性的分析方法,常用分析方法包括短時傅里葉變換(short-time Fourier transform, STFT)、時頻分布、Gabor展開、小波分析、分?jǐn)?shù)階傅里葉變換(Fourier transform, FT)等[17-20]。將上述方法與自適應(yīng)陣列處理相結(jié)合,可以提升算法對非平穩(wěn)干擾的抑制能力。已有研究者對此開展了相關(guān)研究工作,文獻[21]將STFT、文獻[22-23]將Cohen類時頻分布、文獻[24-25]將小波分析分別引入了陣列信號處理,并對其性能進行了分析。然而,上述研究所采用的陣列處理方法主要是STAP算法,對于目前導(dǎo)航系統(tǒng)常用的SFAP算法與非平穩(wěn)信號處理方法結(jié)合的有效性,相關(guān)研究成果較少。文獻[26]雖然開展了將STFT與SFAP相結(jié)合的研究工作,根據(jù)干擾時頻特征對采樣數(shù)據(jù)進行分組以提升抗干擾性能,但該算法僅適用于存在單個非平穩(wěn)干擾的場景。當(dāng)環(huán)境中存在多個非平穩(wěn)干擾時,該算法的干擾抑制性能仍然會出現(xiàn)下降。
LFM信號是一種典型的非平穩(wěn)信號,廣泛被應(yīng)用于雷達、聲吶、通信等領(lǐng)域[27-30]。衛(wèi)星導(dǎo)航接收機在實際工作中很容易受到LFM信號的干擾,因此有必要針對該信號進行抗干擾設(shè)計。
基于以上分析,本文提出了一種基于數(shù)據(jù)空時頻三維特征分組的SFAP算法,對常規(guī)SFAP算法進行了改進,提升其對LFM干擾的抑制能力。考慮到STFT與SFAP中的離散傅里葉變換(discrete Fourier transform, DFT)結(jié)構(gòu)相近,首先通過STFT獲取采樣數(shù)據(jù)的時域、頻域聯(lián)合分布,并利用空間相關(guān)系數(shù)分析相同頻率干擾在不同時間的空間相關(guān)性,然后對采樣數(shù)據(jù)進行分組,將不同時間具有相同頻率和到達角度(direction of arrival, DOA)參數(shù)的采樣點分到相同組,最后利用分組數(shù)據(jù)進行協(xié)方差矩陣估計和權(quán)值計算,并利用所得權(quán)值對相同統(tǒng)計特征的數(shù)據(jù)進行抗干擾濾波,從而達到提高干噪比、增加零陷深度的目的,實現(xiàn)抗LFM干擾性能的提升。該算法對存在單個和多個LFM干擾的場景均能適用,仿真結(jié)果證明了該算法的有效性。
假設(shè)陣列包含L個天線單元,且各天線單元各方向響應(yīng)一致。陣列接收到I個導(dǎo)航信號和J個干擾,且滿足J (1) 式中:X(t)=[x1(t),x2(t),…,xL(t)]T為L×1陣列數(shù)據(jù)向量;Si(t)為陣列接收第i個導(dǎo)航信號si(t)形成的L×1向量;Jj(t)為陣列接收第j個干擾jj(t)形成的L×1向量;n(t)為L×1噪聲向量。 將第l個天線單元接收到的信號si(t)相對參考點的時延表示為τil,接收到的干擾jj(t)相對參考點的時延表示為τjl,則第i個信號向量和第j個干擾向量可分別表示為 Si(t)=[si(t-τi1),si(t-τi2),…,si(t-τil)]T (2) Jj(t)=[jj(t-τj1),jj(t-τj2),…,jj(t-τjl)]T (3) SFAP需要在頻域?qū)π盘栠M行處理,將陣列接收數(shù)據(jù)進行FT,得 (4) 式中:X(f)、Si(f)、Jj(f)和n(f)分別為X(t)、Si(t)、Jj(t)和n(t)各組成元素經(jīng)FT后對應(yīng)的頻域L×1向量,且 Si(f)=[si(f)e-j2πfτi1,si(f)e-j2πfτi2,…,si(f)e-j2πfτil]T= [e-j2πfτi1,e-j2πfτi2,…,e-j2πfτil]Tsi(f)=aSi(f)si(f) (5) Jj(f)=[jj(f)e-j2πfτj1,jj(f)e-j2πfτj2,…,jj(f)e-j2πfτjl]T= [e-j2πfτj1,e-j2πfτj2,…,e-j2πfτjl]Tjj(f)=aJj(f)jj(f) (6) 式中:si(f)、jj(f)分別為導(dǎo)航信號si(t)和干擾jj(t)的FT;aSi(f)、aJj(f)分別為導(dǎo)航信號si(t)和干擾信號jj(t)在頻率f處的導(dǎo)向矢量。 由式(4)、式(5)和式(6)可得 AS(f)S(f)+AJ(f)J(f)+n(f) (7) 式中:AS(f)=[aS1(f),aS2(f),…,aSI(f)]和AJ(f)=[aJ1(f),aJ2(f),…,aJJ(f)]分別為信號陣列流形矩陣和干擾陣列流形矩陣;S(f)=[s1(f),s2(f),…,sI(f)]T為信號頻譜向量;J(f)=[j1(f),j2(f),…,jJ(f)]T為干擾頻譜向量。 SFAP是一種寬帶抗干擾算法,可以理解為將寬帶數(shù)據(jù)分解為多個窄帶數(shù)據(jù),分別對各窄帶數(shù)據(jù)進行加權(quán)濾波,最后將濾波后的各窄帶數(shù)據(jù)合成寬帶數(shù)據(jù)輸出。SFAP的主要流程如圖1所示。 圖1 SFAP流程圖Fig.1 Flowchart of SFAP SFAP需要在每個子帶獨立進行自適應(yīng)權(quán)值計算及濾波處理,常規(guī)算法的主要實現(xiàn)過程如下[26]: (1) 將時域采樣數(shù)據(jù)分為M個數(shù)據(jù)塊,每個數(shù)據(jù)塊的長度為N,采樣數(shù)據(jù)記為xlm(n)。其中,l為陣元號,m為數(shù)據(jù)塊號,n為采樣計數(shù)。 (8) 式中:w(n)為長度N的窗函數(shù);WN=e-j2π/N。 (9) (10) 假設(shè)信號、干擾與噪聲不相關(guān),則有 (11) (12) (4) 求取第k個頻帶的濾波權(quán)值 (13) 式中:μ為常數(shù);a為約束矢量,當(dāng)采用自適應(yīng)調(diào)零算法時,取a=[1,0,…,0]T。 (5) 利用所得權(quán)值在各個頻帶進行自適應(yīng)濾波 (14) (6) 對濾波后的頻域數(shù)據(jù)進行IDFT變換,并對變換后的時域數(shù)據(jù)加反窗: (15) 上述步驟為SFAP的主要處理過程,經(jīng)抗干擾濾波后的時域數(shù)據(jù)還可根據(jù)需要進行重疊處理,去除每段數(shù)據(jù)的起始及末尾的部分?jǐn)?shù)據(jù),以消除DFT變換對首末數(shù)據(jù)的不利影響。 LFM信號是一種典型的非平穩(wěn)信號,本文選用LFM信號作為干擾進行算法性能分析。 持續(xù)時間為T的LFM信號可表示為 (16) 式中:A為信號幅度;f0為初始頻率;K為調(diào)頻斜率;φ為初始相位。LFM信號的頻率變化范圍為[f0:f0+KT]。 在進行導(dǎo)航抗干擾測試時,LFM信號的頻率通常在導(dǎo)航信號帶寬內(nèi)周期性變化,構(gòu)成掃頻干擾,可表示為 (17) 當(dāng)采用常規(guī)SFAP對掃頻干擾進行抑制時,考慮到降低運算量的需求,DFT運算長度不能太大,通常導(dǎo)致每個數(shù)據(jù)塊的時間長度小于掃頻周期。假設(shè)接收到的J個干擾均為LFM信號,第j個干擾的掃頻周期對應(yīng)Pj個數(shù)據(jù)塊,即在一個掃頻周期內(nèi)進行Pj次DFT運算,則頻率fk對應(yīng)的干擾頻域采樣非零概率為1/Pj。 (18) 式中:diag(·)表示對角矩陣;I表示L×L單位矩陣。 根據(jù)算法原理,SFAP在某一頻格對干擾的抑制能力與相應(yīng)頻格對應(yīng)的干擾特征值成正比,即干擾特征值越大,零陷越深,干擾抑制能力越強。由式(18)可知,常規(guī)SFAP算法由于未對不同數(shù)據(jù)塊的干擾特性進行區(qū)分,將導(dǎo)致協(xié)方差矩陣中的干擾功率值降低,干擾特征值相應(yīng)降低,從而影響抗干擾能力。 (19) 式中:α1(k)、α2(k)分別是兩個干擾在頻率fk的導(dǎo)向矢量。 對R(k)進行特征值分解,則 (20) 將特征值按照從大到小的順序排列,由式(19)、式(20)可得 (21) (22) 當(dāng)干擾功率遠大于噪聲功率時,則兩個干擾方向在第k個頻帶的陣列響應(yīng)幅值為 (23) 特別地,當(dāng)β=0時 (24) 可見,當(dāng)干擾功率遠大于噪聲功率時,干擾特征值近似與干擾非零數(shù)據(jù)塊的占比成正比,干擾方向的陣列響應(yīng)幅值近似與干擾非零數(shù)據(jù)塊的占比成反比。當(dāng)非零數(shù)據(jù)塊的占比較低時,干擾方向的陣列響應(yīng)幅值將明顯增加,即干擾零陷深度降低,從而影響算法的抗干擾能力。進一步分析式(18)和式(19)可知,當(dāng)干擾個數(shù)增加時,上述結(jié)論依然成立。 針對以上問題,本文將信號時頻分析的時域、頻域特征與陣列處理的空域特征相結(jié)合,按照干擾的時間、頻率以及DOA參數(shù),對SFAP算法的采樣數(shù)據(jù)進行分組。將不同時間具有相同頻率和DOA參數(shù)的采樣點分到相同組,并利用分組后的數(shù)據(jù)進行協(xié)方差矩陣估計、權(quán)值計算和自適應(yīng)濾波,實現(xiàn)提高干擾特征值、增加零陷深度、提升抗干擾能力的目的。 本文的時頻特征分析與空域特征分析方法如下。 (1) 時頻特征分析:選用STFT算法,分塊采樣數(shù)據(jù)xlm(n)的STFT可定義為 (25) 由式(25)和式(8)可知,當(dāng)STFT的窗函數(shù)與DFT窗函數(shù)相同時,有 (26) 可見,當(dāng)STFT的時間分辨率取N時,可以用分塊DFT對其進行替換。 (2) 空域特征分析:為了判別不同時刻相同頻率的干擾是否來自同一干擾源,即其干擾來向是否一致,本文不直接估計干擾來向,而是利用空間相關(guān)系數(shù)進行判斷,空間相關(guān)系數(shù)的定義見式(22)。當(dāng)不同時刻相同頻率的干擾來向相同時,空間相關(guān)系數(shù)模值接近于1,本文取閾值0.985,即等效干擾來向估計誤差約為10°(arccos 0.985≈10°)。 改進后的算法處理過程如下: (1) 將時域采樣數(shù)據(jù)進行連續(xù)分塊,每個數(shù)據(jù)塊的長度為N,采樣數(shù)據(jù)記為xlm(n),其中l(wèi)為陣元號,m為數(shù)據(jù)塊號,n為采樣計數(shù),并對分塊數(shù)據(jù)進行時域加窗和N點DFT運算: (27) (28) 在實際工程應(yīng)用中,陣列各接收通道的頻率響應(yīng)無法完全一致,因此在利用數(shù)據(jù)向量進行干擾DOA參數(shù)估計時,必須考慮通道間幅相誤差的影響。但本算法不需要獲取準(zhǔn)確的干擾DOA參數(shù),不同數(shù)據(jù)塊的接收向量受幅相誤差的影響一致,因而大大降低了幅相誤差對本算法的影響。 (3) 利用分組編號相同的M個頻域樣本來估計第k個頻帶的協(xié)方差矩陣: (29) 式中: (30) 仍然設(shè)信號、干擾與噪聲不相關(guān),且不考慮信號對協(xié)方差矩陣的影響,則 (31) (4) 求取第j組第k個頻帶的濾波權(quán)值: (32) (5) 利用權(quán)值在各個頻帶進行自適應(yīng)濾波: (33) (6) 對濾波后的頻域數(shù)據(jù)進行IDFT變換,并對變換后的時域數(shù)據(jù)加反窗: (34) 為了消除DFT變換對首尾數(shù)據(jù)的不利影響,本文算法在抗干擾過程中對時域數(shù)據(jù)進行重疊處理,重疊數(shù)據(jù)的處理過程相互獨立。設(shè)置本文數(shù)據(jù)重疊率為50%,即抗干擾處理完成后僅保留每段數(shù)據(jù)中間50%的數(shù)據(jù)。 需要說明的是,在步驟2對接收數(shù)據(jù)進行分組時,當(dāng)接收數(shù)據(jù)中無干擾或包含平穩(wěn)干擾時,則判別結(jié)果為:所有數(shù)據(jù)塊在某一頻點不存在干擾或存在干擾且不同數(shù)據(jù)塊之間的空間相關(guān)系數(shù)接近于1。此時,將所有接收數(shù)據(jù)統(tǒng)一處理,改進算法的處理過程與常規(guī)處理算法相同,即在處理無干擾或平穩(wěn)干擾數(shù)據(jù)時,本文所提算法等效為常規(guī)算法。 (35) 式中:αj(k)表示第j個干擾在頻率fk的導(dǎo)向矢量。 對R(j,k)進行特征值分解,則 (36) 將特征值按照從大到小的順序排列,由式(35)和式(36)得 (37) 式中:η1為干擾特征值;v1為干擾特征向量,其余為噪聲特征值和噪聲特征向量。則第j個干擾方向在第k個頻帶的陣列響應(yīng)幅值為 (38) 由式(24)和式(38)可知,本文提出的算法在進行自適應(yīng)濾波后,干擾方向在相應(yīng)頻帶的陣列響應(yīng)值更低,可獲得更深的干擾零陷,因此具有更強的抑制非平穩(wěn)干擾的能力,性能提升增益值(單位為dB)為 (39) 基本仿真參數(shù)設(shè)置為:導(dǎo)航信號選擇北斗導(dǎo)航系統(tǒng)BDS B3I信號,標(biāo)稱載波頻率為1 268.52 MHz,調(diào)制方式為二進制相移鍵控(binary phase shift keying, BPSK),偽碼速率為10.23 Mcps,并將信噪比設(shè)置為-25 dB。天線陣列為四陣元方陣,相鄰陣元間距為半波長。基帶信號采樣率設(shè)置為62/3 MHz,SFAP中的DFT運算數(shù)據(jù)長度為256,窗函數(shù)選用旁瓣抑制100 dB的Chebyshev窗。 仿真中的坐標(biāo)系定義為:陣列放置于水平面,陣列幾何中心與1號陣元連線為方位角0°,逆時針方向為方位角正角度方向,陣列平面法線方向為俯仰角0°。 施加2個由LFM信號構(gòu)成的掃頻干擾,干擾的參數(shù)設(shè)置如表1所示。表中2個LFM干擾的頻率分別周期性地按照從小到大和從大到小的規(guī)律變化。為便于分析,設(shè)置的調(diào)頻斜率使一次頻率變化周期對應(yīng)1 024個采樣點。 表1 干擾參數(shù)設(shè)置 對陣列第一通道的數(shù)據(jù)按照50%重疊進行DFT變換,則此時的DFT運算等效為時間分別率為N/2的STFT。變換后的頻譜圖如圖2所示,圖中的頻率值0對應(yīng)BDS B3頻點的中心頻率??梢?STFT變換可反映LFM干擾頻率隨時間變化的特征,并可根據(jù)頻譜圖識別不同數(shù)據(jù)塊的干擾頻率。 圖2 接收數(shù)據(jù)頻譜圖Fig.2 Spectrum of receiving data 根據(jù)接收數(shù)據(jù)的時頻特征,可分析干擾的時間及頻率分布。以數(shù)據(jù)塊2、數(shù)據(jù)塊6和數(shù)據(jù)塊10為例,上述3個數(shù)據(jù)塊均在-7.5~-2.5 MHz與2.5~7.5 MHz頻段周圍存在干擾。 然而,相同頻率的干擾可能來自不同方向。進一步對數(shù)據(jù)空域特征進行分析,以判斷不同數(shù)據(jù)塊相同頻率的干擾是否對應(yīng)同一干擾源。數(shù)據(jù)塊2與數(shù)據(jù)塊6、數(shù)據(jù)塊10在相應(yīng)頻段的空間相關(guān)系數(shù)如圖3、圖4所示。 可見,數(shù)據(jù)塊2和數(shù)據(jù)塊10在相應(yīng)頻率的空間相關(guān)系數(shù)接近于1,大于閾值0.985,因此可判定數(shù)據(jù)塊2和數(shù)據(jù)塊10相應(yīng)頻率包含的干擾來自同一干擾源,可將其分到相同數(shù)據(jù)組。綜合利用接收數(shù)據(jù)的空時頻三維特征,即可完成數(shù)據(jù)分組。 圖3 空間相關(guān)系數(shù)分析(-7.5~-2.5 MHz)Fig.3 Analysis of spatial correlation coefficient (-7.5~-2.5 MHz) 圖4 空間相關(guān)系數(shù)分析(2.5~7.5 MHz)Fig.4 Analysis of spatial correlation coefficient (2.5~7.5 MHz) 對協(xié)方差矩陣進行特征值分解,由式(21)和式(37)可見,常規(guī)算法的最大和次大特征值為干擾特征值,其余特征值為噪聲特征值,改進算法的最大特征值為對應(yīng)分組的干擾特征值,其余特征值為噪聲特征值。由于噪聲特征值基本穩(wěn)定,因此可以用干擾特征值和噪聲特征值的比值(干噪比)來表征干擾特征值的變化趨勢。同時對干擾方向的陣列響應(yīng)進行分析,干擾方向的陣列響應(yīng)值越小,合成方向圖形成的零陷越深,抗干擾效果越好。 分別對第32、第96、第160和第224個頻點的干噪比和干擾方向的陣列響應(yīng)進行了100次仿真,仿真結(jié)果如圖5、圖6所示。可見,本文所提算法相比原算法的干噪比提升約6 dB,干擾方向的陣列響應(yīng)降低約6 dB,即干擾零陷深度增加約6 dB,仿真結(jié)果與第3.2節(jié)的理論分析一致。 圖5 干擾1特征值及陣列響應(yīng)分析Fig.5 Eigenvalue and array response analysis of Interference 1 圖6 干擾2特征值及陣列響應(yīng)分析Fig.6 Eigenvalue and array response analysis of Interference 2 由式(21)可知,圖5和圖6中常規(guī)算法的最大特征值、次大特征值雖然與干擾1特征值、干擾2特征值不能完全對應(yīng),但由于式(21)中空間相關(guān)系數(shù)模值較小,不影響上述結(jié)論成立。 將干信比提升至110 dB和115 dB,導(dǎo)航星的方向設(shè)置為方位角90°,方位角45°,其他仿真參數(shù)保持不變。 對經(jīng)過抗干擾處理的數(shù)據(jù)與本地導(dǎo)航信號進行互相關(guān)處理,歸一化后的相關(guān)圖如圖7和圖8所示。可見,當(dāng)干信比為110 dB時,本文算法與原算法均能在正確位置形成相關(guān)峰值,但本文算法的底噪更低。當(dāng)干信比進一步提升至115 dB時,本文算法仍然能夠在正確位置形成相關(guān)峰值,但原算法的相關(guān)峰值已不明顯。 圖7 干信比為110 dB時的歸一化相關(guān)峰值Fig.7 Normalized correlation peak value of interference to signal ratio of 110 dB 圖8 干信比為115 dB時的歸一化相關(guān)峰值Fig.8 Normalized correlation peak value of interference to signal ratio of 115 dB 本文針對LFM干擾統(tǒng)計特征時變引起的常規(guī)SFAP抗干擾性能下降問題,提出了一種基于數(shù)據(jù)空時頻三維特征分組的SFAP改進算法,根據(jù)采樣數(shù)據(jù)的空時頻特征對數(shù)據(jù)進行分組,并進一步利用分組數(shù)據(jù)進行協(xié)方差矩陣估計、權(quán)值計算和自適應(yīng)濾波,從而提升了SFAP對LFM干擾的抑制能力,理論分析和仿真實驗均驗證了所提算法的有效性。 在實際工程應(yīng)用中,本文算法的實現(xiàn)復(fù)雜度相比常規(guī)算法有一定提升,主要是需要增加存儲器的容量。一方面,芯片化是當(dāng)前抗干擾算法工程應(yīng)用的趨勢,大容量的存儲器通過芯片化較易獲取,另一方面也可以通過外擴存儲器的方式實現(xiàn)。因此,本算法除具有理論價值,也具有一定的工程應(yīng)用價值。2 常規(guī)SFAP算法及其存在的問題
2.1 常規(guī)SFAP算法概述
2.2 LFM干擾抑制問題分析
3 基于空時頻特征分組的SFAP算法
3.1 算法描述
3.2 改進算法抗干擾性能分析
4 仿真實驗
4.1 干擾空時頻特征分析
4.2 協(xié)方差矩陣特征值及干擾抑制能力分析
4.3 導(dǎo)航相關(guān)性能分析
5 結(jié) 論