时序InSAR大范围地表形变监测的关键问题分析

刘媛媛，晏 霞

(1. 东华理工大学 测绘工程学院, 江西 南昌 330013; 2. 长安大学 地质工程与测绘学院, 陕西 西安 710054)

城市工程建设以及社会经济的快速发展导致人口密集程度越来越高，由此引发的地面沉降波及区域逐年扩大，影响范围高达数万平方千米，并且具有成片分布的空间分布特征(何庆成等, 2006)。

大范围(数万平方千米)地表形变监测要求监测手段既满足给定的精度指标，又具有大范围覆盖和同步测量的能力。传统的地面测量技术(如GNSS、精密水准测量等)布设的监测网络仅能覆盖数百平方千米的范围，而且难以保证监测的实效性。合成孔径雷达干涉测量(Synthetic Aperture Radar Interferometry，InSAR)及其衍生技术的出现为解决这一难题提供了强有力的技术支撑。然而，随着地面沉降区域的不断扩大，受到SAR影像幅宽范围(如Envisat ASAR标准幅宽为100 km×100 km)的限制，单幅SAR影像无法满足整个区域或者整个盆地尺度的地表形变监测需求，给大范围的地质灾害形变监测带来不便。另外，尽管时序InSAR技术(如永久散射体技术、小基线集技术等)已广泛应用于地表形变监测中(朱煜峰等, 2017;董少春等, 2019; 王升等，2019;Liu et al., 2018b, 2019)，却只聚焦于小范围、局部地区的形变分析，同时进行数千乃至数万平方千米的形变监测研究尚不多见(李广宇等, 2018; 张永红等, 2016; Chaussard et al., 2014; Peng et al., 2020; Yang et al., 2014)。

要实现高精度的大范围地表形变监测，采用的SAR影像必须同时兼顾在方位向能够连续获取长条带SAR数据和在距离向覆盖范围广(葛大庆, 2013)。然而，随着数据处理范围的不断扩大，影响形变信号解算精度的干扰因素(如基线误差、大气附加相位、相位解缠误差等)也随之增加，如何有效降低这些不利因素的影响是时序InSAR大范围地表形变监测面临的主要问题之一。同一形变信号在不同入射角下的InSAR监测结果在平面位置上并不能严格的一一对应，而且仅能获取相对给定参考点的形变量，导致相邻轨道的形变解算结果间存在明显的参考基准偏差，如何实现相邻轨道InSAR形变监测结果的坐标系与参考基准的统一是时序InSAR大范围地表形变监测亟待解决的又一主要问题。

针对上述问题，笔者开展时序InSAR大范围地表形变监测的关键技术分析，并将此应用于山西省太原盆地，获取盆地尺度的高精度地表形变监测结果。

1 长条带SAR数据的时序InSAR监测

长条带SAR数据的时序InSAR监测即对连续扫描获取的SAR影像序列进行分析，以获取长条带模式下的大范围地表形变场。

1.1 长条带SAR数据的分块处理法

考虑到计算效率、误差传递等因素对大范围形变监测的影响，尤其是轨道误差、大气附加相位等长波长误差的影响，笔者采用分块处理方式。其基本思路是首先将长条带SAR影像在方位向与距离向分别进行分块，分割成若干个具有一定重叠区域的数据子块，然后在各子块内利用时序InSAR技术进行参数估算，最后对估算出的线性形变速率、形变时间序列等进行拼接处理从而获取大范围的地表形变场。

1.2 轨道趋势性误差改正

由于SAR影像获取时卫星定轨精度可能不够精确，导致干涉组合的空间基线估值产生偏差，因此在长条带时序InSAR形变监测中，轨道残余相位成为影响形变监测精度不可忽略的因素之一。当监测区域范围逐渐扩大时，轨道残余相位随之逐渐累积，最终在形变图中表现为一种趋势性的影响。鉴于此，笔者通过建立多项式模型拟合相位解缠的趋势面，以消除其对形变解算的影响。数学模型表示如下：

φorb=a0+a1x+a2y+a3x2+a4y2+a5xy+a6h

(1)

式中，φorb表示轨道误差相位，(x,y)表示参与多项式模型构建的像元行列坐标，h表示参与多项式拟合运算的像元高程值，a0,a1,a2,a3,a4,a5,a6为模型拟合系数(即待求参数)。

1.3 连续刈幅的参考基准统一

在各数据子块中利用时序InSAR技术估算的线性形变速率、形变时间序列等参数具有各自独立的参考基准，且各参考点的活动性通常未知。由于各数据子块具有相同的轨道编号且已进行精确配准，因此连续刈幅参考基准统一的实质就是补偿不同数据子块间的基准偏差(可看作一个系统偏差)，通过对重叠区域内高相干点目标的形变参数进行统计分析实现。

以线性形变速率为例，连续刈幅参考基准统一的基本思路是首先选定一个子块作为参考，提取出所有相干点目标的形变速率，对其余数据子块的形变监测结果进行空间内插，生成连续的地表形变场；然后根据参考子块上相干点目标的空间位置提取其余数据子块对应的形变监测结果；最后按照公式(2)计算重叠区域内同名高相干点线性形变速率的整体偏差，从而完成连续刈幅线性形变速率参考基准的统一。

(2)

2 相邻轨道InSAR监测成果的拼接

2.1 坐标系的统一

为实现相邻轨道长条带时序InSAR监测成果的有效拼接，首要工作是在空间上统一相邻轨道InSAR监测成果坐标系(空间参考基准)。基本思路是将不同轨道下的InSAR观测值均转换至地理坐标系，然后以某一轨道对应的地理坐标系作为参考基准(主轨道)，将其他地理坐标系下(辅轨道)的InSAR观测值通过坐标转换使其位于统一的坐标系。

2.2 参考基准的统一

与长条带时序InSAR监测中的参考基准统一类似，相邻轨道InSAR监测成果参考基准统一的核心问题同样在于如何有效提取相邻轨道重叠区域内高相干点目标的InSAR形变参数。在此基础上估算相邻轨道监测结果间的整体参考基准偏差。由于不同轨道SAR数据对应的入射角各异，即使是地面上的同一地物在不同轨道SAR影像上的空间位置也会略有差异。因此，根据主轨道上高相干点目标的空间位置，提取辅轨道上对应的高相干点目标的形变参数，利用式(2)估计参考基准整体偏差，进而统一相邻轨道InSAR监测成果参考基准。考虑到相邻轨道各高相干点目标与形变解算参考点的距离对误差传播的影响，故采用反距离加权平均法公式(3)估算重叠区域高相干点目标的形变值。

(3)

3 应用分析

为检验上述时序InSAR方法的可行性，本次选取山西省太原盆地作为实验区，进行时序InSAR大范围地表形变监测试验与分析。

3.1 研究区概况

太原盆地位于山西省中部，东临太行山，西傍吕梁山，覆盖面积约6 195 km2，涉及太原市、晋中市、孝义市、介休市、清徐县、交城县等。该盆地是目前我国地面沉降等地质灾害最为严重的地区之一。太原盆地作为一个典型的内陆新生代断陷盆地，其地貌格局完全受控于活动构造，主要受到交城断裂带和范村-洪山断裂的控制(刘媛媛, 2018)。

3.2 数据源及数据处理

3.2.1 SAR数据

本研究采用的SAR影像主要包括2008年8月至2010年9月的20景Envisat ASAR数据和2007年2月至2011年3月的55景ALOS PALSAR-1数据，数据覆盖范围如图1所示。其中，Envisat ASAR数据包含2个连续的图幅，覆盖范围约为150 km×100 km；ALOS PALSAR-1数据包含3个相邻的轨道(Track455, Track456, Track457)，各轨道包含3个连续的图幅，数据覆盖范围约为210 km×70 km，相邻轨道的重叠度约为30%，地面距离约为21 km。

图1 太原盆地的SAR数据覆盖范围Fig.1 Spatial coverage of the SAR datasets in Taiyuan basin

3.2.2 时序InSAR处理

为实现太原盆地高精度的地表形变监测，采用多级配准策略实现主、从影像间的精确配准与重采样，并采用小基线模式组合干涉对，Envisat ASAR和ALOS PALSAR-1数据的时间阈值分别设为250天和1 000天，空间基线阈值分别设为500 m和1 500 m。利用欧空局提供的DORIS精密轨道数据进行基线精化，并利用研究区30 m分辨率的SRTM DEM数据模拟与去除平地和地形相位的影响。

由于太原盆地内既存在城镇等高相干点目标密集的区域，又包含大片的农田等低相干区域，因此本次主要采用干涉点目标分析技术(刘媛媛, 2018; Werner et al., 2003)进行时序InSAR处理，获取各高相干点目标的形变时间序列和年平均形变速率图。

3.3 结果分析

3.3.1 相邻轨道InSAR监测成果的拼接

为获取太原盆地完整的地表形变场，需要对3个相邻轨道的ALOS PALSAR-1形变监测结果进行拼接。鉴于不同轨道SAR数据获取时刻存在差异，因此仅能对相邻轨道获取的年平均形变速率进行有效拼接。此外，由于相邻轨道SAR数据入射角略有不同，假定研究区地表形变仅发生在垂直方向，利用各轨道的入射角将视线向形变速率投影至垂直向，完成相邻轨道的成果拼接。在本实验中，以Track456作为参考轨道，首先将Track455和Track457的形变速率投影转换至该坐标系，然后统计相邻轨道重叠区域的同名高相干点目标，利用公式(3)进行参考基准偏差校正，得到具有统一参考基准的太原盆地尺度的年平均形变速率(图2)。从图2可以看出，相邻轨道InSAR形变监测结果在重叠区域具有较好的一致性，表明笔者提出的时序InSAR大范围形变监测方法具有一定的可行性。

图2 ALOS PALSAR-1数据的年平均形变速率(2007—2011年)Fig.2 Annual deformation rate derived from ALOS PALSAR-1 from 2007 to 2011

3.3.2 形变结果分析

图3为Envisat ASAR数据获取的2008年8月至2010年9月的太原盆地年平均形变速率。图2和图3均显示，2007—2011年期间整个盆地普遍存在不同程度的地面沉降灾害，沉降区具有明显的差异特征，主要分布在太原盆地的重点城市区域，如太原市中南部地区、清徐县、交城县、祁县、太谷县以及介休市等区域，其年平均沉降速率一般超过90 mm/a，最大沉降中心位于清徐县，沉降速率高达240 mm/a。笔者分析发现形变区与太原盆地地下水源地的空间分布具有一定的关联性，沉降区基本与地下水开采严重的区域对应。Liu等(2018a, 2018b)详细分析了盆地典型区域(太原市、清徐-交城区域)的地面沉降与地下水之间的相关性。此外，从空间分布来看，太原盆地的地面沉降受到断裂带的控制，断裂带两侧具有明显的差异沉降特征。

图3 Envisat ASAR数据的年平均形变速率(2008—2010年)Fig.3 Annual deformation rate derived from Envisat ASAR from 2008 to 2010

3.3.3 精度评定

本次首先对ALOS PALSAR-1相邻轨道重叠区域的形变观测值进行求差统计，检验独立观测值的精度(内符合精度评定)。图4为Track455与Track456、Track457与Track456重叠区内线性形变速率差值的统计直方图，其基本符合正态分布，对应的标准差分别为3.7 mm/a和4.2 mm/a，表明本次提出的时序InSAR大范围地表形变监测方法是可靠的。另外，笔者还对ALOS PALSAR-1与Envisat ASAR数据重叠区域的形变速率进行了求差统计(图5)，对应的差值标准差为5.6 mm/a。经分析发现，差值较大的区域恰好与形变速率较大的区域对应(例如清徐-交城区域)，可能是受水平位移影响所致。

图4 ALOS PALSAR-1相邻轨道形变速率差值统计直方图Fig.4 Histograms of deformation rate difference between the results from adjacent track ALOS PALSAR-1

图5 ALOS PALSAR-1与Envisat ASAR形变速率差值(a)及对应的统计直方图(b)Fig.5 Deformation rate difference between the results from ALOS PALSAR-1 and Envisat ASAR (a) and the corresponding histogram (b)

此外，笔者还利用中国地震局第二监测中心提供的16个GPS监测点位(图3所示)观测值评价InSAR形变监测结果的可靠性(外符合精度评定)。考虑到InSAR监测的是各高相干点目标的形变值，而GPS技术反映的仅是特定监测点的变形，二者在空间位置上可能存在不匹配现象。因此，本次以GPS监测点为中心提取100 m范围内的高相干点目标，采用反距离加权平均方法获取各GPS监测点对应位置的InSAR形变值。InSAR与GPS间的相关性分析(图6)表明，InSAR形变监测结果与GPS监测数据间具有很高的一致性，相关系数分别为0.992 8(ALOS PALSAR-1)和0.988 2(Envisat ASAR)，与GPS间差值的均方根误差(RMSE)分别为2.8 mm/a(ALOS PALSAR-1)和2.5 mm/a(Envisat ASAR)，表明本次监测结果具有较高的精度。

图6 InSAR与GPS之间相关性分析Fig.6 Correlation analysis between InSAR result and GPS observations

4 结语

笔者针对时序InSAR大范围地表形变监测的关键问题开展了有针对性的研究，包括长条带SAR数据分块处理法以及相邻轨道InSAR监测成果的拼接，并以山西省太原盆地为研究对象开展了应用试验分析，获取了山西省太原盆地尺度的高精度地表形变监测结果。结果显示，ALOS PALSAR-1和Envisat ASAR数据与GPS监测数据的相关系数分别为0.992 8和0.988 2，与GPS间差值的均方根误差分别为2.8 mm/a和2.5 mm/a，表明本次提出的时序InSAR技术有效地解决了大范围地表形变监测同步观测的难题，实现了地表形变大范围、高精度的监测，有效地揭示了区域性地表形变的空间分布特征，展现了常规地面监测手段无法比拟的技术优势，对于宏观揭示盆地尺度的地表形变的空间变化特征具有重要意义。