基于SBAS技术的九寨沟震后地表形变监测

胡 波，吴 洋，魏德宏，杨 斌，余俊鹏，陈志谋

(1.广东工业大学，广东 广州 510006；2.伟志股份公司，福建 晋江 362200)

2017-08-08T21：19于四川九寨沟(103.82°E，33.20°N)发生Ms7.0地震，震源深度约为20 m。地震发生后，九寨沟景区遭受到较大的破坏，本次地震造成了25人死亡，520余人受伤，同时7万余房屋受到不同程度的损坏[1-2]。而地震之后一段时间内当地地形的变化成为灾后重建的一项重要监测内容。

合成孔径雷达干涉测量(Interferometry Synthetic Aperture Radar，InSAR)是一门测量地形及地表形变的技术。经过近30年的发展，InSAR技术已经成为大面积区域地面沉降监测的重要手段[3]。相比较常规测量技术而言，InSAR技术具有全天候、宽覆盖、高精度、低成本的优势[4]。但是，传统的差分合成孔径雷达干涉测量(Differential Interferometry SAR，D-InSAR)技术容易受到大气扰动、时间失相关和空间失相关等因素的影响，其测量的精度受到一定限制，并且无法获得监测区域地表形变的时间演化情况[5]。针对D-InSAR中存在的诸多问题，科研工作者相继提出了时序InSAR技术。其中，具有代表性的时序InSAR技术有永久散射体干涉测量[6](Permanent Scatter Interferometry, PS)技术和小基线集[7-8](Small Baseline Subset, SBAS)技术。目前时序InSAR技术已广泛应用于城市地下水开采[9]、地质灾害[10-11](滑坡、火山、地震等)和城市地表[12-13]等形变测量当中。

本次实验选取25景Sentinel卫星数据，使用小基线集合成孔径雷达干涉测量技术对九寨沟地震区及其附近山区进行测量，得到该地区震后近1年内的时间序列形变，为探究该地区震后的沉降规律、灾后重建及地质灾害预防提供技术和数据支持。

1 研究区域及实验数据

1.1 研究区域

本次研究的区域位于四川省阿坝州松潘县北部，九寨沟地震震中附近。中心坐标为103.80°E，33.12°N，其覆盖面积约为6 712 km2。研究区域包含九寨沟景区，区域内地形以山脉为主，如图1所示。

图1 研究区域覆盖范围图

1.2 实验数据

本次研究选用了25景Sentinel-1A的升轨数据(见表1)，采用C波段观测。为提高Sentinel-1A影像数据处理精度，本研究使用卫星的精密轨道数据(Precise Orbit Ephemerides,POD)。实验中所使用的DEM数据是由日本JAXA公司提供的ALOS数字表面模型“ALOS World 3D- 30 m”。干涉影像对的时间与空间基线关系如图2所示。

表1 Sentinel-1A影像信息

图2 干涉影像对时空基线图

2 SBAS技术原理

Berardina和Lanari等人于2002年首次提出了SBAS技术[7]，该技术以多幅影像作为主影像生成干涉对，并利用干涉图中高相干点恢复研究区域的时间序列形变信息。同时，SBAS技术利用多视处理降低差分干涉图的相位噪声，并应用奇异值分解(Singular Value Decomposition, SVD)方法对相位进行求解得到形变相位速率，进而解得整个观测时段的形变时间序列[14]。

SBAS能有效提高研究区域形变的时间和空间分辨率，其步骤如下：

1)针对传感器、观测条件、研究区域的不同设定空间基线阈值和时间基线阈值以及平均相干性阈值，可以排除空间和时间失相干的像对，同时生成M对小基线干涉对，干涉对数量应满足：

(1)

式中：N+1为影像个数。在地形起伏较大的区域，需借助外部DEM数据去除地形起伏导致的地形相位影响，从而得到n幅差分干涉图的形变相位；

2)对于第i景差分干涉图中，方位向坐标A以及距离向坐标R的像元的干涉相位值可以表示为：

δφi(A,R)≈[d(TB,A,R)-d(TA,A,R)]+

(2)

3)利用得到的干涉图进行二次筛选，将干涉图中干涉条件差或是解缠结果中存在大面积整周跳变结果的干涉图剔除；

4)采用选取多个相对稳定的地面点的方法求取干涉图里的大气相位以及噪声项成分。在求取误差项以后，使用高次多项式插值的方法，对研究区域中的大气相位进行插值，达到恢复研究区域中每一景干涉图所包含的大气相位及噪声的目的。

本次实验采用三次多项式插值，为了去除DEM误差相位和大气延迟带来的影响，可假设地表形变的低频部分为：

(3)

(4)

Ax=Δφ.

(5)

5)在去掉残余地形相位、大气相位和相干噪声后，可以将式(2)简化为：

δφi(A,R)=d(TB,A,R)-d(TA,A,R).

(6)

用矩阵形式表达式(6)，可写成：

Bφ=δφ.

(7)

式中含有M个观测向量和N个未知数，计算结果时，可以使用最小二乘法解得结果，然而在差分干涉中，为了抑制时空基线过长带来的去相干，常出现M幅干涉图并不处于同一子集的情况，该情况下用最小二乘得到的结果并不唯一，此时使用SVD方法可以得到未知参数的最小范数解。

3 平均形变速率及其时间序列形变反演

图3 九寨沟震区2017-08-011—2018-07-01平均速率图

本研究利用SBAS技术得到的该地区整体垂直方向平均速率如图3所示。图3显示该区域在2017-08-11—2018-07-01期间整体形变量较小，大部分区域形变速率保持在-50～50 mm/a。同时，研究区域的东北部地区下沉速率较大，其年形变速率达到-70～-50 mm/a。除此之外，在中部地区和南部地区出现3块面积较小，但下沉速率较快的区域，它们的面积小于1 km2，但最快下沉速率已经超过120 mm/a。为了进一步探究研究区域整体形变随时间的变化，选取3块区域共11个采样点进行分析(见图3)，采样区域信息及采样点分布如表2所示。

表2 九寨沟采样区域信息及采样点分布信息

图3中区域1占据了研究区域的大部分面积，该区域主要表现为较为稳定的形变，其累积形变量多在-50～50 mm之间。为探究该区域地表形变规律，本文提取4个采样点，其中P1位于震中(103.82°E，33.20°N)附近，P2为该区域中下沉明显的点，而P3和P4分别位于区域的西北部和西南部。从P1～P4的时间序列形变图(见图4)中可以看出，P3和P4的时序形变较为稳定，波动幅度较小。而P1点由于位于九寨沟地震震中附近，在地震之后仍伴有较大的起伏，其最大抬升量和最大下沉量的差值达到37 mm。而P2虽然表现为较为平稳的下沉，但其累积下沉量仍超过80 mm，达到82 mm。

图4 区域1采样点时间序列形变图

区域2位于研究区东北部，九寨沟县北部。该区域表现为整体下沉，且形变速率多在-40～-70 mm/a之间。在该区域中，本文选取3个点(P5～P7)作为采样点进行时序分析。从图5中可以看出，2017-08-11—2018-01-02，该地区处于快速下沉状态，其累积下沉量分别达到41 mm、49 mm和58 mm。而2018-01-02后，该区域开始趋于稳定，其形变量波动小于±5 mm。

图5 区域2采样点时间序列形变图

图6 区域3采样点位置图

区域3位于研究区域南部(见图6)，其最大下沉区域面积均不超过1 km2，但其形变速率较大，且累计形变超过-100 mm，本研究同样在区域3和4中分别取两个采样点进行时序分析，时序形变图如图7所示，从图7中可以看出，P8～P11始终处于下沉状态中，且下沉量不断累积。其中P11下沉最为明显，其最大下沉量达到118 mm。由于该区域中两个下沉明显的部分沉降速率较快且保持持续下沉，故可作为重点对象进行监测和后续的考察。

图7 区域3采样点时间序列形变图

4 结论

本次研究基于欧空局Sentinel-1A卫星影像，运用SBAS技术对九寨沟地震区及其周边地区进行了时间序列形变的监测，基本查明震后灾区附近的地表形变趋势，研究结果表明：

1)研究区域中大部分区域形变速率均维持在-50～50 mm/a。研究区的东北部表现为大面积下沉，且与中部区域相比下沉幅度更大。另外研究区域中存在3处面积不足1 km2的地区(103°41′46″E，33°13′8″N；103°52′E，32°54′N；103°59′E，32°51′N)出现了速率为-80～130 mm/a的大幅度下沉，在2017-08-11—2018-07-01期间最大累积下沉量也达到118 mm。

2)SBAS技术能有效监测大范围的时序形变，其监测精度可以达到mm级。另外，InSAR技术在监测大范围形变的同时，还可以监测到小面积、大幅度的地表形变，在地表持续监测上可以减少人力物力财力的投入，实现对形变异常区域更精准的定位。

除此之外，可以对当地历史影像进行数据处理，还能进行后续形变监测，从而得到更全面的形变规律，为山区形变监测和地质灾害预警及防治提供更强大的理论基础和数据支持。