罗雪玮， 向喜琼，2， 吕亚东

(1.贵州大学喀斯特地质资源与环境教育部重点实验室,贵阳 550025；2.贵州大学资源与环境工程学院，贵阳 550025)

0 引言

为克服差分合成孔径雷达干涉测量技术(differential interferometric synthetic aperture Radar，D-InSAR)时空失相干、大气效应及噪声等因素的影响，一系列的时序合成孔径雷达干涉测量技术(interferometric synthetic aperture Radar，InSAR)技术应运而生。具有代表性的有Stacking-InSAR、永久散射体(permanent scatter，PS)PS-InSAR、小基线集(small baseline subset，SBAS)SBAS-InSAR以及相干永久目标分析(interferometric persistent target analysis)IPTA-InSAR等技术[1]。Colombo等[2]利用ESA ERS1/2卫星采集的合成孔径雷达(synthetic aperture Radar，SAR)数据，采用D-InSAR和PS-InSAR进行了干涉分析，在监测的工业区内发现了不同的地面沉降模式； Wu等[3]采用SBAS-InSAR方法对延安新区2015—2019年的地表变形进行了调查，发现区内地基沉降主要来源于填方和城市荷载的作用，而建筑物的应力释放是造成地面隆起的主要因素； 刘一霖等[4]基于SBAS-InSAR技术，引入偏移量追踪法、FFT过采样方法、多窗口迭代自适应滤波技术对矿区进行了地表时序监测，分析了研究区开采工作面内地表的时空变化规律； 葛大庆等[5]融合了PS-InSAR和SBAS-InSAR两种方法，以短基线差分干涉生成相位图后利用点目标识别技术提取出相干目标，分析了德州地区地面沉降-回弹的动态变化过程； 李金超等[6]将基于灰色模型和支持向量回归的组合模型(gray model and support vector regression，MG-SVR)与SBAS-InSAR技术相结合，对淮南潘集矿区发生沉降的区域进行了监测及预警。

在使用SBAS方法对地表形变进行监测时，需要在轨道精炼步骤中引入稳定的地面控制点(ground control point，GCP)估算和去除残余的恒定相位和解缠后还存在的相位坡道，从而使形变监测结果更为准确。但在没有精确GCP的情形下，如果研究区内存在较大高差、或者影像数据噪声过大，无法进行消除，人工来选择GCP可能会引起较大误差。

针对此问题，本文提出了一种基于PS的SBAS-InSAR方法。该方法通过设置相干系数的阈值、振幅离差指数的阈值以及地表形变速率的阈值选出稳定的PS点，并将这些点作为SBAS-InSAR轨道精炼中的地面控制点，以期能得到准确度更高的地表形变监测结果。最后，以位于贵州省龙里县洗马镇上石坎村的一处地面塌陷为实例验证了本文方法的有效性。

1 研究方法

1.1 相干系数法

相干系数法的基本思路为： 对于某一像元选定一定大小的移动窗口(m,n)，利用该移动窗口内的像素复数信息来估计相干系数值。计算完毕后，移动该窗口，重复计算[7-8]。相关系数γ计算公式为：

(1)

式中:M为主影像;S为辅影像;*为共轭相乘。

(2)

1.2 振幅离差法

Ferretti等[9]提出可以采用像素中振幅信息的时间序列来度量干涉相位的稳定性。R和I分别为影像的实部和虚部，若存在标准差为σn的高斯噪声，则振幅值A服从Rice分布，公式为：

(3)

(4)

式中：v为干涉相位；σv为相位离差指数；σnI为虚部的标准差；mA和σA为时序振幅的均值和标准差；DA为振幅离差指数。

1.3 相干性、振幅离差指数以及形变速率结合的选点方法

首先采用相干系数阈值法初步识别出具有高相干系数的像元的PS点目标； 然后设定振幅离差值进一步筛选具有稳定强散射性的PS点作为目标； 最后设定一个形变速率区间来进行最终的PS点选取。如图1所示。

图1 方法流程Fig.1 Method and procedure

1.4 SBAS-InSAR技术

SBAS-InSAR技术可以限制时空基线失相干和大气延迟的影响[10]。将N+1景覆盖同一区域的SAR影像序列，影像获取时间为(t0,…，tN)，根据空间基线和时间基线的条件组合干涉对得到M幅差分干涉图，M的数量与影像数量满足以下关系[11]：

(5)

假设干涉图j由tA和tB时刻获取的影像干涉生成，在去除平地效应与地形相位影响后，干涉图j中距离向为r与方位向为x的某一像素的干涉相位表示为：

(6)

式中：φ(tB,x,r)和φ(tA,x,r)分别为tB和tA时刻SAR影像上的相位值；φdef,j(x,r)为tA～tB时刻内卫星视线向(los)的形变相位；φtoppo,j(x,r)为地形相位误差；φatm,j(x,r)为大气相位误差；φnoise,j(x,r)为噪声相位。

对所有干涉图构建差分干涉相位与影像获取时刻形变量之间的线性方程组，其矩阵形式为：

(7)

式中:A为M×N的系数矩阵;φ为参数矩阵，由未知像元(x,y)在N个时刻对应的未知形变相位值构成；δφ为由M个式(6)所示的δφj组成的观测方程组；M为差分干涉图的数量。

为求解研究区域各高相干点的形变速率，可用两景影像间的平均相位速率代替相位值，则式(7)变为：

(8)

式中：B为M×N的系数矩阵；v为形变速率向量，vΤ可表示为：

(9)

当系数矩阵B为满秩(即M≥N)时，可用最小二乘法求解形变速率； 当M

2 实验数据和处理

2.1 数据

本文以龙里县洗马镇的一处地面塌陷为实例(图2)，选取20景由欧洲航天局提供的C波段Sentinel-1A SAR影像数据，时间范围为2019年9月1日—2021年4月11日。选用POD精密定轨星历数据去除轨道相位，采用SRTM-DEM去除地形相位并进行坐标系转换。表1列举了Sentinel-1A影像数据的基本参数。

图2 Sentinel-1A影像位置Fig.2 Sentinel-1A image position

表1 Sentinel-1A数据参数Tab.1 Sentinel-1A parameter

2.2 PS点的选取

首先对20景SAR数据进行配准、干涉、去平处理并计算出数据的时序振幅标准差与平均值。先采用相干系数法识别出具有高相干系数的像元作为PS点选取目标，考虑到研究区为城镇，建筑物较多，故设置相干性阈值为0.75； 然后使用振幅离差指数进一步筛选出具有稳定相位的PS点，根据Ferretti等[9]的研究，通常振幅离差指数只考虑小于0.25的像素点，故设置振幅离差指数阈值为0.25； 最后基于形变速率来确定PS点，因为本文选取的PS点需要作为SBAS-InSAR轨道精炼时的GCP点，而GCP点的基本选取原则为远离形变区域，故本文将形变速率区间设置为[-0.1，0.1] mm/a。基于以上3个限制条件，共选出相位稳定的PS点36个，如图3所示，且将PS点导入地图对比查看后发现大部分PS点都位于稳定的建筑物(道路、屋顶等)上。

图3 基于3个阈值得到的PS点Fig.3 PS points based on the three thresholds

2.3 SBAS-InSAR处理

将得到的PS点由地理坐标转为雷达斜距坐标，然后将其作为GCP点与干涉处理后的像对进行轨道精炼，估算和去除残余的恒定相位和解缠后还存在的相位坡道； 接着通过2次反演估算形变速率和残余地形，并去除大气相位； 最后获取了研究区地表形变速率(图4)和时序形变结果(图5)。可以看出： 区内平均地表形变速率在-6.5～5.7 mm/a之间，年平均最大沉降速率为-6.52 mm，说明该地变形程度轻微，且监测到的形变主要集中在塌陷区边界及道路周围，与实地地质调查的结果相符合。除开植被覆盖区域出现失相干外，塌陷区内大部分农田没有监测到相关形变是由于春耕时节居民对农田进行翻土种地等操作，农田里长出了新作物，所以也导致了该区域的部分失相干情况。以2019年9月1日为第一时相，其余各期影像相对于第一时相的形变如图5所示。

图4 SBAS-InSAR形变速率图Fig.4 Deformation rate diagram of SBAS-InSAR

图5 研究区地表形变时间序列

3 结果评述

为了验证该方法的优越性，本文使用相同时间段的SAR影像数据，基于选择GCP点的重要标准[12]： ①没有残余地形条纹； ②没有形变条纹，远离形变区域； ③没有相位跃变； ④至少20～30个点。通过对比滤波后的干涉图和解缠后的干涉图，在没有相位跃变和残余地形相位的区域，人工手动选择了30个控制点，得到研究区的形变速率。并收集了研究区内设置的3处北斗位移监测站点的监测结果(选取时间分别为2020年11月6日，2020年12月12日，2021年1月17日，北斗监测结果的时间跨度为2020年10月28日—2021年11月25日，与本文InSAR时间序列相对一致)，然后在2种SBAS形变结果中选取了监测站点周围不超过50 m的3个变形点，计算这3个变形点在同一位移监测时间内的相对垂直向形变量，并与位移监测数据进行比较，比较结果如表2，图中点号为北斗监测的站点名。

表2 2种方法与北斗监测点比较情况Tab.2 Comparison between the two methods and the Beidou monitoring station (mm)

由表2看出，2种方法中，结合了PS点的SBAS-InSAR方法的监测结果比人工选取GCP点的传统SBAS-InSAR方法更接近位移监测数据，且2种监测结果与位移监测数据均有一定的误差。通过查阅文献再结合实地考察，分析造成误差的原因主要有以下3点： ①由于客观条件的限制，研究区内北斗位移监测站点的设置与所选择的InSAR形变点的位置并不是完全重合的，故选择了在监测站点50 m以内的3个邻近InSAR形变点来进行对比分析； ②InSAR形变结果的分析处理过程是通过像元而不是单个点来处理的，像元覆盖了一定的范围，所以得到的结果可能会受到周围临近像元及整体区域上的影响而造成误差； ③除了InSAR技术本身存在的误差源，在时序InSAR中还存在雷达影像对精度的影响、干涉相位对函数模型观测量与未知量之间函数关系的影响[13]。InSAR的相位组成可通过函数表现为以下形式：

(10)

式中：ε为失相干随机噪声，E{φ}=0。基于时序InSAR技术，可以将多幅干涉图组成其函数模型观测量，增加多余观测，从而进行误差的解算与分离。

4 结论

1)在进行SBAS-InSAR处理时，为了使选择的GCP点不引起较大误差从而影响到形变结果，本文通过设置相干系数的阈值、振幅离差指数的阈值以及地表形变速率的阈值选取了稳定的PS点，然后将其作为GCP点进行轨道精炼及后续的反演处理得到地表形变信息。并将本文方法与采用人工选取GCP点的方法所得结果和北斗位移监测结果进行对比。由本文方法得到的形变结果比通过人工选择GCP的SBAS方法得到的形变结果更为精准，能更加有效地分析研究区地表时序形变。在实例研究中，研究区内平均地表形变速率在-6.52～5.66 mm/a之间，年平均最大形变速率为-6.52 mm，该地整体变形程度轻微，变形主要集中在塌陷区边界及道路周围。

2)由于本文方法所选取的PS点大部分都位于建筑物上，当建筑物没有发生竖直方向上的变形时，也会被判定为PS点。那么基于这样的PS点得到的时序形变结果也可能会出现较大误差。因此，后续的研究可以通过将升降轨数据进行融合得到三维的形变量，再进一步解算水平与竖直方向的形变，通过共同判断x和y方向上的位移，这样得到的永久散射体就会更为可靠。