朱曉秀, 劉利民, 胡文華, 郭寶鋒, 史 林, 朱瀚神
(1. 陸軍工程大學石家莊校區(qū)電子與光學工程系, 河北 石家莊 050003;2. 中國人民解放軍32398部隊, 北京 100192)
對于傳統(tǒng)的單基逆合成孔徑雷達(inverse synthetic aperture radar, ISAR)成像系統(tǒng)[1-5]而言,圖像的距離分辨率和方位分辨率分別受雷達的發(fā)射信號帶寬和觀測累積角度的限制,因而無法突破其理論分辨率[6-10]。為提高二維分辨率,多視角多頻帶ISAR融合成像技術(shù)利用工作在不同頻段的多部雷達從不同角度觀測目標,通過對觀測回波進行融合處理,等效得到一個更大帶寬和更大視角的回波[11-13]。該方法打破了單基地雷達成像分辨率的約束,可同時提高ISAR成像的二維分辨率。
多視角多頻帶ISAR融合成像方法主要基于譜估計和基于稀疏表示兩大類。譜估計類融合成像方法[14-17]將ISAR融合成像模型轉(zhuǎn)化為二維指數(shù)和模型,并分解為兩個一維矢量,利用譜估計方法分別估計參數(shù)。此類方法在散射點個數(shù)精確已知的情況下可較好地實現(xiàn)信號融合,但實際情況下散射點個數(shù)一般很難準確估計,而且還存在散射點的二維坐標配對問題,影響了算法性能。稀疏表示類融合成像方法[18-20]利用成像場景網(wǎng)格化處理模擬散射點可能出現(xiàn)的位置,將ISAR融合成像問題轉(zhuǎn)化為信號的稀疏表示問題,利用稀疏重構(gòu)方法進行求解。此類方法無需估計散射點個數(shù),也不涉及散射點坐標配對問題,算法性能優(yōu)于譜估計類方法。
在高頻區(qū),目標的電磁散射可近似等效為多個獨立散射點的后向電磁散射之和。在窄帶小角度觀測條件下,散射點的散射系數(shù)可認為是常數(shù),近似為理想散射點模型。文獻[21]利用理想散射點模型建立ISAR成像回波模型,提出了一種基于稀疏表示的多雷達信號二維融合成像方法,利用基追蹤算法求解稀疏表示問題,得到的融合成像結(jié)果要優(yōu)于譜估計方法。然而,在寬帶小角度觀測條件下,散射點的散射系數(shù)是隨頻率變化而變化的。為更準確地反映目標的散射特性,文獻[22]提出了一種基于幾何繞射理論(geometrical theory of diffraction, GTD)模型的多雷達信號二維融合方法,將多雷達信號的二維融合問題轉(zhuǎn)化為稀疏表示問題,利用正則化方法能較好地估計目標的散射參數(shù)。然而,在建立多視角多頻帶ISAR融合成像的稀疏表示模型時,需要將二維回波數(shù)據(jù)矢量化處理,此時涉及到大規(guī)模數(shù)據(jù)重構(gòu),利用正則化方法迭代求解時運算量大,耗時較長,而且需要人工調(diào)整正則參數(shù),影響了算法性能。
基于此,本文提出一種新的基于GTD模型的多視角多頻帶ISAR融合成像方法,在利用GTD模型建立ISAR成像回波模型的基礎(chǔ)上,采用矢量化處理,將多視角多頻帶ISAR融合成像問題轉(zhuǎn)化為信號稀疏重構(gòu)問題。為減少運算復雜度,利用正交匹配追蹤(orthogonal matching pursuit, OMP)算法實現(xiàn)稀疏重構(gòu),在保證成像質(zhì)量的同時提高融合成像效率,仿真實驗驗證了算法的有效性。
假設雷達發(fā)射LFM信號,在成像時間Ta內(nèi)共發(fā)射N個脈沖,脈沖重復時間為Tr,tn=nTr為慢時間,其中n=0,1,…,N-1。若目標共有P個獨立散射點,經(jīng)解線頻調(diào)處理后,基于理想散射點模型目標回波可寫為
(1)
(2)
在寬帶成像時,散射點的散射系數(shù)隨頻率變化而變化,此時利用理想散射點模型不能準確描述目標的散射特性。為考慮散射系數(shù)隨頻率變化的特性,引入頻率依賴因子αp,在式(2)的基礎(chǔ)上,利用GTD模型[23]構(gòu)建目標回波,可表示為
(3)
式中:f0為頻帶的起始頻率;αp為散射點p的頻率依賴因子。在GTD模型中,不同的αp可表征不同的散射點類型,典型的散射點類型及對應的頻率依賴因子如表1所示[24]。
表1 典型的散射點類型及對應的頻率依賴因子
在小角度觀測條件下,有cos Δθ(tn)≈1,sin Δθ(tn)≈Δθn。經(jīng)運動補償后,目標運動模型可近似為轉(zhuǎn)臺模型,假設勻速轉(zhuǎn)動的角速度為ω,則有Δθ(tn)=ωtn,此時ΔRp(tn)可近似寫為ΔRp(tn)≈yp+xpωtn,式(3)可近似表示為
(4)
將距離頻率離散化,令fc+f=f0+mΔf,Δf為頻率采樣間隔,M為頻率采樣點數(shù),m=0,1,…,M-1。此時,式(4)可寫為
(5)
在成像時間較短且目標尺寸較小時,可忽略越分辨單元徙動(migration through resolution cells, MTRC)的影響,若在成像過程中存在MTRC現(xiàn)象,可利用文獻[16]中的方法進行MTRC校正,校正后的目標回波可表示為
(6)
利用網(wǎng)格化處理思想,分別令2Δfyp/c=k/K和2f0xpωTr/c=l/L,其中,k=0,1,…,K-1,l=0,1,…,L-1,則ISAR成像場景可離散化為大小為L×K的二維網(wǎng)格,其中距離維有K(K≥M)個網(wǎng)格,方位維有L(L≥N)個網(wǎng)格,此時,式(6)可表示為
(7)
式中:αk l表示網(wǎng)格點(k,l)上散射點的頻率依賴因子;σk l表示網(wǎng)格點(k,l)上散射點的散射系數(shù)常數(shù)復幅度,若網(wǎng)格點(k,l)上存在真實散射點,則σk l不為零,若網(wǎng)格點上不存在真實散射點,則σk l為零。
從式(6)中可以看出,由于頻率依賴因子的存在,使得回波信號中的幅度與相位是耦合的。為方便對耦合的二維回波信號進行稀疏表示,將目標回波沿視角方向進行堆疊,則式(7)可矢量化[22,25-27]表示為
(8)
式中:s為矢量化處理后的回波矢量,大小為MN×1維,可寫為
s=[S(0,0),…,S(M-1,0),…,S(0,N-1),…,S(M-1,N-1)]T
(9)
(10)
(11)
式中:I表示單位矩陣;?表示矩陣的Kronecker乘積;Ti可表示為
(12)
W為MN×KL維的字典矩陣,可表示為
(13)
其中,
Wn=[Φ(0,0),…,Φ(K-1,0),…,Φ(0,L-1),…,Φ(K-1,L-1)]M×KL
(14)
(15)
本文以工作在不同頻帶且從不同視角觀測目標的兩部雷達為例進行多視角多頻帶ISAR融合成像分析。為保證實現(xiàn)融合成像,兩部雷達觀測目標時獲得同一散射中心的信息差距相差不能太大。假設兩部雷達相近放置,對目標的觀測視角無重合,雷達間的回波信號已經(jīng)過運動補償以及互相干等預處理。若在全視角全頻帶回波中存在因外界干擾或數(shù)據(jù)傳輸錯誤等情況,導致雷達信號中的某一觀測頻段或觀測視角的回波缺失或受到污染,不能用于成像,此時需要利用有效的觀測數(shù)據(jù)實現(xiàn)成像。假設在全視角全頻帶中有效的觀測數(shù)據(jù)為雷達1和雷達2的回波,通過多視角多頻帶ISAR融合成像方法可補全缺失的回波數(shù)據(jù),與單部雷達相比,可等效提高觀測視角和發(fā)射帶寬,從而改善成像的二維分辨率。
假設Δf為頻率采樣間隔,雷達1和雷達2的發(fā)射信號頻帶分別為f=f0,f1,…,fM1-1和f=fM-M2,fM-M2+1,…,fM-1,分別包含M1個和M2個頻率采樣點。全頻帶可表示為f=f0,f1,…,fM-1,共包含M個頻率采樣點,則全頻帶的頻率采樣數(shù)據(jù)可表示為f=f0+mΔf(m=0,1,…,M-1)。假設Δθ為角度采樣間隔,雷達1和雷達2的觀測角度分別為θ=θ0,θ1,…,θN1-1和θ=θN-N2,θN-N2+1,…,θN-1,分別對應N1個和N2個角度。全視角可表示為θ=θ0,θ1,…,θN-1,共有N個角度,全視角的角度采樣數(shù)據(jù)可表示為θ=θ0+nΔθ(n=0,1,…,N-1)。
(16)
(17)
(18)
式中:I和0分別表示單位矩陣和零矩陣。
圖1 基于矢量化處理的多視角多頻帶雙雷達觀測數(shù)據(jù)融合示意圖Fig.1 Multi-angle and multi-band dual radar observation data fusion based on vectorization processing
從式(16)可以看出,矢量化后的融合成像模型涉及到大規(guī)模的數(shù)據(jù)重構(gòu),因此需要尋找快速有效、簡單穩(wěn)定的稀疏重構(gòu)算法實現(xiàn)多視角多頻帶ISAR融合成像。OMP算法[28]原理簡單,易于實現(xiàn),且運算復雜度較低,是一種常用的稀疏重構(gòu)算法,故本文采用OMP算法求解式(16)。OMP算法的主要實現(xiàn)流程如圖2所示。
圖2 OMP算法實現(xiàn)流程Fig.2 Implementation process of OMP algorithm
OMP算法具體的實現(xiàn)步驟如下:
在迭代過程中,OMP算法利用了正交化思想,避免了重復選擇原子。對于終止迭代條件,若信號稀疏度已知,則當?shù)螖?shù)為稀疏度大小時終止迭代,若稀疏度未知,則當殘差小于預設門限值時終止迭代。
基于OMP算法的多視角多頻帶ISAR融合成像方法實現(xiàn)流程如圖3所示,主要步驟如下:
步驟 1基于GTD模型得到各雷達回波,進行運動補償及互相干等預處理,得到相干的距離頻域-方位慢時間域回波信號;
步驟 2構(gòu)造字典矩陣,將回波信號離散化;
圖3 多視角多頻帶ISAR融合成像流程Fig.3 Multi-angle and multi-band ISAR fusion imaging process
本文仿真實驗環(huán)境為Windows 10 64位操作系統(tǒng),Matlab 2016A軟件平臺,仿真所用計算機主要參數(shù)如下:處理器為Intel酷睿i5-8265U,主頻為1.60 GHz和1.80 GHz,內(nèi)存為16 GB。本節(jié)分別利用簡單散射點模型和復雜飛機模型的融合成像仿真實驗驗證所提融合成像算法的有效性。
本仿真實驗利用簡單散射點模型為成像目標進行多視角多頻帶ISAR融合成像,目標中包含6個散射點,仿真模型如圖4(a)所示,其中散射點A和B的坐標分別為(-0.1,-0.1)和(0.1,0.1),對應的頻率依賴因子分別為0.5和-1,其余散射點對應的頻率依賴因子均為0。雷達工作頻帶為17~18.5 GHz,采樣頻率為1.92 GHz,發(fā)射脈沖寬度為6 μs,PRF為50 Hz,共發(fā)射90個脈沖。在成像時間內(nèi),假設目標以0.04 rad/s的角速度勻速轉(zhuǎn)動,累積轉(zhuǎn)角為4.08°,此時雷達的距離分辨率和方位分辨率理論值分別為0.10 m和0.12 m。每個脈沖回波內(nèi)采樣75個距離采樣單元的回波數(shù)據(jù),經(jīng)運動補償后作為全視角全頻帶回波數(shù)據(jù),大小為75×90。利用距離-多普勒(range-Doppler, RD)算法得到的全視角全頻帶ISAR成像結(jié)果如圖4(b)所示,由于分辨率足夠,從成像結(jié)果中可以完全區(qū)分所有散射點。
圖4 目標模型及其全視角全頻帶RD成像結(jié)果Fig.4 Target model and the RD imaging result of the full-angle and full-band echo data
在全視角全頻帶回波中取左上角大小為M1×N1的數(shù)據(jù)作為雷達1的觀測數(shù)據(jù),取右下角大小為M2×N2的數(shù)據(jù)作為雷達2的觀測數(shù)據(jù)。令M1=M2=20,N1=N2=25,此時雷達1和雷達2的帶寬均為400 MHz,對應的距離分辨率理論值均為0.375 m;雷達1和雷達2的觀測角度均為1.13°,但由于各自的中心頻率分別為17.2 GHz和18.3 GHz,故兩部雷達的方位分辨率理論值分別為0.44 m和0.41 m。當信噪比(signal-to-noise ratio, SNR)為20 dB時,雷達1和雷達2的RD算法成像結(jié)果分別如圖5(a)和圖5(b)所示。從圖5可以看出,由于發(fā)射帶寬和觀測視角有限,單部雷達成像的二維分辨率低,導致無法完全區(qū)分所有的散射點,影響了成像質(zhì)量。
圖5 兩部雷達的RD成像結(jié)果Fig.5 RD imaging results of two radars
直接采用RD算法對兩部雷達的觀測數(shù)據(jù)進行融合成像,融合后的距離分辨率和方位分辨率理論值分別為0.187 5 m和0.217 m,成像結(jié)果如圖6(a)所示。由于頻帶和視角缺失,RD算法在利用快速傅里葉變換進行二維壓縮成像時引起了強烈的副瓣和大量的能量泄露等問題,嚴重影響了成像質(zhì)量,無法有效地實現(xiàn)融合成像。分別采用正則化算法[22]和OMP算法實現(xiàn)多視角多頻帶融合成像,融合后的距離分辨率和方位分辨率理論值分別為0.10 m和0.12 m,與全視角全頻帶回波的理論分辨率一致,得到的融合成像結(jié)果分別如圖6(b)和圖6(c)所示。從圖6(b)可以看出,采用正則化算法可以有效地實現(xiàn)融合成像,與單部雷達成像結(jié)果相比改善了二維分辨率,但從融合成像結(jié)果中無法完全區(qū)分散射點A和散射點B。這是因為正則化算法中需要設置正則化參數(shù),若參數(shù)設置不合適則會影響算法的重構(gòu)性能。從圖6(c)可以看出,采用OMP算法能夠得到干凈清晰的融合成像結(jié)果,正確估計出各個散射點的大小和位置,且從成像結(jié)果中能夠完全區(qū)分所有散射點。另外,采用正則化算法和OMP算法實現(xiàn)稀疏重構(gòu)時所耗的時間分別為1 213.43 s和4.65 s。由于正則化算法在迭代過程中涉及構(gòu)造大規(guī)模矩陣并求逆問題,所需的內(nèi)存較大且運算時間較長,而OMP算法易于實現(xiàn)且運算過程簡單,更適合求解大規(guī)模矢量問題。
圖6 不同算法的融合成像結(jié)果Fig.6 Fusion imaging results of different algorithms
為進一步體現(xiàn)算法對不同頻率依賴因子的散射點的重構(gòu)效果,重點關(guān)注散射點A和散射點B的參數(shù)估計結(jié)果,如表2所示。從表2中的參數(shù)估計結(jié)果可以看出,利用兩種算法均能成功估計出目標散射點的頻率依賴因子,其中OMP算法對散射點位置坐標的估計精度比正則化算法更高。
表2 目標散射點參數(shù)估計結(jié)果
為進一步驗證算法在復雜目標模型條件下的融合成像性能,利用波音727飛機的回波數(shù)據(jù)實現(xiàn)多視角多頻帶ISAR融合成像。飛機的回波數(shù)據(jù)由美國海軍實驗室在其網(wǎng)站上提供[29],雷達的載頻為9 GHz,帶寬為512 MHz,PRF為20 kHz,每個脈沖內(nèi)截取的距離單元為64個,共發(fā)射256個連續(xù)脈沖,但考慮到目標的機動性,只利用其中100個脈沖的回波數(shù)據(jù)作為全視角全頻帶回波,大小為64×100。全視角全頻帶雷達的二維回波及其RD成像結(jié)果分別如圖7(a)和圖7(b)所示,從其RD成像結(jié)果中可以看出飛機的基本輪廓。
圖7 全視角全頻帶二維回波及其RD成像結(jié)果Fig.7 Full-angle and full-band two-dimensional echo data and RD imaging result
4.2.1 不同觀測情形下算法融合成像性能驗證
為比較在不同觀測數(shù)據(jù)條件下的融合成像效果,在SNR為20 dB條件下,改變雷達1和雷達2的回波數(shù)據(jù)大小,分別采用正則化算法和OMP算法實現(xiàn)融合成像。
情形 1令M1=M2=25,N1=N2=40。雷達1和雷達2的觀測二維回波如圖8(a)所示,采用正則化算法和OMP算法得到的融合成像結(jié)果分別如圖8(b)和圖8(c)所示。從圖8(b)可以看出,采用正則化法進行融合成像時,雖然可以估計出基本的強散射點,但存在弱散射點丟失的情況,導致丟失了目標結(jié)構(gòu)的一些細節(jié)信息,影響了目標的基本形狀。從圖8(c)可以看出,采用OMP算法進行融合成像時,可從成像結(jié)果中分辨出目標形狀及其部分細節(jié)結(jié)構(gòu)信息,成像效果較好。
圖8 M1=M2=25,N1=N2=40時的觀測二維回波數(shù)據(jù)及融合成像結(jié)果Fig.8 Two-dimensional observation echo data and fusion imaging results with M1=M2=25 and N1=N2=40
情形 2令M1=M2=20,N1=N2=30。雷達1和雷達2的觀測二維回波如圖9(a)所示,采用正則化算法和OMP算法得到的融合成像結(jié)果分別如圖9(b)和圖9(c)所示。從圖9(b)可以看出,隨著有效觀測回波數(shù)據(jù)的減少,采用正則化算法進行融合成像時容易丟失散射點,特別是機翼和機尾部分散射點丟失較多,導致從成像結(jié)果中分辨出目標的基本輪廓比較困難。從圖9(c)可以看出,隨著觀測回波數(shù)據(jù)的減少,采用OMP算進行融合成像時,仍能得到較為清晰且完整的目標圖像。
圖9 M1=M2=20,N1=N2=30時的觀測二維回波數(shù)據(jù)及融合成像結(jié)果Fig.9 Two-dimensional observation echo data and fusion imaging results with M1=M2=20 and N1=N2=30
情形 3令M1=M2=16,N1=N2=25。雷達1和雷達2的觀測二維回波如圖10(a)所示,采用正則化算法和OMP算法得到的融合成像結(jié)果分別如圖10(b)和圖10(c)所示。從圖10(b)可以看出,隨著有效觀測回波數(shù)據(jù)的進一步減少,采用正則化算法進行融合成像時,成像結(jié)果中不僅容易丟失散射點,還容易引入虛假散射點,導致難以分辨目標的基本輪廓。從圖10(c)可以看出,在觀測回波數(shù)據(jù)較少的情況下,采用OMP算法實現(xiàn)融合成像時,雖然在融合成像結(jié)果中丟失了少量的散射點,但仍能分辨出目標的基本形狀。
圖10 M1=M2=16,N1=N2=25時的觀測二維回波數(shù)據(jù)及融合成像結(jié)果Fig.10 Two-dimensional observation echo data and fusion imaging results with M1=M2=16 and N1=N2=25
為進一步對比不同有效觀測回波數(shù)據(jù)條件下算法的融合成像性能,采用算法運行時間,圖像對比度(image contrast, IC)和目標背景比(target-to-background ratio, TBR)作為算法性能指標[30]。其中,IC可以評價目標圖像的整體聚焦質(zhì)量,值越大表示圖像越聚焦,質(zhì)量越好;TBR能有效表征目標圖像的能量聚焦程度,可以評價目標圖像的噪聲抑制能力和聚焦能力,值越大越好。不同情形下的算法性能衡量指標如表3所示。
表3 不同情形下的算法性能指標對比
從表3可以看出,在相同的觀測回波數(shù)據(jù)條件下,利用OMP算法得到的成像結(jié)果的IC值和TBR值均比正則化算法大,這是因為正則化算法中涉及到正則化參數(shù)的估計和調(diào)整,若參數(shù)設置不合適則會影響算法性能,而OMP算法無需設置參數(shù),可得到較好的融合成像質(zhì)量。另外,在同一觀測條件下,OMP算法的運行時間明顯比正則化算法小,這是因為正則化算法中涉及到大規(guī)模的矩陣求逆運算,導致運算復雜度較高,運行時間長,不適合大規(guī)模數(shù)據(jù)重構(gòu)。而OMP算法原理簡單且易于實現(xiàn),擁有快速重構(gòu)能力,可用于大規(guī)模數(shù)據(jù)稀疏重構(gòu),比正則化算法更適用于多視角多頻帶ISAR融合成像。
4.2.2 不同SNR條件下算法融合成像性能驗證
為比較在不同SNR條件下的融合成像效果,令M1=M2=20,N1=N2=30,改變回波數(shù)據(jù)的SNR,在SNR為20 dB、10 dB和0 dB條件下分別采用正則化算法和OMP算法實現(xiàn)融合成像,成像結(jié)果分別如圖11和圖12所示。從圖11可以看出,隨著SNR減小,利用正則化算法得到的融合成像結(jié)果中目標的輪廓越來越不清晰,特別是當SNR為0 dB時,融合成像結(jié)果中存在大量的噪點沒有被抑制,嚴重影響了成像質(zhì)量,導致無法辨別目標形狀結(jié)構(gòu)。這是因為正則化算法的重構(gòu)性能與噪聲水平大小密切相關(guān),在重構(gòu)過程中需要估計噪聲水平參數(shù),若參數(shù)估計不準確將會直接影響融合成像質(zhì)量。從圖12中可以看出,隨著SNR減小,利用OMP算法均能得到較好的融合成像結(jié)果,即使在SNR為0 dB的條件下,仍能從成像結(jié)果中分辨出目標的基本形狀。實驗結(jié)果說明與正則化算法相比,OMP算法在低SNR條件下具有更強的魯棒性。
圖11 不同SNR條件下正則化算法的融合成像結(jié)果Fig.11 Fusion imaging results of regularization algorithm under different SNRs
圖12 不同SNR條件下OMP算法的融合成像結(jié)果Fig.12 Fusion imaging results of OMP algorithm under different SNRs
不同SNR條件下算法性能衡量指標如表4所示。從表4中可以看出,兩種算法成像結(jié)果的IC值和TBR值均隨著SNR減小而減小。在相同的SNR條件下,利用OMP算法所得圖像的IC值和TBR值均比正則化算法大,說明在低SNR條件下,利用OMP算法能得到更高質(zhì)量的圖像,體現(xiàn)了本文所提算法的優(yōu)越性。
表4 不同SNR條件下的算法性能指標對比
為提高ISAR成像二維分辨率,提出一種基于GTD模型的多視角多頻帶ISAR融合成像方法。此方法基于GTD模型構(gòu)建了更符合實際的目標回波,利用稀疏表示理論建立了多視角多頻帶ISAR融合成像模型,為在保證重構(gòu)質(zhì)量的同時提高成像效率,采用OMP算法進行稀疏重構(gòu),實現(xiàn)融合成像。仿真實驗利用不同目標模型驗證了所提算法的有效性,并通過在不同的有效觀測數(shù)據(jù)和不同SNR條件下算法的融合成像性能對比,體現(xiàn)了所提算法在提高成像效率和抗噪性能等方面的優(yōu)越性。