基于改進OMP算法的稀疏目標微波關聯成像方法

張 瑞, 全英匯, 朱圣棋, 李亞超, 邢孟道

(1. 西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071;2. 西安電子科技大學電子工程學院, 陜西 西安 710071)

0 引 言

由于成像雷達具有全天時、全天候、抗干擾能力強等顯著優勢,在災害監測、對地觀測等方面具有廣泛應用,已逐步成為獲取目標精細特性的遙感設備[1]。基于逆合成孔徑或合成孔徑技術的雷達設備,其成像原理均依賴于目標與雷達間的相對運動關系,通過發射大帶寬信號獲得距離向的高分辨率,方位向分辨率取決于合成孔徑長度[2-4]。然而,在某些凝視成像或前視成像的場合中,由于在方位向上雷達與觀測目標不存在相對運動關系或合成孔徑陣列與視線方向相同,不能實現方位向成像,難以獲得高分辨圖像。

近年來,起源于光學關聯成像技術[5]的雷達微波關聯成像技術受到海內外越來越多學者的關注。雷達發射端通過發射隨機調制的探測信號,使其在空間中疊加產生時間—空間二維隨機波動的輻射場,該輻射場與目標相互作用后,信號回波在接收端與事先存儲的參考輻射場模板進行關聯處理,實現對波束內目標信息的解耦與提取[6-10]。

雷達微波關聯成像技術不依賴于雷達與目標間的相對運動關系,與傳統的高分辨雷達體制形成互補,可工作于前視模式或凝視觀測模式,將其應用在靜止或準靜止平臺上可實現對地的長時間連續監測、資源勘查等功能,具有廣泛的應用前景。

目前,有關微波關聯成像技術的研究主要集中在分辨力理論、空時二維隨機輻射輻射源設計方法[11-14]、關聯重構算法[15-17]、三維成像方法[18]、運動目標關聯成像[19-21]等方面。其中,由于微波關聯成像模型與壓縮感知(compressed sensing, CS)存在天然一致性,對稀疏目標成像時,能夠取得較好的成像效果,因此稀疏重構類成像算法成為研究熱點。國防科技大學、西安電子科技大學、中國科學技術大學等單位對微波關聯成像技術相繼開展研究,并取得一系列研究成果[22-26]。

針對關聯重構算法——正交匹配追蹤(orthogonal matching pursuit, OMP)算法在每次迭代時均需求解最小二乘解而導致計算量大的問題,本文提出一種改進OMP算法的稀疏目標微波關聯成像方法,該方法可獲得良好的成像效果并且降低了計算量,有利于算法的實時化處理。本文首先闡述了雷達微波關聯成像機理,構建了成像數學模型;然后結合共軛梯度法推導了改進OMP算法公式,給出了算法實現步驟;最后通過仿真對比實驗驗證了所提算法的有效性。

1 基于隨機輻射場的成像機理與信號模型

1.1 成像機理描述

基于隨機輻射場的雷達前視微波關聯成像技術由經典熱光源強度關聯成像技術發展而來。經典熱光源強度關聯成像技術通過使用熱光源在特定的成像平面上產生強度隨機起伏的光場來探測目標信息,并且將接收到的探測通道和參考通道中隨機起伏的光場強度進行關聯得到目標物體的像。

雷達微波關聯成像技術采用微波作為信號源,通過發射各種帶有調制的電磁波信號(如隨機相位調制信號),在感興趣的觀測空間中人為地構造空時二維隨機輻射場來模擬經典熱光源強度關聯成像中隨機起伏的光場。微波信號源發射隨機調制的信號,并且該信號在空間中非相干疊加,形成分別在空間維度上和時間維度上均具有非相干特性的場。該場與觀測區域的目標相互作用后,使得回波中攜帶有目標對不同輻射場分布的調制信息,將這些回波信號與參考隨機輻射場進行關聯處理,得到探測平面內目標的成像結果。微波關聯成像原理示意如圖1所示。

圖1 微波關聯成像原理示意圖Fig.1 Schematic diagram of microwave correlation imaging

輻射場可用U(r)來表示,其空間相干性的定義為

(1)

式中:r1,r2表示波束覆蓋區域內不同的兩點;R(r1,r2,t)表示隨機輻射場在時刻t的空間自相關函數。

輻射場U(r)的時間相干性定義為

(2)

式中：t1,t2表示兩個不同時刻;R(r,t1,t2)表示隨機輻射場在位置r處的時間自相關函數。

理想情況下,空時二維隨機輻射場表現為理想的非相關特性,即

(3)

(4)

在波束覆蓋范圍內,同一時刻t時不同分辨單元r1、r2的電磁輻射場U(r1,t)、U(r2,t)具有非相干差異性分布特征,使得不同分辨單元被具有不同空間分布特征的輻射場所標度;同一分辨單元r不同時刻t1、t2的電磁輻射場U(r,t1)和U(r,t2)也具有非相干差異性分布特征,使得目標在不同時刻受到具有不同時間分布特征的輻射場所照射。

傳統的成像雷達如實孔徑雷達,由于發射天線發射單一形式的信號波形,其在探測空間中形成的是起伏規律的電磁輻射場,空間中不同位置處的輻射場相干性較強。一方面,分辨率受實孔徑天線尺寸的限制,另一方面,位于同一天線波束范圍內的不同目標被相干性較強的輻射場所標度,輻射場不能為目標分辨帶來額外的有用信息,因此雷達分辨率較低。此外,發射天線通常發射的是具有周期規律的雷達信號,該信號在觀測空間中形成的電磁輻射場在空間上的分布并不隨時間的變化而變化,導致雷達回波信號在不同時刻也是相干的,增加照射時間導致信噪比(signal to noise ratio, SNR)提高,對雷達分辨率的提高無任何實質性幫助。而微波關聯成像的信號源通過發射隨機調制的信號,形成具有空間和時間二維非相干特性的隨機輻射場,為雷達分辨同一波束內的不同目標提供額外信息。若能充分利用空時二維隨機輻射場的非相干特性,再結合相應的成像算法,則可使得目標超分辨成像成為可能。

1.2 雷達微波關聯成像數學模型

假設雷達陣列有M個發射陣元和N個接收陣元。發射陣元排列成均勻線陣,且發射獨立的隨機信號。在目標成像平面建立直角坐標系,并將成像平面離散劃分為Q個相同大小的網格,沿X軸方向的網格尺寸為ΔX,沿Y軸方向的網格尺寸為ΔY。假設目標散射點位于成像網格中心,并以此時散射點的位置矢量r表示該網格的位置。

將目標的等效后向散射系數矢量表示為

σ=[σ1,σ2,…,σq]T

(5)

式中:σq表示第q個成像網格中心rq處等效后向散射系數,若該成像網格單元沒有散射點,則σq=0。

為簡化起見,現考慮具有M個發射陣元和1個接收陣元的多發單收天線結構。接收天線選定為發射天線中的參考天線(即第1個發射天線)。如此,第m(m=1,2,…,M)個發射陣元在第l(l=1,2,…,L)個脈沖內的發射信號為

(6)

(7)

信號經過調制后,每個陣元的發射信號頻率在脈沖之間隨機捷變。一方面,這種隨機調制的信號在空間中非相干疊加,形成具有空間和時間二維非相干特性的隨機輻射場,為后續微波關聯成像處理提供重要基礎。另一方面,這種形式信號的抗干擾能力較強,特別是對基于數字射頻存儲器(digital radio frequency memory, DRFM) 的轉發式干擾抑制方面效果顯著[28]。一般而言,對方干擾機轉發的是上一個或上幾個脈沖周期的信號波形,對于當前脈沖周期的信號來不及截獲、存儲、復制和轉發。而我方雷達接收機以當前脈沖周期的信號作為參考進行脈沖壓縮。由于發射信號在每個脈沖重復周期內都是隨機捷變的,并且這種變化只由信號發射方掌握,干擾方只能被動轉發,因此干擾信號經脈沖壓縮時由于失配而被成功抑制,只有真目標才能被正確檢測。

為方便公式推導,假設空間中存在1個點目標,其等效后向散射系數記為σ0,當存在多個點目標時,只需在公式上稍作修改即可,不影響公式推導。

在第l個脈沖內,接收陣元接收到的該單個點目標的回波信號為

(8)

(9)

式中：c表示光速;r0表示等效后向散射系數為σ0的點目標在成像平面坐標系中的位置矢量;Rt,m和Rr,1分別表示第m個發射陣元和接收陣元在成像平面坐標系中的位置矢量。

記第l個脈沖內,M個發射陣元的頻率分別為

(10)

將第l個脈沖內接收到的目標回波信號進行下變頻及脈沖壓縮處理后的信號可以寫作:

(11)

式中:B為信號帶寬。

在成像平面內,圍繞目標劃分出Q個大小相同的成像網格,網格的位置矢量由其中心點的坐標矢量表示。由于已知各個陣元的位置坐標(xarray(m),yarray(m),zarray(m))和各個成像網格的位置矢量(xq,yq,zq),于是可計算得到第m個陣元到第q個網格的距離為

(12)

因此,由第m個陣元發射的信號經第q個網格(σq=1)反射后,被接收陣元接收而產生的時延可表示為

(13)

這樣,第l個脈沖內第q個網格處(即rq處)的參考信號可寫為

(14)

因此,雷達接收端接收到來自第l個脈沖內第q個網格的回波信號可表示為

(15)

將成像模型寫為矩陣形式:

y=Sσ

(16)

式中：y=[y(t1),y(t2),…,y(tL)]T為回波信號向量;σ=[σ1,σ2,…,σQ]T為待恢復的目標散射系數向量;S為采樣矩陣,其具體形式為

(17)

從式(17)可以看出,觀測矩陣不僅與雷達發射、接收陣列幾何構型有關,還與雷達信號波形、成像網格數、脈沖數有關。此外,由式(16)還可以看出,解雷達微波關聯成像方程是一個典型的線性逆問題求解過程,利用方程的解可以恢復目標、重構目標圖像。

如果從另外一個角度來解釋待恢復目標與空時二維隨機輻射場間的關系,可以認為空時二維非相干隨機輻射場的差異性分布特征相當于空間波前的隨機編碼,隨機輻射場與目標相互作用后,目標的信息就蘊含在散射回波信號中。這種相互作用等于對感興趣的目標進行采樣或測量。從這個方面來看,式(17)描述的采樣矩陣S與稀疏重構理論里的采樣矩陣或者觀測矩陣類似。然后,利用本地預先存儲的采樣矩陣和接收到的目標采樣數據進行恢復和反演就可以得到原始目標圖像。這種關于目標與隨機輻射場間關系的解釋與CS理論不謀而合。

對于微波關聯成像,空時二維完全非相干隨機輻射場天然地滿足CS中對觀測矩陣的要求:限制等距特性。此外,雷達目標回波信號可以看作是少數幾個強散射點回波信號疊加的結果[29]。目標的稀疏特性是雷達成像中比較常用的先驗信息。通常情況下,所觀測的目標本身就具有稀疏性,此時CS中字典矩陣是一個單位陣,這樣采樣矩陣就與CS中的測量矩陣對應起來。因此,式(16)所示的成像模型就是CS模型。另一方面,即使當前目標不具有稀疏特性,通常利用各種稀疏變換手段(例如小波變換)便可快速得到其稀疏表示。結合本文中的應用場景,目標飛機在空域中可以看作是由K個主要散射點組成的,在構造成像模型時,需要將目標飛機所在成像平面劃分成Q個成像網格,通常K?Q。因此,飛機在空域上滿足稀疏條件。稀疏度的定義可以沿用壓縮感知理論中對于稀疏度的描述。如向量σ=[σ1,σ2,…,σQ]T中只有少于K(K?Q)個非零元素,稱K為信號的稀疏度。因此,可以將CS相關理論用于微波關聯成像處理中重構恢復目標圖像。

由于對成像平面進行網格劃分時,網格數量多且尺寸小,接收信號相干性較強。因此,求解式(16)是一個欠定問題,且解通常不唯一。為了得到符合要求的唯一解,這里采用稀疏恢復算法求解。

已有文獻證明[30],若σ稀疏,求解式(16)可轉化為以下優化問題:

(18)

式中：η表示誤差允許量；‖·‖k表示向量的lk范數。但是,求解式(18)是一個非確定性多項式難題(nondeterministic polynomially problem， NP)。最常見的一種解決方法是用范數代替l0范數:

(19)

可以證明,若σ是稀疏的且測量矩陣S滿足一定條件時,式(18)和式(19)具有相同的解。

2 微波關聯成像算法

從前面的分析可以看出,求解式(16)是一個病態的逆問題。導致這個逆問題病態的原因是測量數據信息量不足,實際目標場景中蘊含的信息量遠遠多于有限的測量數據量,僅僅采用諸如最小二乘法是無法求解或者無法獲得穩定解的。因此,考慮利用目標稀疏的先驗信息,結合稀疏重構類算法以獲得方程穩定的解,從而恢復目標圖像。

2.1 OMP算法

考慮采用OMP算法求解上述問題,其基本思想是從字典矩陣中,選擇與輸入信號相關性最強的原子,形成一個稀疏逼近,同時求出匹配殘差,然后繼續選擇與殘差最匹配的原子,反復迭代直到滿足循環終止條件為止。需要稀疏恢復的信號可以由迭代中選擇相關性最強原子的線性和再加上最后的殘差值來表示。如果殘差值在可以接受范圍內,則信號就是這些原子的線性組合。

對于壓縮觀測:

y=Φx

(20)

式中:yM×1為觀測所得向量;ΦM×N為測量矩陣;xN×1為原信號,其在某個變換域ΨN×N是稀疏的,即

x=Ψθ

(21)

式中:θN×1是K稀疏的,一般有K?M?N。

令傳感矩陣A為

A=ΦΨ

(22)

則

y=Aθ

(23)

OMP重構算法的具體步驟如下。

輸入:傳感矩陣A=ΦΨ,觀測向量y,信號的稀疏度K。

步驟 1初始化Λ0=?,A0=?,r0=y,i=1。

步驟 3令Λi=Λi-1∪{λi},Ai=Ai-1∪aλi。

步驟 4計算y=Aiθi的最小二乘解:

步驟 5更新殘差

步驟 6迭代次數i=i+1,并判斷i≤K是否滿足,若滿足則返回步驟2,否則終止迭代。

其中:i表示迭代數目;ri表示殘差值;Λi表示第i次迭代索引構成的集合;λi表示第i次迭代中得到的索引;θi為i×1維列向量;aj表示矩陣A的第j列構成的向量;Ai為依照集合Λi選出的矩陣A的列構成的集合。

在利用OMP算法求解式(19)過程中,涉及大量向量與向量相乘、矩陣與向量相乘的運算,而且每迭代一次就要求解一次最小二乘問題,利用直接矩陣求逆法求最小二乘解時,非常耗費硬件空間、時間資源,這降低了算法的實時性,增加了算法的實現難度。因此,優化OMP算法以利于實時處理是非常有必要的。

2.2 改進OMP算法的稀疏恢復算法

觀察第2.1節中OMP算法步驟4可知,每一次迭代都要利用直接矩陣求逆法來求解下式的最小二乘解:

y=Aiθi

(24)

式中:Ai為M×i維矩陣,表示按索引Λi選出的矩陣A的列集合;θi為i×1維的待重建信號。式(24)是一個超定方程,這里只能求得其最優解。定義最優解為使殘差y-Aiθi的l2范數取得極小值時的解為

(25)

對式(25)關于θi求梯度并令梯度值為0可得

(26)

進一步化簡可得

(27)

Gθ=b

(28)

當矩陣G維數較大時,直接矩陣求逆法的計算量會相當大,這將耗費巨大的時間和空間資源。因此,本文采用共軛梯度法來求解式(28),從而避免了矩陣求逆運算,改進了原始的OMP算法。共軛梯度法只涉及少量的矩陣與向量相乘、向量與向量相乘運算,從而大大降低計算數據量,節約計算資源和空間資源,增強整個成像恢復算法的實用性。

共軛梯度法是用來解決無約束凸二次規劃問題的一種方法。其基本思想是對于一個目標函數,找到一組方向向量,依次按此方向向量中的方向對迭代點進行更新,使得函數值在該方向上取得最小值。在每一個新的更新方向對迭代點進行更新時,不會影響在之前方向上的更新結果。如果能找到這樣一組方向向量,那么就可保證在多次迭代后找到函數的全局極小點。

共軛梯度法是介于最速下降法與牛頓法之間的方法,既克服了最速下降法收斂慢的缺點,又避免了存儲和計算牛頓法所需要的二階導數信息。該算法具有收斂速度快、所需存儲空間少以及二次終止等優點,因此在解決實際問題中得到了廣泛的應用。

將式(28)轉化為如下問題:

(29)

要應用共軛梯度法,最關鍵的是要確定算法中的兩個參量:優化方向p和優化步長ρ,下面給出其推導過程。

(30)

在上式兩邊同時左乘p(k)TG,得

(31)

解得

(32)

已知迭代點θ(k)和優化方向p(k),利用一維搜索確定優化步長ρk,也就是求解下列問題:

(33)

為了求解式(31),這里記:

φ(ρ)=f(θ(k)+ρp(k))

(34)

[G(θ(k)+ρp(k))+b]Tp(k)=0

(35)

(36)

求解式(36)可獲得第k次的優化步長ρk,其表達式為

(37)

利用共軛梯度法求解式(28)的詳細步驟如下所示。

算法 1 共軛梯度法輸入:對稱正定矩陣G,等式右端項b,初始值θ0;輸出:方程的解θ;計算殘差:r0=b-Gθ0;初始化:p0=r0;repeat for k=0,1,… do ifpk=0then return θ0 else αk=rTkrkpTkGpk; θk+1=θk+αkpk; rk+1=rk-αkGpk; βk=rTk+1rk+1rTkrk; pk+1=rk+1+βkpk; endend

2.3 算法計算復雜度分析與對比

在OMP算法迭代過程中,相比于直接矩陣求逆,采用共軛梯度法求最小二乘解能夠使得問題求解復雜度明顯降低,大大節省計算時間,提高算法的有效性和可靠性。下面詳細分析對比利用直接矩陣求逆法和利用共軛梯度法求最小二乘解的計算復雜度。

在第2.1節OMP重構算法的步驟4中,利用直接矩陣求逆法求得y=Aiθi的最小二乘解為

(38)

通過以上分析可以看出,由于利用直接矩陣求逆來求解最小二乘解時,不僅要頻繁進行矩陣與矩陣相乘操作,還要進行矩陣求逆運算,其計算量隨信號稀疏度的增加而快速增長。此外,在硬件實現上需要對矩陣、向量頻繁地進行存取,耗費了大量的時間資源,不利于硬件電路性能的提升和整個成像稀疏恢復算法的實時化處理。

共軛梯度法只涉及少量的矩陣與向量相乘、向量與向量相乘、標量間運算,數據量小且數據間的一致性較強,這樣就節約了計算時間資源,利于算法實現和實時性處理。

3 仿真結果及分析

根據式(14),可將常用的微波關聯成像方法分為匹配濾波法、最小二乘法、稀疏重構法。本小節通過數值仿真對比這三種方法的性能。表1給出了具體的仿真參數取值。

表1 參數設置

仿真場景如圖2所示,圖2(a)是一個由3個散射點組成的簡單目標,圖2(b)是一個由8個散射點組成的飛機目標,兩種目標都是稀疏的,但稀疏度不同。

圖2 原始目標Fig.2 Original target

針對圖2(a)的目標場景,分別利用匹配濾波法、最小二乘法、傳統OMP算法和改進OMP稀疏重構算法進行仿真,其中SNR=5 dB。

由于匹配濾波法是將回波信號y與采樣矩陣S進行相關處理,該方法利用了空時二維隨機輻射場的非相干特性,因此回波信號中與采樣矩陣相關程度較高的列會在圖中產生尖峰。但實際中形成的隨機輻射場不能滿足理想的時空完全非相干特性,這會大大影響成像質量。如圖3所示,圖中紅色圓圈標出的點即為目標散射點,目標勉強可見。

圖3 匹配濾波法仿真結果(3個點目標)Fig.3 Simulation result of matched filtering method (three point targets)

圖4 最小二乘法仿真結果(3個點目標)Fig.4 Simulation result of least square method (three point targets)

從現有的仿真結果質量上來看,傳統的OMP算法和改進的OMP算法在重構結果的質量方面差異性不大，如圖5和圖6所示，但從算法處理用時、硬件實現難度以及算法實時性上來看,本文所提方法在運算量上具有巨大優勢,有利于算法地實現和實時性處理，4種算法進行目標恢復成像的用時如表2所示。

圖5 傳統OMP稀疏重構結果(3個點目標)Fig.5 Sparse reconstruction result of traditional OMP method (three point targets)

圖6 改進OMP稀疏重構結果(3個點目標)Fig.6 Sparse reconstruction result of improved OMP method (three point targets)

表2 4種算法運行時間對比

由于利用傳統OMP算法求解成像方程時,需要涉及大量的向量與向量相乘、矩陣與向量相乘運算,而且每迭代一次就要求解一次最小二乘問題,非常耗費硬件空間、時間資源,降低了算法的實時性,增加了算法的實現難度。此外,由于改進的OMP算法增加了稀疏先驗信息的約束,所以削弱了逆問題求解過程中的病態性,降低了采樣矩陣S病態性造成的成像不穩定,提高了雷達成像質量。

針對圖2(b)的飛機目標,分別利用匹配濾波法、最小二乘法、傳統OMP算法和改進OMP稀疏重構算法進行仿真。其中SNR=5 dB。

如圖7所示,由仿真結果很難分辨出飛機目標。這是由于匹配濾波法是將回波信號y與采樣矩陣S進行相關處理,其本質還是利用隨機輻射場的時空非相干特性。若采樣矩陣S中各列的相關性很強,則經相關處理后,圖上會出現很多幅度值較大的點,導致原始目標被其淹沒,難以分辨。

圖7 匹配濾波法仿真結果(飛機目標)Fig.7 Simulation result of matched filtering method (aircraft target)

最小二乘法的仿真結果如圖8所示。由于最小二乘法對噪聲和誤差極其敏感,且輻射場的時空相干性較強,所以該算法不能正確得到飛機目標圖像,成像失敗。

圖8 最小二乘法仿真結果(飛機目標)Fig.8 Simulation result of least square method (aircraft target)

圖9和圖10分別給出了利用傳統OMP算法和所提的改進OMP稀疏重構算法的仿真結果。兩種算法在有噪聲的情況下能夠準確恢復目標,得到正確的目標圖像。但所提算法優化了傳統OMP算法以利于實時處理,改進后的OMP算法數據量小且數據間的一致性較強,節約了計算時間資源,有利于算法硬件實現和實時性處理。

圖9 傳統OMP稀疏重構結果(飛機目標)Fig.9 Sparse reconstruction result of traditional OMP method (aircraft target)

圖10 改進OMP稀疏重構結果(飛機目標)Fig.10 Sparse reconstruction result of improved OMP method (aircraft target)

由于實現匹配濾波法只需進行矩陣的相關運算,因此耗時最短。但如前面所述,該方法的性能受眾多因素影響,其成像效果差。改進OMP算法耗時與匹配濾波法相當,但不論是在有噪聲的情況下還是當采樣矩陣S的非相干性不理想時,其成像結果均令人滿意,為3種算法中最優異的一種。最小二乘法耗時最長,原因在于該方法計算過程中需要進行大型矩陣與矩陣相乘操作以及矩陣求逆運算,計算量隨著網格劃分數目、脈沖數、信號稀疏度的增加而迅速增加,因此非常耗時。此外,由于受到噪聲和采樣矩陣S非相干性不理想等因素的影響,該方法的成像效果非常差。傳統OMP算法步驟涉及大量的向量與向量相乘、矩陣與向量相乘運算,而且每迭代一次就要求解一次最小二乘問題,利用直接矩陣求逆法求最小二乘解時,非常耗費硬件空間與時間資源。

4 結 論

針對傳統OMP算法在每次迭代時,均需利用直接矩陣求逆法求解目標函數的最小二乘解而導致計算量大的問題,本文提出一種改進OMP算法的稀疏目標微波關聯成像方法。首先,闡明了微波關聯成像原理并構建了微波關聯成像數學模型;然后,利用共軛梯度法對OMP算法中的直接矩陣求逆步驟進行了改進,并分析了改進后算法的計算量;最后,通過仿真實驗驗證了所提算法的有效性。實驗結果表明,與匹配濾波法、最小二乘法、傳統OMP算法相比,所提算法能夠獲得良好的成像效果,降低了算法計算量,有利于算法的實時化處理。