多时相SAR振幅信息的大量级沉陷区识别研究

褚夫玉，胡海峰，廉旭刚，朱利超

(1.太原理工大学矿业工程学院，山西 太原 030024;2.山东省水利勘测设计院，山东 济南 250000)

煤层开采所引起的地表沉陷是采空区周围应力场平衡-破坏-再平衡的力学现象，是我国主要的地质灾害之一，长期影响着人民生命财产安全和生态环境[1-2]。山西省作为煤矿采出大省，地质灾害频发，开采沉陷区分布广泛，农田、建筑物损毁、植被退化、地质环境恶化问题时有发生，矿区可持续发展受到严重威胁[3]。因此，利用先进的监测手段，研究开采沉陷区的运动特征及规律，识别并监测沉陷区破坏程度，对促进开采沉陷区损害的早期预警、矿区地质灾害治理和可持续发展有重要意义。

开采沉陷是煤矿开采工程中常见的地质灾害类型。传统监测矿区地表形变的方法有水准测量、全球定位系统(global positioning system，GPS)等，测点精度高，但在地表沟壑纵横、地形起伏较大的地区点位布设困难且不能维持监测点的完整性，难以有效开展大范围、长时间的监测活动[4]。合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)是监测和评估地质灾害的新技术，该技术无需地面控制点，可进行全天候、高精度、大范围监测，目前已应用于我国煤矿开采沉陷的监测研究中[5]。CARNEC等[6]在1996年首次将差分干涉测量技术(differential interferometry synthetic aperture radar，DInSAR)应用到法国Gardanne煤矿开采沉陷研究中，获取了单个重访周期内的最大下沉量，但受大气延迟、地形误差和时空失相干的影响，监测结果的有效性和可靠性较低；FERRETTI等[7]在1999年针对常规DInSAR受失相干和大气影响的缺陷，率先提出永久散射体干涉测量技术(persistent scatterer interferometric synthetic aperture radar，PS-InSAR)，并利用ERS-1/2数据对意大利Camaiore镇地表沉降进行监测，获得了较为可靠的形变信息；WALTER等[8]在2004年利用PS技术对德国Prosper-Haniel煤矿进行开采沉陷监测，证明了该技术在矿区沉陷监测的可行性；2005年，CASU等[9]利用短基线集方法(small baseline subset，SBAS)进行地表沉降监测，获得了与水准及GPS测量一致的结果，但忽略了InSAR形变数据与实测离散点时空尺度上的差异。

目前，国内学者先后利用DInSAR、PS-InSAR、SBAS等多种InSAR技术在我国煤矿开采沉陷区开展了一系列研究[10-16]。其中，PS-InSAR数据量要求较大(一般为20景以上)、PS点的分布呈随机性且受时空失相干的影响较大，在植被覆盖密集、地形变化复杂的矿区，难以找到足够的相位稳定点，导致观测结果不能客观反映沉陷盆地整体变化[17]；SBAS技术通过选取短时空基线的干涉对，减小了时空失相干的影响，提高了形变监测点的密度，但在地形变化较大的区域目标的信噪比较低，进而影响形变监测结果的可靠性[18]。DInSAR、PS-InSAR、SBAS等是基于电磁波的相位信息提取微小形变量的技术，在形变剧烈的矿区，相邻两像元间的形变量大于雷达波长的1/4时，会产生相位混叠，难以解缠到大量级形变。基于SAR影像振幅信息偏移量追踪技术的迅速发展，GRAY等[19]在1998年首次采用SAR影像的振幅信息解译得到研究区地表形变；ZHAO等[20]在2013年采用该技术处理PALSAR影像获取了内蒙古某煤田4.5 m的垂直位移。 利用SAR影像振幅信息通过多级配准获取地表形变信息，不需要解缠，不受矿区相干性制约，能够获取大量级形变，但对配准精度要求极高，获取的沉陷盆地边缘形变精度较低。欧空局(ESA)发射的Sentinel-1A/B SAR系列卫星，具有开放数据源、短重访周期、超幅宽、多极化等特点[21]，大大缩短了时间基线，在矿区开采沉陷准动态监测中具有很好的应用前景。

本文基于SAR影像振幅信息的动态时序开采沉陷监测方法，利用时序振幅数据提取开采沉陷区地表后向散射系数，分析开采沉陷过程中地表不同区域的后向散射系数与沉降速率的关系，较好地克服了常规SAR技术难以识别大量级沉陷区域的问题，可以提高开采沉陷监测的时空分辨率。以阳煤新元矿某工作面为研究区，利用该方法和16景Sentinel-1A数据，获取研究区振幅时间序列图和地表后向散射系数曲线图，并使用地表后向散射数据分析了沉陷区不同区域时序变化幅度，对地表后向散射系数与沉降速率关系进行研究，发现二者具有极强的相关性。

1 研究区与数据

1.1 研究区概况

阳煤新元矿位于沁水煤田西北部，太原至阳泉的中间地带，寿阳-阳泉构造堆积盆地地区，海拔1 058.3～1 110.9 m，属黄土高原温带季风气候区，天气寒冷干燥，温差变化大，年均降水量约500 mm，多集中在6～9月。 研究区地表大部分被第三系、第四系松散地层覆盖，经长期风化剥蚀，沟谷纵横，地形复杂。矿区南北长4～9 km，东西长15～16 km，地处东经112°58′51″～113°09′33″，北纬37°49′54″～37°52′09″。 工作面倾向长240 m，走向长1 207.08 m，平均采深为579.05 m，煤层厚度为3.40 m，平均倾角为3°，采用走向长壁全陷综合机械化一次采全高方法开采。研究区上方布设观测站进行地表岩移监测，点间距为19～38 m，监测日期为2017年5月23日至2017年10月30日。研究区域位置见图1。

图1 研究区域地理位置Fig.1 Geographical location of the study area

1.2 影像数据

本文选用欧空局Sentinel-1A卫星自2017年5月13日至2018年2月25日获取的16景C波段升轨单视复数影像(SLC)，数据重访周期12 d，空间分辨率为5 m(R)×20 m(A)，极化方式为VV，成像模式为Interferometric Wide Swath (IW)，原始数据宽总幅超过250 km，使用分辨率为90 m的SRTM3 Version 4数据来消除地形相位误差。 SAR影像信息见表1。

表1 所用Sentinel-1A数据信息Table 1 Sentinel-1A used for data information

2 基于振幅的开采沉陷时序InSAR监测方法

2.1 振幅反演地表形变特征基本原理

煤矿开采引发地表破坏，造成地物目标的物理结构、表面粗糙度或地物目标类型发生变化，向散射能量也会发生相应的变化，这将导致SAR图像的亮度和色泽发生变化[22]。通过雷达回波信号所携带的振幅数据，获取不同地表散射目标的后向散射系数信息，数字化处理后得到像元亮度值(DN)，利用像元值表示开采沉陷过程中地表目标的后向散射变化来反演地表形变特征[23]。像元DN值与地表的后向散射系数σ0关系如图2所示。

图2 像元DN值与地表后向散射系数关系Fig.2 Relation between pixel DN value andsurface backscattering coefficient

像元(i,j)的强度或DN值与相应地面单元的后向散射的能量Pr的关系见式(1)[24]。

(1)

式中，a、b为比例常数。

雷达所接收的能量Pr与地面单元雷达后向散射系数σ0(i,j)成比例关系见式(2)。

Pr(i,j)=k(i,j)σ0(i,j)

(2)

联合式(1)和式(2)得DN(i,j)与σ0(i,j)的关系式，见式(3)。

(3)

结合分布目标的雷达方程见式(4)。

(4)

由式(4)可得式(5)。

(5)

从式(5)可以看出k(i,j)与发射功率Pt、到地表散射目标的距离R、地面分辨率单元A、波长λ及天线的增益函数G有关。当k(i,j)在绝对或相对尺度上为已知值时，可以生成绝对尺度或相对尺度上的雷达回波图象。

2.2 雷达影像回波强度的影响因素

1) 传感器因素。雷达图像表示地面雷达后向散射的估算值，在高后向散射区域大部分能量反射回雷达系统中表现浅色调，反之，为深色调。从雷达方程可知，雷达回波强度与入射波长直接相关。同时，雷达遥感系统所选择的波长长短，决定了表面粗糙度的大小以及入射波穿透深度的能力，也直接影响到雷达回波的强弱；不同的极化方式也会导致目标对电磁波有不同的响应，使雷达回波强度不同，需要根据地面不同特点选择极化方式[25]。除此之外，俯角、照射带宽度也会影响回波强度。

2) 地面因素。影响回波强度的地面因素主要包括地表特性和负介电常数。开采沉陷过程中地表特性(坡度、坡向等)发生随机变化，雷达图像中会出现透视收缩、叠掩、阴影等现象[26]。同时，回波路径的改变使原有地物的散射机制(面散射、体散射等)发生改变，进而影响回波强度。负介电常数由物质组成及温度决定，其值越大，回波强度越强，矿区裸露的岩石、建(构)筑物反射能力强，在影像中表现为高亮度区域，植被覆盖密集区表现灰暗[27]。

3 时序振幅变化及地表沉降速率分析

3.1 研究区时序后向散射系数分析

利用多时相影像处理技术，将研究区影像进行多视、配准、滤波、地理编码和辐射定标后，获得2017年5月13日至2018年2月25日期间新元矿某工作面时序后向散射系数图(图3)。从图3可以看出，工作面附近的房屋等建(构)物相干性最高，反射微波能力强，后向散射系数值接近5，在影像中表现较为密集；其次，地表及山坡裸露的岩石相干性较高，反射微波能力较强，数量最多，在影像中分布不均且呈点状分布，后向散射系数值为3～4。工作面内地物整体后向散射系数值较小，很难选取到密集高的相干点，但仍有部分高相干地物随机出现。这是因为开采活动的剧烈进行，造成地表破坏，裸露的岩石或高相干地物发生出现-掩埋-再出现的现象，同时受卫星重访周期和植被季节性变化的影响，无法获取短暂暴露的高相干地物回波信号。利用开采活动影响范围内地表后向散射系数的变化，反应整个沉陷盆地内地表沉降速率变化是本文的研究重点。

图3 后向散射系数时间序列图Fig.3 Backscattering coefficient time series diagram

3.2 矿区分析模型的建立与地表沉降速率分析

3.2.1 建立分析模型

煤矿开采活动诱发地表发生移动，形成移动盆地，各个区域的移动和变形的性质及大小不尽相同，分为中性区域、压缩区域和拉伸区域[28]。为分析整个沉陷盆地的沉降特征，依据开采活动过程中工作面及附近地表变化造成其后向散射系数不同，结合地表达到充分采动时采空区的长度和宽度均要达到并超过1.2～1.4H0(H0为平均开采深度)[28]，考虑走向(C线)和倾向(A线)主断面的监测点实际布设情况，将研究区整体划分为工作面中心、工作面内边缘与工作面外边缘(图4)。工作面内受开采沉陷影响严重，不同区域地表剧烈变化程度不同，为更加准确地研究工作面内地表后向散射系数与沉降速率的关系，利用ArcGIS软件将工作面划分为九块相同大小的分区[29]。其中，工作面外边缘监测点包含点A1～A20、B1～B14、C1～C19；工作面内边缘1区包含点B30～B32，2区包含点B22～B29、点C20～C39，3区包含点B15～B21，4区包含点A35～A38，6区包含点A21～A25；工作面中心5区包含点A26～A34。在实际分析过程中，依据开采沉陷盆地的对称性，选取走向和倾向线部分监测点的像元后向散射系数进行分析，进而反演整个沉陷盆地的地表后向散射系数变化。

图4 开采沉陷分析模型图Fig.4 Mining subsidence analysis model

3.2.2 矿区地表沉降速率分析

为了分析研究区的开采活动与后向散射系数的关系，在走向和倾向观测线上分别以4个监测点为间隔提取出21个监测点内像元的后向散射系数。其中，走向线上的监测点主要集中在工作面外边缘(C1、C5、C10、C15、C20)与内边缘(C25、C30、C39)，为保证工作面走向信息的完整性，增加内、外边缘分界点C19与中心点A30。倾向观测线上选取工作面外边缘点(A1、A5、A10、A15、A20)、内边缘点(A25、A35、A38)与中心点A30,为保证沉降特征区域后向散射系数的特殊性，增加内、外边缘区分界点A21和内边缘、中心分界点A26进行研究。 利用简单的折线绘制出时序后向散射系数图(图5和图6)。

图5 走向监测点时序后向散射系数Fig.5 Backscattering coefficient of strike monitoring points

由图5可知，后向散射系数在0.10～0.36之间，整体呈低-高-低趋势。走向观测线各点在沉陷开始阶段与衰退阶段地表后向散射系数较低，2017年5月至2017年6月，地表后向散射系数较小但整体处于上升趋势，该时间段采矿活动刚开始进行，地表移动变形出现滞后。2017年11月后地表后向散射系数减小并趋于稳定，该时间段采矿活动停止，地表存在微小的残余变形。2017年7月至2017年11月期间地表后向散射系数较高呈起伏变化，该日期与开采活动日期基本吻合。同时，工作面中心A30、内外边缘分界点C19、C20、靠近工作面中心点C39，在采矿活动进行时后向散射系数整体较高，且最大后向散射系数依次降低。该研究区域工作面中心地表沉降速率大，地表活动剧烈，实测沉降大于2.7 m；边界两侧受采矿活动影响，地表活动存在差异。在内外边缘区的其他各点后向散射系数较低，无明显规律，同时这些区域受采矿活动影响小，地表移动变形不剧烈。

在倾向时序后向散射系数图(图6)中各点与走向各点规律相似，在采矿活动开始阶段、活跃阶段、衰退阶段各点后向散射系数呈现低-高-低变化，沉陷活跃阶段后向散射值较高且呈高低起伏变化。2017年5月至2017年6月地表后向散射系数较小但逐渐增加，2017年7月至2017年11月期间地表后向散射系数值较大且变化明显，2017年11月后地表后向散射系数减小并趋于稳定。倾向线上所有监测点后向散射可分为两类：①在工作面中心A30、中心与内边缘分界点A25，后向散射系数整体较高，这两个区域沉陷量极大，地表移动变形也较为剧烈；②其他各区域点后向散射系数整体较低，中心区域的最外侧点A26值较大，最外边缘点A1值最低。

图6 倾向监测点时序后向散射系数Fig.6 Backscattering coefficients ofinclination monitoring points

由走向线和倾向线上监测点的后向散射系数变化情况可知，工作面边缘区域受开采活动影响较小，后向散射系数较低；工作面中心和各分界点处采矿活动影响剧烈，地表后向散射系数值也较大。

进一步分析工作面中心区域(5区)可以发现其后向散射系数分三类(图7)：①高后向散射系数区(A29-A33)；②中后向散射系数区(A26、A28)；③低后向散射系数区(A27、A34)。由沉陷盆地对称特征分析可知，工作面中心区域后向散射系数整体由边缘至中心逐渐增大。

图7 工作面中心监测点时序后向散射系数图Fig.7 Backscatter coefficients of monitoring pointsat the working face center

4 开采沉陷区地表沉降速率与后向散射系数相关性分析

由上述分析可知采矿过程中移动盆地不同区域内地物单元后向散射系数不同。为了验证开采沉陷引起的地表变化速率与后向散射系数的相关性，评估多时相后向散射系数结果的可靠性，结合实测沉降速率，选取走向外边缘点C1、内外边缘分界点C20、内边缘点C39、倾向外边缘点A1、内边缘与中心区域分界点A25、最大下沉点A29，分别进行走向和倾向监测点沉降速率与后向散射系数对比分析。

从图8和图9可知，只有在部分时段沉降速率与后向散射系数有一定关系，2017-05-13至2017-08-21时段，走向地表沉降速率与后向散射系数变化趋势不同，倾向上二者变化趋势一致；2017-08-21至2017-11-29时段，走向地表沉降速率与后向散射系数变化明显且趋势相同，倾向上地表沉降速率较小，对应后向散射系数起伏变化明显，二者相关性较低。由此可知，2017-05-13至2017-08-21时段，走向监测点距推进面较远，开采活动还未影响该区域地表沉降，地物长势变化造成后向散射系数变化，造成地表沉降速率与后向散射系数相关性降低；2017-08-21至2017-11-29时段，倾向地表沉降处于衰退期，沉降速率变小，而地物长势、地物类型变化叠加影响后向散射系数呈伏变化，造成地表沉降速率与后向散射系数相关性降低。由上可得，靠近工作面中心的大量级沉降区域后向散射系数变化主要由地表沉降速率变化引起。

图8 走向监测点沉降速率与后向散射系数关系Fig.8 Relationship between settlement rate andbackscattering coefficient at strikemonitoring point

图9 倾向监测点沉降速率与后向散射系数关系Fig.9 Relationship between settlement rate andbackscattering coefficient of inclinationmonitoring point

利用SPSS软件对工作面特征点后向散射系数与沉降速率作相关性分析见表2。由表2知，走向外边缘点C1、内外边缘分界点C20、内边缘点C39的后向散射系数与沉降速率的相关关系均大于0.850，平均相关性达到0.878，说明存在极强的相关性；倾向方向监测点后向散射系数与沉降速率相关性由外边缘区至中心区域相关性逐渐增加，最大下沉点A29的相关性为0.960，最外边缘点A1的相关系数为0.696，相关性较弱，工作面中心与内边缘分界点A25的相关性达到0.836，相关性较强。倾向主断面平均相关性为0.831略小于走向，说明倾向主断面的点相关性整体较走向低，结合地表沉陷规律，发现煤矿开采引发地表移动变形使不同区域地表后向散射系数发生改变，同时地表后向散射系数可以作为描述整个沉陷盆地变化特征的指标。

表2 工作面特征点后向散射系数与沉降速率相关关系Table 2 Correlation between backscattering coefficient ofcharacteristic points of working face and settlement rate

5 结 论

1) 针对传统InSAR利用相位信息难以解缠矿区大量级沉陷问题，本文提出利用多时相SAR振幅数据提取沉陷区地表后向散射系数分析矿区地表沉降速率规律的方法。利用时序振幅法得到了多时相后向散射系数与地表沉降速率的关系，并发现地表后向散射系数可以描述大量级沉陷区域的地表沉降速率变化。

2) 通过建立研究区后向散射系数分析模型，发现沉陷区地表后向散射系数由工作面边缘至中心增大，与工作面走向和倾向地表沉降速率的平均相关系数r大于0.8；大量级沉陷区地表沉降速率与后向散射系数变化一致，相关性大，边缘区地表后向散射系数会受地物类型、地物长势及临采影响，两者相关性较低。