基于雷达卫星时序分析技术的荆江沿岸堤防形变研究

罗天文

(1.武汉大学 遥感信息工程学院， 武汉 430079；2.长江科学院 空间信息技术应用研究所，武汉 430010;3.贵州省水利水电勘测设计研究院, 贵阳 550002)

1 研究背景

近几十年来，雷达卫星时序分析技术已应用于分析各类与不同现象相关的地面变形[1]。基于雷达卫星的时序分析技术可以区分为2大类：①永久散射体干涉测量(Persistent Scatterer-Interferometric Synthetic Aperture Radar,PS-InSAR)，基于单个主图像的使用以及时间序列中稳定和高度相干散射体的选择[2-3]。②小基线子集(Small Baseline Subset,SBAS)干涉测量，其组合多个计算的干涉图，选择以轨道(基线)之间的小空间间隔为特征对，以限制空间去相关效应[4-5]。Berardino等[6]、Ferretti等[7]提出原始PS-InSAR方法，随后几年又提出了几种改进方法。这些方法已被用于欧洲遥感(ERS)和环境卫星(Environmental Satellite，Envisat)数据的分析，证明永久散射体干涉测量(PS-InSAR)是相对低成本高精度的陆地变形监测的强大工具。该技术在大范围区域高精度形变监测领域的巨大优势，使得大量国内学者参与到雷达卫星时序分析技术的研究中，也取得了飞速的进展与卓越成就。2013年裴媛媛等[8]利用Envisat ASAR影像得到了上海市长江口南岸和杭州湾北侧堤坝的沉降速率；2019年李陶等[9]利用COSMO-SkyMed数据对公明水库路堤边坡的沉降过程进行了监测，都取得了很好的效果。得益于Sentinel-1A(欧洲航天局哥白尼计划中载有C波段合成孔径雷达的地球观测卫星)的开放存取数据，每12 d提供一次10 m分辨率的全球观测雷达数据，可以利用PS-InSAR技术高效快捷地对地表进行大面积形变监测。本文展示了如何通过永久散射干涉测量法和Sentinel-1A数据监测荆江大堤的形变，对雷达图像进行处理，提取出持久散射体，并对图像处理得到的数据进行分析建模，获得荆江大堤形变情况。

2 PS-InSAR监测形变算法

保持较高时间相干性的相干点目标为永久散射体PS(Permanent Scatterer)，PS-InSAR形变监测的基本思路是构建一系列基于相同主影像的累计差分干涉图，然后分析时序影像的幅度与相位变化，识别选取保持较高时间相干性的永久散射体，进而建立基于永久散射体的时序差分干涉模型，根据PS点的相位特征与点之间的关系，分离出有效相位，并转换为堤防形变数据。

首先在研究区域获取的雷达卫星时序数据中选择一幅合适的雷达图像作为主图像，然后将其他图像作为辅图像生成干涉图。再利用DEM(Digital Elevation Model)去除地形相位，从而得到差分干涉图。差分干涉相位Δφi可以表示为

Δφi=φtopo+φdef+φatm+φnoi。

(1)

式中：φdef为地表形变相位，φdef由线性形变相位Δφlin(x,y)和非线性形变相位Δφnon-lin(x,y)相加而成；φtopo为DEM误差相位；φatm为大气延迟相位；φnoi为噪声相位(散射体变化，热噪声，误配准造成的误差)。在雷达时序数据中选择振幅和相位稳定的点，作为候选PS点，通过设置一定的阈值，将候选PS点连成空间上相互连通的稀疏网格，其中2相邻PS点x和y的相位差Δφdiff(x,y)可以表示为

Δφdiff(x,y)=Δφtopo(x,y)+Δφlin(x,y)+

Δφnon-lin(x,y)+Δφatm(x,y)+Δφnoi(x,y)。(2)

运用干涉获取形变最重要的是获取形变相位，对于相邻的PS点来说，Δφnon-lin(x,y)很小，近似为噪声忽略不计，Δφlin(x,y)为有效形变相位，式(2)可以转换为

Δφatm(x,y)+Δφnoi(x,y)=

Δφdiff(x,y)-Δφtopo(x,y)-Δφlin(x,y)。

(3)

令Δω(x,y)=Δφatm(x,y)+Δφnoi(x,y)，则有

Δω(x,y)=Δφdiff(x,y)-(Δφtopo(x,y)+

Δφlin(x,y)) 。

(4)

求解目标函数式(4)可获得有效形变相位和DEM误差相位。选取某一PS点候选点作为参考点，在所有时序差分干涉图中迭代搜索，当Δφtopo(x,y)和Δφlin(x,y)满足所有干涉图的整体相干性最高时，则可分别作为DEM误差相位与线性形变相位的最优估计。定义PS点的相干性指标为

(5)

式中：M为干涉图的个数;xi和yi表示第i幅干涉图中相邻的2个候选PS点；j为复数的虚数单位;γs的取值范围在[0，1]，候选PS点在所有干涉图上的相位离散度越低，γs越接近1，因此可通过最大化γs来求得线性形变速率和DEM误差。

获得DEM误差相位最优估计后，再用差分干涉相位减去大气相位，利用时序差分干涉模型，对所有的像素点进行分析，当像素点的γs值满足设定的阈值时则可最终确定为PS点，然后重新构建稀疏网格，对PS点进行相位解缠，利用测量平差方法和最小二乘平差法估计所有 PS 点的线性形变速率从而获得堤防的形变量。

3 试验数据与方法

3.1 研究区域

本文选择的研究区域为下荆江石首至监利段荆江大堤。荆江位于江汉平原的南麓,是长江在湖北境内的部分河段，上起枝城,下迄城陵矶,全长约400 km，它是世界上典型的自由河曲之一[10]，河道蜿蜒曲折,特别是下荆江段。荆江大堤被列为长江防洪重点确保堤，保护了江汉平原大量城镇，被称为江汉平原的“生命之堤”，大堤的形变监测工作非常必要。

3.2 试验数据

本文利用Sentinel-1A数据进行堤防形变研究，Sentinel-1A是一颗C波段雷达卫星，于2014年4月3日发射，具有12 d重复周期，空间分辨率为10 m，中心频率为5.405 GHz。数据由欧洲航天局免费提供(https:∥scihub.copernicus.eu/dhus/)。Sentinel-1A有4种模式:条带模式(Strip Map，SM)、干涉宽幅模式(Interferometric Wide，IW)、极宽幅模式(Extra-Wide Swath，EW)和波谱模式(Wave)。本研究利用单视斜距复数影像(Sentinel-1A IW Level 1 (L1) SLC)产品。收集了2018-03-26—2019-04-02观测到的Sentinel-1A IW模式(包含VV极化方式和VH极化方式)数据共32景，数据参数如表1所示。

表1 本研究中Sentinel-1A的数据参数

3.3 试验方法

PS-InSAR处理技术分为5个阶段，具体步骤如下。

(1)配准和生成干涉图。实现PS-InSAR算法首先在研究区域获取的32幅SAR影像中选择一幅SAR影像作为主图像，然后将其他图像作为辅图像生成31幅干涉图，于2018年10月4日获得的Sentinel-1A图像具有最合适的垂直基线和一致性。因此，该图像被认为是创建干涉图堆栈的主图像，图1为干涉相对基线分布情况。

图1 干涉相对基线分布

(2)生成差分干涉图。在生成干涉图之前辅影像需要被配准到主影像。首先对图像进行粗配准，然后对高分辨率图像和相对于主图像的亚像素图像进行精配准。2个配准SAR图像的组合生成干涉图。干涉图是辅影像和主影像各值的乘积，每个像素的输出结果是新因子的关联。干涉图采用二维相位差图的形式，在本文中设置相关性阈值为0.75。在生成干涉图后，本文中利用空间分辨率为30 m的SRTM version-4数据从生成的干涉图中去除地形相位，得到没有地形相位影响的差分干涉图。

(3)选择永久散射体。得到干涉图以后需要选择候选的持久散射点，候选持久散射点的选择一般分为2个步骤。第1步根据后向散射振幅选择初始的持续散射点;第2步根据校正后的相位选择主要的持续散射点，这些点上的相位在时间上应该是稳定的。在PS点的选择上，由于地形、轨道和地面变形等因素没有考虑干涉相位，因而不能利用干涉相位的相干信息来选择持久散射点，所以采用基于振幅的方法来选择初始点。传统的基于振幅的PS点选择方法有基于振幅色散指数的振幅色散法[11]和信杂波比法。本研究采用振幅色散指数，理论上来说随着时间的推移，相位变化越小的像素，振幅色散指数也越低。通过在振幅色散指数上定义特定的阈值限制来选择初始候选像素。

(4)相位解缠。2幅SAR复图像获得的干涉相位差值是被周期折叠后位于(-π，π)之间的相位主值，想要得到正确的形变信息，必须进行相位解缠,恢复被模糊掉的相位周期，获得目标的真实相位差。关于相位解缠算法国内外进行了深入的研究，如区域增长法[12]、最小费用流法[13]、最小二乘法[14]等。在本研究中，采用最小费用流法进行相位解缠。该方法适用于由于存在大面积低相干性区域或因生长限制因素而导致解缠困难的情况，在这种情况下相对于其他方法最小费用流法的性能更好。

(5)地理编码。最后将展开相位转换为相应的形变并进行地理编码。地理编码的结果导入ArcGIS©10.5中，对PS点进行适当的插值，以便更好地展示PS-InSAR结果。通过对研究区域的地质统计，并对基于点空间分布的数据进行正态性检验，选择逆距离加权(IDW)算法，利用已知PS点单元计算未知单元，生成栅格图，作为荆江堤防的位移模型。在插值图制作完成后，利用Sentinel-1A图像处理生成的位移值，生成时间序列形变图，以便更好地比较位移随时间的变化。

4 试验结果

利用PS-InSAR技术计算的形变速率结果如图2所示。 对试验结果进行分析， 从图2的结果中可以看出， 沿荆江堤防主干区域出现不均匀的缓慢形变现象(17～30 mm/a)， 形变导致堤顶高程相对变化， 堤防的形变速率在不同区域表现较大差异， 最大形变可以达到30 mm/a。 形变在软土基础上的堤防段表现较为严重， 而混凝土基础的堤防形变程度则较轻。 荆江堤防有分布广、 岸线长的特征， 不同区域的土质不同， 形变情况也不同， 而不均匀形变会对堤防的安全性将构成危险。 基于时序SAR影像监测形变的技术可以高精度、 大面积地监测堤防形变， 为堤坝监管和防汛抗风暴潮提供了有效信息。

图2 最终形变速率结果

为了更精细地监测形变结果，分析在1 a时间里堤防的形变变化，本文在目标沿线选择采样点进行时间序列的形变分析，采样点选择的依据为：采样点必须在荆江堤防沿线，并且周围的PS点数量足够，具体表现为该点100 m范围内至少存在5个PS点，通过邻域分析可减小形变速率误差。这里选取了8个采样点进行分析，分别包含软土为基础的堤防和混凝土为基础的堤防，各类采样点的分布情况见图3。雷达数据时间序列跨度为2018-03-26—2019-04-02，平均每12 d有一次时间序列的结果，由此得到的结果如图4所示。

图3 采样点分布(A-H共8个点)

图4 采样点的时间序列形变

从图4可以发现，所有采样点都呈现出了比较明显的共同规律：①2018年9月前后，形变速率变化规律明显不同，从明显的地表沉降转变为地表的抬升；②形变速率在2018年9月左右达到最大，在2018年12月左右达到最大抬升。分析原因是江汉平原地区上半年雨水丰富，特别在7月份的梅雨季节，荆州地区经历了持续降雨，根据国家气象中心的数据，降水量达到了100～180 mm，持续的降水造成了堤防的沉降，而到了秋冬季雨水减少、气温降低、土壤凝结，从而出现了抬升的现象，到了1月份开始出现降雪，雨水增多，又出现沉降。

5 结 语

由于越来越多地使用高分辨率雷达图像，近年来在基础设施空间监测方面取得了重大进展。利用雷达卫星时序分析技术对堤防的形变进行监测，其消耗成本远比地面测量手段的消耗成本低。本文利用Sentinel-1A数据得到了石首至监利段在2018年4月至2019年4月间的平均沉降速率，结果表明，堤防在不同的时间，其形变速率不同，特别在雨季和旱季形变情况尤为不同。试验证明利用雷达卫星时序分析技术可实现荆江堤防的沉降情况高效大范围监测，表明该技术在长江堤防形变监测上具有很大的应用潜力。