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

劉成洲, 李 斌, 許藝騰

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

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

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

1 算法的基本原理

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

Xu等[15]針對Snake模型無法收斂至凹型邊界以及對初始化輪廓較為敏感等問題，提出了GVF-Snake模型。該模型重新定義了外部能量場，即外部能量場用全局梯度向量場來表示。具體實現(xiàn)流程：首先利用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表示沿著水平不連續(xù)參數一階鄰近，vij是沿著豎直不連續(xù)參數一階鄰近，hij、vij∈{0,1}，當hij、vij=0時，信號是不連續(xù)的，像元不連續(xù)。

式(6)各分量計算式為

(7)

(8)

(9)

(10)

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

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

算法的實現(xiàn)包括以下步驟：

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

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

(3)對“不為真”的情況，邊界探測的相位解纏已經實現(xiàn)了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為判別率，規(guī)定R=0時，探測邊界線“為真”，否則，“為假”；τ為最小距離閾值，如圖2所示。

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

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

1.4 高通插值濾波

圖像粗解纏后會在其邊界形成很多孤立點，這些孤立點相對高頻，因此可以使用高頻濾波將這些邊界點標記出來，進行去除，然后再用雙三次內插法將去除的點進行插值。實現(xiàn)步驟如圖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)進行雙三次內插法插值填補空區(qū)域，得到連續(xù)的解纏圖像。

2 實驗結果與分析

2.1 研究區(qū)域概況

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

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

2.2 數據選取

選擇sentinel-1A Interferometric Wide (IW) 模式，這種模式下SAR衛(wèi)星的空間分辨率平均13 m(空間分辨率5～20 m),SAR衛(wèi)星的波長為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)所示?？梢钥闯?，由于干涉條紋內部不連續(xù)性情況嚴重，致使探測邊界“不為真”，使得解纏失敗。通過本文方法的改進后，能夠正確進行解纏，得到粗解纏圖像。這些粗解纏圖像條紋邊界上分布少量孤立的散點,可以通過高通濾波插值法進行消除，得到最終的解纏圖像。

圖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中可看出各種相位解纏算法都能夠解纏出相位的變化趨勢，但部分算法缺點表現(xiàn)得很明顯。其中最小費用流容易產生局部最優(yōu)解，致使條紋不連續(xù)情況嚴重時，局部地區(qū)無法正確解纏；Lp最小范數法，對于整體解纏最優(yōu)，忽略了細節(jié)的表達;高金斯枝切法由于受到間斷相位的干擾，在圖像下部出現(xiàn)解纏錯誤。解纏效果最好的三種算法，分別為本文方法、最小不連續(xù)法和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可知，最小不連續(xù)法是僅次于本文方法的解纏方法，從各項指標上看均優(yōu)于其他四種；Lp最小范數法雖然直觀上看上去解纏效果不佳，但其面分析指標上要比高金斯枝切法和最小費用流法，這也說明了其整體參數最優(yōu)的解纏特性；其次比較好

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

的是Snaphu相位解纏方法，各項指標均優(yōu)于最小費用流法和高金斯直切法；最好的解纏方法是本文方法，除平方根誤差未達到最佳，其他指標均優(yōu)于其他五種；最差勁的是高金斯枝切法。按照面性精度排名：本文方法>最小不連續(xù)法>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最小范數法獲得的絕對相位各向指標均為最大，因此可以確定此法不適合礦區(qū)變形監(jiān)測。其他五種算法都能很好地反映線性變形情況，其中本文方法要明顯優(yōu)于另外四種算法，Snaphu解纏線性性能要高于最小不連續(xù)算法，總體的線性精度排名為：本文方法>Snaphu>最小不連續(xù)>高金斯枝切法>最小費用流>Lp最小范數法。

3 結論

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

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

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

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

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