基于SBAS技术的金属矿山沉陷规律研究

2015-03-28雷广渊

测绘工程 2015年3期

雷广渊，周辉

（1．广州市城市规划勘测设计研究院，广东广州510060；2．广东交通职业技术学院公路学院，广东广州510650）

矿产资源是国民经济中重要组成部分，也是工业发展的重要物质保障。矿产资源的开采造成了大量地表塌陷、建筑物变形、道路破坏等，严重危害人类生命财产安全。目前煤炭开采的沉陷规律的研究已有显著成果，但是金属矿山开采沉陷研究国内外仍然处于探索阶段，尚且没有形成成熟的理论，主要原因在于金属矿山地质采矿条件与煤炭相比要复杂很多。研究金属矿山开采沉陷规律可以为指导和管理矿山开采提供依据。

地下有用矿物被开采之后，原有的应力平衡被破坏，岩层移动波及到地表，损害矿区地表的建筑物、铁路等。文献［1］对构造应力下的地表移动规律做了研究，总结出了急倾斜崩落法的地表移动规律。文献［2－3］对铁矿的地表移动规律做了探索性研究，总结出崩落法开采铁矿石可采用概率积分法实现，并将其应用在程潮、北洺河铁矿中。但是目前对矿区开采沉陷的研究都是基于传统测量方式为基础进行分析研究，本文将采用小基线集（Smallbaseline Subset，SBAS）技术对某金属矿山的沉陷规律进行探索研究。针对复杂的时空沉陷过程，研究其沉陷规律和下沉机理，SBAS技术避免了基线的限制，可以生成时序的地表形变。本文针对湘中地区某大型的锑金属矿山，利用SBAS技术综合分析该金属矿山的岩层与地表移动规律，指导矿山开采并作为类似金属矿山的工程类比依据。

1 SBAS技术原理

SBAS技术将影像通过基线的长短分成多个小集合，对于每个基线差异较小的集合，利用最小二乘方法计算，得到地表形变时间序列。在解算的过程中，由于影像之间存在相关性以及小基线集的不同组合会使方程产生秩亏，对于秩亏问题，一般采用奇异值分解附加最小范数的方法解算。

假设针对N＋1幅影像，时间向量t＝［t0，…，tN］T，选取合适的基线进行干涉，生成G个差分干涉图，其中

忽略大气误差、去相关以及高程误差的影响，对于时间tA和tB模型简化为

其中：λ为雷达波长，d（t，i，j）为t时刻雷达视线方向相对参考时刻的形变。

SBAS是按像元逐个计算时间序列的形变，所以不失一般性，以某一像元建立关系。设SAR影像中该像元的相位表示为 φT＝［φ（t0），φ（t1），…，φ（tN）］，为待求量，相干图中的观测相位表示为δφT＝［δφ1，δφ2，…，δφM］，为观测量。其中φT表示像元在不同时刻的相位，δφT表示干涉图中相对于主参考点解缠后的相位。若对主辅影像按时间进行排序，则主影像的时间大于辅影像的时间，对于干涉图中的相位可表示为

δφi＝φ（tMaster（i））－φ（tSlave（i）），i＝1，2，…，M．（3）用矩阵表示为δφ＝Aφ

A矩阵每一行有两个非零元素（1和－1），代表一个干涉图，矩阵的列代表对应时间的一景SAR影像。若M＞N，可以利用最小二乘方法获得唯一解，若M与N接近或者小于N，需采用奇异值分解以及最小范数求解。公式可表示为

带入公式δφ＝Aφ，可得

表示第j幅干涉图的相位值可以用两幅主辅影像的获取时间区间内各时间段的沉降速率的积分求得，写成矩阵表示为

其中B为M×N的矩阵，在主辅影像获取时间之间，B的系数不全为0，否则B＝0。与式（5）同样需要奇异值分解，得到速度的最小范数解。时间段上的形变可以利用速度在时间内的积分计算。

2 工程概况

研究矿山位于湘中地区，该矿山属于大型锑矿，锑储量非常丰富。矿山分南北两个矿区，本文针对北矿区的某矿床进行分析。该矿区地质构造复杂，处于多种构造体系的交界部位，而且受到不同性质的构造应力作用。矿体主要受背斜和F3断层控制。矿田赋存Ⅰ、Ⅱ、Ⅲ号3个矿体，Ⅰ、Ⅱ号为主要开采矿体。矿体呈不规则的似层状和透镜状，倾角15°左右，矿体厚度为0．56～8．81 m，平均2．43 m。顶板岩层为页岩、硅化灰岩、破碎的灰、页岩。主要采矿方法有普通房柱法和胶结充填法。目前矿床共有8个中段，每个中段间隔距离30 m，一中段标高为＋567 m，八中段标高为＋357 m。四中段以上矿体已开采完毕，但是对于部分矿体残留，仍进行残采。五、六中段为主要开采中段开采矿体为Ⅰ、Ⅱ号矿体，八中段目前仍在开拓阶段。

3 实验分析

3．1 实验数据

实验采用的数据包括15景日本空间局ALOS卫星的PALSAR数据以及NASA的SRT M（Shuttle Radar Topography Mission，航天飞机雷达地形测绘任务）数据。ALOS卫星获取SAR影像的时间分别为：2007－01－25、2007－06－12、2007－07－28、2007－12－13、2008－08－14、2008－12－15、2009－01－30、2009－08－02、2009－09－17、2009－12－18、2010－02－02、2010－03－20、2010－06－20、2010－09－20、2011－02－05。SRT M数据在网上公开发布，其DEM水平分辨率为90 m，数据可以从JPL提供的地址：htt p：／／dds．cr．usgs．gov／srt m／下载。覆盖实验区域以及周边区域的DEM 有 9 块，分别为 N26E110、N26E111、N26E112 N27E110 N27E111 N27E112 N28E110、N28E111、N28E112。

3．2 SBAS监测金属矿山

采用瑞士Gamma软件作为实验平台，对15景PALSAR数据使用SBAS技术进行处理。假设2007－01－25作为参考影像，即形变为0，那么可以得到14个时间段的形变序列，地表形变序列如图1所示。

图1 地理编码后形变序列图

由图1可知，影像的中间区域在监测时间段内发生较大的形变，根据资料分析得知，该区域有大量煤矿，煤炭的大量开采导致沉降范围较大，形变也较为显著，在监测周期内发生最大变形约为60 c m。研究区域位于该影像的左下角，该区域在监测周期内发生地表形变达16 c m。该区域为某大型金属矿山，本文主要针对这一局域进行分析，从图像中提取该矿地表沉降序列，见图2。

图2 地理编码后地表形变序列

选取2011－02－05的形变进行分析，如图3所示，该区域为矿床地表的沉降情况，该区域在2007－01－25至 2011－02－05 监测期间地表最大沉降达16 c m，沉降走向为北偏东50°左右。

图3 2011－02－05累积沉降等值线图（单位：m）

3．3 矿山地表移动规律探讨

为了深入了解该矿床地表的沉陷规律，本文对沉降序列进行分析。该矿床的矿体走向为北东东，与地表沉降的发展方向一致。分别选取矿体走向方向和倾向方向的两个断面进行分析，图4、图5分别为矿体走向方向上和倾向方向上的断面沉降序列。

图4 矿体走向方向上的地表下沉曲线

图5 矿体倾向方向上的地表下沉曲线

由图4、图5所示，该矿床的地表塌陷在空间上与高斯曲线相似，这说明在一定范围内，该金属矿床的地表移动规律宏观上符合概率积分法的原理，经实地调查发现，沉降较大的区域主要是受该矿床第三中段开采的影响，四中段和五中段在不同程度上也影响上覆地表。由于无法获得地下完整的采空区数据，本文采用高斯曲线拟合沉降值，从曲线的形态上分析下沉盆地的规律。拟合的结果见图6、图7，使用的高斯拟合函数为

通过曲线拟合求得的倾向方向和走向方向上的高斯曲线结果可知，曲线在宏观上与实际数据拟合效果良好，但在局部区域拟合结果较差，可能由于SBAS观测精度的问题或者地表变化不连续。在曲线的边缘拟合效果较差，这一问题可能的主要原因有：①概率积分法模型误差，较多工程实例表明，概率积分法在边缘拟合效果较差。②实验区域为金属矿山，采矿方法也是采用房柱法等，上覆岩层并不是理想的随机游动模型，部分区域可能存在断层、构造应力等的多重影响。

图6 倾向方向的高斯曲线拟合以及残差分析

走向方向和倾向方向求得的高斯曲线的参数和残差中误差见表1。

图7 走向方向的高斯曲线拟合以及残差分析

表1 高斯拟合函数相关指标

表1中均方根误差

3．4 地质采矿因素与地表下沉联系

矿山开采导致了地表的下沉，大范围的地表下沉不仅对地表建筑物和道路造成损害，而且使得在开采过程中，由于上覆岩层的破裂，坍塌等，造成采场、运输巷道、通风系统等受到影响，无法进行安全生产，从而危害人类的生命财产安全。为了深入了解矿床的地表移动规律，进一步分析矿山开采与地表下沉之间的联系。图8反映了该矿床35号勘探线3个矿体空间分布特征以及相应的地质条件。由图中可知，该矿床主要受背斜以及F3（右第3条断层）断层控制，即控矿因素为“背斜加一刀”，I和II号矿体为主要开采矿体。

为了分析矿体开采与地表沉降的因果关系，将InSAR获取的沉降场与地形图配准，配准结果见图9。

图8 35勘探线剖面图

图9 InSAR沉降值与矿区开采图配准结果

分析图中A和B两个区域的变形特点，图中所示的时间是2011年2月（最后一次观测）的累积形变，从图中可知，在14次的监测中，A区域没有明显的沉降，B区域有较大的变形。为分析A、B区域的形变原因，对以上区域进行的采矿活动进行统计。通过2012采选厂储量计算图、各中段采矿验收报告中得知：

A区域是第五中段，采场编号为5100～5127以及采场5211～5217等，采深大约在150 m左右，在2008年至2009年期间，该区域进行了采矿活动，主要对I号薄矿体和II号厚矿体进行了开采。开采的方法为充填法开采，由于开采离地表较远，而且对大部分采空区进行了管理，地表在此期间没有发生大范围、剧烈的形变。

B区域为34线到80线之间的区域，该区域矿体深度主要分布在90～150 m，即矿体主要分布在三、四、五中段。由于矿区开采的年限较长，三中段和四中段采空区未进行处理，留有较大的采空区。采空区由于仍有矿体遗留（矿柱、边角等矿体），在2007年至2011年期间，该区域三中段和四中段一直存在残采，矿柱和边角矿体被开采之后，对采空区未进行处理，采空区离地表又较为接近，所以直接导致了地表产生了下沉。五中段开采主要由于离地表较远，部分采场也进行了充填，对上覆岩体进行了支撑，故可认为对地表影响较小。所以该区域主要的地表下沉一方面是由于存在大量空区所致，另一方面，矿柱和边角矿体的回收，导致上覆岩体失去支撑，重复采动也对地表影响加剧，最终导致了该区域地表下沉较大。表2为对该矿床空区的统计，从表中可知，三中段采空区体积最大，四中段次之，从图9中可知，下沉较大的区域主要范围是三中段的巷道之上，这与地表下沉较大的区域较为符合。

表2 矿区空区统计

为更加清楚地认识矿山开采与地表下沉的联系，本文对66号勘探线上的沉降值进行分析。66号勘探线位于矿床的中上部，在勘探线上每10 m取一个下沉点，下沉序列见图10。64～70勘探线在2008年之后对矿床进行残采工作，15号点附近为下沉较为剧烈的区域，该区域为三中段采场249～259残采范围。随着点号的增加，地表下沉变缓，点号23附近，在2008－2010年期间五中段5217～5243采场进行了采矿活动，由于离地表相对较远，地表变化变缓。

图10 66号勘探线地表下沉序列

倾斜曲线表示地表的倾斜以及地表的下沉速度大小，在获得了66号勘探线上沉降曲线之后，利用下沉值计算倾斜曲线，对最后一期数据进行计算获得的倾斜曲线如图11所示。由于影像受到大气、失相关、影像分辨率等因素影响，虽然下沉曲线在趋势上呈U型，但是在局部区域，相邻地表沉降值出现相等或者相近，所以在倾斜计算中，使得曲线出现跳动，与实际的物理意义不符合。为了削弱这种由于上述因素影响下的相邻沉降值出现相等或者相近而引起的倾斜曲线跳动，本文使用一种滤波去噪算法，既要去掉大部分干扰噪声，又不失原始下沉曲线的沉降特点，并且保证修正后的数据与原始数据差距在一定范围内，拟合的数据能够反应地表下沉的主要特征。本文采用滑动平均对原始数据进行平滑，选用的窗口为5，目的是消除相邻两个点沉降值相等或者相近以及噪声等因素的影响，实验对原始数据分别进行了一次、二次、三次平滑，结果见图11、图12。

图11 原始沉降进行修正

图12 修正后的倾斜曲线

平滑3次后，结果见图13、图14。

图13 利用滑动平均修正沉降数据图

图14 利用滑动平均修正沉降数据后倾斜曲线

由图13、图14知：地表最大下沉值在点号15附近，此处计算的倾斜接近0，倾斜在15号点之前，有一个增大和下降的过程，与开采沉陷的理论中倾斜曲线和下沉曲线的关系符合。点号15之后倾斜仍然在下降，下降达到最小（负极值）后开始递增，这里倾斜并没有继续增大，形成一个对称的曲线，原因是由于20点之后，矿山开采的深度和采矿方法以及顶板管理方法等因素的制约，导致曲线变缓。