CN112184902A - 一种面向越界开采识别的地下开采面反演方法 - Google Patents

Abstract

本发明公开了一种面向越界开采识别的地下开采面反演方法，在描述和构建矿山地表开采沉陷监测模型的主要类及其相互关系的基础上，结合InSAR技术和开采沉陷预计方法，提出了根据地表形变信息反演地下地下开采范围和深度的方法，通过揭示以采动形变为核心的地表采动信息与地下开采区域的关联机理，可根据获得地表采动信息、已知或可掌握的影响因素信息，推演、发现深藏在地表以下的非法开采区域。

Description

一种面向越界开采识别的地下开采面反演方法

技术领域

本发明涉及越界开采识别方法，尤其涉及一种面向越界开采识别的地下开采面反演方法。

背景技术

矿区越界开采是一种普遍存在的违法行为，往往隐蔽性较强且容易造成重大安全事故，矿区的地表沉降在一定程度上影响正常的采矿秩序，当沉降十分严重时，会影响采矿活动的安全进行。为了矿区的安全生产，为了矿区的合法利益，监管部门非常需要具有高精度且监测时间短的设备用以对矿区进行实时地监测，以便及时有效地发现越界开采行为。

传统识别非法开采手段主要是实地调研，这往往受到矿区相关人员的阻挠，且花费大。因此如何能够高效、快速、自动识別越界开采区是遏制非法开采问题的关键。

确定地下开采范围和深度是进行地下越界开采识别的关健环节，在不能进入井下的条件下快速、自动圈定地下开采范围的现实选择就是利用地表采动信息，随着各类SAR数据源类型的多样化以及InSAR技术的迅猛发展，使得InSAR技术具备了大范围监测矿区地表形变的能力，但为了推断出地下开采范围，仅掌握地表形变信息还是不够的，还需建立地表形变信息与地下开采区域的映射关系，找出影响这些关系的主控因素，如开采时间、开采厚度与深度、矿层倾角、推进速度、上覆岩层特性、松散层厚度、开采及顶板管理方法、地形地貌等等资料，这些信息对开发监督管理等部门而言是很容易获取的。通过揭示以采动形变为核心的地表采动信息与地下开采区域的关联机理，可根据获得地表采动信息、已知或可掌握的影响因素信息，推演、发现深藏在地表以下的非法开采区域。

发明内容

针对上述存在的问题，本发明旨在描述和构建矿山地表开采沉陷监测模型及其相互关系的基础上，结合InSAR技术和开采沉陷预计方法，提出了根据地表形变信息反演地下地下开采范围和深度的方法。

为了实现上述目的，本发明所采用的技术方案如下：

一种面向越界开采识别的地下开采面反演方法，其特征在于，包括以下步骤：

步骤S1：通过矿山地表工作面上方的时序SAR影像，利用D-InSAR技术获取时序相邻两景SAR数据间开采沉陷盆地的地表形变信息，并将获取到的下沉图依次叠加，得到地表移动盆地的沉降图。

步骤S2：对获取的沉降图进行处理并提取出地表形变信息和开采沉陷特征，得到地表移动盆地主断面上的地表最大下沉点、盆地边界点参数；

步骤S3：根据获得的地表移动盆地主断面上的参数，根据随机介质理论得到地表下沉曲线与曲率曲线；

步骤S4：根据地表下沉曲线和曲率曲线计算得出边界点和拐点以及拐拐点偏移距和边界角；

步骤S5：根据所述边界点和拐点得出采空区范围，并且根据拐点偏移距和边界角进行反演采深；

步骤S6：将计算出的采空区和采深进行空间上的叠加反演出地下开采面，并基于采矿权范围数据判断是否存在越界开采行为。

进一步地，步骤S2中提取出地表形变信息的具体步骤包括：

S21：采用重复轨道干涉测量模式，从获取的干涉相位图中定义干涉相位的综合贡献的组成为：

其中ω表示缠绕算子，

为卫星视线方向的地表形变相位；

为参考面相位，

为地形相位，

为大气相位，

为噪声相位。通过两幅或多幅干涉影像进行差分处理去除平地效应，并逐一将上式中后四种相位消除，即可分离出地表形变信息；

S22：通过对基线进行精确估算去除参考面相位，通过模拟外部DEM去除地形相位，通过提高信噪比来降低大气相位，通过低通滤波的方式抑制噪声相位；

S23：最终分离出地表形变信息。

进一步地，步骤S2中从干涉相位图中得出的开采沉陷特征包括空间特征、几何特征以及形变特征；

空间特征是地下开采会引起相应的采矿活动区域上方出现地表下沉现象，最大的沉降量主要发生在采矿区域的地表中心，从中心到边缘下沉幅度逐渐减小，最终在该采矿区域表面形成一个空间漏斗，则漏斗中各点的沉降量为：

f_z(x_n,y_n)＜f_z(x_m,y_m)＜0且df_z(x,y)/dr＞0 (2)，

其中，f_z(x,y)是下沉区z的表达式，(x_n,y_n)和(x_m,y_m)代表下沉范围内的点，且(x_m,y_m)离下沉中心更近；

几何特征基于空间特征得出，即所述空间漏斗的沉降中心由四周的斜坡包围，形变中心点的梯度绝对值大于没有形变的区域的梯度绝对值，梯度的大小代表地表下沉的幅度信息，方向代表的是相位信息，且相位信息的计算公式为：

形变特征为地下开采引起的底板沉陷区域，在干涉相位图中通常为一系列闭合的圆形，这些圆被近似为干涉图坐标系上的一组小椭圆，表示为：

其中，∑(x_n,y_n)＜∑_t，且n＝1,…,N；x和y是在轮廓中分别沿着距离向和方位向上点的坐标；x_n和y_n是椭球体的中心点坐标；a和b分别是长半轴和短半轴；θ为主轴上的倾角；∑(·)代表这是椭球体的大小，且都小于常数∑_t。

进一步地，步骤S3中所述的地表下沉曲线表示开采活动导致的地表下沉的分布规律，即在地表最大下沉点O处下沉值最大，自盆地中心至盆地边缘地表下沉值逐渐减小，在盆地边界点A和B处下沉值为零，且地表下沉曲线以采区中央对称。

进一步地，步骤S3中所述的曲率曲线是表示地表移动盆地内曲率的变化规律，曲率曲线可表示为下沉的二阶导数，曲率曲线的分布规律表现为：盆地边界点A和B处和拐点E处曲率为零，且盆地边缘区为正曲率区，盆地中部为负曲率区。

进一步地，步骤S5中所述的计算采空区范围的具体操作步骤包括：

S51:在三维空间坐标系中，基于概率积分法得出单元开采引起地表任意点的下沉最终值公式为：

W_g(x，y)＝1/r²·exp(-π(x-x_i)/r²)·exp(-π(y-y_i+l_i)/r²) (5)，

其中，r＝H₀/tgβ，l_i＝H_iCtgθ，r为主要影响半径，H₀为平均采深，β为主要影响角，l_i为采动影响长度，H_i为单位的深度，C为采动影响系数，_θ为最大下沉角，B(x_i，y_i)为单元中心点的平面坐标，A(x，y)为地表上任意一点的坐标；

S52：设工作面的开采范围为0：D₁和0：D₂组成的矩形采空区，即工作面走向方向的开采长度为D₁，沿工作面倾向方向的开采宽度为D₂，则整个工作面开采引起地表任意点下沉的概率积分法计算公式为：

其中，W₀为该地质采矿条件下的最大下沉值，且W₀＝mqcosα，m为煤层开采厚度，q为下沉系统，α为煤层倾角；

公式(3)可表示为：

其中，W₀为走向和倾向均达到充分采动时地表的最大下沉值，W⁰(x)为倾向方向达到充分采动时走向主断面横坐标为x的点的下沉值，W⁰(y)为走向方向达到充分采动时倾向主断面横坐标为y的点的下沉值；

S53：根据获得的地表下沉的公式，计算出地表上任意一点A(x，y)的倾斜、曲率、水平移动、水平变形值。

进一步地，步骤S53中所述的任意一点的倾斜、曲率、水平移动、水平变形值具体计算公式分别为：

a.倾斜：对于坐标为(x,y)的点沿

方向的倾斜i

为下沉W(x,y)在

方向上单位距离的变化率，其计算公式为：

公式(5)可简化为：

b.曲率：对于坐标为(x,y)的点沿

方向的曲率k

为倾斜i

在

方向上单位距离的变化率，其计算公式为：

公式(6)可简化为

c.水平移动：对于坐标为(x,y)的点沿

方向的水平移动U

的计算公式为：

d.水平变形：对于坐标为(x,y)的点沿

方向的水平变形ε

的计算公式为：

其中，i⁰(x)、k⁰(x)、U⁰(x)、ε⁰(x)分别别表示倾向方向达到充分采动时走向主断面上横坐标为x的点的倾斜、曲率、水平移动和水平变形值，i⁰(y)、k⁰(y)、U⁰(y)，、ε⁰(y)分别表示走向方向达到充分采动时倾向主断面上横坐标为y的点的倾斜、曲率、水平移动和水平变形值。

进一步地，步骤S5中所述的计算采深的具体操作步骤包括：

S51：根据经验值，给出边界角δ₀的初始值，并且边界角δ₀取值和岩层性质的关系为：

δ₀＝g(k) (14)，

其中，k＝f(k₁，k₂,k₃，L,k_n)，(k1，k₂，k₃，L，k_n)k为地下工程上覆岩层岩性系数；

S52：根据覆岩层岩性系数k得出采深H的计算公式为：

H＝h(k) (15)，

S53：根据采深H得知覆岩岩性，由岩性计算k，并根据公式(14)得出边界角δ₀；

S54：判断采深H、边界角δ₀、下沉盆地边界点到拐点的距离D以及拐点偏移距D₀是否满足以下公式：

其中r为主要影响半径，

如果满足公式(16)则得到采深H，如果不满足公式(16)则重新执行步骤S51-53。

本发明的有益效果是：

本发明通过介绍开采沉陷基本规律和开采沉陷预计方法，说明了概率积分法的基础理论和基本原理，探析了边界角与覆(围)岩岩性、开挖深度等地质环境因素存在的函数关系，并以此反演出地下开釆平面范围以及深度等参数，该方法为矿区地下越界非法开采识别提供了一种高效低成本的识别方法，也为采空区定位提供了新的思路。

附图说明

图1基于InSAR和概率积分法越界开采识别模型主要对象的交互图；

图2为本发明中的地表移动盆地的形成过程；

图3为本发明中的动态和静态移动盆地示意图；

图4为本发明中的开采沉陷规律示意图；

图5为本发明中的颗粒体介质的理论模型；

图6为本发明中的开采三维空间坐标系；

图7为本发明中的InSAR相位提取示意图；

图8为本发明中的不同结构的地下开釆引起的地表形变示意图；

图9为本发明中的地下工作面和地表形变的关系示意图；

图10为本发明中的地表移动盆地下沉曲线与曲率曲线分布规律图；

图11为本发明中的地表移动盆地动态下沉示意图；

图12为本发明的面向越界开采识别的地下开采面反演流程图；

图13为实施例中模拟下沉盆地，模拟采空区及反算采空区位置对比图；

图14为实施例中模拟下沉盆地在倾向方向的倾斜与曲率。

具体实施方式

为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。

为了对矿山地表的开采沉陷信息以及对地下开采面的反演进行建模与应用，因而将矿山地表抽象为地质对象，监测矿山地表形变的星载SAR传感器和地下开采面分别抽象成相应的图层对象。星载SAR传感器观测到地表的沉陷信息、地下开采面的采掘进度平面图以及反演出的地下开采面范围都作为矿山地表对象的状态。矿山地表的实时数据来源主要依赖于对地观测的星载SAR传感器，且矿山地表对象是衔接矿山地下开采程度与星载SAR传感器的重要纽带，它是将星载SAR传感器数据和地下开采面接入到矿山开采时空过程的必要操作。当地下有新的开采事件或者开采面推进事件时，开采面对象则会根据地下开采位置和范围，生成对应等级的地质事件，并且将生成的地质事件发送到开采沉陷的地质时空过程中，该时空过程又将此事件发送到影响区域内的矿山地表，矿山地表根据InSAR技术监测地表开采沉陷信息的约束规则和接收到的矿权范围对象信息，决定是否响应地下越界开采事件的驱动，基于InSAR和概率积分法越界开采识别模型主要地质对象和事件的交互作用如附图1所示。

由附图1可知，矿山地表对象在地表形变状态满足相应条件时才产生对应级别的事件，随着地下开采面在不断移动和推进，地质时空过程需不断地通知区域内矿山地表所发生的地表形变移动事件，并生成开采沉陷事件。矿山地表通过接收到InSAR技术获取的形变信息后，能动态地提取出各开采沉陷事件的沉陷信息。根据开采沉陷信息，结合几何沉陷预计理论和概率积分法，反演出地下开采面的范围，再通过对照井上下对照图采掘进度情况，可实时地掌握地下开采事件是否严格按照计划的采掘进度安排来推进。故当矿权范围对象接收到反演出的开采面事件后，就可判断其是否存在越界开采事件，如果越界开采事件得以确认后，则对该事件进行预警预报响应。综述所述，在搭建的地下非法采矿识别平台上，要实现基于InSAR和概率积分法越界开采事件的快速识别，还需解决好InSAR矿区地表形变监测、沉陷信息提取和地下开采面反演等关键问题。

当地下局部矿体被采出以后，采区岩体内部形成一个采空区，在采空区上覆岩层重力的作用下，导致周围岩体应力平衡状态受到了破坏，从而使岩体内部产生移动变形和破坏，直至岩体达到新的平衡。当开采工作面自开切眼开始向前推进的距离约为平均采深的0.25～0.5倍时，岩层移动即发展到地表，使地表产生移动和变形，并引起地表下沉。随着地下采空区范围的继续推进，地表受影响范围也不断扩大，逐渐在地表形成一个比地下开采范围大得多的下沉盆地。

参考附图2可以看出地表移动盆地的形成过程，当工作面由开切眼推进到位置1时，在地表形成了一个小盆地；工作面继续推荐到位置2时，在移动盆地的范围内，地表继续下沉，同时在工作面前方原来尚未移动区域的地表点，先后进入移动，从而使下沉盆地扩大而形成移动盆地，这种移动盆地是在工作面推进过程中形成的，故称为动态移动盆地，图2中的W₁、W₂、W₃和W₄为形成的动态移动盆地。当工作面回采结束后，地步移动不会立即停止，移动盆地的边界还将继续向工作面方向扩展，移动首先在开切眼一侧稳定，而后在停采线一侧逐渐形成最终的移动盆地，又称为静态移动盆地。在工作面推进过程中，如附图3所示的工作面停在1、2、3、4的位置上，待地表稳定后，其对应的每一个位置都会有一个静态的移动盆地，而W₀₁、W₀₂、W₀₃和W₀₄为地表移动稳定后最终形成的静态移动盆地。在地表移动动态发展过程中，开采工作面后方的地表点仍在继续移动，但其移动的激烈程度会随着工作面的远离而逐渐减弱，直至稳定，一般最先进入移动的点，也最先达到稳定。

开采沉陷规律就是探究由地下开采引起地表形变值和空间分布特征及其与地质采矿条件之间的关联机理。参考附图3可以看出在通常的地质采矿条件下，由单一地下开采面引起地表移动盆地的特征。对于水平煤层，地表移动盆地位于采空区的正上方，它的形状与采空区对称，且移动盆地内外边缘区的分界点，大致位于采空区边界的正上方。对于倾斜煤层，在倾斜方向上，移动盆地的最大下沉点偏向采空区的下山方向，盆地与采空区的相对位置在倾斜方向上不对称，且盆地的上山方向较陡，移动范围小，而下山方向较平缓，移动范围大。

再参考附图4的开采沉陷规律示意图可以看出其中的红色曲线、黄色曲线、蓝色曲线、绿色曲线和紫色曲线分别表示了地表移动盆地内垂直向下沉W(x)、水平移动U(x)、倾斜i(x)、曲率K(c)、水平变形ε(x)五项指标的变化规律。

对于水平煤层，沉陷最大下沉值在盆地中央，且下沉曲线凹凸分界的拐点处，下沉值约为最大值的一半；边界点和采空区中点的水平移动为零，边界点和采空区中点之间有极值；盆地边界点至拐点间倾斜渐增，拐点至最大下沉点间倾斜渐减，在最大下沉点处倾斜为零，且倾斜和水平移动的曲线变化趋势同步；盆地边缘区为正曲率，盆地中部为负曲率，且盆地边界点和拐点处的曲率为零；盆地边缘区为拉伸区，盆地中部为压缩区，且盆地边界点和拐点处的水平变形为零。对于倾斜煤层，下沉曲线将失去对称性，随着煤层倾角的增大，指向上山方向的水平移动值逐渐增大，而指向下山方向的水平移动值逐渐减小；最大拉伸变形在下山方面，最大压缩变形在上山方向，水平变形为零的点与最大水平移动点重合，且水平移动曲线和倾斜曲线，水平变形曲线和曲率曲线变化趋势不再相似。

从19世纪工业革命开始，人类对各种能源利用的深度和广度在不断增长，地面表层煤矿的开采已经无法满足社会需求，煤炭资源的开采不断向地下延伸，地下采动容易造成地面建筑物、自然景观和生态环境的破坏。对于一个正在计划进行的地下开采，在进行开采之前，可以根据相关地质采矿条件和选用的预计函数、参数，预先计算出岩层和地表的形变，因此，利用开采沉陷预计方法能为合理安排矿区的生产和生态环境的保护提供技术支持。开采沉陷预计的方法有很多种，我国采用较多的主要有剖面函数法、典型曲线法和影响函数法。

1、剖面函数法：利用某些函数来表示各种开采条件下地表下沉盆地主断面内移动，我国常用的剖面函数有皮尔森函数、双曲线函数和负指数函数等，这些函数是典型曲线的解析表示形式。剖面函数可以根据实测资料快速确定函数的参数值，比较直观，便于数学处理和计算机解算。由于剖面函数和实际下沉盆地形状有一定差异，特别是在地表形变特征点处差异明显，故求得的参数与实际值会存在一定的出入。但对于相同地质条件下，剖析函数法比较适合矩形工作面开采的地表形变预计。

2、典型曲线法：通过对研究区域设立观测站，观测工作面主断面上方的地表移动情况，把地表形变观测值绘制成无因次曲线，通过无因次曲线来描述地表的下沉、倾斜、曲率、水平移动和水平变形的分布规律，适用于规则采空区上方的地表移动和变形预计，在我国使用典型曲线法较多的是峰峰矿区和平顶山矿区。典型曲线参数的确实是直接根据实测资料，预计精度较高，但由于典型曲线是针对特定矿区来建立的，其他矿区由于地质条件和实测资料的不同，无法套用该方法，故局限性大，难以大范围地进行推广和使用。

3、影响函数法：该方法是介于经验法和理论模型法之间，把井下煤层划分为若干个小单元，先计算每个小单位对地表形变的影响，然后再叠加所有小单元对地表的影响值，叠加后的总值就是地表形变值。影响函数法符合先微分再积分的思维过程，与剖面函数相比，影响函数法能推广到三维空间。它是由德国学者凯因霍尔斯提出，源于波兰学者J.Litwiniszyn创立的岩层移动随机介质理论，在此基础上，我国学者刘宝琛将其实用化，发展为概率积分法。概率积分法将采空区上方覆岩岩石看成松散介质，把地表移动看成服从统计规律的随机过程，符合影响函数的叠加原理。概率积分法自引入国内以来，众多学者对其开展研究，并提出了适应中国地质采矿条件的修正模型，由于概率积分法偏于计算机编程实现，在我国开采沉陷预计中得到广泛的应用，也是本发明开展地下开采面反演研究所采用的方法。

在开采沉陷理论研究中，常用连续介质模型和非连续介质模型来模拟岩体，而开采引起的岩体和地表移动规律与作为非连续介质的颗粒体介质模型所描述的运动规律在宏观上相似，概率积分法就事基于非连续介质理论。颗粒体介质的理论模型如图附图5所示，在理论模型中，假设介质颗粒为一些大小相同、质量均一的小球，并被装在大小相同、排列均匀的方格内，若在岩体中，这些方格即为地表至地下开采面的矿体空间。若第1层方格中的小球(煤矿)被移走后，由于重力作用，第2层上的两个相邻方格中的小球滚入第1层方格的概率应均是1/2。由此向上每一层类推，就可以得到附图5(b)中的颗粒移动概率分布图，从中选取一个坐标系，则介质内任意一个水平的概率分布可以绘制成图5(c)所示的概率分布直方图(虚线)。若方格和颗粒无限小，则该直方图就趋近于一条光滑的曲线，如图5(c)中的实线所示。

根据数学推导，可得到单元开采引起下沉盆地的表达式为：

其中，A为反映方格尺寸的一个常数，r_z为主要影响半径。

单元开采引起地表水平移动可表示如下：

对于等深开采，r_z和k_z均为常数，则公式(1)可简化为：

同理，公式(3)可简化为：

本发明提出的一种面向越界开采识别的地下开采面反演方法，包括以下步骤：

步骤S1：获取矿山地表工作面上方的同一地区两幅InSAR影像，将其中一幅作为主影像，另一幅作为辅影像，并且将两幅InSAR影像对应像素的相位值相减得到干涉相位图；

影像中的每一个像素不但包含了地面分辨元的雷达后向散射强度信息，还包含了传感器到目标距离间有关的相位信息，如附图7所示，星载InSAR提取的这些信息包含了大气延迟、地形起伏、参考椭球面以及地表形变等因素的影响；

步骤S2：对获取的干涉相位图中进行处理并提取出地表形变信息和开采沉陷特征，得到地表移动盆地和其主断面的地表最大下沉点、盆地边界点参数；

步骤S3：根据获得的地表移动盆地主断面的参数，根据随机介质理论得到地表下沉曲线与曲率曲线；

步骤S4：根据地表下沉曲线和曲率曲线计算得出边界点和拐点以及拐拐点偏移距和边界角；

步骤S5：根据所述边界点和拐点得出采空区范围，并且根据拐点偏移距和边界角进行反演采深；

步骤S6：将计算出的采空区和采深进行空间上的叠加反演出地下开采面，判断反演出来的采空区范围是否大于地下采矿权范围来判断是否存在越界开采行为；

可以借助于矿产管理部门采矿批准和采矿计划的详细信息，如合法采矿权、采深和开采时间表等数据，结合反演出的开采区范围，可初步确定越界开采的事件，这个非法采矿识别的过程可以用以下的一个数据集来判断，该数据集表示为：

其中，D_x为沿工作面走向方向的开采长度，D_y为沿工作面倾向方向的开采宽度，H为采深，T为开采时间，数据集g_n(D_x，D_Y，H，T)为采矿权边界范围和采深等数据。

进一步地，为了从获取的干涉相位图中提取出形变信息，需从中去除大气延迟等其他因素引起的相位值，从而步骤S2中提取出地表形变信息的具体步骤包括：

S21：采用重复轨道干涉测量模式，从获取的干涉相位图中定义干涉相位的综合贡献的组成为：

其中ω表示缠绕算子，

为卫星视线方向的地表形变相位；

为参考面相位，

为地形相位，

为大气相位，

为噪声相位。通过两幅或多幅干涉影像进行差分处理去除平地效应，并逐一将上式中后四种相位消除，即可分离出地表形变信息；

S22：通过对基线进行精确估算去除参考面相位，通过模拟外部DEM去除地形相位，通过提高信噪比来降低大气相位，通过低通滤波的方式抑制噪声相位；

S23：最终分离出地表形变信息。

进一步地，步骤S2中从干涉相位图中得出的开采沉陷特征包括空间特征、几何特征以及形变特征；

在一幅差分干涉图中，地下开采引起地表沉陷区域具有一些唯其独有的典型特征，总结这些特征对于通过InSAR技术获取地表形变值，并以此计算开采沉陷盆地的模型参数，反演出地下采空范围来说有着至关重要的意义，其中：

空间特征，即地表的最大沉降量主要发生在地下采矿区所对应的地表中心，从中心到边缘下沉幅度逐渐减小，最终在该区域表面形成一个沉降漏斗，定义一个指向漏斗中心的径向矢量，则漏斗中各点的沉降量为：

f_z(x_n,y_n)＜f_z(x_m,y_m)＜0且df_z(x,y)/dr＞0，

式中，f_z(x,y)是下沉区z的表达式，(x_n,y_n)和(x_m,y_m)代表下沉范围内的点，但(x_m,y_m)离下沉中心更近；

几何特征，是由第一个特征推出的，这个明显的特点就是，由于地下开采引起的地表沉陷通常会出现一个沉降漏斗，所以沉陷中心应该是由四周的斜坡包围着的。它也意味着形变中心点的梯度绝对值大于没有形变的区域的梯度绝对值。另一方面，梯度方向大约呈现出由边缘指向沉降中心的反向扩展格局。关于梯度，二维梯度算子是个复杂的表达式，在此用来表达煤矿引起的地表沉陷情况，在梯度算子中，梯度的大小代表地表下沉的幅度信息，方向代表的是相位信息，公式如下：

形变特征，就是地下开采引起的地表沉陷区域，在差分干涉图中一般是以典型的圆形或椭圆形的形状来呈现的。由于陆地表面一般可以被认为是弹性的，也意味着开采沉陷是由中心向外稳定传播的，因此，下沉边缘轮廓在雷达视线方向也可以看作是一系列闭合的圆形。从数学上讲，这些圆可以近似为干涉图坐标系上的一组小椭圆，可以用下式来表示：

式中，∑(x_n,y_n)＜∑_t，且n＝1,…,N；x和y是在轮廓中分别沿着距离向和方位向上点的坐标；x_n和y_n是椭球体的中心点坐标；a和b分别是长半轴和短半轴；θ为主轴上的倾角；∑(·)代表这是椭球体的大小，且都小于常数∑_t；

对地表采动信息与地下开采区域的时空关系研究较多的当属地表形变(沉陷)，开采形变(沉陷)的分布规律取决于地质和采矿因素的综合影响，而不同结构的地下开釆引起的地表形变分布和特征也判然不同，如附图8所示，理论与实践表明地下资源开采诱发的形变(沉陷)是一个时间和空间过程，随着工作面的推进，不同时间的回采工作面与地表点的相对位置不同，开采对地表点的影响也不同，地下开釆诱发的地表形变特征主要表现在以下两个方面：一是地表形变范围与地下开釆范围存在定量关系，据沉陷学理论，该量值随埋深一定幅度波动；二是地表形变轮廓与地下开釆轮廓存在上下空间对应关系。根据以上这两种对应关系，结合深部岩土力学和开采沉陷学等理论知识和相关技术支撑，使得根据地表形变信息和特征反演出地下开釆面的范围、规模、轮廓和走向成为可能。

在主断面上，根据随机介质理论，地下水平煤层开采诱发的地表下沉曲线与曲率曲线分布规律如附图10所示，图中红色曲线为下沉曲线，蓝色曲线为曲率曲线，O为最大下沉点，在水平煤层开采时，在采区中央正上方；A和B为盆地边界点，此处地表下沉值为零；δ₀为边界角，即盆地边界点A和B至采空区边界的连线与水平线在矿柱一侧的夹角；E为拐点，即下沉曲线凹凸的分界点，在地表达充分采动条件下，拐点处的下沉值约为最大下沉值W_m的一半，E₀为曲率曲线的零点，H为采深。

进一步地，步骤S3中所述的地表下沉曲线表示开采活动导致的地表下沉的分布规律，即在地表最大下沉点O处下沉值最大，自盆地中心至盆地边缘地表下沉值逐渐减小，在盆地边界点A和B处下沉值为零，且地表下沉曲线以采区中央对称。

进一步地，步骤S3中所述的曲率曲线是表示地表移动盆地内曲率的变化规律，曲率曲线可表示为下沉的二阶导数，曲率曲线的分布规律表现为：盆地边界点A和B处和拐点E处曲率为零，且盆地边缘区为正曲率区，盆地中部为负曲率区。

附图10中D₀为拐点偏移距，即自下沉曲线拐点在地表面上投影点按影响传播角作直线与煤层相交，该交点与采空区边界沿煤层方向的距离。拐点偏移距的变化规律会受到岩层特性、地质构造和采动程度等众多因素的影响，但在一般情况下，同一开采区域的地质采矿因素、覆岩岩性、采厚、采深等参数相差不会太大，因此，在开采沉陷预计工作中将拐点偏移距认为是一个常数来采用，虽然会带来一些误差，但也易于计算和消除。

进一步地，步骤S5中所述的计算采空区范围的具体操作步骤包括：

S51：在三维空间坐标系中，基于概率积分法得出单元开采引起地表任意点的下沉最终值公式为：

W_g(x，y)＝1/r²·exp(-π(x-x_i)/r²)·exp(-π(y-y_i+l_i)/r²) (2)，

其中，r＝H₀/tgβ，l_i＝H_iCtgθ，r为主要影响半径，H₀为平均采深，β为主要影响角，l_i为采动影响长度，H_i为单位的深度，C为采动影响系数，θ为最大下沉角，B(x_i，y_i)为单元中心点的平面坐标，A(x，y)为地表上任意一点的坐标；

S52：设工作面的开采范围为0：D₁和0：D₂组成的矩形采空区，即工作面走向方向的开采长度为D₁，沿工作面倾向方向的开采宽度为D₂，则整个工作面开采引起地表任意点下沉的概率积分法计算公式为：

其中，W₀为该地质采矿条件下的最大下沉值，且W₀＝mqcosα，m为煤层开采厚度，q为下沉系统，α为煤层倾角；

公式(3)可表示为：

其中，W₀为走向和倾向均达到充分采动时地表的最大下沉值，W⁰(x)为倾向方向达到充分采动时走向主断面横坐标为x的点的下沉值，W⁰(y)为走向方向达到充分采动时倾向主断面横坐标为y的点的下沉值；

S53：根据获得的地表下沉的公式，计算出地表上任意一点A(x，y)的倾斜、曲率、水平移动、水平变形值。

进一步地，步骤S53中所述的任意一点的倾斜、曲率、水平移动、水平变形值具体计算公式分别为：

a.倾斜：对于坐标为(x，y)的点沿

方向的倾斜i

为下沉W(x，y)在

方向上单位距离的变化率，其计算公式为：

公式(5)可简化为：

b.曲率：对于坐标为(x，y)的点沿

方向的曲率k

为倾斜i

在

方向上单位距离的变化率，其计算公式为：

公式(6)可简化为

c.水平移动：对于坐标为(x，y)的点沿

方向的水平移动U

的计算公式为：

d.水平变形：对于坐标为(x，y)的点沿

方向的水平变形ε

的计算公式为：

其中，i⁰(x)、k⁰(x)、U⁰(x)、ε⁰(x)分别别表示倾向方向达到充分采动时走向主断面上横坐标为x的点的倾斜、曲率、水平移动和水平变形值，i⁰(y)、k⁰(y)、U⁰(y)、ε⁰(y)分别表示走向方向达到充分采动时倾向主断面上横坐标为y的点的倾斜、曲率、水平移动和水平变形值。

进一步地，步骤S5中所述的计算采深的具体操作步骤包括：

S51：根据经验值，给出边界角δ₀的初始值，并且边界角δ₀取值和岩层性质的关系为：

δ₀＝g(k) (11)，

其中，k＝f(k₁，k₂，k₃，L，k_n)，(k₁，k₂，k₃，L，k_n)_k为地下工程上覆岩层岩性系数；

由于地表移动过程和相关参数跟地质采矿条件是紧密相关的，尤其是覆岩分层岩性及其赋存条件对各移动参数的影响很大，覆岩分层岩性评价系数作为反映岩性总体状况的重要指标，它的值取决于各分层的岩性组成和厚度，且随采动程度而变化，在我国各矿区的一般情况，不包括特殊地质采矿条件，覆岩分层岩性评价系数的取值如表1所示。

表1覆岩分层岩性评价系数表

而对于地下工程上覆岩层岩性系数，取值方法见表2：

表2地下工程上覆岩层岩性系数取值表

S52：确定了覆岩层岩性，而拐点偏移距D₀可在《三下采煤规程》中寻找获取，下沉盆地边界点到拐点的距离D可根据下沉盆地求得，根据以上数据，便可确定采深H的计算公式为：

H＝h(k) (12)，

S53：根据采深H得知覆岩岩性，由岩性计算k，并根据公式(11)得出边界角δ₀；

S54：判断采深H、边界角δ₀、下沉盆地边界点到拐点的距离D以及拐点偏移距D₀是否满足以下公式：

其中r为主要影响半径，

如果满足公式(13)则得到采深H，如果不满足公式(13)则重新执行步骤S51-53。

利用获取地表移动盆地的四个边界点都可以用来计算采深，但利用每个边界点计算得出的采矿深度并不都一致，这主要是由于掘进边界的误差较大，倾向方向的两边界可能未达到充分采动，因此使用掘进方向后侧的边界进行采深计算会得到较为理想的效果。

实施例：

1、模拟实验

为了验证所本发明所提出的方法，首先使用FLAC3D软件模拟了地表下沉，然后利用所提出的方法进行采空区位置的计算，并与模拟预设参数进行对比。FLAC3D采用了显式拉格朗日算法和混合离散分区技术，能够准确模拟材料的塑性破坏，尤其对于本次实验的大变形非线性问题，其运算速度与精度都优于有限元程序，具体的操作步骤包括：

(1)模型的构建

本次模型共搭建了9层岩层，包括煤层与底板，且在每层覆岩间设置了介面，以保证岩层间可以发生滑动或分离。所有岩层均采用摩尔库伦本构模型。模型块体的水平尺寸为10m*10m，在保证精度的同时提高运算效率。当最大不平衡力低于50N时认为模型达到稳定。模拟的采空区位置如图13与表2所示。为评估非充分采动的影响，令采空区宽度大于长度。煤层厚度设为2m以使下沉盆地形状更为明显。采深与煤层厚度的比值远大于30，以保证地表形变连续，无裂缝产生。

(2)模拟结果

根据上述建立的模型可以得到如图13所示的下沉盆地与图14所示的盆地倾向方向的倾斜与曲率。图14中(a)表示的是倾斜，(b)表示曲率，且图14中(a)底部绿色直线与(b)中底部黄色直线分别为倾斜与曲率为0的位置；从图13中可以看到由于煤层倾斜，下沉盆地中心与采空区中心并不重合，而是向下山方向移动。从图14中可以看到倾斜煤层的倾斜不对称。相比于下山方向，上山方向的倾斜变化更为剧烈。图14(a)中下沉盆地倾斜波峰与波谷之间的零值为一条线，表明此时下沉盆地并未达到充分采动。由于本次模拟采空区宽度大于长度，因此可以得知下沉盆地在走向和倾向两个方向均未达到充分采动。

2、反算结果及分析

根据上述方法得到反演结果及误差如表2所示。由表2计算可得到平均误差为4.71％。由于模拟实验下沉盆地不受InSAR监测精度影响，因此可以认为表2中误差均为模型误差。除倾向长度与上山方向采深外，其余参数反算误差均在2％以下。倾向长度直接由下沉盆地拐点及拐点偏移距计算得到。由于走向长度误差较小，因此可以认为造成倾向长度误差较大的原因为煤层倾斜导致的下沉盆地拐点位置发生变化。上山方向采深由主要影响半径计算得到。与下山方向采深相比，误差明显增大。同样的，产生该现象的原因是煤层倾斜导致上山方向主要影响半径与采深之间的关系发生变化，综上可以看出，本申请提出的方法在原理上可行，能够得到工作面上的反演结果。

表2工作面位置反演结果及误差

以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (8)

1.一种面向越界开采识别的地下开采面反演方法，其特征在于，包括以下步骤：

步骤S1：通过矿山地表工作面上方的时序SAR影像，利用D-InSAR技术获取时序相邻两景SAR数据间开采沉陷盆地的地表形变信息，并将获取到的下沉图依次叠加，得到地表移动盆地的沉降图。

步骤S2：对获取的沉降图进行处理并提取出地表形变信息和开采沉陷特征，得到地表移动盆地主断面上的地表最大下沉点、盆地边界点参数；

步骤S3：根据获得的地表移动盆地主断面上的参数，根据随机介质理论得到地表下沉曲线与曲率曲线；

步骤S4：根据地表下沉曲线和曲率曲线计算得出边界点和拐点以及拐拐点偏移距和边界角；

步骤S5：根据所述边界点和拐点得出采空区范围，并且根据拐点偏移距和边界角进行反演采深；

步骤S6：将计算出的采空区和采深进行空间上的叠加反演出地下开采面，并基于采矿权范围数据判断是否存在越界开采行为。

2.根据权利要求1所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S2中提取出地表形变信息的具体步骤包括：

S21：采用重复轨道干涉测量模式，从获取的干涉相位图中定义干涉相位的综合贡献的组成为：

其中ω表示缠绕算子，

为卫星视线方向的地表形变相位；

为参考面相位，

为地形相位，

为大气相位，

为噪声相位。通过两幅或多幅干涉影像进行差分处理去除平地效应，并逐一将上式中后四种相位消除，即可分离出地表形变信息；

S22：通过对基线进行精确估算去除参考面相位，通过模拟外部DEM去除地形相位，通过提高信噪比来降低大气相位，通过低通滤波的方式抑制噪声相位；

S23：最终分离出地表形变信息。

3.根据权利要求1所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S2中从干涉相位图中得出的开采沉陷特征包括空间特征、几何特征以及形变特征；

空间特征是地下开采会引起相应的采矿活动区域上方出现地表下沉现象，最大的沉降量主要发生在采矿区域的地表中心，从中心到边缘下沉幅度逐渐减小，最终在该采矿区域表面形成一个空间漏斗，则漏斗中各点的沉降量为：

f_z(x_n,y_n)＜f_z(x_m,y_m)＜0且df_z(x,y)/dr＞0 (2)，

其中，f_z(x,y)是下沉区z的表达式，(x_n,y_n)和(x_m,y_m)代表下沉范围内的点，且(x_m,y_m)离下沉中心更近；

几何特征基于空间特征得出，即所述空间漏斗的沉降中心由四周的斜坡包围，形变中心点的梯度绝对值大于没有形变的区域的梯度绝对值，梯度的大小代表地表下沉的幅度信息，方向代表的是相位信息，且相位信息的计算公式为：

形变特征为地下开采引起的底板沉陷区域，在干涉相位图中通常为一系列闭合的圆形，这些圆被近似为干涉图坐标系上的一组小椭圆，表示为：

其中，∑(x_n,y_n)＜∑_t，且n＝1,…,N；x和y是在轮廓中分别沿着距离向和方位向上点的坐标；x_n和y_n是椭球体的中心点坐标；a和b分别是长半轴和短半轴；θ为主轴上的倾角；∑(·)代表这是椭球体的大小，且都小于常数∑_t。

4.根据权利要求1所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S3中所述的地表下沉曲线表示开采活动导致的地表下沉的分布规律，即在地表最大下沉点O处下沉值最大，自盆地中心至盆地边缘地表下沉值逐渐减小，在盆地边界点A和B处下沉值为零，且地表下沉曲线以采区中央对称。

5.根据权利要求1所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S3中所述的曲率曲线是表示地表移动盆地内曲率的变化规律，曲率曲线可表示为下沉的二阶导数，曲率曲线的分布规律表现为：盆地边界点A和B处和拐点E处曲率为零，且盆地边缘区为正曲率区，盆地中部为负曲率区。

6.根据权利要求1所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S5中所述的计算采空区范围的具体操作步骤包括：

S51:在三维空间坐标系中，基于概率积分法得出单元开采引起地表任意点的下沉最终值公式为：

W_g(x，y)＝1/r²·exp(-π(x-x_i)/r²)·exp(-π(y-y_i+l_i)/r²) (5)，

其中，r＝H₀/tgβ，l_i＝H_iCtgθ，r为主要影响半径，H₀为平均采深，β为主要影响角，l_i为采动影响长度，H_i为单位的深度，C为采动影响系数，θ为最大下沉角，B(x_i，y_i)为单元中心点的平面坐标，A(x，y)为地表上任意一点的坐标；

S52：设工作面的开采范围为0：D₁和0：D₂组成的矩形采空区，即工作面走向方向的开采长度为D₁，沿工作面倾向方向的开采宽度为D₂，则整个工作面开采引起地表任意点下沉的概率积分法计算公式为：

其中，W₀为该地质采矿条件下的最大下沉值，且W₀＝mqcosα，m为煤层开采厚度，q为下沉系统，α为煤层倾角；

公式(3)可表示为：

其中，W₀为走向和倾向均达到充分采动时地表的最大下沉值，W⁰(x)为倾向方向达到充分采动时走向主断面横坐标为x的点的下沉值，W⁰(y)为走向方向达到充分采动时倾向主断面横坐标为y的点的下沉值；

S53：根据获得的地表下沉的公式，计算出地表上任意一点A(x，y)的倾斜、曲率、水平移动、水平变形值。

7.根据权利要求5所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S53中所述的任意一点的倾斜、曲率、水平移动、水平变形值具体计算公式分别为：

a.倾斜：对于坐标为(x,y)的点沿

方向的倾斜i

为下沉W(x,y)在

方向上单位距离的变化率，其计算公式为：

公式(5)可简化为：

b.曲率：对于坐标为(x,y)的点沿

方向的曲率k

为倾斜i

在

方向上单位距离的变化率，其计算公式为：

公式(6)可简化为

c.水平移动：对于坐标为(x,y)的点沿

方向的水平移动U

的计算公式为：

d.水平变形：对于坐标为(x,y)的点沿

方向的水平变形ε

的计算公式为：

其中，i⁰(x)、k⁰(x)、U⁰(x)、ε⁰(x)分别别表示倾向方向达到充分采动时走向主断面上横坐标为x的点的倾斜、曲率、水平移动和水平变形值，i⁰(y)、k⁰(y)、U⁰(y)、ε⁰(y)分别表示走向方向达到充分采动时倾向主断面上横坐标为y的点的倾斜、曲率、水平移动和水平变形值。

8.根据权利要求1所述的一种面向越界开采识别的地下开采面反演方法，其特征在于，步骤S5中所述的计算采深的具体操作步骤包括：

S51：根据经验值，给出边界角δ₀的初始值，并且边界角δ₀取值和岩层性质的关系为：

δ₀＝g(k) (14)，

其中，k＝f(k₁，k₂，k₃，L，k_n)，(k₁，k₂，k₃，L，k_n)k为地下工程上覆岩层岩性系数；

S52：根据覆岩层岩性系数k得出采深H的计算公式为：

H＝h(k) (15)，

S53：根据采深H得知覆岩岩性，由岩性计算k，并根据公式(14)得出边界角δ₀；

S54：判断采深H、边界角δ₀、下沉盆地边界点到拐点的距离D以及拐点偏移距D₀是否满足以下公式：

其中r为主要影响半径，

如果满足公式(16)则得到采深H，如果不满足公式(16)则重新执行步骤S51-53。

Images