融合梯度向量流條紋探測與圖切法相位解纏

劉成洲, 李 斌, 許藝騰

(中交天津港灣工程研究院有限公司巖土工程研究所， 天津 300202)

相位解纏是一個將纏繞相位恢復到絕對相位的過程，是合成孔徑雷達干涉測量(interferometric synthetic aperture radar, InSAR)處理過程中至關重要的一環。解纏結果的正確與否直接影響InSAR在礦區變形監測的可靠性和精度[1]。大多數的相位解纏算法是建立在理想的纏繞條件下，即鄰近像元的相位差分絕對值小于π，而實際的相位表面往往是不連續，或是纏繞相位含有噪聲，這類算法很難正確解纏出真實相位。目前，相位解纏算法大致可以分為四類：路徑跟蹤[2]、最小Lp范數[3]、貝葉斯正則化[4]和參數模型[5]。路徑跟蹤法，是通過設置合理的積分路徑，以積分相鄰像元上的纏繞相位梯度的方式進行相位解纏。比較常見的路徑跟蹤法有Goldstein枝切法[6]和質量圖引導的枝切法[7]。最小Lp范數法，是一種滿足全局最優條件的相位解纏算法，該算法使解纏圖像的誤差達到相對最小，其解纏的結果具有連續性，不會出現孤島現象。比較常見最小Lp范數法有最小二乘算法[8]、圖切法[9]。貝葉斯正則化是依靠數據觀察機制，以及相位先驗知識建立的模型，常見的貝葉斯正則化方法有隨機非線性相位重構法、基于局部平滑的自適應正則化相位估計。參數模型算法將解纏相位約束到一個參數表面。在文獻[10]中，對圖像進行分塊，每一個塊應用不同的參數模型進行解纏，最后將擬合后的圖像融合成解纏相位。此外，還有遺傳算法[11](利用遺傳算法優化相位解纏參數，使得纏繞和解纏繞相位梯度差別位置總數最小)、蟻群算法[12](通過蟻群算法求取“枝切線”上殘差點的最短路徑，使得枝切線不易形成封閉路徑，避免了部分區域無法解纏)、卡爾曼濾波[13](其將相位解纏問題轉化為狀態估計問題，建立相位的動態方程和觀測方程來求解真實相位，在解纏的同時也消除了相位噪聲)、區域生長法[14](利用區域增長法對纏繞相位進行編碼標記，通過疊加的方式求取解纏相位)等應用到相位解纏。解纏方法很多，很難說哪種方法具有絕對的優勢，因此，需要根據不同的數據特點選取不同的解纏方法。

針對礦區在沉降變形是出現的漏斗狀變形，同時礦區地表由于覆蓋農田、水域和樹木等散射性強的地物，使其在干涉相位上表現為同心環狀，不連續，高噪聲的特點，利用本文提出的融合GVF-Snake條紋探測與馬爾科夫隨機場(Markov random field,MRF)圖切法能夠實現礦區變形精細化解纏。

1 算法的基本原理

1.1 基于GVF-Snake模型邊界探測

Xu等[15]針對Snake模型無法收斂至凹型邊界以及對初始化輪廓較為敏感等問題，提出了GVF-Snake模型。該模型重新定義了外部能量場，即外部能量場用全局梯度向量場來表示。具體實現流程：首先利用Sobel算子提取邊界圖f(x,y),然后通過極小化能量方程，得到整個圖像域的梯度矢量流場V(x,y)=[u(x,y),v(x,y)]。則能量E方程可表示為

(1)

GVF場即式(1)中V(x,y)達到最小時E的值，并用它來代替Snake中Eext[Φ(s)]。為求能量函數的極小值，通常采用變分法，式(1)的解滿足

(2)

(3)

引入時間變量t,把u、v看成時間t的函數，采用有限差分離散化迭代求解，解得u、v為

(4)

(5)

1.2 MRF圖切法基本原理

Bioucas-Dias等[9]將纏繞相位當做目標函數，并用一階馬爾科夫隨機場來表示。通過求取目標函數的最大后驗估計來估計邊界，并用圖論中的最大流最小割對圖像進行分割達到解纏的目的。

定義能量函數為

(6)

式(6)中：i、j表示像元的行列數，(i,j)∈G0≡{(k,l):k=1,2，…,M,l=1,2，…,N}(G0是圖像二維格網行列索引號)；V(·)是由多個相互毗鄰集合定義出的集合函數；hij表示沿著水平不連續參數一階鄰近，vij是沿著豎直不連續參數一階鄰近，hij、vij∈{0,1}，當hij、vij=0時，信號是不連續的，像元不連續。

式(6)各分量計算式為

(7)

(8)

(9)

(10)

式中：k≡{kij∈Z:(i,j)∈G0}是標記場；ψ≡{ψij∈[-π,π):(i,j)∈G0}是觀測得到的纏繞相位； 上標h、v表示水平向和豎直向。

1.3 融合條紋探測與MRF圖切法相位解纏

算法的實現包括以下步驟：

(1)利用基于GVF-Snake模型邊界探測的相位解纏算法進行分步解纏，假設解纏到第n步。

(2)判斷n步模型檢測邊界是否為真。

(3)對“不為真”的情況，邊界探測的相位解纏已經實現了n-1步解纏，定義n-1步探測條紋外側為解纏外部塊，內側為纏繞內部塊。由于第n步條紋探測失敗，需要對纏繞內部塊進行MRF圖切法相位解纏，將其解纏結果按照解纏準則增加或減少2π(n-1)后，再與纏繞外部塊進行融合(具體的解纏準則為由外到內，相位從π向-π變化，為加，否則為減)，得到含有邊界孤立點的粗解纏圖像。

(4)高通濾波插值法消除邊界上的孤立點，獲得最終的解纏圖像。

其流程如圖1所示。

圖1 融合條紋探測與MRF圖切法相位解纏Fig.1 Fusion fringe detection and phase unwrapping by MRF graph cuts

流程中通過距離判別條紋是否為真，首先目視法手動選取初始邊界點(與真實邊界相差很小，可看作真實邊界)，尋找每一個初始邊界點到與之對應的GVF-Snake模型探測邊界的最小距離的點，即

(11)

式(11)中：Xm、Ym為GVF-Snake模型探測邊界的任意點；xm、ym為第m個初始化邊界點坐標；m是小于等于初始化邊界點個數的任意正整數，m=1,2,3,…。

通過多次試驗獲得的最小距離閾值，判定依據是通過真實條紋與探測邊界的對比，探測條紋真假很容易分辨。詳情如式(12)和圖2所示。

(12)

式(12)中：R為判別率，規定R=0時，探測邊界線“為真”，否則，“為假”；τ為最小距離閾值，如圖2所示。

圖2 判別探測邊界是否為真Fig.2 Discriminate whether the detection boundary is true or not

圖2是判斷GVF-Snake探測邊界是否為真的示意圖。圖2中紅線表示其實際探測曲線，黑色線表示通過多次試驗獲得最小距離閾值曲線，藍色點表示手動選取得初始化邊界點，黑色箭頭表示初始化邊界點到探測邊界得最小距離，當黑色箭頭線在兩條閾值邊界之間，表示條紋探測線與真實邊界相符，否則條紋探測失敗。

1.4 高通插值濾波

圖像粗解纏后會在其邊界形成很多孤立點，這些孤立點相對高頻，因此可以使用高頻濾波將這些邊界點標記出來，進行去除，然后再用雙三次內插法將去除的點進行插值。實現步驟如圖3所示。

圖3 高通濾波插值法流程圖Fig.3 Flow chart of interpolation method for high-pass filtering

從圖3可知，該濾波算法的核心是將圖像轉換到頻率域，通過設置高通濾波器來檢測高頻點，然后通過傅里葉逆變換轉換到圖像域。傅里葉變換與反變換，可以表示為

(13)

(14)

式中：M、N分別是圖像的行數和列數；u=x=1,2,…,M,v=y=1,2,…,N，傅里葉變換后的相位譜為

(15)

將高頻變換濾波器設置為

(16)

式(16)中：D0為濾波器的阻帶半徑；D(u,v)是點到濾波器中央的距離。

最后，將檢測出的高頻點進行剔除，并用雙三次內插法對空值進行插值達到濾波目的。該算法的缺點是對不是邊界的高頻值也具有濾波作用，濾波后整幅圖會更加平滑,濾波效果如圖4所示。

圖4 高通插值濾波效果圖Fig.4 High-pass interpolation filtering effect map

圖4(b)反映了高通濾波剔除高頻噪聲后的圖像，白色表示空值，為高頻噪聲點，通過對圖4(b)進行雙三次內插法插值填補空區域，得到連續的解纏圖像。

2 實驗結果與分析

2.1 研究區域概況

研究區域為菏澤市巨野礦區郭屯煤礦某工作面，工作面走向全長(平距)1 795 m，面長(平距)125 m，開采標高-727～-815 m，表土層厚度約572.1 m，地面標高+42.9～+45.4 m。研究區域及布線情況如圖5所示。

圖5 研究區域與布線情況Fig.5 Research area and location

2.2 數據選取

選擇sentinel-1A Interferometric Wide (IW) 模式，這種模式下SAR衛星的空間分辨率平均13 m(空間分辨率5～20 m),SAR衛星的波長為5.546 cm，選取的兩組數據的基本信息如表1所示。

表1 選取的SAR數據的基本信息Table 1 Basic information of selected SAR data

2.3 數據處理預處理

實驗預處理數據通過SNAP ESA軟件進行干涉生成、去平地效應、去地形相位，然后經過地理編碼得到，如圖6所示。

圖6 干涉數據獲取流程圖Fig.6 Flow chart of interference data acquisition

2.4 解纏算法分析

GVF-Snake模型解纏過程如圖7(a)～圖7(c)及其解纏結果如圖7(d)與本文的解纏結果如圖7(e)所示。可以看出，由于干涉條紋內部不連續性情況嚴重，致使探測邊界“不為真”，使得解纏失敗。通過本文方法的改進后，能夠正確進行解纏，得到粗解纏圖像。這些粗解纏圖像條紋邊界上分布少量孤立的散點,可以通過高通濾波插值法進行消除，得到最終的解纏圖像。

圖7 GVF-Snake模型解纏過程及融合前后解纏結果Fig.7 The unwrapping process of GVF-Snake model and the result of unwrapping before and after fusion

Bioucas-Dias等[16]將PEARLS算法用于相位重構，利用重構后絕對相位與真實地形進行比對，其誤差眾數集中在零上，且誤差上限為+0.6 mm,誤差下限-0.8 mm,在±0.3 mm之間的誤差占比為0.936，這說明該方法對地形變形重構是很有效的，可以作為檢核其他相位解纏算法的依據。在面分析上，本文列出了原始的相位圖像[圖8(a)]，選取最小費用流、高金斯枝切法等5種經典的相位解纏方法[圖8(b)～圖8(f)]與本文方法[圖8(h)]進行對比分析，其中以圖8(g)PEARLS相位重構為參照相位。

圖8 多種相位解纏圖像進行對比Fig.8 Comparisons of various phase unwrapping images

從圖8中可看出各種相位解纏算法都能夠解纏出相位的變化趨勢，但部分算法缺點表現得很明顯。其中最小費用流容易產生局部最優解，致使條紋不連續情況嚴重時，局部地區無法正確解纏；Lp最小范數法，對于整體解纏最優，忽略了細節的表達;高金斯枝切法由于受到間斷相位的干擾，在圖像下部出現解纏錯誤。解纏效果最好的三種算法，分別為本文方法、最小不連續法和Snaphu法相位解纏。

為定量評價本文相位解纏的精度，使用了5種評價指標來綜合評定解纏相位的精度。

(1)平均絕對誤差(average absolute error，MAE)。

(17)

(2)均方誤差(mean square error,MSE)和均方根誤差(root mean square error,RMSE)。

(18)

(19)

(3)絕對誤差中值(median absolute error different,MED)。

(20)

(4)相位的最大絕對誤差(max absolute error,MAX)。

(21)

式中：M是相位的行數；N是相位的列數；A為標準相位；B為待評價相位，median為取中值函數；max為取最大值函數。以PEARLS相位重構作為參照相位，將6種解纏方法的解纏圖像作為待評價相位，計算出面評價精度指標，如表2所示。

由表2可知，最小不連續法是僅次于本文方法的解纏方法，從各項指標上看均優于其他四種；Lp最小范數法雖然直觀上看上去解纏效果不佳，但其面分析指標上要比高金斯枝切法和最小費用流法，這也說明了其整體參數最優的解纏特性；其次比較好

表2 面分析評價指標Table 2 Evaluation index of surface analysis

的是Snaphu相位解纏方法，各項指標均優于最小費用流法和高金斯直切法；最好的解纏方法是本文方法，除平方根誤差未達到最佳，其他指標均優于其他五種；最差勁的是高金斯枝切法。按照面性精度排名：本文方法>最小不連續法>Snaphu>Lp最小范數>最小費用流>高金斯枝切法。

也對解纏算法進行了線分析。主要步驟如下：

(1)將解纏后的絕對相位變形方向與視線向(LOS方向)轉換到豎直變形，并去除大氣相位，使絕對相位只含有地面變形信息。

(2)將相位圖像上提取工作面以及布點信息。

(3)與實際水準測量數據進行對比。

絕對相位由視線向變形轉換到豎直向變形，即

(22)

式(22)中：v為豎直形變；d為視線向形變；λ為哨兵1?的波長；θ為入射角。

去除大氣相位，選取Yu等[17-18]提出的迭代對流層分解模型(iterative tropospheric decomposition, ITD),此法是對海拔與涌動對流層延遲部分進行耦

合所得到的模型。將ITD修正相位轉化到絕對相位下，然后通過解纏后的圖像與之相減，獲得大氣改正圖形。圖9是大氣相位改正以及改正后的相位評估。由圖9可知，大氣相位去除本質上是剔除解纏相位中的大氣差分相位，從圖9(d)中可以提取絕對相位的工作面信息和工作面布點的值。從而可以將實際的水準測量數據與相位解纏數據進行對比，如圖10所示。

圖9 大氣相位改正及改正后的相位評估Fig.9 Atmospheric phase correction and phase assessment after correction

圖10 水準數據與各解纏相位提取的變形對比Fig.10 Comparisons between leveling data and deformations extracted from unwrapping phases

由圖10可知，6種解纏方式獲取的變形，都能夠定性地反映出工作面變形情況，但是在變形量級上存在差距。其中，Lp最小范數法大部分不符合實際變形，將A、B線上所有的點作為定量分析的目標，以實際水準測量值為標準，可評價出6種相位解纏的精度指標，如表3所示。

表3 線分析指標Table 3 Line analysis index

由表3可以看出，Lp最小范數法獲得的絕對相位各向指標均為最大，因此可以確定此法不適合礦區變形監測。其他五種算法都能很好地反映線性變形情況，其中本文方法要明顯優于另外四種算法，Snaphu解纏線性性能要高于最小不連續算法，總體的線性精度排名為：本文方法>Snaphu>最小不連續>高金斯枝切法>最小費用流>Lp最小范數法。

3 結論

通過上述闡述，可以得到以下結論。

(1)融合GVF-Snake條紋探測和MRF圖切兩種方法對相位進行解纏，提出了邊界探測“為真”的判別方法，將邊界探測錯誤的內部纏繞塊用MRF圖切法相位解纏，最后按照解纏規則將外部解纏塊與內部解纏塊相融合，得到粗解纏相位(邊界處含有孤立點)。

(2)將高通濾波插值法應用于解纏邊界孤立點去除，得到了可觀的濾波結果。

(3)將ITD大氣相位模型應用到InSAR的大氣去噪，通過分析顯示，大氣去除效果比較顯著。

(4)解纏精度分析上應用了面分析指標和線分析指標，充分說明了各個解纏方法在礦區變形相位解纏方面精度性能方面的優劣性，并就其精度性能進行了排名。在面分析上，其精度性能排名為：本文方法>最小不連續法>Snaphu>Lp最小范數>最小費用流>高金斯枝切法；在線分析上，其精度排名為：本文方法>Snaphu>最小不連續>高金斯枝切法>最小費用流>Lp最小范數法。