开采煤矿所造成的地面塌陷、裂缝以及滑坡等地质环境灾害, 不仅降低了土地的使用价值, 而且会造成沉降区建筑物、道路及相关基础建设的损坏, 从而制约矿业城市的可持续发展。新采空区地面将快速、非线性沉陷; 老采空区则会在相当长一段时间内产生相对缓慢的残余沉降。为了有效地控制矿区沉降, 减少损失, 就必须通过一定的监测手段及时掌握沉降情况。合成孔径雷达差分干涉测量 (DInSAR) 是近年来发展起来的一种监测地表形变的有效技术手段, 能够获取研究区域短周期内空间连续的地表形变信息, 监测精度可达mm级^[1]。国内外相关学者陆续开展了应用DInSAR技术监测矿区沉降的研究: 1996年Carnec等^[2]首先用DlnSAR技术监测了法国Gardanne附近煤矿的沉降, 2003年Raucoules等^[3]用DlnSAR技术对法国Vauvert地区的一个盐矿进行了地面沉陷监测; 国内研究者应用该技术分别在武安、济宁及平朔等煤矿的监测中取得了很好的研究成果^[4]。为了克服时间、空间失相干和大气效应等因素的影响, 目前又发展了永久散射体(permanent scatter, PS)和小基线集(small baseline subset, SBAS)等时间序列形变分析方法, 这些方法能够监测长时间范围内地表缓慢形变, 获得研究区域的沉降发展和演化情况。2008年Baek等^[5]利用SBAS方法获取了1992— 1998年韩国Gangwon-do地区煤矿的形变结果; 2012年张学东等^[6]选用2004— 2010年27景ENVISAT ASAR影像进行分析, 查明了唐山城区的地面沉降量及其空间分布特征。

山西省长治地区分布着潞安矿区众多不同规模的煤矿。由于几十年的煤矿开采, 已有地质灾害实地调查表明该地区的沉降区面积大、分布集中、潜在塌陷区大。小基线集InSAR时序分析方法可以很好地解决地表缓慢形变监测问题, 但在矿区监测中难以获取快速大形变以及总体变形。为此, 本文针对开挖矿区变形区域, 应用DInSAR方法进行快速大形变监测分析; 用SBAS方法监测微小缓慢形变信息, 研究该区沉降形变分布特征及时间演化规律, 并对具体地区沉降情况进行相关分析。

2001年, Berardino等^[7]提出了小基线集(SBAS)方法。该方法将所有获得的SAR数据组合成若干个集合, 集合内SAR影像对要求的时间和空间基线距较小。利用最小二乘方法求解每个小集合的地表形变时间序列。为了增大时间采样率, 可以通过奇异值分解(singular value decomposition, SVD)方法将多个小基线集联合起来求解^{[7, 8]}。此后, 随着各种相干点目标提取算法、3D相位解缠算法和各误差项去除算法的提出, 以及时间序列形变模型的改善等关键技术的进步^{[9, 10]}, SBAS-InSAR方法也越来越完善, 测量精度也越来越高。

假设同一地区有按照时间顺序(t₀, …, t_N)排列的N+1幅雷达影像。主影像(IE)和从影像(IS)按时间顺序排列。则

δ ϕ _j=ϕ ( tIEj)-ϕ ( tISj), ∀ j=1, …, M(1)

写成矩阵形式为

Aϕ =δ ϕ , (2)

式中: ϕ 为各SAR影像中高相干点对应的相位组成向量; δ ϕ 为各干涉图解缠后相位组成向量; A为M× N阶矩阵, ∀ j=1, …, M, 若IS_j≠ 0, A(j, IS_j)=-1; 若IS_j=0, A(j, IS_j)=0; A(j, IE_j)=+1, 其他值为0。由式(2)可见, 矩阵A依赖于数据集生成的干涉图, 如果所有的数据属于单一的小基线子集, 那么M≥ N, A是一个N秩矩阵。可由最小二乘法来求解, 计算式为

ϕ˙=A^#δ ϕ , A^#=(A^TA)^-1A^T 。 (3)

然而, 所有的数据属于一个单一的子集是不常见的。所以, 为了提高形变信号的时间采样率, 必须要面对数据属于不同子集的情况。对于式(2)中A表现出秩亏时, 式(3)中的A^TA是奇异矩阵。假设有L个不同的小基线子集, A的秩则为N-L+1, 那么这个系统将有无穷多个解。引入奇异值分解算法, 可以估算矩阵A的伪逆, 得到式(2)的最小二乘解( ϕ˙)。

采用SVD方法直接对相位进行求解, 得到的累积形变在时间上表现出不连续性, 不符合沉降的物理规律。为了获得具有物理意义的沉降序列, 将相位表示为2幅SAR影像获取时间之间的平均相位速度和时间的乘积。利用SVD分解可以得到速度矢量的广义逆解。根据各个时间区间的沉降速度, 对各时段速度在时间域上进行积分即可得到各个时间段的形变量。

SAR影像相位中除了形变相位外, 还有由区域地形误差带来的相位、由获取影像时大气不均质带来的大气相位及由基线、时间失相干现象和热噪声等引起的相位。所以, 在处理过程中需要分离出各误差项, 最终获取形变相位信息, 具体技术流程如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 SBAS-InSAR技术流程Fig.1 SBAS-InSAR processing flow

受开采的影响, 研究区内一些矿区在30 d或者更短的时间内发生了较大形变, 以至于在SBAS方法中未能提取出相干点, 本文利用DInSAR的监测结果进行分析说明。对于县(市)区以及矿区周围相干性保持较好的居民区, 本文利用SBAS方法监测形变结果, 并对其中形变速率较大的地区进行时间序列分析。

3.1 DInSAR结果分析

由影像20050120与20050224干涉处理生成的差分干涉图, 得出研究区内分布着多个形变区。以图2中王庄矿区中部部分区域为例, 图4(a)为该区域差分干涉图, 图中每个形变区均表现出不同大小的干涉周期, 干涉周期最大达4个。由于所用数据为C波段, 波长约为5.6 cm, 故形变区最大沉降约为11 cm。图4(b)为对应区域Google Earth光学影像, 图中矿区位置与干涉图中形变区位置分布较一致。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 DInSAR监测的矿区形变结果(部分)Fig.4 Deformation results by DInSAR(part)

3.2 SBAS结果分析

研究区内整体表现为沉降趋势, 不同矿区沉降状况与其地形、开采时间及开采方式等有关。城区沉降状况与矿区影响有关, 也与城市建设、城市用水有关。由图5可知, 该研究区表现为沉降形变, 平均沉降速率在5~10 mm/a左右, 最大沉降速率达到15 mm/a。其中沉降较严重的区域在长治市西部(图5C区)、屯留县东部(图5B区)西北部(图5A区)以及襄垣县东部和南部地区。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 研究区平均形变速率图(A, B, C分别为余吾矿、常村矿及三元矿居民地形变区)Fig.5 Mean velocity map of research area

3.2.1 形变结果分析

实验处理过程中均以第一幅影像(20030710)为基准, 共生成22幅时序累积形变图, 图6为其中5幅。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 研究区视线向累积形变图Fig.6 Accumulated deformation maps in LOS direction of research area

从图6可以看出, 自 2005年7月起, 屯留县东部和西北部地面开始产生较大的沉降形变, 且随着时间的增加, 形变量逐渐增大, 特别是2009年12月以后, 屯留县有些区域的沉降量达到60 mm。长治市自2007年以后沉降明显, 且西部和南部沉降量较大, 这可能与三元矿和司马矿的开采有关。潞城市在2007年5月开始表现出明显的沉降趋势。襄垣县自2005年7月起开始沉降, 且随着时间的增加, 沉降量日益增大。

图7为图5中A, B, C 3个区域的时间序列形变图。结合图2中各矿区范围分析可得, 余吾矿区在2003年7月— 2010年7月期间沉降量较大, 达到60 mm, 平均形变速率约10 mm/a, 最大可达13 mm/a。常村矿区在2003年7月— 2006年2月期间沉降量为10 mm, 平均形变速率为4 mm/a, 此后形变速率增大, 2006年4月— 2007年5月形变速率为15 mm/a, 至2010年7月累积形变量近40 mm。三元矿区2003年7月— 2006年4月期间平均形变速率约为10 mm/a, 2010年1月— 2010年7月期间平均形变速率达到12 mm/a。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 矿区居民地代表区域时间序列形变图Fig.7 Time series deformation maps of residential area around mining area with large deformation rate

野外调查发现, 常村煤矿和漳村矿采空区的塌陷使208国道屯留路段和襄垣县闫村段均匀下沉^[13]。本次除了对居民区沉降进行监测外, 也对研究区内国道、省道形变进行了研究分析。图8为208国道南渔泽村— 上村镇段平均形变速率图, 研究发现, 该路段沉降较严重, 最大沉降速率达5~7 mm/a。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 道路形变监测结果Fig.8 Road deformation results

3.2.2 形变特征分析

矿区沉降一般不具有明显的线性特征, 开采塌陷常常表现出非线性特征。矿区开采一定时间后, 随着开采深度的增加, 造成地表发生移动的时间加长, 下沉速度减缓。残余沉降减少, 但会持续较长时间。由于王庄、五阳及漳村矿区开采时间较早, 残余形变逐渐减小。常村矿区、余吾矿区可采储量很大, 且开采时间不长, 故形变量比其他矿区的要大。司马矿区开采后, 也开始有了较大沉降形变。

对于王庄矿区, 根据已有的地表移动测量结果及山西煤矿开采的实践经验和观测资料分析结果, 可知区内采空塌陷地表移动正处于一个残余移动变形过程中, 最大残余移动变形期达44 a, 最晚残余移动变形终结年度预计至2040年; 2002年以后最大剩余下沉量为175~439 mm, 年平均极限下沉量为7~12 mm^[14]。本文监测到的王庄矿区最大沉降速率为7~10 mm/a, 符合残余移动变形规律。

在地形平坦地区, 当开采深度大、黄土层厚时, 沉陷盆地内就会出现积水。潞安矿区内的积水主要有季节性积水和外来水, 前者积水来源主要是大气降水, 如煤矿的季节性积水塌陷地; 后者则主要是排水管管道遭受破坏发生泄漏, 水流注入沉陷盆地所形成的积水。本文实验发现, 地表形变趋势有时出现周期性特点, 某些时刻矿区会出现少量抬升形变, 分析认为, 这可能与矿区积水有关。

1)本文采用DInSAR方法对山西省长治地区一些矿区发生的快速较大形变进行了监测分析, 矿区最大形变量在30 d内达11 cm; 利用SBAS-InSAR方法监测长治矿区微小缓慢形变, 获取了2003年7月— 2010年7月期间研究区地面沉降的空间展布及时间序列相对形变量。分析认为, 该区整体形变表现为沉降趋势, 沉降面积较大, 沉降速率为5~15 mm/a, 最大累计沉降为90 mm。

2)对发生形变较大的余吾矿、常村矿及三元矿3个矿区周围居民区区域进行了详细的时间序列形变分析, 结果发现, 矿区沉降常表现出非线性特征, 地表形变趋势有时出现周期性特点。不同矿区因开采时间、开采方式、采储量及地形等因素的不同而呈现出不同的沉降结果。

3)本文对王庄矿区残余沉降的监测结果与前人的研究分析结果相一致; 对国道形变的监测也与野外调查报告和已有报道资料较为一致。

需要说明的是, 由于研究区无2007年5月— 2009年12月期间的存档数据, 给时序形变监测带来了一定挑战。但是在此期间, 王庄、五阳及漳村等老矿区处于残余形变期, 司马、余吾等新矿区处于刚开始形变期, 所以, 应用SBAS方法依然可测整体时间序列形变, 只是未能获取2007年5月— 2009年12月期间时序形变值。另外, 矿区水准资料、GPS监测资料较少, 与InSAR时序监测结果未能进行同步比较验证。然而, 实地调查显示, 本文的结果较好地吻合了矿区开采范围和地表及建筑物的破坏情况。说明本文方法可以很好地揭示老采空区残余沉降形变趋势, 其结果对预防矿山地质灾害具有重要的意义。