侯育星 徐 剛(陜西黃河集團(tuán)有限公司 西安 710043)
②(東南大學(xué)信息科學(xué)與工程學(xué)院毫米波國家重點(diǎn)實(shí)驗(yàn)室 南京 210096)
基于傳統(tǒng)合成孔徑雷達(dá)(Synthetic Aperture Radar, SAR) 2維成像,干涉合成孔徑雷達(dá)(Interferometric SAR, InSAR)通過在高度向配置基線,利用多個天線相位中心在高度上的差異能夠獲取目標(biāo)場景多視角觀測,具備對觀測場景3維測繪的能力[1-3]。相應(yīng)地,InSAR在數(shù)字高程測繪(Digital Elevation Model, DEM)生成、場景高度估計和地形變化檢測等方面已經(jīng)得到了廣泛應(yīng)用[4,5]。在InSAR處理中,獲取高精度的干涉相位是其中一項(xiàng)關(guān)鍵技術(shù)。然而,干涉相位不可避免的存在噪聲是InSAR固有問題,這是因?yàn)橥ǖ乐g視角差異引入SAR復(fù)圖像之間的去相干問題,對應(yīng)干涉相位噪聲[6]。在幅度方面,相干斑噪聲是SAR圖像典型性特征,這是由SAR成像的相干積累特性所決定的。直觀地講,相干斑噪聲造成相同散射介質(zhì)的物體在SAR圖像幅度表現(xiàn)出“椒鹽”噪聲特性,從而影響目標(biāo)的分類和識別。因此,InSAR干涉圖像幅度和相位同時存在噪聲問題。為了滿足后續(xù)圖像應(yīng)用的要求,InSAR圖像的相干斑抑制和相位噪聲濾波處理是必要過程。
近十幾年來,許多學(xué)者較為系統(tǒng)的研究了干涉相位噪聲濾波和相干斑抑制算法[7-10]。在現(xiàn)有大多數(shù)研究工作中,干涉相位噪聲建模為加性噪聲模型[11,12],而相干斑則基于乘性噪聲模型假設(shè)。絕大多數(shù)算法都是基于場景目標(biāo)局部平穩(wěn)假設(shè),即同一散射介質(zhì)目標(biāo)相鄰散射系數(shù)具有平穩(wěn)特性,并且目標(biāo)高度上連續(xù)變化(除了邊緣區(qū)域)。近幾年來,越來越多的學(xué)者研究利用結(jié)合幅度和相位提高兩者噪聲抑制的性能。文獻(xiàn)[13,14]提出了一種聯(lián)合幅度和干涉相位的稀疏正則化算法,通過聯(lián)合構(gòu)造幅相差分(Total Variation, TV)字典進(jìn)行聯(lián)合稀疏約束,以提升在城區(qū)場景的降噪能力。在文獻(xiàn)[15]中,作者的前期工作研究了利用自適應(yīng)表征字典進(jìn)行InSAR圖像干涉相位和幅度的聯(lián)合稀疏表征,在稀疏成像中能夠?qū)崿F(xiàn)幅相濾波的特征化增強(qiáng)。文獻(xiàn)[16]利用InSAR干涉處理,構(gòu)造稀疏字典進(jìn)行InSAR圖像稀疏化表征,實(shí)現(xiàn)在降數(shù)據(jù)率采樣下的圖像重構(gòu)。
在稀疏信號處理中,貝葉斯算法作為一種典型性的算法理論,相比其他算法,具有更高的穩(wěn)健性和精度。相應(yīng)的,稀疏貝葉斯算法在雷達(dá)成像中得到了較為廣泛的研究和應(yīng)用。進(jìn)一步,結(jié)構(gòu)化稀疏貝葉斯成像被提出用于提升雷達(dá)成像性能。文獻(xiàn)[17]研究了聯(lián)合稀疏模型在稀疏ISAR(Inverse SAR,ISAR)成像中的應(yīng)用。文獻(xiàn)[18]研究了結(jié)構(gòu)化統(tǒng)計模型在ISAR成像中的應(yīng)用。在文獻(xiàn)[19]中,作者的前期工作研究了結(jié)構(gòu)化稀疏貝葉斯模型在ISAR成像中的應(yīng)用。文獻(xiàn)[17-19]都是利用ISAR圖像中散射目標(biāo)相鄰結(jié)構(gòu)相關(guān)性以提升成像性能,其區(qū)別在于統(tǒng)計建模和算法求解的差異。對于InSAR而言,通道之間具有相干性,利用其相干性,通過結(jié)構(gòu)化貝葉斯方法聯(lián)合通道圖像處理,能夠進(jìn)一步提升干涉相位濾波和相干斑抑制的性能,對應(yīng)本文的主要研究內(nèi)容。
本文提出了一種基于貝葉斯聯(lián)合通道稀疏的InSAR成像算法,通過貝葉斯結(jié)構(gòu)化稀疏化表征,提升干涉相位噪聲濾波和相干斑抑制性能。首先,基于貝葉斯準(zhǔn)則建立InSAR成像模型,通過多層級統(tǒng)計建模,將InSAR復(fù)圖像聯(lián)合幅度和干涉相位,建立為結(jié)構(gòu)化稀疏統(tǒng)計模型。然后,利用最大期望(EM, Expectation Maximization)算法進(jìn)行InSAR圖像重構(gòu)和多層級統(tǒng)計參數(shù)估計。相比較文獻(xiàn)[15]的方法,本文方法由于利用貝葉斯結(jié)構(gòu)化稀疏模型,能夠進(jìn)一步提升InSAR圖像幅度和相位噪聲濾波的能力。最后,通過實(shí)驗(yàn)分析進(jìn)一步驗(yàn)證本文算法的有效性。
圖1所示InSAR成像幾何,在SAR成像基礎(chǔ)上,InSAR通過在高度向增加通道數(shù),即通道(本文主要考慮兩個通道情況,此時),可實(shí)現(xiàn)對目標(biāo)場景不同高度視角的觀測。InSAR與SAR 2維成像幾何相同,平臺飛行方向?yàn)榉轿幌?,對?yīng)圖1中X軸,與X軸垂直的波束指向?yàn)榫嚯x向,對應(yīng)圖1中的Y軸。如圖1所示,在方位時間時刻,通道與目標(biāo)P的距離可以表示為。假設(shè)已經(jīng)利用SAR算法進(jìn)行成像處理,此時數(shù)據(jù)在距離時域和方位多普勒域,第()個通道的信號可以表示為
在第2節(jié)中,式(2)實(shí)質(zhì)是一個求逆以及降噪問題。在文獻(xiàn)[15]中,作者的前期工作已經(jīng)研究了利用自適應(yīng)表征字典實(shí)現(xiàn)InSAR圖像稀疏化表征,其稀疏成像模型為
基于式(5)和式(8),由貝葉斯準(zhǔn)則求得最大后驗(yàn)概率密度函數(shù)
至此,基于貝葉斯建模,本文已經(jīng)建立了聯(lián)合通道稀疏成像模型。假設(shè)給定一組估計參數(shù),的最大后驗(yàn)概率估計為。在實(shí)際中,由于InSAR數(shù)據(jù)量通常較大,因此式(10)中的矩陣求逆運(yùn)算量巨大,為了避免矩陣求逆,可以利用求解線性方程求得
式(11)的問題可以通過共軛梯度算法進(jìn)行求解,參考文獻(xiàn)[20],利用句柄操作直接對2維圖像或者進(jìn)行距離和方位分維度處理,而不拉成向量操作,能有效提高算法的實(shí)際可行性。貝葉斯方法利用EM算法將問題求解分解為圖像重構(gòu)和參數(shù)的估計兩個部分,兩個部分交替迭代,構(gòu)成一個整體。
運(yùn)算量分析:通過分析可知,本文算法的主要運(yùn)算量在于式(10)和式(11)中和的計算,即矩陣的求逆操作。為了避免矩陣求逆,前面已經(jīng)分析,可以利用式(11)通過共軛梯度算法求解。對于式(11)而言,其運(yùn)算復(fù)雜度主要依賴于和操作,分別對應(yīng)于SAR投影算子以及其逆過程操作,在本文中為方位傅立葉變換矩陣及其逆變換矩陣,其運(yùn)算復(fù)雜度為。假設(shè)EM算法迭代次數(shù)為,在每次迭代中,求解式(11)的共軛梯度迭代次數(shù)為,那么本文算法的運(yùn)算量為。其中,運(yùn)算量與通道數(shù)成正比。
表1 本文算法流程框圖Tab. 1 Algorithm flow chart in this paper
下面通過實(shí)驗(yàn)分析以驗(yàn)證本文算法的有效性,并通過與現(xiàn)有方法比較,揭示本文算法的優(yōu)勢。在本節(jié)中,將實(shí)驗(yàn)分為2個部分。第1個實(shí)驗(yàn)部分為仿真數(shù)據(jù)實(shí)驗(yàn),第2個部分為實(shí)測數(shù)據(jù)實(shí)驗(yàn),包含山區(qū)數(shù)據(jù)和城區(qū)數(shù)據(jù)。下面所有實(shí)驗(yàn)均在個人電腦進(jìn)行運(yùn)行,電腦配置為CPU 3.50 GHz雙核處理器,編程軟件為R2015b版本MATLAB。
在仿真數(shù)據(jù)實(shí)驗(yàn)中,理想情況下圖像幅度和干涉相位值如圖2(a)和圖2(b)所示。其中,幅度值對應(yīng)平地和山區(qū)兩個區(qū)域,每個區(qū)域?qū)?yīng)相同的散射系數(shù)。圖2(c)和圖2(d)對應(yīng)存在相干斑和相位噪聲的圖像。為了抑制噪聲,利用本文算法進(jìn)行稀疏化特征成像處理,其成像結(jié)果如圖2(e)和圖2(f)所示??梢姡鹊南喔砂咴肼暫透缮嫦辔辉肼暥嫉玫搅撕芎玫囊种?。為了更好地評估本文算法的性能,此處本文利用小波變換相位濾波[21]算法進(jìn)行比較,其結(jié)果如圖2(g)所示。通過比較圖2(e)和圖2(g)可知,本文算法在地形平坦區(qū)域和高度陡變區(qū)域具有更好的相位濾波效果。圖2(h)和圖2(i)進(jìn)一步給出兩種算法結(jié)果與真實(shí)相位的差異圖。通過比較可知,本文算法的相位估計誤差更小。同時,計算了兩種算法相位估計的均方誤差,本文算法和文獻(xiàn)[21]算法結(jié)果分別為0.0082 rad2和0.0405 rad2,與圖2(h)和圖2(i)一致。在運(yùn)算量方面,EM算法的迭代次數(shù)設(shè)置為20,算法可近似達(dá)到收斂,此時兩個通道圖像處理時間約為118 s。
在實(shí)測數(shù)據(jù)實(shí)驗(yàn)中,我們利用RADARSAT-2重復(fù)航過干涉數(shù)據(jù)進(jìn)行實(shí)驗(yàn),每個通道為單視復(fù)(Single Look Complex, SLC)圖像。兩個通道的錄取回波數(shù)據(jù)間隔時間為24 h,并被應(yīng)用于DEM生成,可以用來高程估計。在第一個實(shí)驗(yàn)中,選取山區(qū)場景的部分?jǐn)?shù)據(jù)。通道1幅度的圖像如圖3(a)所示,兩個通道的原始干涉相位如圖3(b)所示。明顯可見,干涉相位存在噪聲,SAR幅度圖像存在相干斑。然后,利用稀疏特征化增強(qiáng)的方法,在稀疏成像中進(jìn)行噪聲濾波處理。其中,本文利用文獻(xiàn)[15]中作者前期工作的算法作為對比。在稀疏化表征中,都選擇雙樹小波變換作為稀疏表征字典。圖3(d)-圖3(i)分別為文獻(xiàn)[15]和本文方法的成像結(jié)果。相比傳統(tǒng)SAR成像,兩種稀疏方法都能夠較為有效的進(jìn)行干涉相位噪聲濾波和相干斑抑制。同時,我們給出了傳統(tǒng)SAR成像、文獻(xiàn)[15]和本文方法結(jié)果的相位“殘點(diǎn)”,如圖3(c)、圖3(f)和圖3(i)所示??梢?,文獻(xiàn)[15]和本文方法在山區(qū)地形陡變部分,能夠保留一定的細(xì)節(jié)信息。相比于文獻(xiàn)[15],本文方法在地形平坦局域具有更好的濾波效果。需要指出的是,文獻(xiàn)[15]和本文方法由于稀疏字典構(gòu)造等局限性因素,其濾波效果具有一定的低通特性,這是后續(xù)工作需要改進(jìn)的地方。
第2個實(shí)驗(yàn)部分,本文選取城區(qū)場景進(jìn)行實(shí)驗(yàn)。文獻(xiàn)[15]和本文的成像結(jié)果如圖4所示。由圖4可知,本文的方法在在建筑物區(qū)域具有更好的噪聲抑制能力,同時在能夠較好地避免相位的不連續(xù)性,相應(yīng)地,干涉相位圖像的突變區(qū)域?qū)?yīng)幅度圖像的突變部分,表現(xiàn)為觀測場景建筑物的邊緣區(qū)域。為了定性分析,本文通過計算圖像的等效視數(shù)(Equivalent Number of Looks, ENL)評估相干斑抑制的均勻場景區(qū)域強(qiáng)度圖像的均值和方差。本文選取2個不同場景區(qū)域,以目標(biāo)像素為中心,通過加窗(大小)選取相鄰窗口像素計算ENL,其結(jié)果如表2所示。由表2可知,傳統(tǒng)SAR成像的ENL結(jié)果近似為1,與使用的單視復(fù)圖像一致??梢姡疚牡乃惴軌蚋佑行?shí)現(xiàn)在同一散射介質(zhì)區(qū)域的相干斑噪聲抑制。相比于文獻(xiàn)[15]的方法,本文算法可以取得更大的ENL值,對應(yīng)較好的相干斑抑制性能,這得益于算法的結(jié)構(gòu)化稀疏貝葉斯模型。為了評估算法對于獨(dú)立散射點(diǎn)的分辨率保持能力,利用孤立散射點(diǎn)進(jìn)行分辨率測量。傳統(tǒng)SAR成像2維分辨率為2.67 m×4.13 m(理論值為1.92 m×3.27 m)。文獻(xiàn)[15]和本文算法分辨率與傳統(tǒng)成像比較接近,圖5給出了方位脈沖響應(yīng)形式,圖中藍(lán)色、紅色和黑色分別對應(yīng)傳統(tǒng)成像、文獻(xiàn)[15]和本文方法結(jié)果。
表2 ENL估計Tab. 2 ENL estimates
本文針對InSAR稀疏成像,提出了一種通道聯(lián)合的結(jié)構(gòu)化稀疏貝葉斯成像算法?;谪惾~斯準(zhǔn)則建立了聯(lián)合通道結(jié)構(gòu)化稀疏模型,并通過EM算法進(jìn)行了圖像重構(gòu)和參數(shù)估計,相比于傳統(tǒng)稀疏算法,本文算法利用了通道之間聯(lián)合稀疏的相關(guān)特性,能夠更加有效地提高InSAR成像算法中的噪聲濾波特性,降低幅度噪聲和相位噪聲對成像結(jié)果的影響,實(shí)驗(yàn)結(jié)果證明了本文算法的有效性。