質量引導下的最小二乘相位解纏算法

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

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

0 引言

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

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

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

1 數據源

1.1 模擬數據

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

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

1.2 真實數據

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

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

采用的真實干涉SAR影像對為Sentinel-1A數據, 產品類型為SLC,成像模式為IW,入射角為39.08°,軌道方向為升軌。獲取2019年11月26日和12月 8 日的InSAR干涉影像對,影像大小為1 000像素×1 000像素,時間基線為12 d,空間基線為24.89 m。

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

圖3 干涉影像對流程化處理

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

(a) 纏繞相位 (b) 相關系數 (c)殘差點分布

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

(a) 秦嶺地區DEM 3D視圖 (b) DEM平面視圖 (c) 相位導數方差

2 研究方法

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

相位導數方差Zm,n是利用干涉相位局部統計特性,判斷相位數據的好壞程度,其計算公式為:

(1)

,

(2)

,

(3)

相位導數方差圖是利用k×k范圍內相位導數求取的質量圖,該質量圖表征的是受到干擾后干涉相位的連續性,即相位導數越大對應的干涉相位連續性越差,可靠性越低。

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

,

(4)

,

(5)

式中:W為最小值;ψ為絕對相位;wi,j為第i行、第j列的權值。

式(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為加權拉普拉斯算子。

由于引入權值后,無法直接采用傅里葉變換求解解纏相位,所以在實際運算中,常通過某種數學變換將纏繞相位和絕對相位之間的關系線性化,再采用傅里葉變換求解ψ的初始值,并進行迭代循環,從而求得解纏相位。

本文實現的質量引導與最小二乘相結合的相位解纏算法,其假設干涉相位是二維的連續函數,則只需要從r0開始沿著某一條路徑C進行積分,可得到任意一點r的相位,即

,

(8)

。

(9)

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

,

(10)

,

(11)

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

3 結果與分析

3.1 模擬數據實驗結果

分別采用枝切法、基于相位導數方差圖的質量引導法、無權最小二乘法、加權最小二乘法和本文提出的質量引導與最小二乘法相結合的相位解纏算法對模擬數據進行解纏處理。

圖6給出實驗中使用的干涉相位相關圖像,其中圖6(a)為根據圖1(c)計算得到的殘差點分布圖,其中±1為殘差點,0為非殘差點; 圖6(b)為枝切法中平衡殘差點極性的枝切線布置圖; 圖6(c)為計算所得的圖1(c)的相位導數方差圖。根據圖1(c)中干涉相位在黑色方框中的密集分布,結合圖6,可以看出相位噪聲使對應位置的殘差點數量增多,導致該區域內枝切法為平衡殘差點極性布置了大量的枝切線,同時相位噪聲使該位置相位連續性降低,質量引導法使用的相位導數方差的數值增加,即該位置相位質量較差,對應的解纏次序也靠后。

(a) 殘差點分布 (b) 枝切線分布 (c) 相位導數方差

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

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

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

表1 不同算法的解纏結果RMSE

3.2 真實數據實驗結果

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

表2 真實數據各算法解纏結果RMSE統計

(a) 質量引導法解纏結果圖 (b) 質量引導法反纏繞圖 (c) 質量引導法解纏誤差分布圖

4 討論與結論

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