1. 引言
卫星合成孔径雷达干涉(Interferometric Synthetic Aperture Radar, InSAR)是一种重要的空间对地观测技术。相比于传统测量技术,这种技术具有观测范围大、空间分辨率高、自动化程度高等优点 [1] [2] [3] [4] 。InSAR技术起初主要应用于生成数字高程模型和制图。后来,随着差分雷达干涉技术(Differential InSAR,DInSAR)的出现,其应用领域扩展到测量地表形变 [5] [6] [7] 。但是DInSAR的应用受到时空失相关和大气延迟的负面影响,限制了其在应用中的精度和可靠性 [4] 。有研究表明,在二轨差分干涉测量中,20%的湿度变化会导致10 cm的形变误差,会影响毫米到厘米级的InSAR形变监测结果的可靠性 [8] 。
针对常规DInSAR中的应用局限性,国内外诸多学者开展了基于多时相卫星SAR影像探测形变时空演变的研究,提出了相应的时序差分雷达干涉理论与方法。其中最典型的两种方法是:Ferretti等人率先提出的永久散射体雷达干涉(PSI)和Berardino等率先提出的短基线集(SBAS)干涉 [9] [10] [11] 。在PS和SBAS方法的基础上,Hooper提出一种提高可用信号空间采样率的方法,优化并融合两种方法提取的高相干点,构建模型解算形变信息 [12] [13] [14] 。虽然构建形变解算模型时有顾及大气延迟的影响,但是大气延迟产生的误差可能达到十几厘米。误差仍会通过模型传递给形变解算结果。因此,有必要在解算前对大气延迟进行校正,提高多时相InSAR形变监测的可靠性。
较早的用于缓解大气延迟的方法是层叠法 [8] 。该方法基本原理是大气信号在时间尺度上超过一天的间隔是不相关的,因此可以通过平均独立的干涉图的方法来缓解大气延迟问题。但是该方法局限是通过大量干涉图进行平均,且不能获得时序形变。目前在时序上缓解大气延迟的方法是应用滤波技术,例如Berardino等人在2002年提出一种通过时间域和空间域滤波提取大气延迟相位的方法 [11] 。该方法根据大气延迟具有高空间相关性和低时间相关性,在先移除模型评估出的形变信号的前提下,通过进行空间域低通滤波和时间域高通滤波的处理,提取出大气延迟相位。但是这种滤波方法在对噪声进行平滑处理的同时也对形变信号进行了平滑处理,而且去除的大气延迟相位受主观因素影响。此外,还有许多直接估计并去除大气延迟的方法,例如基于GPS数据或雷达数据辅助InSAR大气校正的方法和基于气象数据的大气模型校正方法 [15] [16] [17] [18] 。但是,基于GPS数据的方法受站点的位置和密度的局限,MODIS和MERIS则不能在夜间进行观测且受云雾影响很大,而基于数值大气模型的校正方法的稳定性相对较差。
大气延迟可以分为与空间异质分布有关的延迟和与大气垂直特性有关的延迟 [19] 。其中,与大气垂直特性有关的大气延迟称为对流层垂直分层延迟,与高程有关。大气延迟受水汽含量影响最大,而水汽含量随高程呈指数级的下降 [20] 。在此基础上,Lin等用泰勒级数将指数表达式展开并保留一次项,得到大气延迟与高程的线性关系式,提出了一种使用多尺度的方法线性评估大气延迟 [21] 。该方法侧重于减轻对流层垂直分层延迟的影响,能够在不借助外部先验数据的条件下,实现精准的多时相合成孔径雷达干涉测量大气校正。
本文以Sentinel-1A SAR影像为数据源,以成都市主城区为实验区域,在MT-InSAR技术提取出高相干点的基础上,引入多尺度稳健估计模型优化MT-InSAR大气延迟分量的解算和校正,大气校正后相位均方根误差总体上减少。探索2016年2月至2017年9月地表形变情况,并用同时段CORS数据进行验证,以证明基于多尺度的MT-InSAR大气延迟校正方法是有效和可靠的。
2. 技术原理
Hooper在PSI方法提取PS点的基础上,提出可以利用短基线干涉对里的滤波后相位微失相干(slowly-decorrelating filtered phase, SDFP)点。这些点在短时间域内滤波处理后仍保持较高相干性,可以与PS点一起通过优化筛选,融合成高相干点集,用于SBAS模型解算形变信息。在此基础上,本文引入多尺度稳健估计模型优化MT-InSAR大气延迟分量的解算,再对校正后的高相干点集进行相位解缠。经完善后的MTI形变提取流程图如图1所示。
假设有覆盖研究区域的p+1幅SAR影像,选择一幅作为公共主影像,将其余影像分别与主影像配准并重采样到主影像空间,与单一主影像干涉得到PS差分干涉对,同时选取合适的时空基线阈值,构建稳健且无孤立子集的短基线集。经过去除参考椭球面和地形相位处理,分别得到两种二次差分干涉图集(PS二次差分干涉图集和SBAS二次差分干涉图集)。通过设置振幅离差(Amplitude Dispersion, AD)和振幅差分离差(Amplitude Difference Dispersion, ADD)阈值分别提取PS候选点集和SDFP候选点集。振幅离差指数的计算公式为
(1)
其中,
为像元在p幅PS二次差分干涉图中的振幅标准差,
为像元在p+1幅干涉图中的振幅平均值。振幅差分离差的计算公式为
(2)
其中,
为像元在各SBAS二次差分干涉图的振幅差标准差;
为像元在各SBAS二次差分干涉图的振幅差平均值。
![](//html.hanspub.org/file/6-2840170x16_hanspub.png)
Figure 1. Flow chart of deformation extraction
图1. 形变提取流程
对干涉图低通滤波后,将信噪比低于阈值的PS候选点和SDFP候选点去除,得到高相干点集。把信噪比作为定权的依据,在空间域上按照加权平均的方式计算各候选点的相位,用于之后的相位解缠。
在相位解缠后,我们引入线性模型校正与高程相关的大气相位。模型根据大范围区域大气的全局属性建立,由整幅差分干涉图线性回归分析得到,所以不会影响局部形变与高程的关系。假设干涉图仅受静态大气延迟的影响,那么差分干涉图相位表达式如下
(3)
其中,
为大气延迟;b是模型的偏移常量;K是模型的系数,描述对流层垂直分层延迟与高程的线性关系;h是高程。
然而,实际情况下差分干涉图受其他干扰因素的影响,很难定义大气相位与高程的线性关系,所以需要一种稳健的方法去评估出模型参数。有关文献中提到,不同的空间尺度对不同的噪声敏感性不同 [21] 。在一定空间尺度范围内,差分能很大程度上抑制噪声的影响,而对流层垂直分层延迟与高程的线性关系不会变。
基于此,本文采用多尺度的方法稳健地评估模型的参数。首先,用重采样和滤波的方法将高程图和差分干涉图分解到不同尺度的空间,生成不同空间尺度下高程图和差分干涉图像,然后分别计算相邻两个尺度间的差值,最后运用选出的高相干点评估模型参数。假设通过不同空间尺度l的滤波和差分处理后,得到高程图和差分干涉图各n组,式3可以写为
(4)
其中,
、
分别为第i(i ≤ n)组范围尺度下的相位差值和高程差值。
对于同一地区m幅差分干涉图,用影像获取时间间隔去定义一组模型参数分量,每一个时间段都有对应的转换系数因子分量和偏移量分量。p + 1幅影像对应p段时间间隔,可以得到线性模型新的表达式
(5)
其中,
为第q(q ≤ m)幅差分干涉图在第i组范围尺度下的相位;bj和kj分别为第j(j ≤ p)个时段对应的模型参数分量;aj为差分干涉图跨越第j个时段的判定系数,取值为0或1。式5的矩阵表达式为
(6)
考虑到解缠相位中噪声的影响,采用最小二乘方法评估出模型参数分量的最优解,根据干涉对的时间跨度可以计算得到每幅干涉对的大气校正模型参数,进而求得大气相位。对高相干点集进行大气相位校正后,构建SBAS形变解算模型,即可求解得到高相干点集雷达视(LOS)方向的年平均形变速率和形变时间序列。
3. 实验分析
成都市地处四川盆地西缘,毗邻龙门山断裂带,多阴雨天气,水汽的不均匀分布和对流层垂直分层效应都对MT-InSAR构成挑战,是理想的实验对象。本文选取实验区域范围大小约为35 × 24 km2,位于成都市区周边(详见图2)。
观测数据选用欧洲空间局哥白尼计划发射的工作波长为5.6cm的Sentinel-1A卫星雷达数据,观测时段从2016年2月6日至2017年9月16日,获取影像共计14景,方位向和斜距向分辨率分别为13.98 m和2.33 m。为对雷达影像坐标转换,实验采用SRTM提供的空间分辨率为30 m的DEM数据,用于模拟地形相位和坐标转换。对于用振幅离差法选取PS候选点时,选取以获取时间为2016年11月20日的影
像为主影像,得到13幅单一主影像差分干涉图。对于用振幅差分离差值法选取SDFP点时,通过时空基线选取,得到65幅短基线干涉对,SBAS时空基线如图3所示。
实验通过融合振幅离差法和振幅差分离差值法提取到的PS点和SDFP点,得到1,099,478个高相干点,高相干点密度为1340个/平方公里。提取过程中,ADI阈值为0.4,ADD阈值为0.6。高相干点地物除了包括人工建筑物等优势散射体,也包括在一段时间内保持较高相干性的散射体。在相位解缠后,按照本文提出的方法,根据实验区域的覆盖范围,将高程图和短基线干涉图通过滤波分解到不同的空间尺度。根据实验区域范围大小,选取2000 m、4000 m、8000 m、16,000 m为空间尺度,相邻空间尺度进行差分后得到3个尺度范围下的高程图和短基线干涉图,将3个尺度范围下65幅干涉对的高相干点代入式6求解模型参数,即可得到线性大气校正模型。用线性模型进行大气校正前后相位对比图如图4所示。
图4中,a、b、c分别为65个SBAS干涉对中的3个干涉对大气校正前的相位图;d、e、f分别为三幅干涉图用线性模型模拟的大气相位图;g、h、i分别为三幅干涉图大气校正后的相位图。从d、e、f可以看到,线性模型模拟的大气相位与高程相关,对于工作波长为C波段(5.6 cm)的Sentinel-1A卫星来说,其改正值约在0~30.2 mm之间。模型参数K(单位:rad/km)由多尺度的高程差分图和相位差分图评估得到,它能反映该区域对流层垂直分层相位延迟的平均值。对比大气校正前后的相位图,可以较明显的看出相位变化梯度有一定的改善。大气校正前后相位的均方根误差(表1)对比也表明,大气校正后相位均方根误差有所减小,说明干涉图相位中确实存在对流层垂直分层延迟相位,证明本文提出的方法是有效的。值得一提的是,大气校正后有的干涉对均方根误差增加,说明在这些干涉对中,与空间异质分布有关的延迟占主导地位。
在进行大气校正后,构建SBAS形变解算模型,对大气校正后的高相干点集形变解算,得到实验区域在该时段内年平均沉降速率图如图5所示。由图5可以看出,实验区域的年平均形变速率在6 mm/y以内,整个城区无相对位移趋势,几乎没有沉降。其中,CHDU为一个CORS站点,位于城区形变稳定的区域;五角星R为形变解算的参考点;a、b、c分别为农田、湖泊和施工用地,相干性有所降低。
![](//html.hanspub.org/file/6-2840170x26_hanspub.png)
Figure 3. Temporal and spatial baselines of SBAS
图3. SBAS时空基线
![](//html.hanspub.org/file/6-2840170x27_hanspub.png)
Figure 4. Comparison of phase before and after atmospheric correction
图4. 大气校正前后相位对比
![](Images/Table_Tmp.jpg)
Table 1. Comparison of RMSE of phase before and after atmospheric correction
表1. 大气校正前后相位均方根误差对比
实验通过CORS的形变观测数据对形变解算结果进行精度验证。由于CORS站点数量有限,实验假设在监测起始时间,CORS和对应高相干点累积形变量为0 mm,将LOS向形变监测结果投影到竖直方向,计算两类监测结果差值的均方根误差来验证解算结果的可靠性。两类监测结果差值RMSE如表2所示。
验证结果表明,两类形变监测结果差值的标准差在4 mm以内,证明了本文所提出模型与方法是有效而可靠的。
4. 结论
本文以Sentinel-1A SAR影像为数据源,以成都市主城区为实验区域,基于MT-InSAR技术提取出高相干点集,引入多尺度稳健估计模型优化MT-InSAR大气延迟分量的解算和校正。模型评估的大气延迟在0~30.2 mm,大气校正后相位均方根误差总体上减少,说明多尺度方法评估出的线性模型应用于大气校正是有效的。大气校正后有的干涉对均方根误差增加,说明该模型局限于校正大气延迟中的垂直分层
![](//html.hanspub.org/file/6-2840170x28_hanspub.png)
Figure 5. Average deformation rate of LOS direction
图5. LOS向平均形变速率
延迟。利用MT-InSAR技术提取出的地表形变结果分析表明,整个实验区域较为稳定,年平均形变速率多在6 mm/y以内,用同时段CORS数据进行验证,两类形变监测结果差值的均方根误差在4 mm以内,证明了基于多尺度的MT-InSAR大气延迟校正方法是有效和可靠的。
致谢
本文获得了国家重点研发计划“地球观测与导航”领域重点专项课题(2017YFB0502704),国家自然科学基金青年科学基金项目(41601503),西南交通大学理工类科技创新项目(2682016CX087)和西南交通大学“雏鹰学者”人才计划项目(2682016CY19)的联合资助。此外,美国地质调查局和欧洲空间局分别为本研究提供了数字高程模型(SRTM DEM)和Sentinel-1A卫星影像数据,成都市勘察测绘研究院提供了CORS监测结果,在此一并表示感谢。
NOTES
*通讯作者。