基于SBAS-InSAR技术的关闭矿井地表多维形变时空监测与分析方法

张连蓬 梁 亮 陈炳乾 胡晋山 于 洋 秦 璐 余 昊 杨家乐 杨 宇

(1.江苏师范大学地理测绘与城乡规划学院,江苏 徐州 221116;2.长安大学地质工程和测绘学院,陕西 西安 710054;3.苏州工业园区测绘地理信息有限公司,江苏 苏州 215000)

煤炭是支撑我国国民经济发展最重要的基础能源,对于我国经济发展具有重大的战略意义。近年来随着我国去产能政策的出台和大量矿井的煤炭资源枯竭,导致大量矿井被关闭。矿井关闭后,煤岩体在应力、地下水等多种因素的作用下发生风化劣化、强度降低,将改变废弃采空区内破裂岩体的应力和承载能力,有可能导致采空区地表发生二次形变或多次形变,同时由于地下水的作用,煤柱逐渐软化脱落甚至垮塌,引起不同程度的地面塌陷等地质灾害[1-2]。合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar,InSAR)技术是近年来兴起的一种新的空间对地监测手段,相比于传统的测量方法,具有全天时、全天候、高精度和大范围等优点[3-4]。传统的In-SAR技术在地面微小形变监测中取得了重大突破[5],但其仍受时空失相干和大气延迟的影响[6]。为了克服该局限,国内外学者先后提出了时间序列InSAR(Time-series InSAR,TS-InSAR),主要包括小基线集差分雷达干涉测量(Small Baseline Subset In-SAR,SBAS-InSAR)[7]、永久散射体雷达干涉测量(Permanent Scaterer InSAR,PS-InSAR)[8],临时相干点雷达干涉测量(Temporally Coherent Point InSAR,TCP-InSAR)[9]等。

最初的时序InSAR技术主要应用于城市地表形变监测,随着该技术应用领域的拓展,逐渐也被广泛应用于矿区形变监测[10-12]。尹宏杰[13]利用L波段的PALSAR数据并基于SBAS-InSAR技术获取了湖南某矿区的沉降序列图;邢学敏等[14]提出了一种角反射器InSAR与PS-InSAR的联合解算方法,并将其应用于区域线性沉降探测中;FAN等[15]利用TCP-In-SAR监测了中国陕西榆阳矿区的地表沉降,得到最大沉降速率为162 mm/a,累计沉降量为86 mm;潘光永等[16]利用SBAS-InSAR技术对济阳井田矿区40景C波段Sentinel-1A升轨数据进行处理,获取了2017年5月20日—2018年10月18日研究区内地面沉降的年平均沉降速率和累积沉降量。结果显示,区内年平均沉降速率最大达到320 mm/a,累积沉降量最大为447 mm。李达等[17]针对传统InSAR技术易受时空失相关、大气相位延迟等因素的影响,采用SBASInSAR技术对13景TerraSAR-X数据进行时序处理,分析了研究区域2012—2013年沉降速率;李路等[18]针对矿区地表发生严重时空失相干现象影响小基线集技术的应用,从提高监测点的数量与质量角度出发对该技术进行了改进,并应用于实地矿区,结果显示:区内共有两处区域发生大量级沉降现象,其中北侧采空区最大沉降速率可达-50.948 mm/a,南侧采空区虽较缓和,但也对周围村庄造成了影响,应引起注意。KUMAR等[19]利用改进的PS-InSAR技术分析了Jharia煤田地下采矿活动引起的缓慢地表变形,发现区内最大缓慢变形速率为29 mm/a,累计沉降值为90 mm,并绘制了沉降区域图;DANG等[20]基于2015—2019年的56景Sentinel-1A影像,利用PS-In-SAR技术,监测了越南Ha Long和Cam Pha 地区的地面沉降,分析了地面沉降主要是由煤炭开采引起,同时讨论了该地区地表形变的空间分布特征,确定了沉降速率随时间的变化关系。

现有研究大多集中在矿区开采过程中产生的地表形变,但近年来针对关闭矿井后地表形变监测的研究也日益受到关注。本研究以徐州东部矿区为例,基于SBAS-InSAR技术,提出了关闭矿井地表形变时空监测与分析方法。利用2015年7月6日—2021年11月19日的171景Sentinel-1A SAR数据,监测并分析了徐州东部矿区闭矿后6 a内的地表多维变形时空演化规律,分析结果对于安全、有效、合理地利用废弃矿井资源和保证各种建(构)筑物以及人民生命财产安全具有重要的理论意义和实用价值。

1 SBAS-InSAR技术原理

SBAS-InSAR技术是一种基于多主影像的时间序列InSAR方法,将同一地区的多景SAR影像数据按照一定规则组成若干个小基线干涉对,运用最小二乘法或者奇异值分解方法进行形变参数估计[21-22]。

假设N+1幅SAR影像按照时间序列t0,t1,t2,…,tn排列,选取其中一幅主影像对其他影像进行配准,基于时空基线阈值进行干涉对组合生成M幅多视差分干涉图[23]。假定主影像获取时刻为tB,辅影像获取时刻为tA(tB＞tA),若不考虑大气延迟相位、残余地形相位和噪声相位,两者生成的第j幅差分干涉图对应的任一像元(x,r)的差分干涉相位Δφj(x,r)可以表示为

式中,j∈[ 1,2,…,M];λ为雷达的中心波长;(x,r)分别为像元在方位向和距离向的坐标;φA(x,r)、φB(x,r)分别为tA、tB时刻相对于初始时刻t0的形变相位;d(tB,x,r)和d(tA,x,r)为tB和tA时刻相对于d(tB,x,r)=0的雷达视线方向的累积形变量。

若IE={IE1,…,IEM}和IS={IS1,…,ISM}分别为干涉数据处理时按时间顺序排列的主影像和辅影像序列,并且满足IEj＞ISj,j=1,2,…,M,则所有的差分干涉图相位可以组成如下观测方程:

将其转换为线性方程为

式中,A为M×N的系数矩阵,其中矩阵内主、辅影像对应的系数分别为1和-1,其余值为0;φ为待求的相位矩阵;δφ为差分干涉图的相位矩阵;以t0时刻获取的影像值为参考值,式(3)中未知数为N个。

当矩阵A为满秩,即M≥N时,即可利用最小二乘法[24]求解出方程(4)的最优估计值

式中,若矩阵A秩亏,即M＜N,对应的ATA为奇异矩阵,可利用奇异值分解法求其最小范数解,进而获取研究区的年平均形变速率[25]。获取形变速率后,在其时间维度上进行积分运算即可获得研究区相应时间段内的累积形变量。

2 研究区域与数据

徐州市位于江苏省北部(图1(a)),处于陇海、津浦铁路干线交点,京杭运河流经矿区,交通发达,煤炭储量较丰、品种较全。研究区域位于贾汪区境内,主要有青山泉、权台、旗山、董庄、大黄山、韩桥等矿井,开采历史超过百年,由于煤炭资源枯竭,目前徐州东部大小煤矿均已关闭。研究区域内包含了青山泉矿、韩桥矿、权台矿、旗山矿、董庄矿和大黄山矿6座矿井,其中大黄山矿关闭时间为2000年1月,董庄矿于2001年11月关闭,青山泉矿和韩桥矿都于2008年12月关闭,权台矿和旗山矿在2016年9月关闭。研究区位置如图1所示,图中3处圆点为二等水准测量点。

图1 研究区概况Fig.1 Overview of the study area

试验过程中共使用171景Sentinel-1A IW模式的SAR卫星影像数据,时间跨度为2015年7月6日—2021年11月19日,影像的距离向和方位向空间分辨率分别为5 m和20 m,极化方式为VV,中心入射角约39.5°,SAR影像参数具体取值见表1。获取的Sentinel-1A SAR影像重访周期大部分为12 d,部分重访周期大于12 d的影像已在表1中列出。

表1 SAR影像参数Table 1 Parameters of SAR image

试验过程中采用NASA提供的空间分辨率为30 m×30 m的STRM 1 DEM数据(http:∥gdex.cr.usgs.gov/gdex/)来消除地形误差影响。另外,利用欧洲航天局提供的定位精度在5 cm以内的POD精密轨道数据(https:∥qc.sentinel1.eo.esa.int/aux_poeorb/)消除轨道误差。

3 试 验

3.1 形变速率与累计形变量分析

试验中进行干涉对组合的垂直基线和时间基线阈值分别设置为120 m与60 d,最终生成可用干涉组合691个,干涉组合的时空基线组合分布如图2所示。首先对生成的691个干涉对进行差分干涉处理得到对应的差分干涉图;然后利用最小费用流方法[26]对其进行相位解缠,并根据干涉图的相干性以及解缠结果筛选出高质量的差分干涉图,而后选取地表稳定的控制点进行轨道精炼和重新去除平地效应,之后去除大气相位误差;最后利用最小二乘方法得到地表年平均形变速率,如图3所示。

图2 SAR影像时空基线组合Fig.2 Spatiotemporal baseline combination of SAR images

图3 地表沿视线向年平均形变速率Fig.3 Average annual rate of surface deformation alone line of sight

由图3可知:研究区域2015年7月6日—2021年11月19日,地面最大抬升速率为36 mm/a,位于图3权台矿点A′处;最大沉降速率为-33 mm/a,位于图3董庄矿点B′处,由光学影像可知,点A′大致位于潘安湖周边的裸地处,点B′位于耕地周围。其中,权台矿以抬升为主(图3A区域),最大抬升为36 mm/a,最大沉降速率相对较小,为-20 mm/a;同时旗山矿也有部分区域(图3B区域)表现为抬升,最大抬升速率为30 mm/a,但还是以沉降为主(图3C区域),最大沉降速率为-29 mm/a。其余4个矿区都以沉降为主,其中董庄矿的最大沉降速率居于4矿之首,为-33 mm/a;大黄山矿次之,为-28 mm/a;青山泉矿和韩桥矿相对来说沉降较小,最大沉降速率分别为-18 mm/a和-23 mm/a。

为了进一步分析矿井关闭后地表形变的发展过程和变化规律,在获取形变速率的基础上对其进行时间维度上的积分运算,获取了各个矿区时间序列的地表累积形变量的变化特征,如图4所示。由于研究时间跨度较大,同时研究区总体形变较小,故将时间序列间隔设置为6个月,共获得了14幅累计形变图。

图4(a)～图4(c)为权台矿和旗山矿关闭前获取的影像。在矿井关闭之前,旗山矿已开始出现沉降区域,至2016年7月24日(图4(c)中椭圆C),沉降量最大达到-130 mm,而权台矿呈微量抬升(图4(c)中椭圆A),至2016年7月24日,最大抬升量为81 mm。这是因为权台矿在2014年4月以前就已停采,当时部分排水设施已关停,导致地下水位上升,作用在采动破裂岩体、第四系松散层内的孔隙压力增大,有效应力减小,使其产生弹性恢复变形,膨胀岩体遇水膨胀也可能导致覆岩及地表上升。为了防止突水事故的发生,同时也是保证旗山矿的安全开采,权台矿保留了部分排水设施待旗山矿关闭后,权台矿-600 m水平及以上排水设施才停止运行[27]。

在所有矿区关闭之后,旗山矿的D区域也开始出现沉降(图4(f)中椭圆D);而旗山矿的C区域出现了很明显的先下沉后抬升的现象,至2019年7月9日,如图4(i)所示,旗山矿达到最大沉降值-216 mm,自此C区域开始逐渐抬升;位于旗山矿西南部的B区域也于2020年7月15日(图4(k)中椭圆B)开始出现抬升,至2021年11月19日,B区域整体呈抬升状态,C区域虽仍表现为沉降,但在闭矿的6 a里,存在一个先下沉后抬升的一个过程,这与文献[28]所得结论一致。

图4 2015年7月6日—2021年11月19日地表LOS方向时序形变特征Fig.4 Time series deformation characteristics in the LOS direction between 6 July 2015 and 19 November 2021

如图4(i)中位于大黄山矿的椭圆E,自2018年7月2日开始出现沉降,至2021年11月19日,大黄山矿的最大沉降量为-200 mm。董庄矿自2017年7月19日起(图4(e)),沉降区域由董庄矿边界向内扩张,至2021年11月19日,董庄矿的最大沉降量达到-218 mm。相对来说,青山泉矿和韩桥矿的沉降较小,最大沉降量分别为-146 mm和-101 mm,值得注意的是大黄山矿和董庄矿已经关闭了十余年,但是沉降量仍然相对较大,推测可能是矿井关闭后人类活动(如地表兴建设施)造成采空区“活化”产生了二次沉降。

3.2 倾斜与曲率形变分析

矿井关闭后地表产生的残余变形除了下沉外,还有倾斜、曲率以及水平移动和水平变形。由于现场观测发现地表目前无明显水平移动,因此本研究重点评估地表倾斜和曲率变形对地表建筑物的影响程度。为获取地表倾斜和曲率的连续形变场,对图4中的LOS方向形变值先转成垂直形变wi再进行插值,之后根据开采沉陷理论公式[29]计算出地表各点时间序列的倾斜变化值(图5)与曲率变化值(图6)。倾斜值计算公式为

图5 2015年7月6日—2021年11月19日各个矿区倾斜值分布Fig.5 Distribution of tilt values in each mining area between 6 July 2015 and 19 November 2021

图6 2015年7月6日—2021年11月19日各矿区曲率值分布Fig.6 Distribution of curvature values in each mining area between 6 July 2015 and 19 November 2021

式中,i2-3地表点2和3之间的倾斜值,mm/m;w2、w3为地表点2和3的下沉值,mm;l2-3为地表点2和3之间的水平距离,m。

由于曲率是倾斜变形的一阶导数,因此曲率可进行如下计算:

式中,k2-3-4为地表点2、3和4之间的曲率值,mm/m2;i3-4为地表点3和4之间的平均斜率,mm/m;l3-4为地表点3和4之间的水平距离,m。

由图5可知:2015年7月6日—2018年1月15日,各个矿区的倾斜值均在《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》的地表允许变形值范围内(允许变形值i=3 mm/m,K=0.2 mm/m2)[30],地表处于相对安全状态(图5(a)至图5(f));2018年1月15日—2021年11月19日,旗山矿出现一处地表倾斜值超过3 mm/m,达到3.35 mm/a(图5(g)至图5(n)中五角星符号标注区),其他5座矿区没有超过允许形变的点;2019年7月9日—2021年11月19,董庄矿出现一区域地表倾斜值超过3 mm/m,该区域内的最大倾斜达到3.5 mm/a(图5(k)至图5(n)中圆圈符号标注区)。

为进一步定量分析地表形变演化规律,分别在图4、图5和图6上提取 2015年7月6日—2021年11月19日期间的地表最大形变点,分析其动态形变规律,结果如图7所示。

图7 各矿区的最大形变点的时间序列变化特征Fig.7 Time series variation characteristics of maximum deformation points in each mining area

由图7(a)可知:旗山矿区的最大沉降点在闭矿后先下沉后再抬升,旗山矿最大累计沉降量为-216 mm,董庄矿最大沉降量为-215 mm,大黄山矿最大沉降量为-200 mm,另外3个沉降略小的韩桥矿、青山泉矿和权台矿的最大沉降量分别为-146、-102、-148 mm。自2018年7月2日之后(如图7(a)中虚线处),6个矿区的形变速率均有所降低,形变逐渐趋于平稳。这可以在一定程度上反映出在闭矿后随着时间的增长地表也趋于稳定,不再有剧烈的形变产生。由图7(b)可知:除了权台矿和旗山矿,其余4个矿区在2017年7月6日之后抬升速率变缓,权台矿的最大抬升量达到230.8 mm,旗山矿的最大抬升量为168 mm。相对来说,大黄山矿、董庄矿、韩桥矿和青山泉矿抬升量较小,分别为96、86、116、49 mm。

图7(c)、图7(d)为各个矿区的最大倾斜曲线,由于旗山矿的闭矿时间为2016年9月,所以在2018年1月15日之后(图7(c)、图7(d)中虚线处)所有矿区的倾斜才基本呈现平缓的变化特征。大黄山矿、董庄矿、韩桥矿、旗山矿、青山泉矿、权台矿的最大倾斜分别为-2.9、3.6、-2.9、-3.4、2.0、-2.9 mm/a,最大倾斜曲线反映出在闭矿后地表的不均匀变化导致倾斜值大于《建筑物、水体、铁路及主要井巷媒柱留设与压煤开采规范》中的允许形变值,需要开展持续监测来预防地质灾害发生。

图7(e)、图7(f)为各个矿区的最大曲率曲线,大黄山矿、董庄矿、韩桥矿、旗山矿、青山泉矿、权台矿最大曲率分别为0.15、-0.19、0.16、0.18、0.08、-0.18 mm/m2,在2017年7月6日之后(图7(e)、图7(f)中虚线处),矿区的曲率变化也变缓,且最大曲率都在允许形变值范围之内,处于相对稳定状态。

3.3 沉降面积统计分析

为分析试验区域内地表形变的年变化情况,对获得的形变值进行叠加,得到以年为间隔的形变值,通过插值得到研究区内完整的形变面积,同时对沉降面积进行了统计分析,得到各矿区沉降面积随时间变化柱状图,如图8所示。

图8 各矿区2015—2021年形变面积变化柱状图Fig.8 Histogram of deformation area change in each mining area from 2015 to 2021

由图8可知:2015年7月6日—2021年11月19日,权台矿整体呈抬升趋势,抬升面积从11.29 km2上升至21.1 km2,并有持续上升的趋势,下沉面积由17.19 km2降至7.38 km2。旗山矿地表呈现先下沉后抬升的趋势,2015年下沉面积为32.78 km2,2017年上升为33.66 km2,此后下沉面积开始减少,2021

年下沉面积减少至25.96 km2,而该矿在2015年的抬升面积为8.44 km2,2017年减少至7.21 km2,此后抬升面积开始上升,到2021年上升至15.27 km2,该矿于2016年9月闭矿之后,推测矿区地表后续会继续抬升。青山泉矿在6 a多内总体较平稳,下沉面积从2015年的9.26 km2上升至2021年的10.35 km2,抬升面积从2015年的2.78 km2下降为2021年的1.69 km2。相较于青山泉矿,韩桥矿下沉趋势较明显,2015年下沉面积为17.6 km2,到2021年上升至20.15 km2,抬升面积由14.96 km2降为12.42 km2。董庄矿整体呈下沉趋势,下沉面积由2015年的10.92 km2升为2021年的12.942 km2,在2015年有略微抬升,抬升面积为2.44 km2,此后抬升面积呈减少趋势,至2021年仅为0.43 km2。大黄山矿整体较为稳定,2015年下沉面积为17.76 km2,2021年为18.17 km2,抬升面积由2015年的5.21 km2下降为2021年的4.8 km2。

同时值得注意的是,虽然大黄山矿与董庄矿已经关闭了20 a以上,但是沉降区域面积仍然不小,甚至还有继续上升的趋势,因此推断是由于矿区关闭后,频繁的人类活动导致采空区的“活化”引起了地表下沉。

3.4 精度验证

在上述分析的基础上,本研究对SBAS-InSAR方法的可靠性进行进一步分析,结果如图9所示。

图9 SBAS-InSAR、PS-InSAR以及水准测量结果对比Fig.9 Comparison of the results of SBAS-InSAR and PS-InSAR and leveling survey

图9(a)是采用PS-InSAR得到的LOS向年平均形变速率结果,分析可知:PS-InSAR所获结果的年平均形变速率为-30～-29 mm/a,而SBAS-InSAR的年平均速率为-33～-36 mm/a,两者最大值差值为7 mm/a,具有高度一致性。同时分别选取SBAS-InSAR与PS-InSAR形变速率图上的31 174个同名点进行拟合,拟合结果如图9(b)所示。由图9(b)可知:SBAS-InSAR与PS-InSAR两者的均方根误差(RMSE)为1.6 mm/a,相关系数为0.83,说明两者在监测结果上具有高度一致性,从监测点空间密度上考虑SBAS-InSAR方法更具优越性。

为了进一步验证SBAS-InSAR监测结果的准确性,将SBAS-InSAR获取的LOS向形变值转为垂直沉降,并提取与图3中水准点相同位置的形变值与对应水准测量数据进行对比,结果如图9(c)所示(水准数据获取时间分别为2021年10月26日和2021年11月19日)。由图9(c)可知:水准数据与SBAS-InSAR监测结果在形变趋势上保持了一致性,同时最大绝对误差为1.3 mm;由图9(d)可知:水准和InSAR拟合结果的空间相关系数为0.81,RMSE为0.3 mm/a。以上分析表明:SBAS-InSAR技术获取的形变结果精度满足实际需求。

4 结 论

(1)基于SBAS-InSAR技术,提出了关闭矿井地表形变时空监测与分析方法,利用2015年7月6日—2021年11月19日的171景Sentinel-1A SAR数据,监测并分析了徐州东部矿区闭矿后6 a内的地表多维变形时空演化规律,发现矿井关闭后早期地表大多呈下沉趋势,而后会出现抬升趋势。

(2)徐州东部关闭矿区内地表最大沉降速率为-33 mm/a,累计最大沉降量为-219 mm;最大抬升速率为36 mm/a,累计最大抬升量为233 mm。将其结果与水准数据对比,表明SBAS-InSAR技术获取的形变结果精度能满足实际需求。此外,监测发现该区域地表最大倾斜与曲率值分别达到3.6 mm/m和-0.19 mm/m2,已经超过了建筑物允许变形值,需要对其进行加固和保持持续监测。

(3)本研究仅使用一种SAR数据来对关闭矿区进行地表形变监测,将多源遥感数据进行充分融合,对矿区开采引发的多方面灾害进行研究和预测是后续关注的重点。