SAR影像中叠掩与阴影区域的识别
——以湖北巴东为例
2019-12-03张同同杨红磊李东明李永杰刘俊男
张同同,杨红磊,2,李东明,李永杰,刘俊男
(1.中国地质大学(北京),北京 100083; 2. 资源环境与灾害监测山西省重点实验室,山西 晋中 030600; 3. 水利部水工金属结构质量检验测试中心,河南 郑州 450044)
合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)是一种新型主动式地表变形监测技术,广泛应用于地质环境灾害监测方面[1-3],其通过微波成像的方式获取影像,与可见光成像方式相比,具有全天候、无接触、大面积、高空间分辨率、高精度等优势[4-5]。
滑坡是坡体在自然因素或人为活动影响下无法维持自身的力学极限平衡状态而发生的一种地质现象[6],滑坡一旦发生会造成巨大的人员伤亡和经济损失,因此滑坡监测在国家经济发展中十分重要。由于地形和雷达卫星运行轨迹的限制,在地形起伏较大的区域会使雷达波束从斜坡顶部返回的时间比从斜坡底部返回的时间短,导致在雷达影像中斜坡的顶部和底部颠倒显示,即为叠掩现象[7-8]。而由于山体信号遮挡,导致雷达影像上阴影区域内没有数据,即为阴影现象。这两种现象都严重破坏了干涉相位的连续性,导致InSAR技术在相位处理工作中无法准确滤波和解缠。
近年来许多专家学者对叠掩和阴影区域的识别与提取做了相应的研究。文献[9]利用SAR图像中的频谱搬移特性识别叠掩和阴影区域;文献[10]针对单幅SAR影像中在瞬时照射到目标区域时会受到最小SAR天线面积的限制,使用多基线和多波段的方法解决叠掩问题;文献[11]利用极大似然估计方法,融合相位图与DEM提取叠掩和阴影区域;文献[12]采用基于干涉信号相关矩阵的特征值提取叠掩和阴影区域;文献[13]分析叠掩区域的干涉相位特性和幅度特性,提出了一种结合SAR幅度图像和干涉相位对叠掩区域进行识别的算法;文献[14]对机载SAR数据在相位域对阴影区域进行干涉处理,利用高质量DEM数据识别阴影区域;文献[15]提出一种结合SAR图像干涉幅度和相关系数来识别叠掩和阴影区域的算法;文献[16]基于Tsai方法,运用连续局部阈值分割的方法对阴影区域进行检测;文献[17]提出一种基于局部频率估计,幅度阈值分割以及结合形态学识别SAR影像中的叠掩和阴影区域的方法。本文利用组合坐标系统提取雷达监测过程中的几何模型,并结合形态学方法实现叠掩和阴影区域的识别,取得了较好的试验效果。
1 原 理
1.1 叠掩和阴影区域的检测
根据雷达采集数据时的几何模型获得该时刻的瞬时局部入射角,试验通过瞬时局部入射角判断该区域是否为叠掩和阴影区域,如图1所示。
局部入射角θ的具体表达式为
θ=β-α
(1)
式中,β为雷达波侧视角;α为目标点的坡度角。当θ<0°时该目标区域产生顶底位移,为叠掩区域;当0°<θ<90°时该目标区域产生透视收缩;当θ>90°时该目标区域产生信号缺失,为阴影区域。
根据数学模型分析可得,在雷达影像成像时,根据信号传输的距离判定目标在SAR影像上的位置,进而判别是否为叠掩和阴影区域。叠掩和阴影区域分为主动区域和被动区域,如图2所示,在检测过程中将两者区分检测,首先对主动叠掩区域进行检测,在获得第一个主动叠掩点后,向近距端方向检测被动叠掩区域,直到确定无被动叠掩区域。在检测出一个完整的叠掩区域之后向远距端方向寻找与叠掩区域初始位置距离一样的位置,该位置与主动叠掩区域之间同样为被动叠掩区域,同理可得,在检测出主动阴影区域之后向该区域远距端方向检测被动阴影区域,直到无被动阴影区域出现。
1.2 坐标系统的组合
常用的坐标系为大地坐标系、高斯平面坐标系、以及WGS-84坐标系,在这几种坐标系的基础上,本文根据实际需要将高斯平面坐标系与大地高程坐标系构成组合坐标系统,获取成像时的几何模型。
试验数据中卫星轨道坐标系为WGS-84坐标系,卫星轨道的参数文件中包含11组位置参数(X(t),Y(t),Z(t),t)。试验利用这些离散参数运用三项式拟合的方法,拟合出雷达卫星的运行轨迹,其拟合公式为
(2)
式中,a0、a1、a2、a3、b0、b1、b2、b3、c0、c1、c2、c3为待定参数。
卫星轨道坐标系为WGS-84坐标系,试验需要统一坐标系,因此将卫星轨迹坐标(xyz)转换到大地经纬度(BLH)坐标系下,其具体表达式为
(3)
式中,N为卯酉圈半径,即
(4)
式中,a=6 378 137 m;e2=0.006 694 379 990 13。
试验中大地坐标系为过渡坐标系,主要是将已知目标点的(B,L)转换到高斯平面直角坐标(x,y)下,从而解算出雷达卫星与照射目标的几何关系。
1.3 数学形态学
图像腐蚀和膨胀算法是形态学图像处理的基础,众多形态学算法也以这两种操作为基础[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。
根据形态学开运算能够去除孤立的小点和毛刺以及滤波灵活的特性,试验选择开运算去除图像中微小错误杂点,即先进行腐蚀运算再进行膨胀运算。
2 数据处理与试验分析
2.1 研究区概况
试验区域为湖北省巴东地区,地形西高东低、南北起伏,且地质条件复杂。在新构造运动及库区水位变化等外力作用下,使得该地区滑坡灾害频发。本文使用的试验数据是ALOS-2卫星采集的单视复数SAR图像和ALOS-1卫星生成的12.5 m分辨率的DEM数据。巴东地区高程介于600 m与1200 m之间,地形纵横交错,地貌环境复杂,使得SAR影像中容易产生叠掩和阴影现象,如图3所示。
2.2 方法流程
本文试验利用瞬时局部入射角的数学模型反演出SAR影像成像时的几何模型,然后进行叠掩与阴影区域的提取,并结合形态学上的开运算去除误差点,其具体流程如图4所示。
2.3 试验分析
根据强度图像特性,在叠掩区域中强度值较大,图像较亮;在阴影区域中强度值较小,图像较暗。本文运用强度影像这一特性作试验结果的对比性分析,其原始图像如图5所示。
由图5可以看出,图中叠掩和阴影区域有明显亮度差异,因此试验将所得结果叠加到强度图上作为对比性试验分析,其试验结果如图6所示。
根据图5与图6对比可知,本文试验检测的结果和强度图像大致相同,但有部分微小地区检测结果和强度图像不一致。根据实际地形和高程数据分析主要有以下两种原因:①有部分微小地区地形起伏较大,与周边整体起伏变化有较大差异,而在强度图像上信号值被被动地区掩盖,不能被分辨出来;②试验数据使用的是ALOS-1卫星获取的DEM数据,时间距今相差较大,部分地区会有植被覆盖或地形改变等因素发生。因此针对此类问题,试验运用形态学上的开运算方法消除误差点。检测结果如图7所示。
对比图6与图7可得,经过开运算处理的结果明显去除错误的检测区域,同时将检测结果更加细化,完善了检测结果的整体性,得到了良好的试验效果。
从检测结果可以看出,叠掩相较于阴影区域更容易产生。对比试验几何模型可知,试验使用的数据日期为2017年10月29日,中心雷达波侧视角为32.411 7°。通过粗略计算可得在12.5 m分辨率的DEM数据中所发生阴影现象必须保证相近两像元之间的高差在20 m以上,而在相同情况下发生叠掩现象只需高差在8 m以上,因此产生两种现象的高差相差巨大。依据本文试验,在今后的雷达影像数据采集工作中根据研究区域地形地貌的不同设计相应的雷达入射角,可以提高数据采集的质量。
3 结 论
本文利用雷达卫星成像时的几何模型与图像处理中形态学方法综合识别叠掩和阴影区域,并以巴东地区为例,验证了本文试验方法在叠掩和阴影识别领域中拥有很好的检测能力。同时随着DEM精度的更新与提高,会更加完善本文的试验结果。