無畸變共形陣列抗干擾算法的降秩處理方法

李 晗, 何 迪, 陳 新, 曲家慶

(1.北斗導(dǎo)航與位置服務(wù)上海市重點實驗室,上海交通大學(xué), 上海 200240;2.上海無線電設(shè)備研究所, 上海 201109)

0 引言

全球衛(wèi)星導(dǎo)航系統(tǒng)(global navigation satellite system,GNSS)如今已經(jīng)成為一重要空間基礎(chǔ)設(shè)施,可為全球用戶提供定位、導(dǎo)航和授時服務(wù)[1-2]。北斗衛(wèi)星導(dǎo)航系統(tǒng)(Beidou navigation satellite system,BDS)是完全由中國自主設(shè)計和建造的多功能全球衛(wèi)星導(dǎo)航系統(tǒng)。與現(xiàn)有的其他GNSS 類似,BDS 也被設(shè)計為能夠承受一定水平的射頻干擾[3],這是通過直接序列擴頻技術(shù)[4]的應(yīng)用來實現(xiàn)的,該技術(shù)能夠抑制常見的噪聲信號。然而,除了接收機熱噪聲之外,在許多情況下還存在無意和有意的干擾,破壞GNSS信號接收的可靠性。為了解決這些問題,提高GNSS 接收機的魯棒性,在過去的幾十年中涌現(xiàn)了許多GNSS干擾抑制方法。實踐證明,空時自適應(yīng)處理(spacetime adaptive processing,STAP)是最有效的方法之一。STAP技術(shù)將天線陣列與有限長度沖激響應(yīng)(FIR)濾波器相結(jié)合,以加權(quán)和的形式自適應(yīng)地重新排列天線的接收信號,以引導(dǎo)天線陣列響應(yīng)波束指向所需方向。將零點指向干擾方向,使GNSS信號原樣通過,同時有效抑制干擾。該方法在抑制高功率干擾方面效果很好[5]。然而,STAP 中引入FIR 濾波器會引起非線性相位響應(yīng),進而導(dǎo)致所需信號的失真。以往的實驗表明,如果不特別注意這個問題,STAP 帶來的定位誤差可以達到幾十米[6]。因此,在高精度GNSS應(yīng)用中,必須考慮該問題,并對信號失真加以補償[7-8]。

共形陣列是一種符合其承載表面的天線陣列,旨在實現(xiàn)更好的空氣動力學(xué)和流體動力學(xué)性能,或出于美學(xué)考慮。因此,共形陣列在許多場景中得到廣泛應(yīng)用,例如用于高速飛機的雷達或BDS 接收機。傳統(tǒng)的干擾抑制技術(shù)不能直接應(yīng)用于共形陣列,因為通常假設(shè)所有共形陣列天線陣元都是全向的,并且具有相同的響應(yīng)模式。然而,由于承載面的非平面曲率,共形陣列天線設(shè)計時應(yīng)采用定向陣元模型,不同的天線陣元通常指向不同的方向,進而導(dǎo)致不同天線陣元對同一信號的響應(yīng)不同。這意味著不能將單陣元響應(yīng)函數(shù)作為平面陣列分析中常用的陣列流形的公因子進行分離。因此,必須修正信號接收過程的數(shù)學(xué)模型以適應(yīng)上述非理想條件。目前,共形陣列信號處理的研究主要集中在模式合成和波達方向(direction of arrival,DOA)估計兩個方面,對共形陣列中的干擾抑制技術(shù)關(guān)注較少[9]。

空時自適應(yīng)處理技術(shù)帶來的另一副作用是,當(dāng)時延抽頭數(shù)增大時,STAP 觀測向量的協(xié)方差矩陣維數(shù)會急速增大,矩陣求逆的速度會影響到算法的實時性。而事實上在許多應(yīng)用中,使用一個維數(shù)較小的子空間對觀測空間進行近似時,抗干擾性能下降并不明顯,但卻大大減少了計算量,因而對基于STAP 的干擾抑制方法進行降秩簡化有重要意義。目前常用的降秩方法包括主成分(principal components,PC)方法、互譜度量(cross spectral metric,CSM)方法和多級維納濾波(multistage Wiener filter,MWF)方法。降秩后處理器穩(wěn)態(tài)性能不可能超過全維處理器,但更小的運算量和更快的收斂速度使其瞬態(tài)性能可能略優(yōu)于全維處理器[10]。降秩處理的關(guān)鍵在于選擇合適的轉(zhuǎn)換矩陣,使處理器的維數(shù)盡可能低,同時盡可能提高瞬態(tài)性能,并使穩(wěn)態(tài)性能接近最佳。

通過對傳統(tǒng)MWF 處理器結(jié)構(gòu)的改進,可以進一步降低計算量,并在保證抗干擾能力的同時,不影響算法的線性相位特性,減小處理過程引起的衛(wèi)星信號畸變。

1 算法原理

1.1 凸共形陣列的信號接收模型

基于共形陣列的北斗接收機不能沿用傳統(tǒng)的陣列信號處理模型。即使采用輻射特性相同的陣元,載體表面曲率也會導(dǎo)致各個陣元的指向不同,進而導(dǎo)致在陣列天線的全局坐標系內(nèi)各陣元的方向圖函數(shù)不同?？紤]到實際天線方向圖的性能需求與模型簡化,使用余弦天線模型進行討論,其方向圖函數(shù)

式中:θ,φ分別為信號在極坐標系下的入射俯仰角與方位角。

首先假設(shè)任意M元凸共形陣列的所有陣元輻射特性相同,均為余弦天線;然后根據(jù)陣列具體結(jié)構(gòu)確定一個全局直角坐標系oxyz。設(shè)pm為第m個陣元在其中的空間坐標矢量,且陣元指向與當(dāng)前位置載體表面的法線方向相同。定義共形因子矩陣G[9]來表征陣元方向圖及陣列的共形帶來的影響,G表示為

式中:g（θd,φd）表示期望信號方向的共形因子,其中θd,φd為期望信號在全局極坐標系下的入射俯仰角及方位角;g（θjk,φjk） 表示第k個干擾信號方向的共形因子,其中θjk,φjk為第k個干擾信號在全局極坐標系下的入射俯仰角及方位角,k∈{1,2,…,K},K為干擾個數(shù);θi(θ,φ)表示第i個陣元處的本地極坐標系的入射俯仰角;T表示矩陣轉(zhuǎn)置運算。結(jié)合歐拉旋轉(zhuǎn)定理可以得到凸共形陣列的觀測向量

式中:A為陣列流形陣列;S(n)為入射信號矢量;d(n)和jk(n)分別為期望信號和第k個干擾信號;V(n)為噪聲信號矢量。定義共形陣列流形C=G☉A,則

1.2 抗干擾算法

在干擾場景復(fù)雜、干擾數(shù)多、陣元數(shù)量有限制的情況下,STAP是替代純空域濾波處理的較好方法。此外,STAP除了調(diào)整天線陣的空域響應(yīng),還能調(diào)整時域響應(yīng),用于補償射頻、中頻信號失真以及天線失配等問題。

當(dāng)時延抽頭數(shù)為P時,基于式(4)定義STAP信號觀測向量

得到基于共形陣列的無失真空時最小方差無失真響應(yīng)(distortionless conformal minimum variance distortionless response,Distortionlessc MVDR)方法[9]。優(yōu)化問題及約束條件表述為

式中:min(·)為取最小值函數(shù);E(·)為求數(shù)學(xué)期望函數(shù);y(n)為輸出信號;|·|為取絕對值運算符;w為抽頭系數(shù)矩陣;wmp為第m個陣元的第p個抽頭的權(quán)系數(shù),其中m∈{1,2,…,M},p∈{1,2,…,P},P為奇數(shù);H 為矩陣共軛轉(zhuǎn)置運算符;R=E（ZZH）為信號觀測向量的自相關(guān)矩陣;c為約束向量;cd為信號對應(yīng)的共形導(dǎo)向矢量。最終抽頭系數(shù)的最優(yōu)解

STAP的傳遞函數(shù)

式中:ω為歸一化頻率;Ts為采樣周期;為cd的第m個元素;* 為共軛運算符;wp=[w1pw2p…wMp]為M個陣元的第p個抽頭的權(quán)系數(shù)組成的向量;為第p個抽頭的系數(shù)。從式(12)的形式可以看出,整個系統(tǒng)的傳輸函數(shù)可以由一個抽頭數(shù)為P的FIR 濾波器來等效,其中第p個抽頭的系數(shù)hp由其權(quán)系數(shù)和信號的共形導(dǎo)向矢量共同來決定。等效FIR濾波器的結(jié)構(gòu)如圖1所示。

圖1 STAP的等效FIR 濾波器

由式(11)和式(12)可以得到等效FIR 濾波器的系數(shù)向量h= [h1h2…h(huán)P]T。h關(guān)于中間系數(shù)共軛對稱,則系統(tǒng)可以等價為一個線性相位濾波器,不會引入載波相位偏差,且群延遲為一個與信號來向無關(guān)的固定常數(shù)(P-1)Ts/2[9]。因此,針對任意方向的北斗信號,經(jīng)過此系統(tǒng)都會引入一個相同且等于此群延遲的偽碼相位偏差,故可以直接補償?shù)?不會影響最終的定位結(jié)果。補償后系統(tǒng)的傳遞函數(shù)可以表示為

式中:f為頻率;‖·‖為矩陣取模運算符。

1.3 降秩處理

使用一個合適的低維子空間對觀測空間進行近似,即降秩處理,可以在保證抗干擾算法性能的同時降低計算量。降秩的廣義旁瓣相消架構(gòu)處理如圖2所示。

圖2 降秩的廣義旁瓣相消架構(gòu)處理

輸入信號Z(n)通過上支路后得到含有北斗導(dǎo)航信號的期望信號d0(n)。下支路通過阻塞矩陣B0阻塞上支路信號得到Z0(n),使其僅包含干擾信號。B0=null(c)是通過求解c的零空間得到的阻塞矩陣,其中c是Distortionlessc MVDR 方法中的約束向量。然后Z0(n)通過維納濾波器W0估計得到對d0(n)的估計0(n)。維納濾波器W0通過適當(dāng)?shù)腗P×D維降秩轉(zhuǎn)換矩陣T,使輸入數(shù)據(jù)從MP維降低到D維的Z1(n),再通過WD進行自適應(yīng)加權(quán)處理,得到對d0(n)的估計0(n)。最后上下支路相減得到誤差信號ε0(n)。WD表示在D維空間內(nèi)處理的最優(yōu)權(quán)值向量。降秩的MWF則是選擇合適的D值,將圖2所示分解過程重復(fù)D次,此時,降維的觀測向量中幾乎不含干擾分量,協(xié)方差矩陣已趨于白化。對應(yīng)的降秩矩陣T和WD分別表示為

式中:hi為第i級分解的約束向量;Bi為第i級分解的阻塞矩陣;wi為最優(yōu)權(quán)的第i項。

MWF的過程中沒有出現(xiàn)矩陣求逆和特征分解運算,進一步降低了運算復(fù)雜度,但仍需計算協(xié)方差矩陣和阻塞矩陣。基于此考慮,使用相關(guān)相減(correlation subtraction architecture,CSA)的方法產(chǎn)生阻塞矩陣,對算法進行簡化,以避開協(xié)方差矩陣的計算,同時在分解的過程中通過采樣均方誤差(SMSE)的比較實現(xiàn)在線秩選。通過對逆向遞歸的綜合濾波器的改進,實現(xiàn)只用一次前向迭代同時完成秩選和濾波器系數(shù)求解,降低計算量。上述濾波器稱為迭代相關(guān)相減多級維納濾波器(iterative-CSA-MWF,ICSA-MWF),其整體結(jié)構(gòu)如圖3所示。

圖3 迭代實現(xiàn)的相關(guān)相減多級維納濾波器結(jié)構(gòu)

迭代實現(xiàn)的相關(guān)相減多級維納濾波器算法流程如圖4 所示,其中size(·)為獲取維數(shù)函數(shù),σSMSE為采樣均方誤差。

圖4 迭代實現(xiàn)的相關(guān)相減多級維納濾波器算法流程

2 仿真分析

2.1 參數(shù)設(shè)置

基于前述信號接收模型,以一種實際工程所使用的帶中心陣元的圓柱共形陣列為例,實現(xiàn)抗干擾算法。圓柱共形陣列結(jié)構(gòu)如圖5 所示。其中陣元數(shù)M=5,圓柱載體半徑R=0.1 m,中心陣元p1指向天頂方向(z向),外圍陣元p2～陣元p5相對于中心陣元繞圓柱軸線旋轉(zhuǎn)角α=30°,圓柱半長L=0.102 5 m。o1x1y1z1是以p2為原點,圓柱徑向為z1向,圓柱軸向為x1向的局部坐標系。

圖5 共形陣列的結(jié)構(gòu)圖

設(shè)時延抽頭數(shù)P=5。為了更好地展示降秩維數(shù)D不同帶來的性能差異,使用3個寬帶干擾和1 個窄帶干擾的配置。其中北斗信號來向為(30°,30°),信噪比為-30 d B;阻塞干擾來向分別為(40°,90°),(60°,180°)和(50°,60°),干噪比均設(shè)為70 d B;窄帶干擾為單頻干擾,頻率與北斗信號中心頻率相同,來向為(70°,240°),干噪比為70 dB。

2.2 協(xié)方差矩陣的特征值

在進行降秩方法性能分析前,應(yīng)首先對觀測向量構(gòu)成的協(xié)方差矩陣的特征值分布情況進行分析,仿真結(jié)果如圖6所示。其中橫軸為按特征值的絕對值大小排序后的特征值編號i,縱軸為特征值的絕對值?？梢钥吹教卣髦档慕^對值大小明顯分為兩段,曲線段的特征值即為干擾分量所對應(yīng)的特征值,水平線部分則為噪聲對應(yīng)的特征值。其中,寬帶干擾和窄帶干擾的特征值分布特征有較大區(qū)別:寬帶干擾對應(yīng)的特征值整體較大,窄帶干擾則相反,可以從窄帶干擾能量更集中的角度來考慮。同時根據(jù)多次不同條件下的實驗,得出如下經(jīng)驗結(jié)論:STAP 觀測向量協(xié)方差矩陣對應(yīng)大特征值數(shù)目?！諲JP,其中NJ即為總的干擾數(shù)目。

圖6 特征值分布情況

2.3 降秩維數(shù)D 對抗干擾算法的影響

設(shè)定不同的D值進行實驗。其中用于對比的PC和CSM 方法,按對應(yīng)準則取D個特征值和特征向量獲得降秩矩陣T;對于ICSA-MWF方法則是設(shè)定固定的迭代次數(shù),不再使用SMSE進行秩選。分別從抗干擾性能和保持等效FIR濾波器的線性相位性的能力兩方面來進行分析,并取全秩處理時的結(jié)果作為對照,標記為Fullrank。4種處理方法在不同快拍數(shù)下的抗干擾性能對比結(jié)果如圖7所示。

圖7 降秩維數(shù)D 對抗干擾性能的影響

由圖7可明顯看出,隨著降秩維數(shù)D的增大,各降秩方法的輸出信噪比均逐漸增大,其中D<10的部分因為維數(shù)過低導(dǎo)致抗干擾算法幾乎失效,參考意義不大。此外,可以看出PC方法的性能相對于另外兩種方法較差,較高快拍條件下達到全秩處理時的輸出信噪比水平所需降秩維數(shù)仍等于觀測向量協(xié)方差矩陣的大特征值數(shù)目Γ(此處Γ=20),明顯高于另兩種方法。CSM 方法得益于自身的特點,即采用遍歷法確定使輸出信噪比最大化的D個特征向量,其抗干擾性能優(yōu)異,且隨著快拍數(shù)增加,達到全秩處理的輸出信噪比水平所需的D值逐漸減小。在快拍數(shù)為10 000時,雖然Γ=20,但CSM 方法在D=16時即可達到全秩處理的輸出信噪比水平。而ICSAMWF 方法的性能介于PC 方法和CSM 方法之間,且在低快拍時性能仍較好。

降秩處理的本意就是為了降低運算復(fù)雜度,提高算法實時性。CSM 方法雖然性能較好,但其引入了過多的額外計算量,與降秩處理的本意相悖。ICSA-MWF 方法省去了協(xié)方差矩陣的計算及其特征分解,且可以通過一次前向迭代實現(xiàn)自適應(yīng)秩選、降秩矩陣計算和降秩權(quán)值計算,大大降低了STAP 算法實現(xiàn)的復(fù)雜度。因此本文認為ICSA-MWF 是更為合適的降秩方法。

選取幾個關(guān)鍵D值,通過仿真信號群延時,分析上述各降秩方法對于Distortionless-c MVDR算法的線性相位特性的影響,仿真結(jié)果如圖8所示。

圖8 降秩維數(shù)D 對群延遲的影響

由圖8可知,當(dāng)D=10時,降秩維數(shù)過低,不僅抗干擾性能差,抗干擾算法的群延遲波動也較大,線性相位特性不理想。當(dāng)D=16 時,CSM、ICSA-MWF和Full rank 曲線幾乎重合,CSM方法的輸出信噪比已經(jīng)接近全秩處理水平,算法也接近線性相位系統(tǒng),ICSA-MWF 方法抗干擾性能雖未達到最佳,但已可以保持較理想的線性相位特性。當(dāng)D=Γ=20時,各方法均接近線性相位系統(tǒng),其中ICSA-MWF和Full rank 曲線幾乎重合。

經(jīng)過本小節(jié)的分析,說明了上述降秩方法取合適的D值時,既可保證輸出信噪比的水平,也可以保證Distortionless-c MVDR 算法的線性相位特性。再綜合各算法的計算復(fù)雜度,ICSAMWF 是最為適合的降秩方法。

3 結(jié)論

本文提出了一種針對無失真空時共形MVDR 算法的降秩處理方法。該方法結(jié)合了相關(guān)相減方法和對逆向遞歸的綜合濾波器的改進,規(guī)避了協(xié)方差矩陣的計算,實現(xiàn)了在線秩選以及只用一次前向迭代同時完成秩選和濾波器系數(shù)求解,極大地降低了計算量。仿真結(jié)果表明,當(dāng)取合適的D值時,既可保證輸出信噪比的水平,也可以保證Distortionless-c MVDR 算法的線性相位特性,減小處理過程對衛(wèi)星信號產(chǎn)生的畸變影響。