基于轨道精炼控制点精选的极艰险区域时序InSAR地表形变监测

潘建平，邓福江，徐正宣，向淇文，涂文丽，付占宝

（1.重庆交通大学土木工程学院，重庆 400074；2.中铁二院工程集团有限责任公司地勘院，四川 成都 610031）

0 引言

合成孔径雷达干涉测量（Interferometry Synthetic Aperture Radar，InSAR）作为近年来新发展起来的一种地表监测技术，具有全天时、全天候、监测精度高、监测范围广等优点[1]，由其发展而来的合成孔径雷达差分干涉测量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）通过对同一研究区的两景SAR影像进行差分干涉，获取精度为厘米级甚至毫米级的地表形变信息，极大得提升了InSAR技术在地表形变监测领域的应用前景[2]。而以D-InSAR为基础发展起来的时序InSAR技术又解决了长期监测中时间失相干、空间失相干、大气延迟相位等因素的影响，成为了一种有效探测长时间缓慢地表形变的监测技术[3−5]。

常见的时序InSAR处理技术有永久散射体合成孔径雷达干涉测量（Persistent Scatters InSAR,PS-InSAR）技术[6]和小基线集合成孔径雷达干涉测量（Small Baseline Subset InSAR，SBAS-InSAR）技术[7]。PS-InSAR技术利用时间序列影像中的永久性散射体（Persistent Scatters，PS）点来进行时序分析，可以有效地减少时间失相干、空间失相干以及大气延迟相位对地表形变监测结果的影响。2017年白泽朝等[8]利用PS-InSAR技术成功获取了天津市2015年至2016年地面沉降速率。SBASInSAR技术通过设置时空基线阈值，从而利用短时间内时序影像的空间相干性，有效减小了时间失相干、空间失相干以及大气延迟相位对地表形变监测结果的影响，同时还增加了时间采样率[9−10]。2019年刘晓杰等[11]利用SBAS-InSAR技术和Sentinel-1A数据得到了厦门新机场海洋填海区地面沉降特征。潘光永等[12]利用SBAS-InSAR技术对济南井田矿区进行了地面沉降监测。PS-InSAR技术和SBAS-InSAR技术各有优点，因此两种技术在不同研究区域下的结合与交叉应用是当前的研究热点。2019年王舜瑶等[13]利用一种顾及永久性散射体的SBAS-InSAR方法获取了郑州市区的地表形变信息，发现形变结果与水准测量吻合较好，说明了该方法的可靠性。

文章以折多山区域为研究对象。极艰险区域永久性散射体分布稀疏，PS-InSAR技术无法获得较好的地表形变结果，同时该区域高差极大，地形复杂，不利于SBAS-InSAR技术选取稳定的地面控制点（GCP），从而降低了地表形变监测的精度。为此，采用了一种改进的SBAS-InSAR技术，通过筛选合适的PS点作为GCP点引入到SBAS-InSAR技术处理流程，从而有效提高该区域地表形变监测精度。

1 技术方法

文章采用SBAS-InSAR技术处理流程为主体，引入顾及研究区实际概况的GCP点筛选研究，提高SBASInSAR技术地表形变监测精度，技术流程如图1所示。改进的SBAS-InSAR技术原理重点涉及GCP点筛选和SBAS-InSAR技术原理两个部分。

图1 改进的SBAS-InSAR技术的基本流程图Fig.1 Basic flow chart of improved SBAS-InSAR Technology

1.1 GCP点筛选

GCP点被引入SBAS-InSAR轨道精炼步骤来进行残余相位估计和重去平处理，从而提高地表形变结果的精度。GCP点应当处于地形平缓，没有相位跃变且远离形变区域的位置，根据GCP点选取特点，可以采用由PS-InSAR技术处理后提取的PS点作为GCP候选点。

PS点是长时间序列影像中具有较高后向散射特性以及稳定性的点，其在不同时期的影像干涉对中均表现出高一致性与高相干性。振幅离差指数法[14]、相干系数阈值法[15]均是常见的PS点选取方法。

（1）振幅离差指数法

PS点较高的后向散射特性和稳定性体现在其回波的相位信息的长时间序列上具有一定的统计特性。根据这一特性，可以用振幅离差指数DA来 定量表示，具体表达示如下：

式中：DA——振幅离差指数；

σA——振幅标准差；

µA——振幅均值。

当某像素点的振幅离差指数DA的大小在0.25～0.4时，该点可以被视为PS点，但为了保证PS点的质量和最终形变监测结果的可靠性，至少需要25幅SAR影像参与计算。

（2）相干系数阈值法

相干系数是衡量干涉影像像对的干涉质量的重要指标，它主要用于描述干涉对中主、副影像同一区域的相似程度。相干系数数值分布区间为[0,1]，0表示完全不相干，1表示完全一致。根据PS点的特性，在相干性好的区域对应的相干系数也相对较高[16]，因此可以采用相干系数为PS点选取的判断阈值，其表达式为：

式中：m、n—主、辅影像内需要计算相干性的数据块大小；

i、j——数据块内的行列号；

M(i,j)、S(i,j)— 主、辅影像数据在 (i,j)位置处的复数值；

*——复数的共轭算子；

|∼|2——数据的二阶范数。

通过以任一像元为中心的i×j窗口大小进行计算，即可得到该像元的相关系数值γ 。γ值越高，像元越稳定，受噪声影响较小，干涉相位的质量越高，高于设定的相干系数阈值的像元点，即可判定为PS点。

实际处理中，联合多种PS点选取方法，通过多层阈值过滤的方式，减少各PS点选取方法考虑特性单一带来的影响，能够有助于提高PS点判断的准确性。选取出的PS点作为GCP候选点引入后续得处理中。

高山区域地形起伏较大，部分裸露岩石位于非平坦区域，此类因雷达散射特性稳定的永久性散射体，依然有被提取并作为GCP候选点的可能。然而此类GCP候选点不符合GCP点的要求，若作为GCP点参与后续流程解算，会降低形变结果精度。因此可以根据光学影像来去除此类GCP候选点，最终得到GCP点并引入SBASInSAR技术处理流程。

1.2 SBAS-InSAR技术原理

覆盖研究区的所有时间序列的SAR影像，假设有N+1幅SAR影像，其时间序列为：

通过选取其中一幅SAR影像作为超级主影像，将其他N幅SAR影像都配准到该超级主影像的雷达坐标系上来，像按照一定的时间、空间基线阈值来组合SAR影像，得到M幅干涉对，并对其进行差分干涉处理。假设第i景差分干涉图由t1，t2(t1

式中：i∈[1,2,···,M]；

dt— —t时刻相对于t0时刻的雷达视线向累积形变量；

式中AV为[M×N]的矩阵，每一行对应每个干涉对。当干涉对属于同一子基线集时，有r(A)=N，由最小二乘法可得：

式中δφ为解缠后干涉相位。然而，受时空基线的限制并非任意干涉对同属于某一小基线集，A为秩亏阵，导致最小二乘法的解不唯一。因此，采用奇异值分解（singular value decomposition，SVD）来求解A的广义逆矩阵，因而可得到最小二乘范数解：

式（7）中U为[M×N]的 正交矩阵且的对角元素为奇异值 σi(i=1,···,N)，将求相位解转换为求相位变换速率，则可得：

式（8）中B为[M×N]的 矩阵，对B进行奇异值分解，可以解出各时间段内相位变化速率v，最后计算得到最终平均形变速率以及形变速率时间序列。

2 实验过程

2.1 研究区概况

研究区位于四川省甘孜自治州康定市西面的折多山区域，面积约为1 000 km2。研究区域地形高差大，沟壑密布，山岭纵横，活动断裂带构造大量发育，属于极艰险区域。研究区地理位置如图2所示。该区域复杂的地理环境，使得常规工程测量施测难以开展，InSAR技术所具备的全天时、全天候、大范围、高精度监测的优势能在该区域得到较好的体现。

图2 研究区地理位置图Fig.2 Geographical location of the study area

2.2 数据源

本文采用由欧空局提供的Sentinel-1A雷达卫星影像升轨数据，成像模式为干涉宽视场（IW）模式，极化方式为VV同极化，重访周期为12天，时间跨度为2017年1月至2019年7月，共计76景。卫星精密轨道数据采用欧空局提供的POD Precise Obrit Ephemerides（POD精密定轨星历数据），共计76个。DEM数据采用美国航空航天局（NASA）和国防部国家测绘局（NIMA）联合测量的SRTM1，分辨率为30 m。

2.3 数据处理

实验数据处理采用了SARscape软件，按照改进后的SBAS-InSAR处理流程，得到研究区时间序列形变结果和地表平均形变速率，并采用GMT软件绘制成年均地表形变速率图。

（1）数据预处理。为减少计算冗余和处理时间，需要按研究区范围进行预先的裁剪。

（2）连接图生成。对预处理后的研究区影像进行基线估算，设置时空基线阈值，并生成连接图。本文空间基线设置为临界基线的3%，时间基线设置为48天，最终形成了289对干涉对，时空基线分布如图3、图4所示。

图3 时间基线连接图Fig.3 Time baseline connection diagram

图4 空间基线连接图Fig.4 Spatial baseline connection diagram

（3）干涉处理。对所有的干涉像对进行干涉处理，并通过查看差分干涉结果图，剔除了49对相干性差、解缠不理想的影像干涉对，为轨道精炼和重去平，以及SBAS反演估算做好数据准备。

（4）GCP点筛选。首先通过PS-InSAR技术获取PS点来作为GCP候选点，根据GCP点的特性，综合考虑振幅离差，相干性以及形变速率，来获取PS点信息。PS点振幅离差阈值在0.25至0.4之间；相干性阈值不能设置过高，否则处于高山山地的极艰险区域很难获取一定数量的PS点，需要反复试验和调整，使获取的PS点数量满足GCP点数量要求；PS点必须是稳定点，形变速率需要小于1 mm/a。因此，本文设置振幅离差阈值Dv为3.5，相干性阈值Tγ为 0.87，形变速率阈值Tv为±1 mm/a，共获取了51个PS点，也就是51个GCP候选点。然后结合光学影像查看51个GCP候选点的位置分布，发现其中9个点没有位于地势较缓的平坦区域，例如图5所示GCP点处于高山半坡，将此类GCP候选点剔除后最终得到42个稳定的GCP点，其结果与人工选取GCP点的结果对比如图6所示。

图5 结合Google earth光学影像筛选剔除的GCP点Fig.5 GCP points filtered out by Google Earth optical image

图6 GCP点选取结果对比（左为人工选取，右为本文方法选取）Fig.6 Comparison of GCP point selection results (left is manual selection, right is method selection in this paper)

（5）轨道精炼与重去平。将筛选得到的GCP点转换到SAR影像雷达斜坐标系中，采用3次轨道精炼多项式估算轨道误差和相位偏移量，消除斜坡相位，最后基于控制点对数据进行重去平。

（6）SBAS反演估算。先是通过第一次反演估算形变速率和残余地形，并设置相应的小波分解等级去除残余地形的影响；再是通过第二次反演，进行时间域高通滤波和空间域低通滤波计算，去除大气相位影响。最终得到时间序列形变结果。

（7）地理编码。将SBAS反演估算得到的形变结果转换到参考DEM坐标系下，最后通过GTM软件绘制成年均地表形变速率图(图7)。

图7 折多山区域年均地表形变速率图Fig.7 Annual surface deformation rate of zheduoshan area

3 结果与分析

对比PS-InSAR技术、SBAS-InSAR技术和改进的SBAS-InSAR技术的地表形变监测结果，可以发现研究区域基于SBAS-InSAR技术提取的地表形变点密度明显高于PS-InSAR技术获取的地表形变点密度，SBASInSAR技术监测结果中该区域中疑似滑移区域体现得更加明显和清晰。这是由于位于高山山地地带的极艰险区域永久性散体较为稀少，PS-InSAR获取的PS点信息也相对较少。区域内山体海拔较高，多存在冰雪覆盖区域，受气温影响，地表土壤呈现出较为明显的季节性反复冻融，地表形变量往往较大，在采用单一主影像的PS-InSAR技术的长时间基线干涉图中，这些区域常出现失相干现象，难以获得一定数量的PS点。通过对比分析可知，SBAS-InSAR在该研究区域的监测效果更好。改进的SBAS-InSAR技术的地表形变结果很好的保留了SBAS-InSAR技术的优势，得到了同样密度较高的地表形变点，整体的形变速率结果与SBAS-InSAR技术处理结果基本一致（图8）。

图8 不同方法的折多山区域年均地表形变速率图Fig.8 Annual surface deformation rate of zheduoshan area with different methods

进一步的分析改进SBAS-InSAR技术与SBAS-InSAR技术的处理结果，对比两者的形变速率、形变速率精度，如表1、表2所示。

表1 两种时序InSAR技术的形变速率统计Table 1 Deformation rate statistics of two time series InSAR techniques

表2 两种时序InSAR技术的形变速率精度统计Table 2 Deformation rate accuracy statistics of two time series InSAR techniques

通过对比改进前后的监测结果，该区域的平均形变速率由−1.208 922变为−1.099 689 1，标准差从9.224 748变为9.172 168。该区域形变结果表明，人工选取的GCP点处于轻微沉降区域，而GCP点代表稳定的位置，这使得使用人工选取GCP点的监测结果中形变点沉降速率相对偏大；而使用精选后GCP点的监测结果更接近形变点的真实沉降速率。并且，在形变速率精度量级上，改进后的形变速率精度平均值从改进前的3.851 148变为3.091 360，标准差从1.501 882变为1.279 466。因此结果表明，改进后的SBAS-InSAR技术的监测结果在形变速率、形变速率精度上均有较好的改善。高山山地地带的地理环境条件使得传统SBAS-InSAR技术中人工选取GCP点较为困难，而利用筛选得到的永久性散体点作为GCP点来代替人工选取GCP点，有助于提高地表形变监测精度。这说明基于轨道精炼控制点精选的SBAS-InSAR技术在极艰险区域的地表形变监测中具有更好的适用性。

4 结论

本文对比了PS-InSAR技术和SBAS-InSAR技术在折多山区域的地表形变监测结果，发现SBAS-InSAR技术在该研究区域获得的地表形变点密度更高，整体效果更好；又考虑到极艰险区域的地理条件，设计了GCP点精选方案来代替GCP人工选点，对SBAS-InSAR技术进行了改进。通过分析，本文得到以下结论：

（1）在极艰险区域的地表形变监测中，SBAS-InSAR技术能够获取较高地表形变点密度的监测结果，而PSInSAR技术的应用效果不佳。SBAS-InSAR技术在极艰险区域比PS-InSAR技术更具有应用优势。

（2）本文提出的基于轨道精炼控制点精选的SBASInSAR技术在极艰险区域的应用中，保留了SBASInSAR技术在极艰险区域地表形变监测中的优势，取得较高密度地表形变监测结果的同时，避免了GCP人工选点引起的误差，提高了地表形变监测精度，具有良好的应用价值。