張同同，楊紅磊,2，李東明，李永杰，劉俊男

(1.中國地質(zhì)大學(xué)(北京)，北京 100083； 2. 資源環(huán)境與災(zāi)害監(jiān)測山西省重點(diǎn)實(shí)驗(yàn)室，山西 晉中 030600; 3. 水利部水工金屬結(jié)構(gòu)質(zhì)量檢驗(yàn)測試中心，河南 鄭州 450044)

合成孔徑雷達(dá)干涉測量(interferometric synthetic aperture radar，InSAR)是一種新型主動式地表變形監(jiān)測技術(shù)，廣泛應(yīng)用于地質(zhì)環(huán)境災(zāi)害監(jiān)測方面[1-3]，其通過微波成像的方式獲取影像，與可見光成像方式相比，具有全天候、無接觸、大面積、高空間分辨率、高精度等優(yōu)勢[4-5]。

滑坡是坡體在自然因素或人為活動影響下無法維持自身的力學(xué)極限平衡狀態(tài)而發(fā)生的一種地質(zhì)現(xiàn)象[6],滑坡一旦發(fā)生會造成巨大的人員傷亡和經(jīng)濟(jì)損失，因此滑坡監(jiān)測在國家經(jīng)濟(jì)發(fā)展中十分重要。由于地形和雷達(dá)衛(wèi)星運(yùn)行軌跡的限制，在地形起伏較大的區(qū)域會使雷達(dá)波束從斜坡頂部返回的時(shí)間比從斜坡底部返回的時(shí)間短，導(dǎo)致在雷達(dá)影像中斜坡的頂部和底部顛倒顯示，即為疊掩現(xiàn)象[7-8]。而由于山體信號遮擋，導(dǎo)致雷達(dá)影像上陰影區(qū)域內(nèi)沒有數(shù)據(jù)，即為陰影現(xiàn)象。這兩種現(xiàn)象都嚴(yán)重破壞了干涉相位的連續(xù)性，導(dǎo)致InSAR技術(shù)在相位處理工作中無法準(zhǔn)確濾波和解纏。

近年來許多專家學(xué)者對疊掩和陰影區(qū)域的識別與提取做了相應(yīng)的研究。文獻(xiàn)[9]利用SAR圖像中的頻譜搬移特性識別疊掩和陰影區(qū)域；文獻(xiàn)[10]針對單幅SAR影像中在瞬時(shí)照射到目標(biāo)區(qū)域時(shí)會受到最小SAR天線面積的限制，使用多基線和多波段的方法解決疊掩問題；文獻(xiàn)[11]利用極大似然估計(jì)方法，融合相位圖與DEM提取疊掩和陰影區(qū)域；文獻(xiàn)[12]采用基于干涉信號相關(guān)矩陣的特征值提取疊掩和陰影區(qū)域；文獻(xiàn)[13]分析疊掩區(qū)域的干涉相位特性和幅度特性，提出了一種結(jié)合SAR幅度圖像和干涉相位對疊掩區(qū)域進(jìn)行識別的算法；文獻(xiàn)[14]對機(jī)載SAR數(shù)據(jù)在相位域?qū)﹃幱皡^(qū)域進(jìn)行干涉處理，利用高質(zhì)量DEM數(shù)據(jù)識別陰影區(qū)域；文獻(xiàn)[15]提出一種結(jié)合SAR圖像干涉幅度和相關(guān)系數(shù)來識別疊掩和陰影區(qū)域的算法；文獻(xiàn)[16]基于Tsai方法，運(yùn)用連續(xù)局部閾值分割的方法對陰影區(qū)域進(jìn)行檢測；文獻(xiàn)[17]提出一種基于局部頻率估計(jì)，幅度閾值分割以及結(jié)合形態(tài)學(xué)識別SAR影像中的疊掩和陰影區(qū)域的方法。本文利用組合坐標(biāo)系統(tǒng)提取雷達(dá)監(jiān)測過程中的幾何模型，并結(jié)合形態(tài)學(xué)方法實(shí)現(xiàn)疊掩和陰影區(qū)域的識別，取得了較好的試驗(yàn)效果。

1 原 理

1.1 疊掩和陰影區(qū)域的檢測

根據(jù)雷達(dá)采集數(shù)據(jù)時(shí)的幾何模型獲得該時(shí)刻的瞬時(shí)局部入射角，試驗(yàn)通過瞬時(shí)局部入射角判斷該區(qū)域是否為疊掩和陰影區(qū)域，如圖1所示。

局部入射角θ的具體表達(dá)式為

θ=β-α

(1)

式中，β為雷達(dá)波側(cè)視角；α為目標(biāo)點(diǎn)的坡度角。當(dāng)θ<0°時(shí)該目標(biāo)區(qū)域產(chǎn)生頂?shù)孜灰疲瑸榀B掩區(qū)域；當(dāng)0°<θ<90°時(shí)該目標(biāo)區(qū)域產(chǎn)生透視收縮；當(dāng)θ>90°時(shí)該目標(biāo)區(qū)域產(chǎn)生信號缺失，為陰影區(qū)域。

根據(jù)數(shù)學(xué)模型分析可得，在雷達(dá)影像成像時(shí)，根據(jù)信號傳輸?shù)木嚯x判定目標(biāo)在SAR影像上的位置，進(jìn)而判別是否為疊掩和陰影區(qū)域。疊掩和陰影區(qū)域分為主動區(qū)域和被動區(qū)域，如圖2所示，在檢測過程中將兩者區(qū)分檢測，首先對主動疊掩區(qū)域進(jìn)行檢測，在獲得第一個(gè)主動疊掩點(diǎn)后，向近距端方向檢測被動疊掩區(qū)域，直到確定無被動疊掩區(qū)域。在檢測出一個(gè)完整的疊掩區(qū)域之后向遠(yuǎn)距端方向?qū)ふ遗c疊掩區(qū)域初始位置距離一樣的位置，該位置與主動疊掩區(qū)域之間同樣為被動疊掩區(qū)域，同理可得，在檢測出主動陰影區(qū)域之后向該區(qū)域遠(yuǎn)距端方向檢測被動陰影區(qū)域，直到無被動陰影區(qū)域出現(xiàn)。

1.2 坐標(biāo)系統(tǒng)的組合

常用的坐標(biāo)系為大地坐標(biāo)系、高斯平面坐標(biāo)系、以及WGS-84坐標(biāo)系，在這幾種坐標(biāo)系的基礎(chǔ)上，本文根據(jù)實(shí)際需要將高斯平面坐標(biāo)系與大地高程坐標(biāo)系構(gòu)成組合坐標(biāo)系統(tǒng)，獲取成像時(shí)的幾何模型。

試驗(yàn)數(shù)據(jù)中衛(wèi)星軌道坐標(biāo)系為WGS-84坐標(biāo)系，衛(wèi)星軌道的參數(shù)文件中包含11組位置參數(shù)(X(t)，Y(t)，Z(t)，t)。試驗(yàn)利用這些離散參數(shù)運(yùn)用三項(xiàng)式擬合的方法，擬合出雷達(dá)衛(wèi)星的運(yùn)行軌跡，其擬合公式為

(2)

式中，a0、a1、a2、a3、b0、b1、b2、b3、c0、c1、c2、c3為待定參數(shù)。

衛(wèi)星軌道坐標(biāo)系為WGS-84坐標(biāo)系，試驗(yàn)需要統(tǒng)一坐標(biāo)系，因此將衛(wèi)星軌跡坐標(biāo)(xyz)轉(zhuǎn)換到大地經(jīng)緯度(BLH)坐標(biāo)系下，其具體表達(dá)式為

(3)

式中，N為卯酉圈半徑，即

(4)

式中，a=6 378 137 m；e2=0.006 694 379 990 13。

試驗(yàn)中大地坐標(biāo)系為過渡坐標(biāo)系，主要是將已知目標(biāo)點(diǎn)的(B,L)轉(zhuǎn)換到高斯平面直角坐標(biāo)(x,y)下，從而解算出雷達(dá)衛(wèi)星與照射目標(biāo)的幾何關(guān)系。

1.3 數(shù)學(xué)形態(tài)學(xué)

圖像腐蝕和膨脹算法是形態(tài)學(xué)圖像處理的基礎(chǔ)，眾多形態(tài)學(xué)算法也以這兩種操作為基礎(chǔ)[18]。

(1) 膨脹運(yùn)算

x⊕B={x|X∩Bx≠φ}={x|Bx↑X}

(5)

式(5)表示數(shù)據(jù)集X用結(jié)構(gòu)元素B進(jìn)行膨脹時(shí)，其結(jié)果為集合x，是BX與X的交集不為空的數(shù)據(jù)集，或者說x是B集中(用↑號表示)而形成的數(shù)據(jù)集。

(2) 腐蝕運(yùn)算

xΘB={x|Bx?X}

(6)

式(6)表示數(shù)據(jù)集X用結(jié)構(gòu)元素B進(jìn)行腐蝕時(shí)，其結(jié)果為集合x。

根據(jù)形態(tài)學(xué)開運(yùn)算能夠去除孤立的小點(diǎn)和毛刺以及濾波靈活的特性，試驗(yàn)選擇開運(yùn)算去除圖像中微小錯(cuò)誤雜點(diǎn)，即先進(jìn)行腐蝕運(yùn)算再進(jìn)行膨脹運(yùn)算。

2 數(shù)據(jù)處理與試驗(yàn)分析

2.1 研究區(qū)概況

試驗(yàn)區(qū)域?yàn)楹笔“蜄|地區(qū)，地形西高東低、南北起伏，且地質(zhì)條件復(fù)雜。在新構(gòu)造運(yùn)動及庫區(qū)水位變化等外力作用下，使得該地區(qū)滑坡災(zāi)害頻發(fā)。本文使用的試驗(yàn)數(shù)據(jù)是ALOS-2衛(wèi)星采集的單視復(fù)數(shù)SAR圖像和ALOS-1衛(wèi)星生成的12.5 m分辨率的DEM數(shù)據(jù)。巴東地區(qū)高程介于600 m與1200 m之間，地形縱橫交錯(cuò)，地貌環(huán)境復(fù)雜，使得SAR影像中容易產(chǎn)生疊掩和陰影現(xiàn)象，如圖3所示。

2.2 方法流程

本文試驗(yàn)利用瞬時(shí)局部入射角的數(shù)學(xué)模型反演出SAR影像成像時(shí)的幾何模型，然后進(jìn)行疊掩與陰影區(qū)域的提取，并結(jié)合形態(tài)學(xué)上的開運(yùn)算去除誤差點(diǎn)，其具體流程如圖4所示。

2.3 試驗(yàn)分析

根據(jù)強(qiáng)度圖像特性，在疊掩區(qū)域中強(qiáng)度值較大，圖像較亮；在陰影區(qū)域中強(qiáng)度值較小，圖像較暗。本文運(yùn)用強(qiáng)度影像這一特性作試驗(yàn)結(jié)果的對比性分析，其原始圖像如圖5所示。

由圖5可以看出，圖中疊掩和陰影區(qū)域有明顯亮度差異，因此試驗(yàn)將所得結(jié)果疊加到強(qiáng)度圖上作為對比性試驗(yàn)分析，其試驗(yàn)結(jié)果如圖6所示。

根據(jù)圖5與圖6對比可知，本文試驗(yàn)檢測的結(jié)果和強(qiáng)度圖像大致相同，但有部分微小地區(qū)檢測結(jié)果和強(qiáng)度圖像不一致。根據(jù)實(shí)際地形和高程數(shù)據(jù)分析主要有以下兩種原因：①有部分微小地區(qū)地形起伏較大，與周邊整體起伏變化有較大差異，而在強(qiáng)度圖像上信號值被被動地區(qū)掩蓋，不能被分辨出來；②試驗(yàn)數(shù)據(jù)使用的是ALOS-1衛(wèi)星獲取的DEM數(shù)據(jù)，時(shí)間距今相差較大，部分地區(qū)會有植被覆蓋或地形改變等因素發(fā)生。因此針對此類問題，試驗(yàn)運(yùn)用形態(tài)學(xué)上的開運(yùn)算方法消除誤差點(diǎn)。檢測結(jié)果如圖7所示。

對比圖6與圖7可得，經(jīng)過開運(yùn)算處理的結(jié)果明顯去除錯(cuò)誤的檢測區(qū)域，同時(shí)將檢測結(jié)果更加細(xì)化，完善了檢測結(jié)果的整體性，得到了良好的試驗(yàn)效果。

從檢測結(jié)果可以看出，疊掩相較于陰影區(qū)域更容易產(chǎn)生。對比試驗(yàn)幾何模型可知，試驗(yàn)使用的數(shù)據(jù)日期為2017年10月29日，中心雷達(dá)波側(cè)視角為32.411 7°。通過粗略計(jì)算可得在12.5 m分辨率的DEM數(shù)據(jù)中所發(fā)生陰影現(xiàn)象必須保證相近兩像元之間的高差在20 m以上，而在相同情況下發(fā)生疊掩現(xiàn)象只需高差在8 m以上，因此產(chǎn)生兩種現(xiàn)象的高差相差巨大。依據(jù)本文試驗(yàn)，在今后的雷達(dá)影像數(shù)據(jù)采集工作中根據(jù)研究區(qū)域地形地貌的不同設(shè)計(jì)相應(yīng)的雷達(dá)入射角，可以提高數(shù)據(jù)采集的質(zhì)量。

3 結(jié) 論

本文利用雷達(dá)衛(wèi)星成像時(shí)的幾何模型與圖像處理中形態(tài)學(xué)方法綜合識別疊掩和陰影區(qū)域，并以巴東地區(qū)為例，驗(yàn)證了本文試驗(yàn)方法在疊掩和陰影識別領(lǐng)域中擁有很好的檢測能力。同時(shí)隨著DEM精度的更新與提高，會更加完善本文的試驗(yàn)結(jié)果。