余 博,李如仁,陈振炜,张 过
(1.沈阳建筑大学 交通工程学院 辽宁 沈阳 110168;2.武汉大学 测绘遥感信息工程国家重点实验室 湖北 武汉 430079)
合成孔径雷达(Synthetic Aperture Radar,SAR)是工作于微波波段的主动式传感器,由于其全天候,全天时获取数据,并能穿透云雾,烟尘和大面积获取地表信息的特点而成为对地观测领域不可或缺的传感器[1]。合成孔径雷达对地观测技术有多种手段,其中合成孔径雷达干涉(Interferometric Synthetic Aperture Radar,InSAR)技术是一种重要的观测手段。虽然InSAR在技术上有着巨大的优势,但在2016年前,仅德国拥有TerraSAR-X高分辨率雷达卫星;加拿大拥有RADARSAT-1和RADASAT-2合成孔径雷达卫星;意大利航天局拥有高分辨率卫星COSMO-SkyMed[2];而国内能够研究的SAR卫星很少,高分三号的发射解决缺少SAR卫星的燃眉之急。
高分三号卫星于2016-8-10在太原卫星发射中心成功发射,作为首颗分辨率达到1 m的合成孔径雷达(SAR)卫星,提升我国对地遥感观测能力[3]。但高分三号作为一颗新星,在卫星对地观测应用方面国内研究少之又少。虽然高分三号卫星产生卫星影像数量多,但是能满足合成孔径干涉处理的数据却很少,所以如何集中挑选出满足干涉条件的数据是一个问题;我国科学家利用高分三号卫星数据生成我国第一幅卫星干涉SAR影像[4],并从影像中提取亚厘米级的地面沉降信息,但是并未就基于高分三号的地面形变信息提取展开研究。因此本文就以上两个问题进行研究,利用河南省登封市周边市区2016-11-29到2016-12-28高分三号数据进行合成孔径雷达干涉测量研究,并与哨兵一号数据在同一地区同一天的实验结果进行分析。
以典型的正侧视SAR成像模式为例,可将InSAR测高几何化简为二维平面模型,如图1所示。
图1 InSAR测高几何模型
图1中P为卫星两次观测的同一目标点,h为P到参考椭球面的高程,S1,S2分别代表卫星两次观测P点时的位置,R1,R2分别代表卫星S1,S2到目标点P的距离(斜距),B为S1,S2之间的空间基线矢量,并且B与水平方向的夹角(基线夹角)为α[5]。设S1为主卫星,则S1到参考椭球面的垂直高度为H,θ为S1观测目标点P时的视角,显然,在不考虑地球曲率的情况下:
h=H-R1cosθ.
微波传播路径为一个波长的距离,相位的变化量为一个周期,即为2π。因此雷达波经过斜距长度R时,相位变化量:
将卫星S1,S2对同一点目标点P所形成的相位差(φ1-φ2)记为干涉相位(φint),同时考虑到雷达波的波程从传感器到目标,再从目标点返回传感器,共经历2倍的斜距距离,在重复轨道观测的模式下,则相位差与斜距差的关系:
图像配准是干涉处理最基础的一步,复数SAR数据的干涉处理是联合两幅单视复数影像(SLC)S1和S2生成干涉图[7],需要把两幅影像配准到亚像元的精度。为提高干涉相关性超过5%,要求配准精度优于0.1个像素。互配准包括组成干涉对的两幅影像之间距离向和方位向偏移值的计算和副影像的重采样,使副影像和主影像完美地配准[8]。计算两个单视复数影像(SLC)之间的偏移值包括两步:搜索影像,整幅影像估计一些小区域局部偏移值和产生偏移多项式,多项式实现副影像匹配参考影像的重采样。
为得到配准关系的各项系数,通过在影像之间的分布的一系列窗口,为主辅影像上每个相对应的窗口估计一个平移量,通过对这一系列的平移量进行最小二乘平差处理以确定这些系数[9-10]。
形变前后获取的两幅SAR图像经过高精度的配准之后,进行相乘并取相角,生成相位干涉图。此时的干涉相位反映地表的形变信息。
U1=u1ejφ1=u1(cosφ1+jsinφ1),
U2=u2ejφ2=u2(cosφ2+jsinφ2),
|u1||u2|×ej(φ1-φ2).
解缠后的形变相位是相对值,此时的形变相位图反映整个观测区的相对形变关系,理论上需已知形变量作为控制点进行绝对形变值解算,通过人工经验指定无地表形变的点作为形变相位零值参考点,来反映整个区域的形变量。
本文以河南省登封市周边地区为例,登封地区地壳活动频繁,岩浆岩活动较为强烈,有利形成多样性矿产地质状况。
文中测区由北至荥阳市,南至新郑市,西至登封市,东至郑州市围成区域。 如图2所示。 该区域有大量的铜矿,地表表面有明显的沉降区域,适合作为实验对象。
图2 实验测量区域图
干涉必须是对同一目标采用同一角度进行观测,高分三号需要29 d才可能回到同一目标上进行相同角度观测,所以在筛选数据时,首先要确定在同一地区的两幅重轨影像数据,其次只有在雷达入射角相同的情况下,雷达接收的电磁波波谱相同,通过距离向和方位向压缩成像才有可能形成干涉对,经过实验验证,高分三号干涉实验数据必须满足两点:数据采集时间为29 d或者29 d的倍数;雷达入射角必须相同。
由于雷达飞行时间和轨道不同,得到两幅主辅影像的宽幅不尽相同,而进行干涉测量要求两幅图像大小相同,并且有重叠区域,两幅图像对于同一点的偏差小于0.1个像素,所以需要对于辅影像进行重采样,与原图像进行配准,达到两幅图像几乎完全一致。图3为主影像,图4为辅影像,辅影像的宽幅远大于主影像,而实验数据为与主图像的重叠区域为下半红色部分,也就是进行干涉测量实验的区域,如图5所示。
图3 主影像
图4 辅影像
图5 重采样后的辅影像
图6 配准后的图像
配准前进行重采样,首先计算局部偏移量,利用局部估计值生成配准多项式。在影像之间的开设256像素×256像素的窗口,然后为主辅影像上每个相应的窗口估计一个平移量,最后通过一系列的平移量进行最小二乘平差处理以确定这些参数。图6为配准后两幅影像对于同一点(图中红色圈区域)进行比较得到的结果,距离向配准精度为0.08个像素,方位向配准精度为0.03个像素,整体的配准精度为0.085个像素,小于0.1个像素,两幅图像几乎完全配准,满足配准要求。将重采样后的辅影像与主图像进行共轭相乘,形成干涉条纹图。
根据InSAR原理,得到未解缠的干涉相位图,将DEM模拟的地形相位从干涉图中去除,将形变相位进行相位滤波和最小费用流法(MCF)进行相位解缠,还原原始相位,分别得到相位解缠图,如图7所示。
图7 高分三号解缠相位图
图8 哨兵一号解缠相位图
图7为高分三号相位解缠图,图8为哨兵一号相位解缠图,由于高分三号图像存在颠倒倒置的,所以哨兵一号与高分三号的图像存在位置上的差异,从图中明显可以看出粉色区域为解缠之后的形变相位,高分三号和哨兵一号都检测出该实验区域的形变量,如图中红色圈中所选的地区,都发生了明显的形变。本文选择了一个形变区域进行研究,对比高分三号与哨兵一号对于山区的形变监测,如图9和图10所示。
图9和图10中可以看出,相同区域下,高分三号对于地形形变有着非常的敏感度,在同样的地区监测出更多地方明显的形变,如图9所示绿色和黄色圆圈区域。哨兵一号监测山区的地表(图10黄色区域)沉降量为±0.05 m,而高分三号监测山区的地表(图10黄色区域)沉降量只有±0.03 m,说明在这个地区,哨兵一号监测的地表形变与高分三号监测的形变基本一致。图11为探测区域卫星图,该研究区域(红色方框区域)为河南省槐树坪村,2015年,国土资源部发现该区域有金矿资源量接近40 t,由于大型矿场不断开采,造成不同地区地表的沉降,这些沉降符合研究结果。
图9 高分三号形变图
图10 哨兵一号形变图
图11 实验区域卫星图
1)本文利用高分三号河南省登封市周围的实验数据,通过进行合成孔径雷达干涉测量实验,获取河南省槐树坪村±0.03 m的地表沉降量,生成该地区的形变图,验证该地区由于金矿开采导致的地表沉降,同时验证高分三号作为SAR卫星具有InSAR的能力。
2)将高分三号得到的实验数据与哨兵一号得到的实验数据对同一地区进行比较,高分三号能够测量多处山区的沉降;哨兵一号作为一颗比较成熟的InSAR卫星,它能够辅助验证通过高分三号数据分析得到的形变。因此高分三号可以与哨兵一号在一定程度上配合使用,也要加大高分三号与其他InSAR卫星的联系,例如COSMO-SkyMed星、TerraSAR-X,这样高分三号可以得到更加广泛的应用。
3)尽管验证高分三号进行一定程度地干涉测量,但是只有观测周期为29 d或者29 d的倍数,并且雷达的入射角相同的数据才有可能进行干涉测量,大部分的高分数据不满足干涉测量条件,1 000多景卫星影像中只能挑出数10对干涉对,这将大大影响高分三号的应用,如何调整高分三号卫星数据质量,得到更多的干涉对,将是未来研究的重要内容。