基于壓縮感知的相關(guān)干涉儀測(cè)向算法

茍曉鳴 孫中森 唐懷玉 單中堯 李永慶

(中國(guó)電波傳播研究所,青島 266107)

引 言

自20世紀(jì)80年代以來(lái),相關(guān)干涉儀測(cè)向技術(shù)被廣泛應(yīng)用于無(wú)線電監(jiān)測(cè)系統(tǒng)中,相比傳統(tǒng)的相位干涉儀,有助于消除設(shè)備有關(guān)的誤差和場(chǎng)地的影響[1-3],其基本原理如下:相關(guān)干涉儀測(cè)量天線陣列中的參考天線單元與其他天線單元收到的信號(hào)之間的相位差異,構(gòu)成復(fù)值矢量,稱為校準(zhǔn)矢量. 在測(cè)向系統(tǒng)校準(zhǔn)過(guò)程中會(huì)記錄不同入射角度和頻率信號(hào)對(duì)應(yīng)的校準(zhǔn)矢量形成數(shù)據(jù)庫(kù). 在進(jìn)行測(cè)向時(shí)通過(guò)尋找測(cè)量矢量與校準(zhǔn)矢量數(shù)據(jù)庫(kù)之間的最佳匹配以確定入射角度. 匹配過(guò)程通常利用計(jì)算矢量相關(guān)系數(shù)來(lái)實(shí)現(xiàn),相關(guān)系數(shù)最大處即為入射角度. 當(dāng)校準(zhǔn)矢量所選取的角度和/或頻率步進(jìn)較大時(shí),還需進(jìn)行插值運(yùn)算以提高估計(jì)精度.

相關(guān)干涉儀測(cè)向信號(hào)處理的數(shù)學(xué)問(wèn)題是利用測(cè)量矢量恢復(fù)出含有各個(gè)入射角度匹配程度的矢量(稱為匹配度矢量),具有如下兩個(gè)特征:1)欠定性(underdetermined):測(cè)量矢量長(zhǎng)度遠(yuǎn)小于匹配度矢量. 2)稀疏性(sparse):在理論上,匹配度矢量只有數(shù)個(gè)非零元素,對(duì)應(yīng)真實(shí)的入射角度. 壓縮感知(compressive sensing, CS)是解決欠定稀疏問(wèn)題的數(shù)學(xué)工具之一,近年來(lái)被廣泛應(yīng)用于各類信號(hào)和圖像處理問(wèn)題[4-6]. 本文基于CS原理提出一種相關(guān)干涉儀測(cè)向信號(hào)處理算法,第1節(jié)將用矩陣代數(shù)的形式介紹相關(guān)干涉儀信號(hào)處理的數(shù)學(xué)模型,第2節(jié)將介紹經(jīng)典相關(guān)干涉儀信號(hào)處理算法,第3節(jié)將介紹本文提出的算法,第4節(jié)則給出計(jì)算機(jī)仿真實(shí)例來(lái)驗(yàn)證算法的性能.

1 數(shù)學(xué)模型

以雙通道相關(guān)干涉儀為例建立信號(hào)處理數(shù)學(xué)模型.假設(shè)測(cè)向天線是由N個(gè)天線單元構(gòu)成的圓陣,其中N號(hào)天線單元作為參考陣元,1～N-1號(hào)天線單元作為選通陣元.在某一頻點(diǎn)上,天線接收到了M個(gè)同頻信號(hào)sm(t)(m=1,2,…,M),令第n(n=1,2,…,N-1)號(hào)天線單元的復(fù)值輸出為xn(t),則有

(1)

式中:vn(t)為第n個(gè)天線單元的噪聲;φn,m為第n個(gè)天線單元接收的第m個(gè)信號(hào)與參考陣元之間的相位差;j表示虛部單位. 本文假設(shè)在該頻點(diǎn)上不存在多徑傳播的信號(hào),即M個(gè)同頻信號(hào)之間相互獨(dú)立;各個(gè)天線單元的噪聲為空間白噪聲(spatially white noise),即各個(gè)噪聲之間相互獨(dú)立;噪聲來(lái)源為元器件噪聲、背景噪聲等等,與信號(hào)之間相互獨(dú)立. 選通陣元與參考陣元輸出的復(fù)相關(guān)為

(2)

y=[y1,…,yN-1]T

(3)

假設(shè)在該頻點(diǎn)上建立的測(cè)向數(shù)據(jù)庫(kù)矩陣為A,即

A=[a1,…,aQ]

(4)

式中:Q為測(cè)向數(shù)據(jù)庫(kù)角度個(gè)數(shù);a1,a2,…,aQ為各個(gè)角度上的校準(zhǔn)矢量;ω為校準(zhǔn)相位差. 當(dāng)不存在誤差且建庫(kù)角度足夠密集時(shí),矩陣B應(yīng)為矩陣A的一個(gè)子集,于是相關(guān)干涉儀測(cè)向信號(hào)處理可以用下面的線性方程來(lái)表示:

y=Ax.

(5)

y=Ax+r.

(6)

式中,r表示殘差矢量,服從零均值分布.

2 經(jīng)典相關(guān)干涉儀信號(hào)處理算法

為了計(jì)算測(cè)量矢量與數(shù)據(jù)庫(kù)之間的匹配度,經(jīng)典方法利用下式計(jì)算相關(guān)度作為x的估計(jì):

c=AHy=AHAx+AHr.

(7)

式(7)是測(cè)量矢量y的空域離散傅里葉變換的結(jié)果,(·)H表示共軛轉(zhuǎn)置.

通過(guò)計(jì)算矢量c的每個(gè)元素的幅度或者幅度的平方,得到譜線,然后視情況進(jìn)行插值操作,譜最大值對(duì)應(yīng)的角度即為信號(hào)入射角度.

因?yàn)樘炀€單元個(gè)數(shù)通常遠(yuǎn)遠(yuǎn)小于測(cè)向數(shù)據(jù)庫(kù)角度個(gè)數(shù)(具體分析詳見(jiàn)第3章),即

rank{A}≤min{N-1,L}=N-1?Q,

(8)

所以,

rank{AHA}≤rank{A}?Q,

(9)

故,

E{c}=AHAx+E{AHr}=AHAx‖/x.

(10)

式中,‖/表示不平行. 因此c既不是x的無(wú)偏估計(jì),也不是無(wú)偏估計(jì)的共線向量.

當(dāng)某一頻點(diǎn)上只有一個(gè)信號(hào)時(shí),矢量x僅有一個(gè)非零元素,不妨令其來(lái)自于第q1個(gè)校準(zhǔn)方向,則有

(11)

于是,

(12)

(13)

這表明,當(dāng)只有一個(gè)信號(hào)入射時(shí),矢量c可以正確指示入射方向.

當(dāng)存在多個(gè)信號(hào)時(shí),不妨令其來(lái)自于第qm個(gè)校準(zhǔn)方向,其中m=1,2,…,M,則有

(14)

由式(14)可知,第qm個(gè)校準(zhǔn)方向的譜值為

(15)

式中,c=[c1,…,cQ]T. 由式(15)可知,不能保證E{cqm}為局部極值,即存在多個(gè)信號(hào)時(shí),矢量c不能正確指示信號(hào)方向.

3 基于CS的相關(guān)干涉儀信號(hào)處理算法

對(duì)于典型的相關(guān)干涉儀系統(tǒng)來(lái)說(shuō),天線陣列通常含有5個(gè)或者9個(gè)天線單元,而校準(zhǔn)入射角度的個(gè)數(shù)則有上百個(gè),例如,當(dāng)步進(jìn)分別為1°和2°時(shí),分別有360和180個(gè). 這意味著,對(duì)于多達(dá)360或者180個(gè)未知數(shù),分別只有4個(gè)或者8個(gè)方程. 這說(shuō)明式(6)是一個(gè)欠定問(wèn)題. 同時(shí),這些未知數(shù)中大部分都是零. 由于具有相同頻率、相同時(shí)隙的信號(hào)只有少數(shù)幾個(gè),矢量x應(yīng)該只有少數(shù)幾個(gè)非零元素(對(duì)應(yīng)信號(hào)的入射角度). 這說(shuō)明式(6)是一個(gè)稀疏問(wèn)題.

鑒于式(6)是一個(gè)欠定、稀疏問(wèn)題,本文提出利用CS數(shù)學(xué)工具對(duì)式(6)進(jìn)行求解. CS的經(jīng)典形式是基于矢量x的l0范數(shù)(即非零元素個(gè)數(shù))約束設(shè)計(jì)的,該范數(shù)是理論上最好的稀疏性約束,其對(duì)應(yīng)的數(shù)學(xué)問(wèn)題如下[4]:

(16)

式中,‖·‖0表示l0范數(shù);‖·‖F(xiàn)表示Frobenius范數(shù);β表示殘差容限;C表示復(fù)數(shù)域. 然而,式(16)是NP-hard問(wèn)題,不易求解. 通常采用凸松弛(convex relaxation)技術(shù),將l0范數(shù)替換為l1范數(shù)(即所有元素的幅度之和)進(jìn)行求解,其對(duì)應(yīng)的數(shù)學(xué)問(wèn)題如下[4]:

(17)

式中,‖·‖1表示l1范數(shù). 由于l1范數(shù)中幅度較大的元素影響偏大,所以可以采用再加權(quán)(reweighted)技術(shù)來(lái)逼近理想的l0范數(shù). 首先將式(17)轉(zhuǎn)化為實(shí)數(shù)域問(wèn)題:

(18)

式中:

(19)

(20)

(21)

Re(·)表示實(shí)部;Im(·)表示虛部;R表示實(shí)數(shù)域.

(22)

式中,w∈R2Q×1為加權(quán)矢量. 設(shè)計(jì)加權(quán)矢量迭代過(guò)程如下[7-8]:

1) 初始化

利用式(7)設(shè)置加權(quán)矢量初始值,即

(23)

式中:ε為小的正常數(shù),以避免分母為零或者過(guò)于接近零;sgn(·)為符號(hào)函數(shù),即

(24)

2) 迭代更新

(25)

式中:(·)(i)為第i次迭代的結(jié)果;

(26)

式(25)中的優(yōu)化問(wèn)題可以利用壓縮感知計(jì)算軟件求解,例如CVX[9].

3) 迭代停止

達(dá)到最大迭代次數(shù),或者兩次迭代之間的代價(jià)函數(shù)差距小于預(yù)設(shè)的門(mén)限,則迭代停止.

在迭代完成之后,通過(guò)計(jì)算矢量x的每個(gè)元素的幅度或者幅度的平方,得到譜線,然后視情況進(jìn)行插值操作,譜極大值對(duì)應(yīng)的角度即為信號(hào)入射角度.

本文以雙通道相關(guān)干涉儀為例,提出了基于壓縮感知的相關(guān)干涉儀信號(hào)處理算法. 因?yàn)閺男盘?hào)處理的角度來(lái)看,多通道和雙通道相關(guān)干涉儀只在獲取復(fù)值相關(guān)矢量y的流程上有區(qū)別,而本文方法的輸入即為復(fù)值相關(guān)矢量y,所以本文方法對(duì)于多通道相關(guān)干涉儀也是適用的,不再贅述.

本文方法能夠同時(shí)估計(jì)的同頻信號(hào)個(gè)數(shù)受到了天線陣列中天線單元個(gè)數(shù)N的限制. 由于復(fù)相關(guān)矢量y的長(zhǎng)度為N-1,所以按照子空間原理,本文方法能夠同時(shí)估計(jì)的同頻信號(hào)個(gè)數(shù)最多為N-2個(gè).

4 計(jì)算機(jī)仿真

假設(shè)相關(guān)干涉儀天線陣列由5個(gè)天線單元組成;接收通道有兩個(gè),其一連接5號(hào)天線單元,其二通過(guò)開(kāi)關(guān)矩陣連接1～4號(hào)天線單元,天線半徑為750 mm,測(cè)向數(shù)據(jù)庫(kù)角度個(gè)數(shù)為360,角度步進(jìn)為1°,譜線未進(jìn)行插值操作.

圖1分別以700 MHz頻點(diǎn)和1 400 MHz頻點(diǎn)為例,給出了當(dāng)存在一個(gè)入射方向?yàn)?0°、信噪比為20 dB的信號(hào)時(shí),經(jīng)典方法(虛線)和本文方法(實(shí)線)的測(cè)向譜線. 從圖中可以看出,當(dāng)只有一個(gè)入射信號(hào)時(shí),經(jīng)典方法和本文方法均可以正確指示信號(hào)方向.

圖2分別以700 MHz頻點(diǎn)和1 400 MHz頻點(diǎn)為例,給出了當(dāng)存在三個(gè)入射方向分別為60°、180°和300°,信噪比均為20 dB的信號(hào)時(shí),經(jīng)典方法(虛線)和本文方法(實(shí)線)的測(cè)向譜線. 從圖中可以看出,當(dāng)只有多個(gè)入射信號(hào)時(shí),本文方法可以正確指示出三個(gè)同頻信號(hào)的入射方向,符合第3節(jié)中關(guān)于最大同頻信號(hào)個(gè)數(shù)的分析,即N-2=3個(gè)同頻信號(hào),而經(jīng)典方法失效了.

(a) 700 MHz頻點(diǎn)(a) 700 MHz

(b) 1 400 MHz頻點(diǎn)(b) 1 400 MHz圖1 存在一個(gè)同頻信號(hào)時(shí),經(jīng)典方法和本文方法譜圖對(duì)比Fig.1 Spectra of conventional and proposed methods in the presence of one signal which has exclusive use of the frequency

(a) 700 MHz頻點(diǎn)(a) 700 MHz

(b) 1 400 MHz頻點(diǎn)(b) 1 400 MHz圖2 存在三個(gè)同頻信號(hào)時(shí),經(jīng)典方法和本文方法譜圖對(duì)比Fig.2 Spectra of conventional and proposed methods in the presence of three signals which share the same frequency

由以上仿真實(shí)驗(yàn)可以看出,本文方法兼具異頻信號(hào)測(cè)向和同頻信號(hào)測(cè)向的能力.

接下來(lái)考察測(cè)向算法的精度,仿真中假設(shè)數(shù)據(jù)接收通道在被校準(zhǔn)后仍存在不超過(guò)±5°的隨機(jī)誤差. 表1給出了3 000次獨(dú)立重復(fù)實(shí)驗(yàn)統(tǒng)計(jì)出的測(cè)向誤差,其中,當(dāng)存在多個(gè)同頻信號(hào)時(shí),測(cè)向誤差給出了多個(gè)測(cè)向誤差的均值.

表1 測(cè)向誤差統(tǒng)計(jì)表Tab.1 Direction-finding accuracy

由表1可以看出,當(dāng)頻點(diǎn)上只存在一個(gè)信號(hào)時(shí),本文方法與經(jīng)典方法的測(cè)向精度類似;當(dāng)頻點(diǎn)上存在多個(gè)信號(hào)時(shí),經(jīng)典方法已經(jīng)失效,而本文方法仍可正確測(cè)向,只不過(guò)測(cè)向精度有所下降.

5 結(jié) 論

本文利用CS數(shù)學(xué)工具進(jìn)行相關(guān)干涉儀測(cè)向的信號(hào)處理,提出了一種高精度測(cè)向算法,解決了經(jīng)典方法不能對(duì)多個(gè)同頻信號(hào)進(jìn)行測(cè)向的問(wèn)題. 需要指出的是,本文方法涉及的數(shù)學(xué)運(yùn)算邏輯復(fù)雜、運(yùn)算量較大,下一步需要從算法工程化的角度做進(jìn)一步的研究.