利用渐进式SBAS-InSAR监测分析北京平原地区地表形变

2023-02-28刘媛媛夏元平赵振宇

大地测量与地球动力学 2023年3期

刘媛媛晏霞夏元平,3 赵振宇

1 东华理工大学江西省数字国土重点实验室，南昌市广兰大道418号，330013 2 东华理工大学测绘工程学院，南昌市广兰大道418号，330013 3 自然资源部环鄱阳湖区域矿山环境监测与治理重点实验室，南昌市广兰大道418号，330013

传统SBAS-InSAR技术通过同时处理覆盖研究区一定时间段内的所有SAR影像来获取该时段内的地表形变结果[1]。若新增SAR影像，则需将所有影像重新处理以获取新增SAR影像时刻的地表形变信息，难以满足高效、持续、动态的监测要求。为解决这一难题，许多学者提出各种解决办法[2-5]，取得了较好效果。本文提出渐进式SBAS-InSAR技术，即在传统SBAS-InSAR技术获取的解算结果基础上，将较新的少部分已有SAR影像结合新增SAR影像进行增量解算，以得到所有影像获取时刻的地表形变结果，并在绝大部分已有数据形变参数和已有数据参数矩阵不变的情况下，实现高效、持续、精确的地表形变监测。以2019-06～2020-04覆盖北京平原地区的26景Sentinel-1A SAR影像为例，利用渐进式SBAS-InSAR进行实验，并从精度和计算效率2个方面与SBAS-InSAR地表形变监测结果进行对比。

1 渐进式SBAS-InSAR原理

如图1所示，渐进式SBAS-InSAR分为2个阶段：1)利用SBAS-InSAR求解已有SAR影像数据集的地表形变时间序列(阶段1)；2)以阶段1获取的已有SAR影像数据集的形变参数为基础，结合已有SAR影像数据集中较新部分SAR影像与新增SAR影像进行高相干点选取、相位解缠、大气延迟、轨道误差和高程误差纠正，获取新增干涉对的解缠相位，并以此作为观测值进行序贯平差解算，得到更新后全部影像时刻的地表形变时间序列(阶段2)。

图1 渐进式SBAS-InSAR流程Fig.1 Flow chart of progressive SBAS-InSAR

首先，采用SBAS-InSAR对已有SAR影像数据集进行处理，获取对应地表形变时间序列结果；在此基础上，将已有SAR影像数据集中较新的部分(tn1时刻以后，共包含n景)SAR影像与新增SAR影像数据组合进行差分干涉处理，若研究区内新增N1景SAR影像(对应时间段为tN+1-tN+N1)，在相同时空基线阈值下形成M1个差分干涉对，以去除大气延迟、轨道误差后的解缠相位作为观测值，则可构建如下观测方程：

(1)

(2)

式中，A、B1、B2、C分别为系数矩阵，由0或相邻两景SAR影像的时间间隔构成(如12、24等)；X1、X2、X3为未知参数，由相邻两景SAR数据间的形变相位速率构成。

若将已有SAR数据和新增SAR数据的解算认为是2次平差过程，则从形式上看，式(1)和式(2)与测量平差中的序贯平差方法一致，故可将SBAS-InSAR求解已有SAR影像集的形变相位速率视为序贯平差的阶段1，将已有SAR影像集中较新的部分SAR影像集与新增SAR影像的形变相位速率视为序贯平差的阶段2，即可利用序贯平差思想更新所有SAR影像获取时刻的形变相位速率。

(3)

式中，δX1、δX2分别为X1、X2的改正数。

假设各观测时刻获取的SAR影像相互独立，且权矩阵为单位矩阵，利用序贯平差思想可得[6]：

(4)

式中，Q为SBAS-InSAR求解中存储的协方差矩阵。

综合式(3)和式(4)即可得到新增N1景SAR影像后各相邻SAR影像获取时刻的形变相位速率，在此基础上，通过数值积分可得到更新后的形变时间序列。

2 实验验证

采用模拟数据和真实数据，从解算精度和计算效率2个方面验证渐进式SBAS-InSAR的可行性与有效性。

2.1 模拟数据实验

以1个像元为例模拟3种形变模型，即线性模型、周期模型和指数模型。模拟时选用覆盖北京平原地区的26景SAR影像参数，假设前16景(2019-06-03～12-24)SAR影像为存档SAR影像，第17～26景(2020-01-05～04-22)为新增SAR影像，设置时空基线阈值为60 d与120 m，共形成112个干涉对，其时空基线组成如图2所示。此外，对3种形变模型分别添加均值为0 mm、方差为1.5 mm的高斯噪声，并进行100次蒙特卡洛随机模拟。分别利用SBAS-InSAR与渐进式SBAS-InSAR求解参数，得到地表形变时间序列，图3为添加高斯噪声与未添加噪声情况下3种形变模型对应的模拟真值、SBAS-InSAR和渐进式SBAS-InSAR获取的形变时间序列。结果显示，2种技术获取的未知参数的标准差小于1 mm，表明渐进式SBAS-InSAR与SBAS-InSAR获得的参数精度相当。

图2 已有SAR影像和已有SAR影像中较新部分SAR影像与新增SAR影像的时空基线分布Fig.2 Spatiotemporal baseline distribution of existing SAR images and newer SAR images in existing SAR images and new SAR images

图3 3种形变模型下SBAS-InSAR和渐进式SBAS-InSAR的形变时间序列Fig.3 Deformation time series derived from SBAS-InSAR and progressive SBAS-InSAR with three different deformation models

2.2 真实数据实验

2.2.1 研究区和数据介绍

选取北京通州区、朝阳区等地面沉降严重的地区进行研究，覆盖范围约38 km×42 km(图4)。研究区位于第四系凹陷区，受地下水水位降低及地下水不断超采的影响存在明显的地表形变[7]。利用SBAS-InSAR和渐进式SBAS-InSAR处理与分析2019-06～2020-04共26景C波段的升轨Sentinel-1A SAR影像。为去除地形和轨道误差对干涉相位的影响，同时使用欧空局提供的30 m分辨率SRTM1-DEM高程数据和精密轨道数据(定位精度优于5 cm)进行校正。

图4 研究区位置Fig.4 Location of the study area

2.2.2 数据处理

假设已有SAR影像集包含16景(2019-06-03～12-24)SAR影像，设置时空基线阈值为60 d与120 m，得到54对干涉对(图2黑色实线)。影像配准和重采样后，结合外部DEM及精密轨道数据去除地形及轨道误差，对高程误差进行建模和估算以将其扣除，采用时空域的高低通滤波去除大气延迟效应，最后利用SBAS-InSAR求解出相应时刻的形变相位值。在此基础上，逐渐新增1～10景(2020-01-05～04-22)SAR影像，并将已有SAR影像中较新部分SAR影像纳入新增SAR影像中(图2红色虚线)，采用渐进式SBAS-InSAR进行处理，获取全部影像的形变时间序列。

为分析渐进式SBAS-InSAR的精度与时间效率，同时利用SBAS-InSAR对所有SAR影像进行解算，获取其地表形变时间序列和所需时间。

3 实验结果与分析

3.1 监测结果分析

3.1.1 形变速率结果分析

图5为新增10景SAR影像时SBAS-InSAR与渐进式SBAS-InSAR分别反演得到的2019-06～2020-04北京平原地区沿雷达视线向的地表形变速率场。可以看出，研究区形变空间分布不均匀，主要集中在朝阳区与通州区交界地带，包括金盏乡、楼梓庄乡、徐辛庄镇、黑庄户乡、台湖乡、梨园、平西府镇、燕丹乡和李桥镇等地，最大形变速率约为-110 mm/a。结合图5(a)和5(b)可知，无论是形变量级还是形变分布区域，2种技术反演得到的形变结果均具有较高的一致性。

图5 年平均形变速率Fig.5 Annual average deformation rate

为更直观地显示本文方法的可靠性，对渐进式SBAS-InSAR与SBAS-InSAR获取的年平均形变速率作差(图6)。可以看出，二者差值的平均值和标准差分别为0.01 mm/a和0.1 mm/a，且2种技术获取的形变速率之差的绝对值基本小于0.5 mm/a，最大不超过1.5 mm/a，表明2种技术获取的形变结果具有较高的一致性。

图6 渐进式SBAS-InSAR与SBAS-InSAR年平均形变速率对比Fig.6 Difference of annual average deformation rate between progressive SBAS-InSAR and SBAS-InSAR

3.1.2 形变时间序列结果分析

图7为利用渐进式SBAS-InSAR获取的形变时间序列。可以看出：1)2019-06～2020-04研究区主要形变位于朝阳区与通州区的交界处，最大累积形变量约为-106 mm。2)研究区西部及东北部表现为轻微的地面抬升现象。3)随着监测时间的累积，金盏乡、楼梓庄乡、徐辛庄镇、黑庄户乡、台湖乡和梨园形变区的累积形变量不断增加，形变范围亦不断扩大；平西府镇、燕丹乡和李桥镇形变区的累积形变量有所增加，形变范围基本不变。

图8为渐进式SBAS-InSAR与SBAS-InSAR获取的累积形变量之间的差异直方图。可以看出，形变差异绝对值基本小于1 mm，最大形变差异绝对值不超过3 mm，表明渐进式SBAS-InSAR获取的累积时间序列与SBAS-InSAR的监测结果具有较好的一致性。

图8 SBAS-InSAR与渐进式SBAS-InSAR时间序列累积形变量之差直方图Fig.8 Histograms of the difference of cumulative deformation time series between SBAS-InSAR and progressive SBAS-InSAR

此外，选取3个特征点(A点位于快速形变区、B点位于较快速形变区、C点位于地表稳定区，位置见图5)对2种技术获取的时间序列进行对比分析，结果如图9所示。可以看出，3个特征点的累积形变量差值均小于2 mm，可忽略不计，表明渐进式SBAS-InSAR与SBAS-InSAR获取的累积形变结果高度吻合。

图9 3个特征点的形变时间序列Fig.9 Deformation time series of three feature points

3.2 时间效率分析

为便于对比分析，采用相同的准则选取渐进式SBAS-InSAR与SBAS-InSAR的高相干点目标，故2种技术处理的高相干点数目相同。以任意一个高相干点为例，分析渐进式SBAS-InSAR与SBAS-InSAR求解所需时间。对于单个高相干点，在相同时空基线阈值下，当新增N1景SAR影像后，SBAS-InSAR需要对M+M1(M为已有数据形成的干涉对数量，M1为新增SAR影像和已有SAR影像中较新部分形成的干涉对数量)幅干涉对去除地形、纠正大气效应和轨道误差等，而渐进式SBAS-InSAR只需对M1幅干涉图进行处理。从理论上讲，当已有SAR影像很多而新增SAR影像不多时，M1≪M，因此渐进式SBAS-InSAR比SBAS-InSAR具有更高的效率。

求解形变相位值时，已知形变相位值的求解主要涉及线性方程组(式(2)和式(3))的最小二乘法计算，计算复杂度与待求解未知数的个数有关，未知数越多，计算越复杂。SBAS-InSAR需要求解N+N1景SAR影像对应的地表形变速率增量，而渐进式SBAS-InSAR则只针对新增的N1景SAR影像对应的形变速率增量进行解算，并更新已有N景SAR影像对应的地表形变速率增量，计算复杂度明显降低，计算所需存储空间更小，时效性更高。

图10为在16景SAR影像的基础上逐渐新增1～10景SAR影像分别进行渐进式SBAS-InSAR与SBAS-InSAR处理的时间成本。数据处理的计算机配置为Intel Core i7(2.60 GHz)，内存32 G。从图10可以看出，执行配准、大气和轨道误差校正等步骤时，SBAS-InSAR的处理时间约为渐进式SBAS-InSAR的2～5倍。在形变参数反演解算阶段，渐进式SBAS-InSAR一方面要引入和使用SBAS-InSAR对已有SAR影像求解的协方差矩阵，另外还要对已解算的结果进行修正与更新，故无法像前一阶段那样明显地缩短运行时间。但随着SAR影像逐渐增加，渐进式SBAS-InSAR的时间效率相较于SBAS-InSAR仍具有明显的改进，在增加第10景SAR影像时，其形变解算时间提高约66%。综上，在处理随时间动态增加的SAR影像时，渐进式SBAS-InSAR明显具有更高的时间效率。