潘光永，陶秋香，陈 洋，王 珂

(山东科技大学测绘科学与工程学院，山东 青岛 266590)

0 引言

随着国家经济和科技水平的飞速发展，矿产资源的需求量越来越大。我国煤炭资源不仅存储量大，而且分布较为广泛，是重要的矿产资源之一，其广泛应用于生产生活的各个领域[1-2]。然而煤炭资源的开采必然会带来矿区及周围地表的沉降与塌陷等一系列地质灾害，为了降低地面沉降带来的安全隐患，预防矿产资源的过度开采，对矿区地面沉降的监测工作尤为重要[3-7]。

现有的地面沉降监测方法主要有水准测量、GPS(Global Position System)测量及合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)等[8]。孟宪纲[9]曾经采用水准测量的方法对京沪高速铁路进行了沉降分析，测量结果精确可靠，但是水准测量的方法在布设点位时难以做到分布均匀，更新周期较长，并且需要很多人力物力支撑[10-11]。张训虎等[12]采用GPS测量的方法对高速铁路进行了沉降监测，得出该方法可快速、高效的完成沉降监测工作的结论。GPS的方法虽然具有较高的定位精度，但对高程信息不敏感，也难以满足地表微小沉降量监测的需求[13-14]。InSAR技术为新型遥感技术，其差分干涉模式D-InSAR技术具有全天时、全天候、大范围等监测地表沉降的特点，其形变监测精度理论上可以达到厘米级甚至毫米级[15]。1989年Gabriel首次采用D-InSAR成功获取了农田区的形变信息，尽管学者们随后对D-InSAR技术做了一系列改进，但其精度易受时空失相干及大气等因素的影响[16]。2001年Berardino等人在D-InSAR基础上克服时间、空间失相干和大气效应的限制，发展了小基线集技术(Small Baseline Subset，SBAS)对长时间缓慢沉降的监测也有很好的效果[17]，因此成为了监测地表沉降的主要手段[18-21]。本文基于SBAS-InSAR技术对覆盖济阳矿区2017年5月20日至2018年10月18日期间40景C波段Sentinel-1A升轨数据进行干涉对的选取、去除大气效应、滤波、相位解缠、去平地效应等处理，得到期间矿区的沉降量及平均沉降速率，结果表明，基于SBAS-InSAR技术可以形象直观地再现矿区沉降变化过程，其沉降监测结果与矿区实际开采进程相符，可作为预防地质灾害的有力依据。

1 研究区概况

济阳井田矿区位于华北平原的南缘，黄河北岸，地势平坦，局部地势低洼，雨季易积水。黄河从井田南部1 km处流过，井田中东部有邢家渡引黄总干渠，中西部有大寺河通过。矿区东西向约9.9 km，南北向约5 km，面积约50 km2。本井田为全隐蔽的华北型石炭二迭系海陆交互相煤田。区内地层自下而上为奥陶系石灰岩、石炭系本溪组和太原组，二叠系山西组和石盒子组，第三系明化镇组，第四系平原组。主要含煤地层为太原组和山西组。具体位置及范围如图1所示(红色矩形框为矿区范围)。研究区主要为四采区7煤层，北靠DF15断层(H=19 m∠47°)，南邻李家断层(H=70～260 m ∠70°)，西至煤7风氧化带边界，东至DF16断层(H=80 m ∠70°)，平均走向长约2 378 m，倾斜长约954 m，面积2.04×106m2。四采区内7煤层位于太原组上部，上距一灰11.50～26.70 m，平均17.97 m，下距二灰30.50～40.30 m，平均34.72 m。下距第10层煤40.00～58.85 m，平均50.87 m。开采上限标高-375.0 m左右，下限标高-581.1 m左右，埋深398.2～604.3 m之间，平均506.3 m，煤厚0～1.81 m，平均1.1 m，倾角5°～12°，平均6°，工业储量为2.596×106t，可采储量为1.37×106t。其结构较简单，属稳定型薄煤层。主要布设有4701、4701东、4701西、4702、4702西、4703、4704上、4704下和4705工作面，其分布如图2所示，开采时间情况如表1所示。

图1 研究区地理位置Fig.1 Geographical location of research area

图2 工作面分布情况Fig.2 Distribution of working face

2 SBAS-InSAR技术原理

常规的D-InSAR的测量方法受制于SAR影像的数量，且易受时空去相关的影响，时间、空间基线过长等都会对差分的结果造成较大影响。由于矿区的开采过程较为缓慢，开采时间较长，因此应用D-InSAR的测量方法监测长时间段内矿区形变结果具有一定局限性。小基线集技术则可通过利用多景SAR数据对长时间的微小形变进行高精度观测，并通过产生的具有时间基线和空间基线特征相对较小的干涉图序列，更好的克服去相关的影响，提高形变监测质量。

假设研究区域在(t0,…,tN)的时间段里共有N+1景雷达影像，根据干涉相对之间的自由组合，便可得到共M幅干涉图，且满足：

(1)

先将tA,tB两景不同时刻的影像去除地形相位，假设tA

δφi(x,y)=φ(tB,x,y)-φ(tA,x,y)

(2)

式中φ(tA,x,y)和φ(tB,x,y)分别表示tA和tB时刻相对应的相位，λ表示雷达波长，d(tA,x,y)和d(tB,x,y)分别表示在tA和tB时刻任意像元(x,y)相对于基准时刻t0在LOS方向上的地表形变。

将公式(2)中的相位φ表示成任意两个获取时间之间的平均相位速度v与时间t的乘积，以此来获得具有物理意义的沉降序列，表达形式为：

(3)

于是，第i幅干涉图的相位可表示为：

(4)

将其改写成矩阵形式：

Bv=δφ

(5)

式中：B——M×N矩阵。最后，利用奇异值分解的方法得到速度矢量v的广义逆解，再根据对各时间段速度在时间域上进行积分便可得到各个时间段的形变量。

3 数据处理

研究采用2017年5月20日至2018年10月18日期间共计40景C波段Sentinel-1A升轨的雷达影像数据，高程数据(Digital Elevation Model，DEM)采用美国地质调查局所提供的SRTM(Shuttle Radar Topography Mission，即航天飞机雷达地形测绘使命)90 m 分辨率DEM，实验通过SARscape软件设置合适的时间阈值和空间阈值，共生成68个差分干涉对。差分干涉对的时空基线组合如图3所示。

图3 时空基线组合图Fig.3 Time-space baseline combination map

3.1 差分干涉处理

利用获取的DEM数据对生成的68个干涉对进行差分处理，生成相干性干涉图，通过Goldstein的滤波方法去除噪声影响，进行增强处理，滤波后的部分差分干涉图如图4所示。采用最小费用流(Minimum Cost Flow)的方法进行相位解缠，保证影像中相对孤立的高相干区域有更好的处理结果，获取研究区域真实相位。

3.2 轨道精炼和重去平

为去除矫正相位偏移，选取位于没有形变或形变较小的区域且具有高相干性的GCP点用于进行修正干涉相位和解缠后的相位，进而利用线性模型对初始位移进行估计，去除残余地形。

3.3 SBAS反演

为了更好的分析形变监测的结果，首先估算形变速率与残余地形，然后对差分干涉图进行二次解缠处理，进而对干涉图进行优化，得到第一次反演结果。在第一次反演得到的形变速率的基础上，利用轨道精炼中所选的GCP点移除恒定相位或斜坡相位，根据大气滤波的定制，估算并去除大气效应，最终得到时间序列上的位移变化。获取该矿区年平均沉降速率图(图5)。

图4 滤波后的干涉图示例Fig.4 Example of filtered interferogram

图5 年平均沉降速率Fig.5 Annual average sedimentation rate

4 结果分析

4.1 与矿区开采进程的对比分析

由图5沉降监测结果分析可知：

(1)监测期间该研究区出现三处沉降盆地，分别位于研究区西北、西南和东南方向的工作面附近，其中4705工作面与研究区西北方向4701和4702工作面附近沉降最为严重，其最大年平均沉降速率可达320 mm/a，矿区的沉降速率由沉降中心到边缘逐渐减小。

(2)为进一步获取监测期间矿区的累积沉降结果，对各时间段沉降速率进行时间域上的积分，获取2017年5月20日至2018年10月18日期间共39幅累积沉降量图。图6给出了研究区部分成像时刻相对于2017年5月20日的累积沉降量图。

图6 累积沉降量Fig.6 Cumulative settlement

由图6结合煤层各个工作面实际开采情况分析得到：

①2017年5月至2017年12月期间，随着开采工作的进行，4702工作面上方沉降持续增大，最大沉降量达到130 mm；4701西工作面随开采工作的开始，上方出现下沉，并且该工作面南邻刚刚结束开采的4701工作面，其上覆岩层尚未达到力学平衡状态，仍在缓慢下沉，从而导致地下岩层结构更趋于不稳定状态，加剧了沉降现象，最大沉降量达220 mm。

②2017年12月至2018年5月期间，4701西、4703及4701东工作面开采工作正在进行。根据图6(c)与图6(d)可看出：随着4701西工作面开采工作继续进行，西南处沉降盆地持续下沉，沉降范围有所扩大，最大沉降量达到390 mm；4703工作面北邻刚刚结束开采的4702工作面，随着其开采工作的进行，加剧了该区域的地面沉降，最大沉降量达300 mm。

③2018年5月至2018年11月期间，4701东工作面因开采工作继续进行，致使该区域沉降量及沉降区域逐渐增大，至2018年10月份开采工作结束，该区域沉降量最大达到300 mm，并且因该区域地下煤层开采导致东北方向地面也发生了沉降现象；4702西工作面从2018年8月份开始开采后便加剧了该地区的沉降，最大沉降量达到447 mm；由于4701西、4702和4703三个工作面的开采，使采区附近岩层结构遭到破坏，原有的力学平衡结构逐步变得不稳固，并使开采区上覆岩层发生移动和变形，从而影响周围地区的岩层结构，导致4705工作面附近在正式开采前就发生沉降现象，加上附近的4704工作面原本就进行过开采工作，地下岩层结构早已破坏，因此受附近三面开采工作的影响使该地区沉降更为严重，在2018年8月4705工作面开采后又加剧了沉降，至2018年10月18日沉降量最大达到447 mm。

(3)分别在沉降较小的4701东工作面附近选取A点，沉降较为严重的4705工作面附近选取B点进行时间序列分析，具体位置如图6(f)所示，结果如图7(a)和图7(b)所示。

图7 时间序列分析Fig.7 Time series analysis

由图7分析可知，两点处随着开采工作的进行，总体上沉降量越来越大且没有减缓的趋势。其中图7(a)可看出，监测初期，由于4701西和4702工作面的开采，A点处开始出现沉降，至2017年8月左右，沉降量达45 mm，之后因为该区域岩层硬度及4701西工作面开采强度减弱等因素致使该处沉降量明显减少。至2018年3月，4701东工作面开采活动的进行导致A点处再次出现较大下沉，随着开采活动的持续进行，沉降量逐渐增大，至2018年10月18日，最大沉降量大于140 mm。

图7(b)可看出在4705工作面附近因受三个方向工作面开采的影响，也在持续发生沉降，并且随着2018年5月4705工作面开采后，沉降量增大，最大大于400 mm。

4.2 与水准监测结果的对比分析

为了验证SBAS-InSAR监测结果的可靠性，分别在沉降盆地边缘及靠近沉降中心处选取了7个水准点的累积沉降结果与SBAS-InSAR监测的结果进行对比分析。选取的水准点布设情况如图8所示，对比结果详见表2。

图8 水准点位置Fig.8 Location of bench marks

表2 SBAS-InSAR形变监测结果与水准测量结果差异比较

根据表2的对比结果可看出，在本次实验中所得到的形变结果与水准测量得到的形变结果吻合较好，两种技术监测结果最小差值为-7.71 mm，最大差值为28.12 mm，经计算平均差值约为5.91 mm。由于SAR影像自身特性的限制、干涉过程中受到各种误差因素的影响，SBAS-InSAR的监测结果与实际水准数据有所差别，但从表2可看出大部分点形变的差值相对较小，表明SBAS-InSAR监测的结果是可靠的。

5 结论

本文基于C波段Sentinel-1A升轨的雷达影像数据，采用SBAS-InSAR技术方法对山东济阳某矿区进行了监测，分别从研究区内地面沉降的年平均沉降速率、累积沉降量以及沉降较为严重区域的时间序列进行了分析。结果表明，2017年5月20日至2018年10月18日期间，研究区沉降最为严重的区域集中在4701西与4702西工作面和4705工作面，年平均沉降速率最大达到320 mm/a，累积沉降量最大为447 mm，且随着时间推移，沉降量逐渐增大。此外，位于矿区周边的104国道与济南绕城高速受开采区影响，也发生了部分沉降，说明了在研究区工作面进行开采工作后，不仅开采地区会有沉降现象发生，周围地区也受其影响而出现下沉。通过与矿区实际开采情况的对比分析，表明SBAS-InSAR的方法在矿区沉降监测上有很大的适用性，能够快速准确地定位矿区沉降的位置及范围，且其监测结果与矿区开采进程相符，并对因矿区开采而导致的地质灾害问题的预防有重要意义。