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