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

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

0 引言

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

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

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

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

1 算法原理

1.1 凸共形陣列的信號(hào)接收模型

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

式中:θ,φ分別為信號(hào)在極坐標(biāo)系下的入射俯仰角與方位角。

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

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

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

1.2 抗干擾算法

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

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

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

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

STAP的傳遞函數(shù)

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

圖1 STAP的等效FIR 濾波器

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

式中:f為頻率;‖·‖為矩陣取模運(yùn)算符。

1.3 降秩處理

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

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

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

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

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

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

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

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

2 仿真分析

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

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

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

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

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

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

圖6 特征值分布情況

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

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

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

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

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

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

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

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

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

3 結(jié)論

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