質(zhì)量引導(dǎo)下的最小二乘相位解纏算法

肖暉, 李惠堂, 顧約翰, 盛慶紅

(1.南京航空航天大學(xué)航天學(xué)院,南京 210016; 2.南京曉莊學(xué)院環(huán)境科學(xué)學(xué)院,南京 211171)

0 引言

干涉合成孔徑雷達(dá)(interferometric synthetic aperture Radar,InSAR[1])是利用2幅或多幅相干合成孔徑雷達(dá)(synthetic aperture Radar,SAR)圖像復(fù)數(shù)數(shù)據(jù)得到干涉相位,并通過成像幾何模型建立相位與地形信息之間的對(duì)應(yīng)關(guān)系,從而獲取地表三維信息[2]和形變信息的對(duì)地觀測(cè)技術(shù)[3-4]。由于InSAR使用的SAR圖像為復(fù)數(shù)數(shù)據(jù),所以在干涉處理時(shí)獲得的干涉相位因三角函數(shù)的周期性被纏繞在(-π,π]之間,導(dǎo)致纏繞相位與絕對(duì)相位之間存在2kπ的差異[5],所以需要通過相位解纏來恢復(fù)相位主值中被模糊掉的整周相位,從而求解地形高程對(duì)應(yīng)的絕對(duì)相位。

目前常見的相位解纏方法主要包括3類[6]: 基于路徑跟蹤、基于最小范數(shù)和基于網(wǎng)絡(luò)規(guī)劃的相位解纏算法?；诼窂礁櫟南辔唤饫p算法通過先驗(yàn)信息制定合適的積分路徑,對(duì)相鄰像元進(jìn)行積分,實(shí)現(xiàn)相位解纏,代表算法有枝切法和質(zhì)量引導(dǎo)法?；诼窂礁櫟南辔唤饫p算法具有解纏速度快、效率高等優(yōu)點(diǎn),但由于積分路徑難以確定,算法精度無法保證[7-9]?；谧钚》稊?shù)的相位解纏方法是將解纏問題轉(zhuǎn)換為纏繞相位與解纏相位差異最小化的數(shù)學(xué)問題,利用相位的離散偏導(dǎo)數(shù)作為最小范數(shù)參量實(shí)現(xiàn)相位解纏,代表算法有無權(quán)最小二乘法[10-11]和加權(quán)最小二乘法[12-13]。基于網(wǎng)絡(luò)規(guī)劃的相位解纏算法通過網(wǎng)絡(luò)流理念將解纏相位和纏繞相位之間的離散偏導(dǎo)數(shù)之差最小化,從而獲取全局最優(yōu)解,代表算法有基于最小費(fèi)用流相位解纏算法和基于不規(guī)則網(wǎng)絡(luò)的相位解纏算法[14]。此類方法大大提高了相位解纏的運(yùn)行效率和解纏精度,但如何確定適當(dāng)?shù)臋?quán)重矩陣,從而獲取更好的積分路徑仍然是研究的重點(diǎn)和難點(diǎn)[15-16]。

針對(duì)復(fù)雜山區(qū)中因斜視成像和地形帶來的相位失真及相位噪聲,本文提出質(zhì)量引導(dǎo)與最小二乘相結(jié)合的相位解纏算法,首先根據(jù)干涉相位質(zhì)量高低在高質(zhì)量區(qū)域利用質(zhì)量引導(dǎo)法獲取可靠的解纏結(jié)果,然后利用最小二乘法的平滑作用,在低質(zhì)量區(qū)域進(jìn)行最小二乘法相位解纏,提高低質(zhì)量區(qū)域解纏的連續(xù)性,從而獲得全局最優(yōu)的解纏結(jié)果。仿真數(shù)據(jù)與真實(shí)數(shù)據(jù)結(jié)果表明該算法很好地克服了質(zhì)量引導(dǎo)法無法阻止相位解纏誤差擴(kuò)散的問題,還可以在最小二乘法的平滑作用下,提高低質(zhì)量區(qū)域解纏的連續(xù)性,避免了低質(zhì)量像元因相位誤差在路徑積分中堆積,導(dǎo)致與高質(zhì)量像元解纏結(jié)果間出現(xiàn)相位跳變的現(xiàn)象。

1 數(shù)據(jù)源

1.1 模擬數(shù)據(jù)

為了盡可能模擬復(fù)雜山區(qū)地形,本文采用Peaks函數(shù)生成255×255的模擬數(shù)字高程模型(digital elevation model,DEM)矩陣,如圖1(a)所示,高程分布在0～63 m,其中包含3個(gè)山峰和2個(gè)山谷區(qū)域。圖1(b)是根據(jù)圖1(a)得到的真實(shí)相位圖,其相位值處在[-34,42] rad,均值為2.622 rad。由于山區(qū)中植物茂盛,所以對(duì)應(yīng)的時(shí)間及空間相干性較低。為了更逼近真實(shí)數(shù)據(jù),對(duì)模擬DEM生成的干涉相位添加最大值為0.65 rad的隨機(jī)噪聲。圖1(c)為在真實(shí)相位加入隨機(jī)噪聲后纏繞在(-π,π]的干涉相位圖,從圖1中可以看出在黑色框選區(qū)域因噪聲影響纏繞相位出現(xiàn)大量不一致情況。

(a) DEM 3D視圖 (b) DEM平面視圖 (c) 纏繞相位

1.2 真實(shí)數(shù)據(jù)

為了驗(yàn)證本文提出的相位解纏算法在真實(shí)數(shù)據(jù)中的有效性,采用位于中國陜西西安的秦嶺山脈作為研究區(qū)域,該區(qū)域中心經(jīng)緯度為E108°50′,N33°35′,相關(guān)數(shù)據(jù)如圖2所示。圖2(b)給出研究區(qū)域的SRTM-1拼接后的影像,可以看出該區(qū)域群山連綿、地形復(fù)雜,經(jīng)過對(duì)該區(qū)域SRTM DEM進(jìn)行統(tǒng)計(jì),可知高程分布在200～3 000 m。

(a) 研究區(qū)SRTM-1影像圖 (b) 研究區(qū)高程

采用的真實(shí)干涉SAR影像對(duì)為Sentinel-1A數(shù)據(jù), 產(chǎn)品類型為SLC,成像模式為IW,入射角為39.08°,軌道方向?yàn)樯?。獲取2019年11月26日和12月 8 日的InSAR干涉影像對(duì),影像大小為1 000像素×1 000像素,時(shí)間基線為12 d,空間基線為24.89 m。

在獲取SAR干涉影像對(duì)后,需要對(duì)獲取的圖像進(jìn)行InSAR流程化處理,處理過程如圖3所示,可以看出該研究區(qū)域基于SAR干涉影像對(duì)生成的干涉相位圖紋理模糊,經(jīng)過干涉圖濾波后出現(xiàn)部分干涉條紋,經(jīng)過去平地效應(yīng)后干涉條紋較明顯,結(jié)合該區(qū)域的SRTM DEM影像可以看出在山谷區(qū)域有較連續(xù)的干涉條紋,而在山峰處干涉相位雜亂無章。

圖3 干涉影像對(duì)流程化處理

在經(jīng)過流程化處理后的干涉相位圖中選取其中的1 000像素×1 000像素進(jìn)行相位解纏研究,圖4給出實(shí)驗(yàn)干涉相位相關(guān)數(shù)據(jù),其中圖4(a)為干涉相位數(shù)據(jù),圖4(b)為根據(jù)SLC影像對(duì)求解的相關(guān)系數(shù)圖,圖4(c)為干涉相位中的殘差點(diǎn)分布圖。通過圖4可以看出,該實(shí)驗(yàn)區(qū)域除相關(guān)性較高區(qū)域有明顯干涉條紋,其他區(qū)域相位雜亂無章。

(a) 纏繞相位 (b) 相關(guān)系數(shù) (c)殘差點(diǎn)分布

由于SAR影像與SRTM DEM數(shù)據(jù)空間分辨率不同,所以利用選取的控制點(diǎn),對(duì)SRTM DEM重采樣,使2張影像在控制點(diǎn)相同的條件下具有相同的尺寸,從而保證相同的分辨率。圖5給出研究區(qū)域?qū)?yīng)的重采樣后DEM數(shù)據(jù)及求解的相位導(dǎo)數(shù)方差圖。通過圖5(a)—(b)中DEM數(shù)據(jù)可以看出,實(shí)驗(yàn)區(qū)域高程分布在400～2 000 m之間,地形多變,主要包含2條山谷和3個(gè)山峰。圖5(c)是根據(jù)纏繞相位獲得的相位導(dǎo)數(shù)方差,從圖中可以看出研究區(qū)域內(nèi)山谷地區(qū)因地形相對(duì)平坦,植被較少導(dǎo)致相位導(dǎo)數(shù)方差較小,相位質(zhì)量高。

(a) 秦嶺地區(qū)DEM 3D視圖 (b) DEM平面視圖 (c) 相位導(dǎo)數(shù)方差

2 研究方法

質(zhì)量引導(dǎo)法利用相位質(zhì)量圖作為輔助信息來布置積分路徑[17]。相位質(zhì)量圖是從干涉影像對(duì)或干涉相位中提取、用來衡量干涉相位質(zhì)量好壞的標(biāo)準(zhǔn)。常用的相位質(zhì)量圖主要有相關(guān)系數(shù)圖、偽相關(guān)系數(shù)圖、相位導(dǎo)數(shù)方差圖和最大相位梯度圖[18]。大量的實(shí)踐證明,在無法獲得相關(guān)系數(shù)圖時(shí),相位導(dǎo)數(shù)方差圖是最可靠的質(zhì)量圖。本文重點(diǎn)介紹相位導(dǎo)數(shù)方差圖的計(jì)算方法。

相位導(dǎo)數(shù)方差Zm,n是利用干涉相位局部統(tǒng)計(jì)特性,判斷相位數(shù)據(jù)的好壞程度,其計(jì)算公式為:

(1)

,

(2)

,

(3)

相位導(dǎo)數(shù)方差圖是利用k×k范圍內(nèi)相位導(dǎo)數(shù)求取的質(zhì)量圖,該質(zhì)量圖表征的是受到干擾后干涉相位的連續(xù)性,即相位導(dǎo)數(shù)越大對(duì)應(yīng)的干涉相位連續(xù)性越差,可靠性越低。

為了減小殘差點(diǎn)引起的全局誤差,獲取更精確的解纏相位值,加權(quán)最小二乘法通過引入質(zhì)量圖、掩模等先驗(yàn)信息,為干涉相位設(shè)置不同權(quán)重,尤其是對(duì)低相干區(qū)域賦予較低權(quán)重。通常情況,權(quán)重wi,j的取值范圍在[0,1],可以由相位質(zhì)量圖或其他先驗(yàn)信息獲得。引入權(quán)重后,最小化梯度差異問題演變成求解下式的最小值,公式為:

,

(4)

,

(5)

式中:W為最小值;ψ為絕對(duì)相位;wi,j為第i行、第j列的權(quán)值。

式(4)中最小二乘解ψi,j滿足:

Ui,j(ψi+1,j-ψi,j)-Ui-1,j(ψi,j-ψi-1,j)+Vi,j(ψi,j+1-ψi,j)-Vi,j-1(ψi,j-ψi,j-1)=ci,j

,

(6)

(7)

式中ci,j為加權(quán)拉普拉斯算子。

由于引入權(quán)值后,無法直接采用傅里葉變換求解解纏相位,所以在實(shí)際運(yùn)算中,常通過某種數(shù)學(xué)變換將纏繞相位和絕對(duì)相位之間的關(guān)系線性化,再采用傅里葉變換求解ψ的初始值,并進(jìn)行迭代循環(huán),從而求得解纏相位。

本文實(shí)現(xiàn)的質(zhì)量引導(dǎo)與最小二乘相結(jié)合的相位解纏算法,其假設(shè)干涉相位是二維的連續(xù)函數(shù),則只需要從r0開始沿著某一條路徑C進(jìn)行積分,可得到任意一點(diǎn)r的相位,即

,

(8)

。

(9)

根據(jù)最小二乘相位解纏的基本思想,最終解纏相位ψi,j將滿足如下要求,即

,

(10)

,

(11)

該算法主要步驟如下: ①根據(jù)干涉相位圖計(jì)算相位導(dǎo)數(shù)方差圖作為質(zhì)量圖; ②搜索相位導(dǎo)數(shù)方差圖中質(zhì)量最好的點(diǎn),并由該點(diǎn)開始,根據(jù)質(zhì)量引導(dǎo)法解纏步驟進(jìn)行相位解纏工作; ③將相位導(dǎo)數(shù)方差數(shù)值歸一化,代入式(5),計(jì)算加權(quán)最小二乘解纏算法中的權(quán)重; ④將質(zhì)量引導(dǎo)法解纏結(jié)果作為新的干涉質(zhì)量圖,結(jié)合步驟③求解的權(quán)重,設(shè)置迭代次數(shù),求解最終解纏結(jié)果。該算法的優(yōu)勢(shì)在于利用質(zhì)量引導(dǎo)法在相位質(zhì)量高的區(qū)域具有較高的解纏精確性,從而保證在進(jìn)行加權(quán)最小二乘運(yùn)算時(shí),該區(qū)域的解纏相位依然具有較高的精確性,同時(shí)在全局算法的作用下,保證了高質(zhì)量區(qū)域與相鄰的低質(zhì)量區(qū)域間的解纏相位連續(xù)性,從而獲取具有更高精確性和一致性的解纏結(jié)果。

3 結(jié)果與分析

3.1 模擬數(shù)據(jù)實(shí)驗(yàn)結(jié)果

分別采用枝切法、基于相位導(dǎo)數(shù)方差圖的質(zhì)量引導(dǎo)法、無權(quán)最小二乘法、加權(quán)最小二乘法和本文提出的質(zhì)量引導(dǎo)與最小二乘法相結(jié)合的相位解纏算法對(duì)模擬數(shù)據(jù)進(jìn)行解纏處理。

圖6給出實(shí)驗(yàn)中使用的干涉相位相關(guān)圖像,其中圖6(a)為根據(jù)圖1(c)計(jì)算得到的殘差點(diǎn)分布圖,其中±1為殘差點(diǎn),0為非殘差點(diǎn); 圖6(b)為枝切法中平衡殘差點(diǎn)極性的枝切線布置圖; 圖6(c)為計(jì)算所得的圖1(c)的相位導(dǎo)數(shù)方差圖。根據(jù)圖1(c)中干涉相位在黑色方框中的密集分布,結(jié)合圖6,可以看出相位噪聲使對(duì)應(yīng)位置的殘差點(diǎn)數(shù)量增多,導(dǎo)致該區(qū)域內(nèi)枝切法為平衡殘差點(diǎn)極性布置了大量的枝切線,同時(shí)相位噪聲使該位置相位連續(xù)性降低,質(zhì)量引導(dǎo)法使用的相位導(dǎo)數(shù)方差的數(shù)值增加,即該位置相位質(zhì)量較差,對(duì)應(yīng)的解纏次序也靠后。

(a) 殘差點(diǎn)分布 (b) 枝切線分布 (c) 相位導(dǎo)數(shù)方差

(a) 枝切法解纏結(jié)果圖 (b) 枝切法反纏繞圖 (c) 枝切法解纏誤差分布圖

(g) 無權(quán)最小二乘法解纏結(jié)果圖 (h) 無權(quán)最小二乘法反纏繞圖 (i) 無權(quán)最小二乘法解纏誤差分布圖

圖7給出各算法針對(duì)圖1(c)干涉相位的解纏結(jié)果相關(guān)圖,包括解纏結(jié)果、反纏繞和誤差分布,表1給出根據(jù)解纏誤差求解的均方根誤差(root mean square error,RMSE)。圖7(a) —(c)為枝切法解纏結(jié)果,其中出現(xiàn)因枝切線布置過于密集導(dǎo)致的無法解纏區(qū)域,在相位積分時(shí)由于缺乏該區(qū)域的積分相位值,導(dǎo)致在積分路徑上的解纏相位整體小于鄰近值,所相位誤差集中在某幾個(gè)范圍內(nèi),其RMSE為7.37 rad。圖7(d) —(f)為基于圖6(c)的質(zhì)量引導(dǎo)法解纏結(jié)果,其中相位誤差在低質(zhì)量區(qū)域中不斷累積,導(dǎo)致解纏結(jié)果與周圍高質(zhì)量像元的解纏結(jié)果有明顯相位跳變,相位誤差分布較集中,其RMSE為2.37 rad。雖然枝切法和質(zhì)量引導(dǎo)法因噪聲區(qū)域誤差較大,但在高質(zhì)量數(shù)據(jù)區(qū)域有良好的解纏效果,同時(shí)反纏繞誤差與原始干涉相位有較好一致性。基于最小范數(shù)算法的無權(quán)最小二乘法(圖7(g) —(i))和加權(quán)最小二乘法(圖7(j) —(l)),解纏結(jié)果表現(xiàn)出良好的相位連續(xù)性,但作為基于全局范圍的解纏算法,使得相位誤差全局平均化,導(dǎo)致反纏繞圖與原始干涉相位吻合度不高,其RMSE分別為7.11 rad和4.07 rad。本文算法解纏結(jié)果(圖7(m) —(o)),在結(jié)果上表現(xiàn)較理想,不僅保持了高質(zhì)量相位的解纏精度,也平滑了高低質(zhì)量像素間的相位跳變,通過誤差分布可以看出,該方法將質(zhì)量引導(dǎo)法中的集中相位誤差平均化,很好地改善了整體相位的一致性,其RMSE為2.00 rad,結(jié)果表明該算法較其他算法在一致性和精確性上都有所提升。

表1 不同算法的解纏結(jié)果RMSE

3.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)結(jié)果

分別采用與模擬數(shù)據(jù)處理過程相同的相位解纏算法對(duì)真實(shí)數(shù)據(jù)進(jìn)行處理,并分別給出解纏結(jié)果、反纏繞和解纏誤差分布,如圖8所示。通過對(duì)解纏誤差進(jìn)行統(tǒng)計(jì),各算法RMSE如表2所示。由于山區(qū)內(nèi)植物茂盛,獲取的2幅干涉影像時(shí)間相關(guān)性很低; 地形起伏大降低了干涉影像的空間相關(guān)性,因此干涉相位中存在大量殘差點(diǎn),枝切線布置過于密集導(dǎo)致枝切法解纏失敗。圖8(a) —(c)為質(zhì)量引導(dǎo)法解纏結(jié)果,從圖8(a)可以看出解纏結(jié)果整體上符合地形趨勢(shì),其中1號(hào)黑色方框中的山峰區(qū)域和兩側(cè)的山谷區(qū)域與地形中對(duì)應(yīng)位置相符,主要是由于山谷區(qū)域相關(guān)系數(shù)稍高,所以解纏結(jié)果在該區(qū)域具有良好的一致性。但圖8(a)中的2號(hào)黑色方框中,結(jié)合圖4(b)可以看出該區(qū)域沒有干涉條紋且相關(guān)系數(shù)低,作為最后解纏的區(qū)域,路徑積分中的誤差積累使本應(yīng)為山峰的區(qū)域被錯(cuò)誤解纏為山谷。質(zhì)量引導(dǎo)法的反纏繞結(jié)果圖8(b)與原始相位圖4(a)保持較高一致性,但因相位噪聲較大,所以解纏結(jié)果誤差圖如圖8(c)分布范圍廣,其RMSE為20.08 rad。圖8(d) —(f)為無權(quán)最小二乘法解纏結(jié)果,從圖8(d)可以看出僅在1號(hào)黑色方框的山谷區(qū)域與地形相符,在其他區(qū)域?yàn)榱双@取平滑的解纏結(jié)果使局部噪聲全局傳播,最終導(dǎo)致解纏相位偏離真實(shí)相位,無法辨識(shí)地形形狀。該算法反纏繞后如圖8(e)所示,與原始相位不符,其解纏結(jié)果誤差分布(圖8(f))較為集中,其RMSE為24.38 rad。圖8(g) —(i)為加權(quán)最小二乘法解纏結(jié)果,由于該方法利用先驗(yàn)信息賦予權(quán)重,但由于先驗(yàn)信息也來自干涉相位數(shù)據(jù),無法避免成像中的相位噪聲,所以無法準(zhǔn)確確定權(quán)重,僅在圖8(g)中1號(hào)黑色方框中山谷處的解纏結(jié)果比無權(quán)最小二乘法理想,其他區(qū)域也是無法辨識(shí)地形,其RMSE為24.04 rad。圖8(j) —(l)為本文算法解纏結(jié)果,通過圖8(j)中1號(hào)和2號(hào)黑色方框中所示,2條地形坡度較小的山谷區(qū)域因具有很好的相干性,解纏結(jié)果具有很好的精確性和一致性,與地形也相吻合。圖8(k)可以看出反纏繞結(jié)果與原始相位圖4(a)高度一致。通過解纏結(jié)果誤差分布圖(圖8(l))可以看出,誤差較大的像素個(gè)數(shù)減少,誤差分布也相對(duì)集中,通過統(tǒng)計(jì)得到RMSE為16.47 rad,較其他傳統(tǒng)方法都低,表明解纏精度有所提升。

表2 真實(shí)數(shù)據(jù)各算法解纏結(jié)果RMSE統(tǒng)計(jì)

(a) 質(zhì)量引導(dǎo)法解纏結(jié)果圖 (b) 質(zhì)量引導(dǎo)法反纏繞圖 (c) 質(zhì)量引導(dǎo)法解纏誤差分布圖

4 討論與結(jié)論

為了提高InSAR在山區(qū)的DEM反演精度,本文針對(duì)山區(qū)嚴(yán)重的相位失相關(guān)和相位噪聲問題,提出質(zhì)量引導(dǎo)與最小二乘法相結(jié)合的相位解纏算法。首先利用質(zhì)量引導(dǎo)法在高質(zhì)量區(qū)域進(jìn)行解纏,從而獲取可靠的解纏結(jié)果。為了避免質(zhì)量引導(dǎo)法在解纏過程中依賴于干涉相位的質(zhì)量高低而出現(xiàn)高質(zhì)量像元被優(yōu)先解纏但低質(zhì)量像元因相位誤差在路徑積分中堆積,導(dǎo)致與高質(zhì)量像元解纏結(jié)果間出現(xiàn)相位跳變的現(xiàn)象,根據(jù)最小二乘法原理,基于解纏相位和纏繞相位之差的平方和最小準(zhǔn)則,使纏繞相位離散偏導(dǎo)數(shù)和真實(shí)相位離散偏導(dǎo)數(shù)之差的平方和最小化,以質(zhì)量引導(dǎo)法解纏結(jié)果作為纏繞相位進(jìn)行最小二乘法相位解纏。實(shí)驗(yàn)結(jié)果表明,本文方法不僅可以保證質(zhì)量引導(dǎo)法獲得的高質(zhì)量區(qū)域解纏相位的精確性,還可以在最小二乘法的平滑作用下,有效地避免低質(zhì)量區(qū)域誤差的傳播,提高山區(qū)低質(zhì)量區(qū)域解纏結(jié)果的一致性和精確性。