張同同,楊紅磊,2,李東明,李永杰,劉俊男
(1.中國地質大學(北京),北京 100083; 2. 資源環(huán)境與災害監(jiān)測山西省重點實驗室,山西 晉中 030600; 3. 水利部水工金屬結構質量檢驗測試中心,河南 鄭州 450044)
合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)是一種新型主動式地表變形監(jiān)測技術,廣泛應用于地質環(huán)境災害監(jiān)測方面[1-3],其通過微波成像的方式獲取影像,與可見光成像方式相比,具有全天候、無接觸、大面積、高空間分辨率、高精度等優(yōu)勢[4-5]。
滑坡是坡體在自然因素或人為活動影響下無法維持自身的力學極限平衡狀態(tài)而發(fā)生的一種地質現象[6],滑坡一旦發(fā)生會造成巨大的人員傷亡和經濟損失,因此滑坡監(jiān)測在國家經濟發(fā)展中十分重要。由于地形和雷達衛(wèi)星運行軌跡的限制,在地形起伏較大的區(qū)域會使雷達波束從斜坡頂部返回的時間比從斜坡底部返回的時間短,導致在雷達影像中斜坡的頂部和底部顛倒顯示,即為疊掩現象[7-8]。而由于山體信號遮擋,導致雷達影像上陰影區(qū)域內沒有數據,即為陰影現象。這兩種現象都嚴重破壞了干涉相位的連續(xù)性,導致InSAR技術在相位處理工作中無法準確濾波和解纏。
近年來許多專家學者對疊掩和陰影區(qū)域的識別與提取做了相應的研究。文獻[9]利用SAR圖像中的頻譜搬移特性識別疊掩和陰影區(qū)域;文獻[10]針對單幅SAR影像中在瞬時照射到目標區(qū)域時會受到最小SAR天線面積的限制,使用多基線和多波段的方法解決疊掩問題;文獻[11]利用極大似然估計方法,融合相位圖與DEM提取疊掩和陰影區(qū)域;文獻[12]采用基于干涉信號相關矩陣的特征值提取疊掩和陰影區(qū)域;文獻[13]分析疊掩區(qū)域的干涉相位特性和幅度特性,提出了一種結合SAR幅度圖像和干涉相位對疊掩區(qū)域進行識別的算法;文獻[14]對機載SAR數據在相位域對陰影區(qū)域進行干涉處理,利用高質量DEM數據識別陰影區(qū)域;文獻[15]提出一種結合SAR圖像干涉幅度和相關系數來識別疊掩和陰影區(qū)域的算法;文獻[16]基于Tsai方法,運用連續(xù)局部閾值分割的方法對陰影區(qū)域進行檢測;文獻[17]提出一種基于局部頻率估計,幅度閾值分割以及結合形態(tài)學識別SAR影像中的疊掩和陰影區(qū)域的方法。本文利用組合坐標系統提取雷達監(jiān)測過程中的幾何模型,并結合形態(tài)學方法實現疊掩和陰影區(qū)域的識別,取得了較好的試驗效果。
根據雷達采集數據時的幾何模型獲得該時刻的瞬時局部入射角,試驗通過瞬時局部入射角判斷該區(qū)域是否為疊掩和陰影區(qū)域,如圖1所示。
局部入射角θ的具體表達式為
θ=β-α
(1)
式中,β為雷達波側視角;α為目標點的坡度角。當θ<0°時該目標區(qū)域產生頂底位移,為疊掩區(qū)域;當0°<θ<90°時該目標區(qū)域產生透視收縮;當θ>90°時該目標區(qū)域產生信號缺失,為陰影區(qū)域。
根據數學模型分析可得,在雷達影像成像時,根據信號傳輸的距離判定目標在SAR影像上的位置,進而判別是否為疊掩和陰影區(qū)域。疊掩和陰影區(qū)域分為主動區(qū)域和被動區(qū)域,如圖2所示,在檢測過程中將兩者區(qū)分檢測,首先對主動疊掩區(qū)域進行檢測,在獲得第一個主動疊掩點后,向近距端方向檢測被動疊掩區(qū)域,直到確定無被動疊掩區(qū)域。在檢測出一個完整的疊掩區(qū)域之后向遠距端方向尋找與疊掩區(qū)域初始位置距離一樣的位置,該位置與主動疊掩區(qū)域之間同樣為被動疊掩區(qū)域,同理可得,在檢測出主動陰影區(qū)域之后向該區(qū)域遠距端方向檢測被動陰影區(qū)域,直到無被動陰影區(qū)域出現。
常用的坐標系為大地坐標系、高斯平面坐標系、以及WGS-84坐標系,在這幾種坐標系的基礎上,本文根據實際需要將高斯平面坐標系與大地高程坐標系構成組合坐標系統,獲取成像時的幾何模型。
試驗數據中衛(wèi)星軌道坐標系為WGS-84坐標系,衛(wèi)星軌道的參數文件中包含11組位置參數(X(t),Y(t),Z(t),t)。試驗利用這些離散參數運用三項式擬合的方法,擬合出雷達衛(wèi)星的運行軌跡,其擬合公式為
(2)
式中,a0、a1、a2、a3、b0、b1、b2、b3、c0、c1、c2、c3為待定參數。
衛(wèi)星軌道坐標系為WGS-84坐標系,試驗需要統一坐標系,因此將衛(wèi)星軌跡坐標(xyz)轉換到大地經緯度(BLH)坐標系下,其具體表達式為
(3)
式中,N為卯酉圈半徑,即
(4)
式中,a=6 378 137 m;e2=0.006 694 379 990 13。
試驗中大地坐標系為過渡坐標系,主要是將已知目標點的(B,L)轉換到高斯平面直角坐標(x,y)下,從而解算出雷達衛(wèi)星與照射目標的幾何關系。
圖像腐蝕和膨脹算法是形態(tài)學圖像處理的基礎,眾多形態(tài)學算法也以這兩種操作為基礎[18]。
(1) 膨脹運算
x⊕B={x|X∩Bx≠φ}={x|Bx↑X}
(5)
式(5)表示數據集X用結構元素B進行膨脹時,其結果為集合x,是BX與X的交集不為空的數據集,或者說x是B集中(用↑號表示)而形成的數據集。
(2) 腐蝕運算
xΘB={x|Bx?X}
(6)
式(6)表示數據集X用結構元素B進行腐蝕時,其結果為集合x。
根據形態(tài)學開運算能夠去除孤立的小點和毛刺以及濾波靈活的特性,試驗選擇開運算去除圖像中微小錯誤雜點,即先進行腐蝕運算再進行膨脹運算。
試驗區(qū)域為湖北省巴東地區(qū),地形西高東低、南北起伏,且地質條件復雜。在新構造運動及庫區(qū)水位變化等外力作用下,使得該地區(qū)滑坡災害頻發(fā)。本文使用的試驗數據是ALOS-2衛(wèi)星采集的單視復數SAR圖像和ALOS-1衛(wèi)星生成的12.5 m分辨率的DEM數據。巴東地區(qū)高程介于600 m與1200 m之間,地形縱橫交錯,地貌環(huán)境復雜,使得SAR影像中容易產生疊掩和陰影現象,如圖3所示。
本文試驗利用瞬時局部入射角的數學模型反演出SAR影像成像時的幾何模型,然后進行疊掩與陰影區(qū)域的提取,并結合形態(tài)學上的開運算去除誤差點,其具體流程如圖4所示。
根據強度圖像特性,在疊掩區(qū)域中強度值較大,圖像較亮;在陰影區(qū)域中強度值較小,圖像較暗。本文運用強度影像這一特性作試驗結果的對比性分析,其原始圖像如圖5所示。
由圖5可以看出,圖中疊掩和陰影區(qū)域有明顯亮度差異,因此試驗將所得結果疊加到強度圖上作為對比性試驗分析,其試驗結果如圖6所示。
根據圖5與圖6對比可知,本文試驗檢測的結果和強度圖像大致相同,但有部分微小地區(qū)檢測結果和強度圖像不一致。根據實際地形和高程數據分析主要有以下兩種原因:①有部分微小地區(qū)地形起伏較大,與周邊整體起伏變化有較大差異,而在強度圖像上信號值被被動地區(qū)掩蓋,不能被分辨出來;②試驗數據使用的是ALOS-1衛(wèi)星獲取的DEM數據,時間距今相差較大,部分地區(qū)會有植被覆蓋或地形改變等因素發(fā)生。因此針對此類問題,試驗運用形態(tài)學上的開運算方法消除誤差點。檢測結果如圖7所示。
對比圖6與圖7可得,經過開運算處理的結果明顯去除錯誤的檢測區(qū)域,同時將檢測結果更加細化,完善了檢測結果的整體性,得到了良好的試驗效果。
從檢測結果可以看出,疊掩相較于陰影區(qū)域更容易產生。對比試驗幾何模型可知,試驗使用的數據日期為2017年10月29日,中心雷達波側視角為32.411 7°。通過粗略計算可得在12.5 m分辨率的DEM數據中所發(fā)生陰影現象必須保證相近兩像元之間的高差在20 m以上,而在相同情況下發(fā)生疊掩現象只需高差在8 m以上,因此產生兩種現象的高差相差巨大。依據本文試驗,在今后的雷達影像數據采集工作中根據研究區(qū)域地形地貌的不同設計相應的雷達入射角,可以提高數據采集的質量。
本文利用雷達衛(wèi)星成像時的幾何模型與圖像處理中形態(tài)學方法綜合識別疊掩和陰影區(qū)域,并以巴東地區(qū)為例,驗證了本文試驗方法在疊掩和陰影識別領域中擁有很好的檢測能力。同時隨著DEM精度的更新與提高,會更加完善本文的試驗結果。