周 柱 張茂軍 周典樂

(1. 國防科技大學(xué)系統(tǒng)工程學(xué)院，湖南長沙 410073；2. 中國人民解放軍31628部隊(duì)，廣東韶關(guān) 512199)

1 引言

地球表面接收到的導(dǎo)航信號(hào)功率很小，以GPS為例，比接收機(jī)熱噪聲低約21 dB。當(dāng)受到人為干擾時(shí)，接收信號(hào)難以恢復(fù)。運(yùn)用天線陣并且進(jìn)行時(shí)延擴(kuò)展，形成空時(shí)二維模式，可以大幅增加自由度從而增強(qiáng)抗干擾能力。設(shè)陣元數(shù)為M，時(shí)間延遲數(shù)為N，則接收數(shù)據(jù)X為MN×1維向量，其協(xié)方差矩陣RX=E{XXH}，為MN×MN維矩陣，其最優(yōu)處理的運(yùn)算量約為O((MN)3)，隨著空時(shí)處理維數(shù)的增加，運(yùn)算量成立方倍增長，大運(yùn)算量對計(jì)算資源消耗和計(jì)算時(shí)長都不可接受。文獻(xiàn)[1]提出一種改進(jìn)側(cè)視陣聲吶的混響抑制方法，該方法利用投影矩陣法設(shè)計(jì)了時(shí)域和空域混響陷波器，分別濾除了與目標(biāo)同頻的旁瓣混響和與目標(biāo)同方位的主瓣混響，但是仍然要用空時(shí)二維聯(lián)合處理濾除剩余的混響，沒有討論空時(shí)二維處理的簡化運(yùn)算。文獻(xiàn)[2]著眼于雷達(dá)運(yùn)動(dòng)平臺(tái)在下視檢測目標(biāo)的過程中，干擾目標(biāo)污染訓(xùn)練樣本引起的功率非均勻問題，首先從受污染樣本與干凈樣本的差異性度量角度入手，采用均值Hausdorff距離度量樣本矩陣相似性，然后結(jié)合凸優(yōu)化包計(jì)算不同樣本的相似度，最后根據(jù)相似度的不同，實(shí)現(xiàn)對受污染樣本的剔除，該方法提高了協(xié)方差矩陣估計(jì)的準(zhǔn)確性，但缺點(diǎn)是需要用凸優(yōu)化包，即需要用最優(yōu)化理論實(shí)現(xiàn)，這一部分難以顯式表述，難于實(shí)現(xiàn)，且該方法仍然包含求空時(shí)二維矩陣特征值及相應(yīng)逆矩陣的計(jì)算，運(yùn)算量大。文獻(xiàn)[3]針對實(shí)際數(shù)字陣列天線陣元間距大于半波長導(dǎo)致大角度估計(jì)時(shí)陣元相位出現(xiàn)模糊問題，提出一種新的基于時(shí)空DOA矩陣的多目標(biāo)高性能二維角度估計(jì)算法。該算法首先計(jì)算時(shí)空DOA矩陣并進(jìn)行特征值分解，接著利用空時(shí)矩陣特征向量與來波導(dǎo)向矢量匹配特性，對估計(jì)得到的各導(dǎo)向矢量分別進(jìn)行相位解模糊，進(jìn)而利用導(dǎo)向矢量聯(lián)立方程組，得到各來波角度信息的最小二乘解。主要解決的是陣列不理想情況下相位模糊問題，運(yùn)用該方法同樣需要求得空時(shí)矩陣的特征向量，有很大的運(yùn)算量。MWF有避免大矩陣分解和降維處理兩個(gè)優(yōu)勢[4]，適用于解決空時(shí)二維抗干擾處理的問題。但是在實(shí)際中運(yùn)用經(jīng)典MWF，噪聲子空間維數(shù)估計(jì)不準(zhǔn)，當(dāng)超過最優(yōu)值時(shí)，將引起信噪比較大降幅。

基于以上分析，本文通過對MWF特性的分析，用MWF輸出的各階期望信號(hào)算得一個(gè)協(xié)方差矩陣，通過對該矩陣的分析，推導(dǎo)出一種滑窗定階的方法。首先用MWF粗略估計(jì)出噪聲子空間維數(shù)，然后從該階開始取一個(gè)小方陣，用冪法算得其最大特征值，與預(yù)設(shè)的門限比較，判定是否為噪聲子空間維數(shù)，沿協(xié)方差矩陣對角線向右下滑動(dòng)，直至估計(jì)出噪聲子空間維數(shù)。在做出估計(jì)之后，可以運(yùn)用經(jīng)典MWF算得最優(yōu)權(quán)值并對接收信號(hào)進(jìn)行濾波抑制干擾。該方法的優(yōu)勢在于界定噪聲和白噪聲子空間分界點(diǎn)時(shí)，相對于經(jīng)典算法具有更大的區(qū)分度。

2 多級維納濾波的基本方法

多級維納濾波是維納濾波器的一種多級等效形式，因此運(yùn)用這一方法首先需要將抗干擾處理轉(zhuǎn)化為維納濾波的問題。通過廣義旁瓣對消器[5]可以從接收信號(hào)矢量分離出期望信號(hào)，形成維納濾波的模式，如圖1所示。接收信號(hào)矢量經(jīng)過上面的支路得到期望信號(hào)d0(n)，d0(n)包含有用信號(hào)和干擾，阻塞矩陣B0=null(h)，其中null表示求零空間，下支路通過B0阻塞期望信號(hào)，則x0(n)僅含干擾。

圖1 廣義旁瓣對消器結(jié)構(gòu)圖

圖2 多級維納濾波器結(jié)構(gòu)圖

由各階標(biāo)量維納濾波器可嵌套組成多級維納濾波器的綜合濾波器部分，求得各階標(biāo)量維納濾波器權(quán)值之后，逐級反推可得到綜合濾波器的權(quán)值

(1)

最終算得的權(quán)值為

wMWF=Tw

(2)

經(jīng)典MWF算法中的迭代次數(shù)P是通過每步迭代后均方誤差ξi=E{|ei(n)|2}與一個(gè)門限值(通常設(shè)為接收機(jī)熱噪聲方差)相比較來確定，只要ξi小于門限值，就讓迭代停止，使P=i。

運(yùn)用經(jīng)典MWF對空時(shí)二維模式抗干擾處理降維，會(huì)使最優(yōu)子空間維數(shù)難以估計(jì)準(zhǔn)確，有兩個(gè)原因：(1)噪聲隨環(huán)境變化，用一個(gè)固定的門限難以準(zhǔn)確估計(jì)干擾子空間維數(shù)；(2)因?yàn)镸WF的每一步按照與接收信號(hào)矢量的最大相關(guān)提取干擾分量，因此前面的迭代提取干擾分量多，后續(xù)提取少，致使后續(xù)接收信號(hào)殘差方差與白噪聲方差區(qū)分不明顯，很難找準(zhǔn)。

3 運(yùn)用滑窗判定噪聲子空間估計(jì)方法

3.1 理論推導(dǎo)

(3)

(4)

(5)

(6)

依此類推，可以得到

(7)

(8)

因RXT的副對角線元素是共軛對稱的[4]，所以只需要計(jì)算δi即可。則協(xié)方差矩陣可表示如下：

(9)

用MWF逐步提取相關(guān)干擾信息過程中，必然存在P

(10)

式中，(·)(P)表示干擾子空間維數(shù)為P的情況，下標(biāo)(·)i表示MWF的第i級，v(n)為白噪聲矢量。

文獻(xiàn)[8]在對噪聲子空間進(jìn)行分解過程中，通過對各階期望信號(hào)矢量生成的協(xié)方差矩陣進(jìn)行分析，類似于本文中的RXT，推導(dǎo)出：當(dāng)階數(shù)超過噪聲子空間維數(shù)時(shí)，RXT的主對角線元素等于白噪聲方差，副對角線元素接近0。本文也用到這一思路，根據(jù)式(4)對矩陣RXT進(jìn)行分析。其副對角線元素可推導(dǎo)如下

(11)

因此在理論上，RXT可以按照下面形式進(jìn)行分塊：

(12)

左上部分仍為P階三對角矩陣，右下部分為對角矩陣，其余位置元素為零。矩陣RXT的特征值為det(RXT-λI)=0的解，有下式

det(RXT-λI)=

(13)

矩陣RXT可以分塊，以噪聲子空間維數(shù)P為界，左上部分為三對角矩陣，直到第P階，右下部分為對角矩陣，對角元素約等于白噪聲方差。該特性可用于估計(jì)噪聲子空間維數(shù)，根據(jù)式(2)，在獲得子空間維數(shù)之后，需要用到降維矩陣T和綜合濾波器的權(quán)值w，即可算得權(quán)值。

3.2 實(shí)現(xiàn)步驟

空時(shí)二維處理的主要環(huán)節(jié)是權(quán)值計(jì)算，圖3為在MWF理論基礎(chǔ)上提出的一種運(yùn)用滑窗定階的噪聲子空間估計(jì)方法。運(yùn)用該方法完成最優(yōu)權(quán)值的計(jì)算，圖1中“w”標(biāo)識(shí)的部分。

運(yùn)用該方法進(jìn)行空時(shí)二維處理的步驟如下：

圖3 運(yùn)用滑窗判定的噪聲子空間分解方法流程圖

Step1將抗干擾處理問題轉(zhuǎn)化為維納濾波模式：將接收信號(hào)矢量通過廣義旁瓣對消器，得到上支路的期望信號(hào)，記為d0(n)，下支路x0(n)為待處理信號(hào)矢量，如圖1，形成維納濾波的模式。

Step3準(zhǔn)確估計(jì)干擾子空間維數(shù)：如圖3所示，在MWF框架下，根據(jù)子空間理論，結(jié)合矩陣最大特征值的計(jì)算，運(yùn)用滑窗逐級搜尋噪聲子空間維數(shù)。這一步分解為以下四步執(zhí)行：

Step3.1協(xié)方差矩陣計(jì)算： 將接收信號(hào)矢量經(jīng)過MWF生成的各階期望信號(hào)di(n)組成矢量形式xT(n)=[d1(n)，d2(n)，…，dMN-1(n)]T，求取該矢量的協(xié)方差矩陣RXT。

Step3.2截取小矩陣：在RXT中用一個(gè)窗口截取一個(gè)小方陣，設(shè)為Rj，不妨設(shè)方陣寬度為Dwin,j為方陣左上角元素在RXT中的階次。首次截取的位置j=P1+1，也就是粗略估計(jì)的噪聲子空間維數(shù)加1，以j=P1+1作為Rj的右下角，從P1向左上取Dwin大小的方陣，如果P1

Step3.3求取最小特征值：用反冪法求得矩陣Rj的最小特征值(Rj為三對角矩陣，且維度小，運(yùn)算量小)，設(shè)為λj。根據(jù)上文分析，RXT的小特征值接近于白噪聲方差，因此用一個(gè)小窗口沿RXT的對角線滑動(dòng)，當(dāng)達(dá)到分界點(diǎn)時(shí)，其最小特征值必定接近于白噪聲方差。

在工程計(jì)算中，可以利用反冪法按模計(jì)算Rj的最小特征值，其步驟如下：

(1)輸入矩陣Rj，初始向量x，誤差限ε，最大迭代次數(shù)Nmax;

(2)將迭代初始值k置為1，初始的最小特征值設(shè)為λ0=0，令矢量y=x/max(abs(x));

(3)對矩陣Rj作三角分解Rj=LU;

(4)解方程組LUx=y;

(5)令μ=max(x)，λ=μ;

(6)若|λ-λ0|<ε，則輸出1/λ，也就算得了最小特征值，否則執(zhí)行下一步;

(7)若k

Step3.4判定噪聲子空間維數(shù)：將λj與1.5倍σ2比較，如果λj大于1.5σ2，則j=j+1，將窗口沿著待分析矩陣RXT的對角線向右下滑動(dòng)一階，返回到S3.2。如果不滿足λj大于1.5σ2，則認(rèn)為已經(jīng)達(dá)到分界點(diǎn)，就是估計(jì)出的噪聲子空間維數(shù)，即為P。

得到最優(yōu)權(quán)值矢量wMWF之后，用該權(quán)值對接收信號(hào)矢量進(jìn)行陣列濾波，就可以達(dá)到抑制干擾的目的，獲得較高的信干噪比輸出。

4 運(yùn)算量分析

經(jīng)典MWF算法包含阻塞矩陣計(jì)算，運(yùn)算量極大[4]，文獻(xiàn)[6-7]提出改進(jìn)方法縮減了經(jīng)典MWF算法前向迭代部分的運(yùn)算量。本文方法在文獻(xiàn)[4,6-7]的基礎(chǔ)上增加了對噪聲子空間維數(shù)的估計(jì)，該方法前向迭代部分的主要運(yùn)算量分析如下:

表1 前向迭代部分主要運(yùn)算量

上表中，第1步使用的乘法器數(shù)量為4(MN-1)，這是一個(gè)(MN-1)階的復(fù)矢量每一元素和一個(gè)復(fù)數(shù)相乘，因此需要4(MN-1)次乘法。第2步中，求取rxi-1di-1的模可以復(fù)用乘法器，只需要8個(gè)足夠，在求取歸一化相關(guān)矢量時(shí)，可以轉(zhuǎn)化為一個(gè)(MN-1)維的復(fù)矢量與一個(gè)實(shí)數(shù)相乘，因此用到2(MN-1)個(gè)乘法器。第4步中，求取期望信號(hào)自相關(guān)是共軛復(fù)數(shù)相乘，只需要2個(gè)乘法器，而求取相鄰兩階期望信號(hào)的互相關(guān)則需要4個(gè)乘法器，而RXT(i-1，i)即是RXT(i，i-1)的共軛，無需另外計(jì)算。第5步中，同理，歸一化相關(guān)矢量與這一階期望信號(hào)矢量相乘也是MN階的復(fù)矢量每一元素和一個(gè)復(fù)數(shù)相乘，所以用到4(MN-1)個(gè)乘法器。

本文方法相對于改進(jìn)的MWF增加的運(yùn)算主要為表1中的第4步，以及對小維度矩陣Rj的運(yùn)算，經(jīng)以上分析得知，其運(yùn)算量相對于大維度矢量運(yùn)算而言很小，因此可以得出結(jié)論，使用本文方法增加的運(yùn)算量很小。

5 仿真試驗(yàn)

仿真主要包括兩個(gè)部分：第一、比較本文方法與傳統(tǒng)方法在噪聲子空間維數(shù)估計(jì)上的性能；第二、比較本文方法與傳統(tǒng)方法的輸出信干噪比；第三、比較本文方法與傳統(tǒng)方法在仿真場景中的抗干擾性能，此項(xiàng)仿真需采用文獻(xiàn)[9]的方法，用導(dǎo)航電文建立GPS衛(wèi)星運(yùn)動(dòng)模型，以獲得貼近于實(shí)際的仿真場景。

采用5元加芯圓陣，成十字型排列，半徑為d=λ/2，λ為接收信號(hào)波長，時(shí)間延遲單元數(shù)N=5。無干擾時(shí)基帶載噪比C/N0=Sr+Ga-10log(kT0)-NF-Ls。其中Sr為GPS信號(hào)接收功率(dBW)，Ga為天線增益(dB)，10log(kT0)=-205 dBW/Hz為接收機(jī)噪聲密度，取常溫T0=300 K，NF為接收機(jī)噪聲系數(shù)，取決于天線和前置放大器噪聲，一般取NF=4 dB，Ls為信號(hào)在通道內(nèi)的傳輸及量化損耗。據(jù)IDC-GPS-200[10]規(guī)定，地球表面上Ga=0 dB增益的GPS接收機(jī)天線接收L1 頻段上C/A碼信號(hào)的能量為-160 dBW。接收機(jī)噪聲功率密度約為-140 dBW，因此無干擾情況下的信噪比為-21 dB。

仿真試驗(yàn)1按照干噪比30 dB設(shè)置兩個(gè)功率相同的部分帶寬干擾和兩個(gè)單頻干擾。進(jìn)行兩組仿真：(1)試驗(yàn)輸出信干噪比(SINR)與迭代次數(shù)的關(guān)系；(2)檢測本發(fā)明對于噪聲子空間估計(jì)的作用。仿真中均重復(fù)多次取平均。

(1)輸出信干噪比(SINR)與迭代次數(shù)的關(guān)系

運(yùn)用經(jīng)典MWF方法進(jìn)行抗干擾濾波，逐漸增加MWF的迭代次數(shù)，觀測輸出SINR的變化。

如圖4所示，隨著迭代次數(shù)增加，SINR在達(dá)到最大值以后將大幅下降，僅僅超過最佳迭代次數(shù)一次，SINR就下降約有5 dB。

圖4 輸出SINR與迭代次數(shù)關(guān)系

(2)本文方法對于噪聲子空間估計(jì)的作用

仿真中將沿RXT對角線截取的小方陣邊長width定為3，即矩陣Rj的階數(shù)設(shè)置為3。

圖5的左圖表示使用經(jīng)典MWF，通過每步迭代運(yùn)算之后均方誤差確定噪聲子空間維數(shù)的情況，從圖5(左)可以看出，在12次迭代之后，干擾信號(hào)提取完，因?yàn)榇藭r(shí)均方誤差與白噪聲方差的比值用分貝表示接近于0，即均方誤差接近于白噪聲方差。而在進(jìn)行第13次迭代之后，均方誤差與白噪聲比值降了3 dB。圖5的右圖表示運(yùn)用本文方法，在粗略判斷最佳迭代次數(shù)(本場景中P1=10)之后，從P1=10階開始，以P1=10為右下角，向左上方截取一個(gè)小矩陣，用反冪法算得該小矩陣最小特征值，求得該特征值與白噪聲方差的比值，然后使窗口向右下方滑動(dòng)，即以P=11為右下角，以同樣的方式算得特征值與白噪聲方差的比值，將所有比值記錄下來，如圖5(右)所示。因?yàn)閳鼍跋嗤遗c運(yùn)用經(jīng)典MWF在同一次仿真中，所以無疑最佳迭代次數(shù)為12。而在圖5(右)中，第12和13階比值的差距有7 dB。這就意味著區(qū)分度大大高于經(jīng)典方法對最佳迭代次數(shù)的區(qū)分度。而且在第13階的位置，方差接近于白噪聲，與理論分析相符。因此可以得出結(jié)論，運(yùn)用本方法可以使得干擾子空間分解更為準(zhǔn)確可靠。

仿真試驗(yàn)2按照一個(gè)部分帶寬、一個(gè)單頻干擾的順序逐步增加干擾，干擾入射角度在入射空間內(nèi)錯(cuò)開，干噪比30 dB。

圖5 經(jīng)典MWF與滑窗定階方法估計(jì)噪聲子空間維數(shù)的判據(jù)

圖6 經(jīng)典MWF與滑窗定階方法輸出信干噪比(SINR)的對比

圖6的表示經(jīng)典MWF和本文方法在抗干擾性能上的對比。由圖可以看出，干擾數(shù)為1和2時(shí)，本文的滑窗定階法獲得的輸出SINR比用經(jīng)典MWF高約1 dB，干擾數(shù)為3到6，滑窗定階法獲得的SINR高約2～3 dB，干擾數(shù)為7時(shí)，滑窗定階法獲得的SINR低約1 dB，干擾數(shù)為8時(shí)，兩種方法獲得的SINR相近。即在可抗干擾數(shù)范圍內(nèi)，滑窗定階法獲得的SINR大部分情況下要更高。需要指出的是：在仿真環(huán)境中，運(yùn)用經(jīng)典MWF時(shí)，是按照算得的白噪聲方差設(shè)置門限，應(yīng)為最佳門限；而在實(shí)際中，往往難以找到這樣一個(gè)最佳門限，而運(yùn)用滑窗定階則大大削弱了這個(gè)問題，因?yàn)檫@種方法的區(qū)分度遠(yuǎn)大于經(jīng)典方法。所以結(jié)論為：在實(shí)際環(huán)境中運(yùn)用滑窗定階方法相對于經(jīng)典方法，可以獲得更高更穩(wěn)定的抗干擾性能。

仿真試驗(yàn)3按照文獻(xiàn)[9]的方法，用RINEX格式的導(dǎo)航電文文件adis0050.13n建立仿真場景。在一段時(shí)間內(nèi)，標(biāo)記出抗干擾處理之后的可用衛(wèi)星。各衛(wèi)星在一段時(shí)間內(nèi)按規(guī)律運(yùn)動(dòng)，發(fā)射信號(hào)，在這一段時(shí)間內(nèi)持續(xù)進(jìn)行抗干擾處理，也就達(dá)到了在可視空域內(nèi)對算法的抗干擾性能進(jìn)行測試的目的。設(shè)接收機(jī)所在位置為東海某處上空，坐標(biāo)為東經(jīng)123° 00′，北緯25° 40′，即25.67°，高程為500 m。模擬產(chǎn)生部分帶寬和單頻兩種干擾。部分帶寬、單頻干擾各2個(gè)。干擾入射方位角分別為：70°，130°，190°，250°，俯仰角分為25°，5°，15°，20°。干噪比均為30 dB。

圖7、圖8中，子圖(a)均表示一段時(shí)間內(nèi)，從接收機(jī)往上空看到的可用衛(wèi)星的運(yùn)動(dòng)軌跡，軌跡表示可用衛(wèi)星在不同時(shí)刻的空間位置；子圖(b)表示同一時(shí)間可用衛(wèi)星數(shù)目。可用衛(wèi)星數(shù)目的多少也就反映了不同抗干擾方法的性能。圖中仿真開始時(shí)間(tow)表示導(dǎo)航電文開始記錄的時(shí)間。由仿真圖可知：(1)使用經(jīng)典多級維納濾波做抗干擾處理，可用衛(wèi)星比較多，能夠滿足定位需求；(2)使用本文方法可以抗干擾性能得到較大提升。圖8(a)中衛(wèi)星的軌跡較為連續(xù)，圖7(a)中部分軌跡有中斷的現(xiàn)象；由圖7(b)和圖8(b)可以看出，在2小時(shí)的測試時(shí)段內(nèi)，圖8(b)中顯示的可用衛(wèi)星數(shù)量明顯要多于圖7(b)中。由此可知，使用本文的滑窗法可以保證接收機(jī)在長時(shí)間內(nèi)有更多數(shù)量的衛(wèi)星可用。此外，在使用經(jīng)典MWF時(shí)，通過反復(fù)比較設(shè)定了相對最優(yōu)的門限，因此能獲得較好性能，但是在實(shí)際環(huán)境中，往往難以進(jìn)行最優(yōu)設(shè)置，因此使用本文方法相對于經(jīng)典MWF可以獲得更好更穩(wěn)定的抗干擾性能。

圖7 用經(jīng)典MWF算法抗干擾之后的可用衛(wèi)星軌跡與數(shù)目

圖8 用滑窗定階方法抗干擾之后的可用衛(wèi)星

6 結(jié)論

本文運(yùn)用子空間理論分析多級維納濾波過程，將接收信號(hào)矢量映射為各階期望信號(hào)并組成矢量型式，根據(jù)此信號(hào)矢量所形成協(xié)方差矩陣的特性，推導(dǎo)出一種運(yùn)用滑窗定階的方法。用此方法判斷噪聲子空間維數(shù)，其區(qū)分度比用經(jīng)典方法明顯要高。仿真證明此方法在噪聲子空間維數(shù)估計(jì)方面顯著優(yōu)于經(jīng)典方法，運(yùn)用此方法能獲得更好的抗干擾性能。