山西省临汾市矿区地表形变InSAR大范围探测与监测

麻学飞， 张双成， 惠文华， 许 强

(1.长安大学地质工程与测绘学院，西安 710054； 2.地理信息工程国家重点实验室，西安 710054)

0 引言

煤炭作为我国的支柱能源，关系到国家经济的命脉，需要不断地进行开采，这势必会形成采空区。采空区上部岩层发生断裂、弯曲等形变，会造成采空区沉陷，破坏当地生态环境，威胁居民的人身财产安全[1-2]。因此对矿区进行大范围探测和动态监测具有重要的研究意义。水准测量、全球定位系统(global positioning system，GPS)等方式是传统的矿区地表形变监测手段，但由于这些方法监测成本高、工作量大等原因，难以进行长期监测，且所获监测结果时空分辨率低，难以反映真实的矿区形变过程[3-5]。

合成孔径雷达干涉测量(interferometry synthetic aperture Radar，InSAR)作为一种高精度的测量手段，利用相位差获取研究区地表形变信息，具有监测范围大、效率高、分辨率高的特点，几乎不受云雾的影响，成为了矿区地表监测的有效手段[6-7]。 差分合成孔径雷达干涉测量方法(differential interferometric synthetic aperture Rader,D-InSAR)可获取到2景合成孔径雷达(synthetic aperture Radar,SAR)影像覆盖期间形变体表面的微小形变信息[7]，从而实现对形变体进行监测，处理数据量小，可以快速获取目标区域的形变信息。Ge等[8]利用D-InSAR技术对Appin，Westcliff和Tower煤矿区地表形变进行监测，所获监测结果精度可达厘米级； 龙四春等[9]利用D-InSAR与GPS集成实现了对资兴唐煤矿区三维形变监测，并对测量结果进行集成与内插，验证了2种监测结果的一致性。但是D-InSAR技术难以获取目标区域时序变化过程，且在实际应用过程中常常受到时空去相干、大气延迟等影响[10-11]。而时序InSAR弥补了D-InSAR技术的不足，在一定程度上解决了时空失相干、数字高程模型(digital elevation model,DEM)误差、大气延迟造成的影响，监测精度可达厘米级甚至毫米级[12-13]。永久散射体测量技术(permanent scatters InSAR，PS-InSAR)利用影像上存在的永久散射点来替代整幅影像上的像素点，通过永久散射点的形变来呈现目标区域的形变[14]，但是只适合小形变速率、小范围的形变，无法很好地适用于沉陷区大范围探测与监测； 小基线集(small baesline subset,SBAS)-InSAR能以更高的时空密度来提取到相干点目标，其形变结果更加明显、直观[15-16]，是目前矿区沉陷监测最主要的方法之一。沙永莲等[17]使用SBAS-InSAR技术对新疆哈密砂墩子煤田矿区进行监测,发现了一个沉降漏斗，时序结果显示沉降漏斗先线性下沉后逐渐趋于稳定； 李达等[18]通过对13景Terra SAR-X数据进行SBAS-InSAR技术处理,监测到了陕西榆林某煤矿的地表形变，并提取多个观测点的时序沉降值进行量化分析，结果表明 SBAS-InSAR技术可为矿区地表形变监测和分析提供新的监测手段。但是时序InSAR技术对影像的数量有一定要求，处理数据量较大，周期长，无法高效率地对大范围的形变区域进行监测。针对矿区地表形变大范围探测与监测的问题，可以将D-InSAR技术和SBAS-InSAR技术结合起来使用： 首先利用D-InSAR技术对目标区域进行大范围探测得到沉陷区位置，再利用时序InSAR技术对目标沉陷区进行时序监测和形变追踪，高效率地完成对目标区域的形变监测。

本文获取了覆盖山西省临汾市全域不同轨道的C波段Sentinel-1A数据，采用联合D-InSAR技术和SBAS-InSAR技术开展临汾市沉陷区的大范围探测和监测工作，总计探测出沉陷区105处，选取其中的2处矿区进行详细分析，并结合光学影像进行验证。

1 研究方法

1.1 D-InSAR基本原理

D-InSAR技术通过采用目标区域形变前后两景影像共轭相乘获得干涉图，与DEM模拟相位差分消除地形相位后获取到地表形变信息。但是在D-InSAR用于具体的地表形变监测过程中，由于大气环境和地表形变处于不断的变化的过程，干涉相位会不可避免地受到影响，进而影响到最终的测量结果[19]。差分相位表达式为：

φ=φflat+φtopo+φdef+φatm+φnoi，

(1)

式中：φ为形变相位；φflat为平地相位，可借助精密轨道数据去除；φtopo为地形相位，可借助DEM数据模拟地形相位去除；φdef为地表形变引起的相位；φatm为大气延迟相位;φnoi为随机噪声相位，可采用滤波的方式进行消除。

1.2 SBAS-InSAR原理

小基线技术利用同一轨道重复观测获取到目标区域的SAR影像集，通过设定时空基线将SAR影像集进行分组，形成多个小基线子集，再分别对每个子集内的SAR影像进行差分干涉生成差分干涉图和相干性图，利用平均相干性从相干性图中选取高相干性的像素点，对差分干涉图依次进行相位解缠得到相位序列，联合多个小基线子集的相位序列来预估形变信息。若在一定时间内有N+1景 SAR 影像覆盖目标区域，设定时空基线后可以组合成M幅干涉图，M需要满足[20]:

(2)

以t0时刻作为初始时刻，则tk时刻的差分干涉相位φtk为观察量，其中1

(3)

式中：d(ta,x)和d(tb,x)为ta和tb时刻雷达视线方向上的形变值;φ(ta,x)，φ(tb,x)分别为d(ta,x),d(tb,x)对应得解缠相位值；λ为波长；Ti为主辅影像时间间隔；i为N幅影像获取时间段内某一时刻；v(x)为像元x处的形变速率。N幅影像可以组成M个方程，其方程组的形式为：

Aφ=δφ，

(4)

式中：A为M×N的系数矩阵，每一行代表一副干涉影像，每一列代表对应时间上的SAR影像。若M>N，且A的秩为N，依据最小二乘法估值为：

(5)

式中ATA为奇异矩阵，存在无数解，可用联合多个小基线采用奇异值分解法求出上式最小范数下的二乘解。

D-InSAR技术和SBAS-InSAR技术2种方法使用的场景不同，根据需求联合2种监测技术实现对于临汾地区矿区的大范围探测和动态监测，技术流程如图1所示。

图1 技术流程图

2 研究区概况及数据源

山西省临汾市位于N35°23′～36°57′,E110°22′～112°34′之间，处于太原、郑州、西安3个省会城市连线的中心点，地理区位优势明显，矿产资源丰富； 地处属温带大陆性气候，四季分明，雨热同期[21-22]； 地表多为黄土，土质疏松，夏季多暴雨，地表植被容易遭到破坏。临汾市包含河东煤田北部、霍西煤田中部和沁水煤田西部，煤炭资源丰富，2019年原煤产量6 201.6万t，占全省原煤产量的6.39%。临汾市矿区多分布于临汾断陷盆地山体之中，地表以植被为主，多年的煤炭开采造成地表沉陷严重，容易引发滑坡、地面塌陷等多种地质灾害，迫切需要有效的沉陷区监测手段。

实验选取432景不同轨道观测获取的Sentinel-1A影像(图2)，其中path=11，frame=111有103景，path=11，frame=116有103景，时间跨度为2017年3月12日—2020年11月27日；path=113，frame=111有113景，path=113，frame=116有113景，时间跨度为2017年3月19日—2020年12月28日。影像均为C波段，周期为12 d，模式为升轨，VV极化，距离向分辨率为5 m，方位向分辨率为20 m。DEM数据为AW3D30 DEM，分辨率优于30 m。卫星轨道参数为欧空局提供的POD(AUX_POEORB，POD Precise Orbit Ephemerides)精确轨道参数。干涉雷达处理软件为瑞士GAMMA软件。

图2 影像和DEM覆盖范围Fig.2 Image and DEM coverage

3 数据处理结果与分析

3.1 矿区大范围探测

矿区形变具有连续性，受非人为因素的干扰较小，只要在矿区连续沉陷过程的某一时间段内的形变过程被监测到，则认为该区域为沉陷区。本文选取2017年末、2019年初和2019年末3个时间段共12景不同轨道的Sentinel-1A影像组成覆盖临汾市全域的6组干涉对共3个组合(组合a、b、c)进行大范围探测。较短的时空基线可以在一定程度上减小时空失相干带来的影响，为了保证干涉对具有较高的相干性，提高观测精度，选取影像时时间基线设置为12 d，空间基线最大值控制在±100 m以内。

在数据处理过程中，将距离向和方位向的多视比设定为4∶1，地形相位利用外部DEM去除，干涉对中存在的噪声采用自适应滤波的方法去除，相位使用奇异值分解法进行解缠，大气相位应用GACOS技术去除，最终获取到高质量的解缠图。由于矿区在短时间内形变量级较大，形变相位在解缠图中特征明显，可直接依靠相位特征获取到沉陷区位置，将沉陷区大致范围用矩形框标识出来。由于沉陷区在地理位置上呈现集群式分布，为了后续的数据方便处理将识别结果划分成6部分，依次编号为A—F，结果如图3所示。

图3 D-InSAR解缠图识别结果

图3(a)—(c)中黄色虚线为获取到的沉降区的大致位置。图3(d)中红色矩形框圈画出来的即为沉陷区，形变相位表现为由青-紫-青变化； 黑色的空洞是因失相干而产生的监测结果缺失。造成失相干的原因主要有2种： ①因沉陷区形变量过大，超出了C波段的监测能力而产生失相干； ②因植被茂密导致地物介电质常数发生变化而导致的失相干。沉陷区的也存在2种状态： ①形变量级较大的沉陷区在解缠结果上表现为中心呈空洞状态，周围有形变相位环绕； ②形变量级较小的沉陷区表现为紫色或者黄色的形变相位环绕。沉陷区主要分布于临汾断陷盆地两侧的山体中，沿盆地东西两侧山体走向分布，呈点状集群分布。分布于盆地中心的附近的沉陷区，地表植被较少，具有较高的相干性，获得结果完整，易于识别沉陷区； 盆地两侧山体因植被影响造成大面积的失相干，给沉陷区识别带来了一定的困难。C区域的沉陷区分布区域面积最大，数量最多，后文将着重以C区域作为重点进行表述。

3.2 矿区时序监测

基于D-InSAR识别到的沉陷区位置，依次对6个区域进行SBAS-InSAR时序处理，数据处理详细参数如表1所示。2个轨道公共主影像时间分别为2019年1月1日和2019年1月8日，将主影像作为统一的基准进行配准并作为时间参考点，设定时间基线不少于90 d，空间基线不超过120 m。 为了提高结果的精度，使用强度互相关算法对主辅影像进行配准，通过设置不同的窗口大小、多项式参数个数等参数提高影像配准精度。通过外部DEM去除地形相位，利用时空滤波方法去除大气误差，选取高相干性点位作为基准进行相位解缠，最终获取到6个区域垂向的形变速率，结果如图4所示。从图4中可以看出，识别所得开采沉陷区范围均在D-InSAR监测结果中，其中灰色线条为手工勾绘的沉陷区边界。本文此次探测出105处沉陷区，总面积达188.75 km2，其中7处与行政界线相交。识别所得沉陷区数量多，单个沉陷区沉陷面积较小，部分沉陷区形变量级较大，多处开采形变速率超过-100 mm/a，可以反映出临汾市开采沉陷情况非常严重，这足以给当地矿区生态环境造成巨大的破坏。

表1 数据处理详细参数Tab.1 Detailed data processing parameters

图4 临汾市2017年3月—2020年12月SBAS-InSAR沉陷区形变速率图Fig.4 Deformation rate map of SBAS-InSARsubsidence area in Linfen City fromMarch 2017 to December 2020

选取形变量级最大的C区域进行详细分析，获取到的累计形变量和形变速率结果如图5所示。C区域存在43处沉陷区，2017年3月19日—2020年12月28日该区域最大累计形变量为1 030 mm，最大形变速率为-381 mm/a。区域右侧植被覆盖较少，影像相干性较高，可以获取到完整沉陷区形变，而左侧区域为山区，植被覆盖率高，白色区域即为影像因植被影响产生失相干所造成的结果缺失。其中C1沉陷区由2个相邻的矿区构成，左侧沉陷区形变量级大于右侧，最大累计形变量分别为922 mm和593 mm。C2沉陷区为单沉陷区，形变轮廓呈现为条带状，沉陷区的轮廓反映出矿区由西南向东北方向不断开采，最大累计形变量为942 mm。

从图6可以看出，C1，C2两个沉陷区均处于不断开采过程中。从形变过程来看，4个点位在矿区开采过程中形变量级在不断增大，形变过程呈现快速沉陷-缓慢沉陷-快速沉陷的周期，以p1点的时序变化表现的最为明显。在2017年5—8月，由于没有获取到有效的干涉对，故时序呈现拟合的直线。p1点整体形变速率最大，形变量总是经历较长时间快速形变到较短时间缓慢形变的周期循环，形变速率的绝对值也是呈大-小-大的周期变化，但是整体的形变趋势逐渐缓和，最大沉降量为922 mm。p2点整体形变趋势较为平稳，周期性形变特征不太明显，最大沉降量为593 mm。p3点位于C2矿区沉陷中心，p3点的整体形变速率要大于p4点，呈现快速下降-稳定-略微抬升-缓慢下降-快速下降的过程，最大沉降量为942 mm。p4点位于矿区开采沉陷方向的中间位置，与p3点为同一矿区，整体形变趋势与p3相似，最大沉降量为404 mm。

将C1，C2矿区分别沿图6中ab，cd剖线截取剖面，结果如图7所示。C1矿区有2个沉降漏斗，左侧的沉降量要大于右侧，2个沉陷区中间的位置已经有100 mm的沉降量，已经成为一个双中心的大型沉陷区。C2矿区沉陷左侧形变倾角大于右侧，说明矿区开采的方向是沿剖线cd方向进行开采，沉降漏斗中心位于起采区，但是随着开采工作的进行，沉降漏斗中心会逐渐向开采的方向移动，沉陷区轮廓会呈现为条带状。

由于没有矿区的实测数据，本文利用Google Earth遥感影像数据检验监测结果的一致性。首先采用阈值分割法将C1，C2沉陷区边界提取出来，经过多次尝试当阈值设置为-10 mm/a时，所获得沉陷区边界效果最好，将矿区边界叠加至光学影像上如图8所示。结合光学影像进行分析，2个矿区均为地下开采，沉陷范围内有建筑物存在，地表分布有大量农田，沉陷面积分别为1.54 km2，2.87 km2。C1矿区南侧地势平坦，已被开垦为农田，北侧为山地，人工地物较少，中间由一条道路横穿沉陷区，左侧100 m附近有一村落，周围发现2处疑似采煤作业区。C2矿区地表起伏较大，但地表已被全部开垦为农田，中间由多条道路穿过，东南侧与小镇相接，东北方向50 m有一条国道，北侧分布有3个疑似采煤作业区。C1,C2矿区沉陷边界直接与农田、道路相接，距离住宅不足100 m，产生的地面沉陷会直接威胁到周围居民的人身财产安全，尤其C2矿区，需要持续进行监测。

4 结语

煤矿开采所造成的沉陷区对地表生态环境和基础设施带来了巨大的破坏。本文利用D-InSAR技术对临汾市沉陷区进行大范围探测和监测结果工作，将探测结果分成了6个区域并依次使用SBAS-InSAR技术进行时序监测，共发现105处沉陷区； 将其中2处典型沉陷区作为重点区域进一步分析，获取到了该沉陷区的形变演化情况，为后续矿区沉陷灾害防治提供重要依据。通过光学影像发现了沉陷区附近存在疑似采煤作业区，验证了InSAR技术大范围探测与监测具有可行性和可靠性。

1)矿区地表形变具有形变速率快、量级大的特点，会对当地的生态环境造成巨大的破坏，难以通过传统的光学遥感进行探测，InSAR技术可以极易捕获沉陷区的大量级地表形变，从而实现对研究区沉陷区的大范围探测。

2)直接利用时序InSAR技术对大范围区域进行沉陷区探测具有巨大的挑战性，处理数据量大，工作效率低，结果不稳定。而基于D-InSAR技术的大范围探测和SBAS-InSAR技术的详细监测的方法可以很快实现对研究区进行层次化的处理，真正做到对大范围区域的高效率监测。

3)本文基于InSAR技术使用Sentinel-1A升轨数据对临汾市矿区地表形变进行了大范围探测和监测，识别和监测结果对临汾市矿区形变监测具有较好的参考价值，但是InSAR技术在植被覆盖率较高的区域适用性不佳，未来还需要借助其他技术来提高大范围探测和监测结果效率。

志谢：感谢欧洲太空局在科学数据中心网站公布Sentinel-1A卫星雷达影像，感谢日本宇宙航空航天局地球观测研究中心提供DEM数据。