劉成洲, 李 斌, 許藝騰
(中交天津港灣工程研究院有限公司巖土工程研究所, 天津 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ū)變形精細化解纏。
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)
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表示水平向和豎直向。
算法的實現(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中紅線表示其實際探測曲線,黑色線表示通過多次試驗獲得最小距離閾值曲線,藍色點表示手動選取得初始化邊界點,黑色箭頭表示初始化邊界點到探測邊界得最小距離,當黑色箭頭線在兩條閾值邊界之間,表示條紋探測線與真實邊界相符,否則條紋探測失敗。
圖像粗解纏后會在其邊界形成很多孤立點,這些孤立點相對高頻,因此可以使用高頻濾波將這些邊界點標記出來,進行去除,然后再用雙三次內插法將去除的點進行插值。實現(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ù)的解纏圖像。
研究區(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
選擇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
實驗預處理數據通過SNAP ESA軟件進行干涉生成、去平地效應、去地形相位,然后經過地理編碼得到,如圖6所示。
圖6 干涉數據獲取流程圖Fig.6 Flow chart of interference data acquisition
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最小范數法。
通過上述闡述,可以得到以下結論。
(1)融合GVF-Snake條紋探測和MRF圖切兩種方法對相位進行解纏,提出了邊界探測“為真”的判別方法,將邊界探測錯誤的內部纏繞塊用MRF圖切法相位解纏,最后按照解纏規(guī)則將外部解纏塊與內部解纏塊相融合,得到粗解纏相位(邊界處含有孤立點)。
(2)將高通濾波插值法應用于解纏邊界孤立點去除,得到了可觀的濾波結果。
(3)將ITD大氣相位模型應用到InSAR的大氣去噪,通過分析顯示,大氣去除效果比較顯著。
(4)解纏精度分析上應用了面分析指標和線分析指標,充分說明了各個解纏方法在礦區(qū)變形相位解纏方面精度性能方面的優(yōu)劣性,并就其精度性能進行了排名。在面分析上,其精度性能排名為:本文方法>最小不連續(xù)法>Snaphu>Lp最小范數>最小費用流>高金斯枝切法;在線分析上,其精度排名為:本文方法>Snaphu>最小不連續(xù)>高金斯枝切法>最小費用流>Lp最小范數法。