CN103091675A - 一种基于InSAR技术的矿区开采监测方法 - Google Patents

Abstract

本发明是一种基于InSAR技术的矿区开采监测方法，首先通过求取与待监测矿区临近的矿区的概率积分法模型系数，通过利用InSAR技术获取的待监测矿区的雷达视线向形变场，将待监测矿区工作面的长度、宽度、厚度、开采深度、走向方位角、中心点坐标作为未知数与待监测矿区临近的矿区的概率积分法模型系数一起带入概率积分法模型，其次利用遗传算法搜索得出待监测矿区工作面的参数值；最后把遗传算法得到的工作面参数值作为模式搜索法的初始值，经过迭代搜索，得出待监测矿区的准确工作面参数值。克服了现有技术中对矿区开采监测时只能得到开采的大致位置，不能精确地获得地下采空区详细参数信息的问题，大大拓宽了InSAR技术在矿区的应用空间，提供了一种低成本、大范围的矿区开采精细化监测的手段。

Description

一种基于InSAR技术的矿区开采监测方法

技术领域

本发明涉及一种基于InSAR技术的矿区开采监测方法。

背景技术

随着对煤炭资源需求的不断增长以及资源日益减少之间的矛盾越来越突出，在利益的驱使下，许多地区出现了没有经过批准和规划的矿区开采现象。这种开采不仅会导致资源的破坏和浪费，妨碍国家对资源的有序管理，同时还会产生许多不明采空区，容易发生突水事故，从而威胁到矿山及工作人员的生命财产安全。此外地下开采会导致地表发生下沉和变形，使采空区上方的建筑物和构筑物发生变形破坏。因此如何准确地获取地下采空区的深度H，中心点坐标(X_c,Y_c)，采空区走向长D₃、倾向长D₁、开采厚度m、走向方位角等信息，实现对地下煤炭资源的开采进行精确监测，显得尤为重要。

为了实现对矿区开采的监测，目前主要的手段有传统的矿政人员定期或不定期巡查、探地雷达探测采空区、微震探测技术、InSAR（Interferometric Synthetic Aperture Radar，简称InSAR）技术监测等，但是这些方法都有很大的局限和不足。传统的矿政人员巡查随机性大，代价较高，效率较低；探地雷达和微震技术由于其工序复杂，监测时间较长，成本较高等限制致使其对国家监管部门而言不适用，InSAR技术以其全天候、全天时、覆盖范围广、成本低、精度高等优点，将其应用于矿区开采监测不仅可以降低监测成本，同时由于其覆盖范围较大，可以实现大范围的监测，提高了监管效率。但是传统的基于InSAR技术的监测方法，是将InSAR监测的形变场与监管部门登记的合法的开采区对比，得出矿区开采区域信息，但是只能得到矿区开采的大致位置，对于地下采空区具体信息则无法获取，从而无法满足国家监管部门对开采矿区定量精细化监测的要求。所以，迫切需要一种既利用InSAR技术的优点，又能获得矿区开采区域地下采空区的详细参数信息的方法。

发明内容

本发明提出了一种基于InSAR技术的矿区开采监测方法，其目的在于提供一种低成本、大范围的矿区开采精细化监测的方法，克服现有的开采监测方法成本较高、效率较低且无法获得开采工作面参数详细信息等缺陷。

本发明采用的技术方案如下：

一种基于InSAR技术的矿区开采监测方法，首先通过求取与待监测矿区临近的矿区的概率积分法模型系数，利用InSAR技术获取的开采矿区雷达视线向形变场，将待监测矿区工作面的长度、宽度、厚度、开采深度、走向方位角、中心点坐标作为未知数与求取的临近矿区的概率积分法模型系数带入概率积分法模型，利用遗传算法搜索得出待监测矿区的工作面的参数值；把遗传算法得到的工作面参数值作为模式搜索法的初始值，经过迭代搜索，最后得出待监测矿区的准确采空区参数值，包括采空区的中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角。

该方法具体操作步骤如下：

步骤1：利用InSAR技术获取待监测的开采矿区经过地理编码后的雷达视线向形变场，即LOS；

所述地理编码，是指将雷达影像坐标系转换到通用横轴墨卡托投影（Universal TransverseMercatol projection，简称UTM）坐标系；

步骤2：利用临近矿区的采空区参数和观测值求取临近矿区的概率积分法模型系数；

其中，临近矿区的采空区参数包括采空区走向长D3、倾向长D1、开采厚度m、开采深度H、走向方位角

，观测值包括利用水准仪、全站仪或GPS测量的临近矿区中任意点的地表垂直下沉值W和水平移动值U；

概率积分法模型中任一点的下沉和水平移动预计公式如下：

其中，(x,y为地表任一点的坐标，erf为概率积分函数，其形式为

u为积分参数，W₀为该地质采矿条件下的最大下沉值，W₀＝m·q·cosα，α为煤层倾角，

为主要影响半径，L为倾向工作面计算长度，

D₁为工作面倾向斜长，l为走向有限开采时的计算长度，l＝D₃-s₃-s₄；

求出临近矿区的概率积分法模型的系数：

包括下沉系数q，取值范围为0.01～1；

下山拐点偏移距、上山拐点偏移距分别为s₁、s₂，取值范围均为0.05H～0.3H；

走向左、右拐点偏移距包括s₃、s₄，取值范围均为0.05H～0.3H；

走向、倾向下山、倾向上山水平移动系数分别为b、b₁、b₂，取值范围均为0.1～0.4；

走向、倾向下山、倾向上山主要影响角正切分别为tanβ、tanβ1、tanβ2，取值范围均为1～3.8；

开采影响传播角θ₀＝90°-kα，α为煤层倾角，可根据实际测量获取，其中k的取值范围为0.5～0.8，煤层倾角α，其取值范围为0～45°；

所述临近地质矿区是指与开采矿区的采煤方法和顶板管理方法相同，工作面煤矿上面的岩石力学性质、岩层分布、开采厚度和深度有70%以上的相同；

步骤3：利用开采矿区的概率积分法模型系数等于临近矿区的概率积分法模型系数，设计搜索算法的适应度函数f，其形式为：f(goaf)＝||LOS-LOS′||；

式中，goaf表示开采矿区的采空区的待测参数，即采空区的中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角，LOS为监测的开采矿区的雷达视线向形变场，从步骤1中获得，LOS′为利用概率积分法模型求得的雷达视线向形变场；

利用概率积分法模型和步骤2中求得的概率积分法模型系数计算获得开采矿区的东西、南北方向水平移动场U_E、U_N及垂直方向下沉值W，按照雷达成像几何条件转换到雷达视线向，得到计算雷达视线向形变场，即LOS′；

开采矿区的地表任一点(x′,y′)在东西、南北方向的水平移动

及下沉W′＝W(x′,y′)；

所述雷达成像几何条件由下式表征：

LOS＝Wcosθ-sinθ[U_Ncos(α_h-3π/2)+U_Esin(α_h-3π/2)]，

其中，θ为雷达卫星入射角，α_h为卫星飞行方位角，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得；

则，LOS′＝W′cosθ-sinθ[U_N′cos(α_h-3π/2)+U_E′sin(α_h-3π/2)]，其中，θ为雷达卫星入射角，α_h为卫星飞行方位角，其值与从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得的值相同；

步骤4：利用遗传算法计算符合迭代次数的种群个体，选取适应度函数最小值对应的个体，即经过遗传算法搜索迭代得到的开采矿区采空区的参数goaf′，其参数包括开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角；

步骤A：设置种群大小为N，迭代次数为M，个体基因包括采空区中心点坐标为(X_c′,Y_c′)，走向长D₃′、倾向长D₁′、开采厚度m′、深度H′，走向方位角

随机生成开采矿区采空区参数goaf的初始种群个体，即初始个体基因，其中goaf是指开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角；

步骤B：按照步骤3根据个体基因计算得到LOS′，同时按照目标函数计算个体的适应度，按照设定的迭代终止条件，即目标函数小于目标设定阈值或者为超过迭代次数，判断迭代结果是否需要继续迭代，若不满足，则对种群个体进行选择、交叉、变异操作，得到新的种群个体，然后再重复步骤B，若满足，则输出种群中适应度函数值最小的个体，即开采矿区采空区的参数值goaf′，即通过遗传算法得出的开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角的值；

步骤5：采用模式搜索法提高经遗传算法得到的开采矿区采空区参数值的精度；

步骤A：设置格网大小、格网扩展因子及格网收缩因子，迭代终止条件即格网大小达到格网大小下限值；令goaf′作为模式搜索法中待求解参数的初始值goaf″_old，即利用模式搜索法得出的开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角的值，初始值包含n个参数，且令f(goaf″_old)为目标函数值A；

步骤B：将一个格网大小作为步长，根据模式搜索法生成2n个模式向量（即搜索方向），将goaf″_old在2n个方向上按照步长移动，得到2n个新的goaf″_new；

按照步骤3分别计算f(goaf″_new)，从2n个f(goaf″_new)选出值最小的f_min(goaf″_new)作为目标函数值B；比较B和A，若B大于A，则将格网按照格网扩展因子增大；否则，则将格网按照格网收缩因子缩小；

判断此时格网大小是否满足迭代终止条件，若不满足，将f_min(goaf″_new)对应的goaf″_new赋值给goaf″_old，返回步骤B；若满足，则将f_min(goaf″_new)对应的goaf″_new作为求解参数的结果，得到开采矿区采空区的参数值，即利用模式搜索法得出的开采矿区的采空区的中心点坐标、走向长、倾向长、开采厚度、开采深度、走向方位角的值，提高经遗传算法得到的开采矿区采空区参数值的精度。

所述步骤4中遗传算法的种群大小为20，种群生成方式为随机生成、种群选择操作的方式为轮盘赌、交叉概率为0.7、交叉方式为两点交叉、变异函数为高斯函数，迭代次数为300代。

所述步骤5中的模式搜索法的初始格网大小为1，扩展因子为2、收缩因子0.5、迭代终止条件即格网尺寸<10^-6。

所述地表下沉W⁰(x)、水平移动U⁰(x)、地表下沉W⁰(y)和水平移动U⁰(y)的预计公式如下：

1）当矿区地下开采工作面倾向充分开采，走向主断面为有限开采时，造成的地表下沉W⁰(x)和水平移动U⁰(x)的预计公式如下：

$\left\{\begin{array}{c}{W}^0\left(x\right)=\frac{W_0}{2}\{\mathrm{erf}\left(\frac{\sqrt{\mathrm{\&pi;}}}{r}x\right)-\mathrm{erf}[\frac{\sqrt{\mathrm{\&pi;}}}{r}(x-l)\left]\right\}\\ U^0\left(x\right)=b\&CenterDot;W_0[e^{-\mathrm{\&pi;}\frac{x^2}{r^2}}-e^{-\mathrm{\&pi;}\frac{{(x-l)}^2}{r^2}}]\\ W_0=m\&CenterDot;q\&CenterDot;\cos \mathrm{\&alpha;}\end{array}\right.---\left(1\right)$

其中，x为地表点的计算横坐标值，横坐标原点，定义为煤壁向采空区或煤柱方向偏移s₃后所对应的地表点，当拐点偏移距s₃>0，则偏向采空区，反之偏向采空区相反方向，横坐标轴沿着工作面走向指向采空区一侧，erf为概率积分函数，其形式为

u为积分参数，W₀为该地质采矿条件下的最大下沉值，

为主要影响半径，b为水平移动系数，α为煤层倾角，q为下沉系数，l为走向有限开采时的计算长度，其值为工作面走向长度D₃与走向左、右边界的拐点偏移距s₃、s₄之差，即：l＝D₃-s₃-s₄；

2）当矿区地下开采工作面走向充分开采，倾向有限开采时，造成的地表下沉W⁰(y)和水平移动U⁰(y)的预计公式为：

$\left\{\begin{array}{c}{W}^0\left(y\right)=\frac{W_0}{2}\{\mathrm{erf}\left(\frac{\sqrt{\mathrm{\&pi;}}}{r_1}y\right)-\mathrm{erf}[\frac{\sqrt{\mathrm{\&pi;}}}{r_2}(y-L)\left]\right\}\\ U^0\left(y\right)=W_0(b_1{e}^{-\mathrm{\&pi;}\frac{x^2}{r_1^2}}-b_2{e}^{-\mathrm{\&pi;}\frac{{(x-L)}^2}{r_2^2}})-\mathrm{ct}{\mathrm{g\&theta;}}_0{W}^0\left(y\right)\end{array}\right.---\left(2\right)$

其中，y为地表点的计算纵坐标，纵坐标原点，定义为倾向下山煤壁向采空区或煤柱方向偏移s1后，按照开采影响传播角投影到地表所对应的点，纵坐标轴沿倾向指向采空区一侧，b₁、b₂为下山、上山水平移动系数；

为下山、上山主要影响半径，

分别为下、上山的采深，s₁、s₂分别为下、上山拐点偏移距，开采影响传播角θ₀＝90°-kα，α为煤层倾角，可根据实际测量获取，其中k的取值范围为0.5～0.8，煤层倾角α，其取值范围为0～45°，L为倾向工作面计算长度，其计算表达式为： $L=(D_1-s_1-s_2)\frac{\sin ({\mathrm{\&theta;}}_0+\mathrm{\&alpha;})}{\sin {\mathrm{\&theta;}}_0},$ D₁为工作面倾向斜长；

当工作面倾向未达到充分开采时，利用公式（1）预计走向地表下沉W⁰(x)和水平移动U⁰(x)，求出的形变值需要乘上一个小于1的倾向采动程度系数W_ym是假设走向已经达到充分开采，利用公式（2）中第一式求出W⁰(y)的最大值W_ym；

当工作面走向未达到充分开采时，倾向地表下沉和水平移动同样利用公式（2）预计，求出的形变值需乘上走向采动程度系数

W_xm是假设倾向已经达到充分开采，利用式（1）中第一式求出W⁰(x)的最大值W_xm；

求取W_xm、W_ym时，x和y的取值范围为矿区采深的正负2-3倍，利用公式（1）和公式（2）求出取值范围中走向和倾向的的最大下沉值，即为W_xm、W_ym。

有益效果

本发明是一种基于InSAR技术的矿区开采监测方法，首先通过求取与待监测矿区临近的矿区的概率积分法模型系数，通过利用InSAR技术获取的待监测矿区的雷达视线向形变场，将待监测矿区工作面的长度、宽度、厚度、开采深度、走向方位角、中心点坐标作为未知数与待监测矿区的临近矿区的概率积分法模型系数一起带入概率积分法模型，利用遗传算法搜索得出待监测矿区的工作面的参数值；把遗传算法得到的工作面参数值作为模式搜索法的初始值，经过迭代搜索，得出待监测矿区的准确参数值。克服了现有技术中对矿区开采监测时只能得到开采的大致位置，不能精确地获得地下采空区的详细参数信息的问题，大大拓宽了InSAR技术在矿区的应用空间，提供了一种低成本、大范围的矿区开采精细化监测的手段。

另外，当分属不同公司的相邻工作面同时对地表建筑物或构筑物造成破坏且无法确定各自的责任大小时，利用本发明的方法,首先分别计算出两个公司相邻的工作面的参数值，然后利用工作面参数值基于概率积分法模型计算出各自造成的地表建筑物和构筑物的移动变形情况，从而可以很好的确定各自的责任；对于公路或者铁路需要穿过老采空区时，同样可以利用本发明的方法在获得该老采空区开采过程中的SAR卫星影像及概率积分法预计参数时，计算获得老采空区的参数信息，从而为公路、铁路穿过老采空区的设计提供依据。

附图说明

图1为本发明的数据处理流程图；

图2为工作面倾向充分开采，走向有限开采时，走向断面下沉计算原理及坐标关系图；

图3为工作面走向充分开采，倾向有限开采时，倾向断面下沉计算原理及坐标关系图；

图4为地下开采导致的地表下沉盆地任意点的空间坐标系图。

具体实施方式

为了使本技术领域的人员能够更好的理解本发明的方法，下面将结合本发明实施例中附图，对本发明的实施方案进行清楚、详细的描述。

如图1所示，为本发明的流程图，其基本原理为：

本发明利用概率积分法模型将矿区采空区与该采空区在开采过程中导致的地表形变建立联系，在得到利用InSAR技术获取该矿区地表形变场以及概率积分法模型系数时，即可利用地表形变场计算出矿区采空区的工作面参数信息，实现利用地表监测数据监测地下采空区参数信息的目的。

概率积分法模型因其所用的移动和变形预计公式中含有概率积分函数（或概率积分法函数导数）而得名，由我国学者刘宝琛、廖国华等从随机介质理论发展得来，目前已成为我国较为成熟的、应用最为广泛的预计方法，众多地区均以该模型为基础，通过地表实际监测得到了该模型的系数，包括下沉系数q；倾向下山、倾向上山、走向左、右拐点偏移距s₁、s₂、s₃、s₄；走向、倾向下山、倾向上山水平移动系数b、b₁、b₂；走向、倾向下山、倾向上山主要影响角正切tanβ、tanβ1、tanβ2；开采影响传播角θ₀。

概率积分法模型的基本原理如下：

矿区地下煤层采出后，地表下沉值达到该地质采矿条件下应有的最大值W₀，此时的采动称为充分开采，反之若没有达到，则为非充分开采；

当矿区煤层从起始开采点一直采到无穷远处，这样的开采叫做半无限开采，反之，则为有限开采；

当矿区地下开采工作面倾向主断面充分开采，走向主断面为有限开采时，造成的地表下沉W⁰(x)，水平移动U⁰(x)预计公式为：

$\left\{\begin{array}{c}{W}^0\left(x\right)=\frac{W_0}{2}\{\mathrm{erf}\left(\frac{\sqrt{\mathrm{\&pi;}}}{r}x\right)-\mathrm{erf}[\frac{\sqrt{\mathrm{\&pi;}}}{r}(x-l)\left]\right\}\\ U^0\left(x\right)=b\&CenterDot;W_0[e^{-\mathrm{\&pi;}\frac{x^2}{r^2}}-e^{-\mathrm{\&pi;}\frac{{(x-l)}^2}{r^2}}]\\ W_0=m\&CenterDot;q\&CenterDot;\cos \mathrm{\&alpha;}\end{array}\right.---\left(1\right)$

式中，x为地表点的计算横坐标值（其坐标原点，定义为煤壁向采空区或煤柱方向偏移s₃后所对应的地表点，当拐点偏移距s₃>0，则偏向采空区，反之偏向采空区相反方向，横坐标轴沿着工作面走向指向采空区一侧，其示意图如图2所示），erf为概率积分函数，其形式为

u为积分参数，W₀为该地质采矿条件下的最大下沉值，

为主要影响半径，b为水平移动系数，α为煤层倾角，q为下沉系数，l为走向有限开采时的计算长度，其值为工作面走向长度D₃与走向左、右边界的拐点偏移距s₃、s₄之差，即：l＝D₃-s₃-s₄。

走向充分开采，倾向有限开采时其地表下沉W⁰(y)和水平移动U⁰(y)的预计公式为：

$\left\{\begin{array}{c}{W}^0\left(y\right)=\frac{W_0}{2}\{\mathrm{erf}\left(\frac{\sqrt{\mathrm{\&pi;}}}{r_1}y\right)-\mathrm{erf}[\frac{\sqrt{\mathrm{\&pi;}}}{r_2}(y-L)\left]\right\}\\ U^0\left(y\right)=W_0(b_1{e}^{-\mathrm{\&pi;}\frac{x^2}{r_1^2}}-b_2{e}^{-\mathrm{\&pi;}\frac{{(x-L)}^2}{r_2^2}})-\mathrm{ct}{\mathrm{g\&theta;}}_0{W}^0\left(y\right)\end{array}\right.---\left(2\right)$

式中，y为地表点的计算纵坐标（其坐标原点，为倾向下山煤壁向采空区或煤柱方向偏移s1后，按照开采影响传播角投影到地表所对应的点，纵坐标轴沿倾向指向采空区一侧，如图3所示，）b₁、b₂为下山、上山水平移动系数；

分别为下、上山的采深）为下山、上山主要影响半径、s₁、s₂为下、上山拐点偏移距，θ₀为开采影响传播角。L为倾向工作面计算长度，其计算表达式为： $L=(D_1-s_1-s_2)\frac{\sin ({\mathrm{\&theta;}}_0+\mathrm{\&alpha;})}{\sin {\mathrm{\&theta;}}_0},$ D₁为工作面倾向斜长。

矿区下沉盆地内，任一点的下沉和水平移动预计公式为：

当工作面倾向未达到充分开采时，同样利用公式（1）预计走向地表下沉和水平移动，但求出的形变值需要乘上一个小于1的倾向采动程度系数

W_ym是假设走向已经达到充分开采，利用式（2）中第一式求出W⁰(y)的最大值；

同样，当工作面走向未达到充分开采时，倾向地表下沉和水平移动同样利用公式（2）预计，求出的形变值需乘上走向采动程度系数

W_xm是假设倾向已经达到充分开采，利用式（1）中第一式求出W⁰(x)的最大值。

本实施例中，求取W_xm、W_ym时，x、y分别从-1000m到1000m，每隔5m取一个x、y，计算出所有坐标对应的走向和倾向的下沉值，找出各自的最大下沉值即为W_xm、W_ym，x和y的取值范围为矿区采深的正负2-3倍。

将非充分开采条件下预计的走向、倾向地表下沉和水平移动带入公式（3），得出非充分开采下地表任一点、任意方向的下沉和水平移动。

本实施例的具体步骤如下：

步骤1：利用InSAR技术获取待监测的开采矿区经过地理编码后的雷达视线向形变场，即LOS；

所述地理编码，是指将雷达影像坐标系转换到通用横轴墨卡托投影（Universal TransverseMercatol projection，简称UTM）坐标系；

步骤2：利用临近矿区的采空区参数和观测值求取临近矿区的概率积分法模型；

其中，临近矿区的采空区参数包括采空区走向长D3、倾向长D1、开采厚度m、开采深度H、走向方位角

，观测值包括利用水准仪、全站仪或GPS测量临近矿区中任意点的地表垂直下沉值W和水平移动值U；

对地下开采导致的地表移动盆地内的任意点的移动和变形，当工作面为矩形或近似矩形时，其模型表达式如下，坐标系如图4所示；

任一点下沉和水平移动预计的概率积分法模型：

其中，erf为概率积分函数，其形式为

u为积分参数，W₀为该地质采矿条件下的最大下沉值，W₀＝m·q·cosα，

为主要影响半径，L为倾向工作面计算长度，

D₁为工作面倾向斜长，l为走向有限开采时的计算长度，l＝D₃-s₃-s₄；

求出临近矿区的概率积分法模型的系数：

包括下沉系数q，取值范围为0.01～1；

倾向下山、倾向上山拐点偏移距s₁、s₂，取值范围均为0.05H～0.3H；

走向左、右方向的拐点偏移距s₃、s₄，其取值范围均为0.05H～0.3H；

走向、倾向下山、倾向上山水平移动系数b、b₁、b₂，其取值范围均为0.1～0.4；

走向、倾向下山、倾向上山主要影响角正切分别为tanβ、tanβ1、tanβ2，其取值范围均为1～3.8；

开采影响传播角θ₀＝90°-kα，α为煤层倾角，可根据实际测量获取，其中k的取值范围为0.5～0.8，煤层倾角α，其取值范围为0～45°；

若能够获得该开采矿区或者地质采矿条件类似的矿区的概率积分法模型系数，将工作面的参数信息作为未知数带入概率积分法估计模型系数，就可得到在该采空区参数和概率积分法模型系数下，导致的地表东西、南北方向的水平移动值U_E、U_N及垂直方向下沉值W，按照雷达成像几何条件（如公式4所示）就可将地表东西、南北方向的水平移动值U_E、U_N及垂直方向下沉值W转换到SAR卫星的雷达视线向，即LOS′，然后通过与InSAR技术监测的开采矿区地表的雷达视线向形变场（LOS）比较，当两者相似度达到终止条件，即目标函数小于目标设定阈值时，得到的采空区参数值，该参数值认为是地下采空区的真实参数信息，即获得了开采矿区的采空区的开采深度H′，中心点坐标(X_c′,Y_c′)，采空区走向长D₃′、倾向长D₁′、开采厚度m′、走向方位角为

利用采空区的参数信息，监管部门可知道开采矿区在哪儿采（由中心点坐标(X_c′,Y_c′)、开采深度H′获得，即采空区的三维坐标）、采多少（由采空区走向长D₃′、倾向长D₁′、开采厚度m′获得采空区的三维边界）以及采空区的分布（由走向方位角为

获得），从而实现基于InSAR技术的矿区开采精细化监测。

本实例中通过计算获得q=0.7，s1=s2=s3=s4=10m，b1=b2=b3=0.3，tanβ=tanβ=tanβ=2.0，θ₀=82°，α=10°，目标阈值设置为0.0005。

由于式（3）较为复杂，直接求解难度较大，同时为了避免求解参数陷入局部最优解，所以，本发明提出了利用遗传算法得出的采空区参数值作为模式搜索法的初始值，然后通过迭代搜索得到精度更高的采空区的参数信息。

所述临近地质矿区是指与开采矿区的采煤方法和顶板管理方法相同，工作面煤矿上面的面的岩石力学性质、岩层分布、开采厚度和深度有70%以上的相同。

步骤3：假设开采矿区的概率积分法模型系数等于临近矿区的概率积分法模型系数，设计搜索算法的适应度函数f，其形式为：f(goaf)＝||LOS-LOS′||；

式中，goaf表示开采矿区采空区的待测参数，LOS为监测的开采矿区的雷达视线向形变场，从步骤1中获得，LOS′为利用概率积分法模型求得的雷达视线向的形变场；

利用概率积分法模型和步骤2中的求得的概率积分法模型系数计算获得开采矿区的东西、南北方向水平移动场U_E、U_N及垂直方向下沉值W，按照雷达成像几何条件合成到雷达视线向，得到计算雷达视线向形变场，即LOS′；

开采矿区的地表任一点(x′,y′)在东西、南北方向的水平移动

及下沉W′＝W(x′,y′)；

所述雷达成像几何条件由下式表征：

LOS＝Wcosθ-sinθ[U_Ncos(α_h-3π/2)+U_Esin(α_h-3π/2)]  （4）

其中，θ为雷达卫星入射角，α_h为卫星飞行方位角，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得；

则，LOS′＝W′cosθ-sinθ[U_N′cos(α_h-3π/2)+U_E′sin(α_h-3π/2)]，其中，θ为雷达卫星入射角，α_h为卫星飞行方位角，其值与从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得的值相同；

步骤4：利用遗传算法计算符合迭代次数的种群个体，选取适应度函数最小值对应的个体，即经过遗传算法搜索迭代得到的开采矿区采空区的参数goaf′；

遗传算法（Genetic Algorithm，简称GA）是一种以遗传学为基础的全局最优寻优搜索算法，遗传算法是将问题的求解参数表示成染色体（chromosome），即基因，利用随机生成的方法生成初始染色体群，即初始种群，并以目标函数作为基因对“环境”（即所求解的问题）适应度的量度工具，根据目标函数的增加或减小的原则选出相对较为适应环境的染色体进行复制，即选择（selection），通过交叉（crossover）、变异（mutation）两种遗传操作生成更加适应环境的种群（population），经过不断的迭代，最后得到收敛于一个最适应环境的个体，即求得该问题的全局最优解；

步骤A：设置种群大小为N为20，迭代次数为M为200代，个体基因包括中心点坐标为(X_c′,Y_c′)，采空区走向长为D₃′、倾向长D₁′、开采厚度m′、深度为H′，走向方位角为

随机生成开采矿区采空区的参数goaf的初始种群个体，即初始个体基因；

步骤B：按照步骤3根据个体基因计算得到LOS′，同时按照目标函数计算个体的适应度，按照设定的迭代终止条件判断迭代结果是否需要继续迭代，若不满足，则对种群个体进行选择、交叉、变异操作，得到新的种群个体，然后再重复步骤B，若满足，则输出种群个体即开采矿区采空区的参数值goaf′；其中种群的选择操作为轮盘赌操作，交叉概率设置为0.7，交叉方式设置为两点交叉，变异函数设置为高斯函数；

步骤5：采用模式搜索法提高遗传算法得到的开采矿区采空区的参数值的精度；

模式搜索法（pattern search）是一种不依赖与问题梯度，对函数连续性和是否可微无要求的数值优化算法，给定求解参数一个初始值x₀，计算该初始值处目标函数值f(x₀)，然后对每个参数按照一定的步长Δx_i生成新的点，比较f(x₀+Δx_i)、f(x₀)、f(x₀+Δx_i)，以其值最小的点作为临时点，继续生成其他参数的增量，并比较目标函数值，当完成全部参数的搜索时，选取当前的临时点作为新的搜索起点，同时搜索格网尺寸增加一倍。若搜索完成后，找不到比当前点更好的点，则搜索格网尺寸减半，如此迭代，直到达到精度要求或其他限制条件，则迭代结束；

步骤A：设置格网大小、格网扩展因子及格网收缩因子，迭代终止条件即格网大小达到格网大小下限值；令goaf′作为模式搜索法中待求解参数的初始值goaf″_old，初始值包含n个参数，且令f(goaf″_old)为目标函数值A；

步骤B：将一个格网大小作为步长，根据模式搜索法生成2n个模式向量（即搜索方向），将goaf″_old在2n个方向上按照步长移动，得到2n个新的goaf″_new；

按照步骤3分别计算f(goaf″_new)，从2n个f(goaf″_new)选出值最小的f_min(goaf″_new)作为目标函数值B；比较B和A，若B大于A，则将格网按照格网扩展因子增大；否则，则将格网按照格网收缩因子缩小；

判断此时格网大小是否满足迭代终止条件，若不满足，将f_min(goaf″_new)对应的goaf″_new赋值给goaf″_old，返回步骤B；若满足，则将f_min(goaf″_new)对应的goaf″_new作为求解参数的结果，得到开采矿区采空区的参数值。

Claims (5)

1.一种基于InSAR技术的矿区开采监测方法，其特征在于，首先通过求取与待监测矿区临近的矿区的概率积分法模型系数，利用InSAR技术获取的开采矿区雷达视线向形变场，将待监测矿区工作面的长度、宽度、厚度、开采深度、走向方位角、中心点坐标作为未知数与求取的临近矿区的概率积分法模型系数带入概率积分法模型，利用遗传算法搜索得出待监测矿区的工作面的参数值；把遗传算法得到的工作面参数值作为模式搜索法的初始值，经过迭代搜索，得出待监测矿区的准确采空区参数值，包括采空区的中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角。

2.根据权利要求1所述的一种基于InSAR技术的矿区开采监测方法，其特征在于，该方法具体操作步骤如下：

步骤1：利用InSAR技术获取待监测的开采矿区经过地理编码后的雷达视线向形变场，即LOS；

所述地理编码，是指将雷达影像坐标系转换到通用横轴墨卡托投影（Universal TransverseMercatol projection，简称UTM）坐标系；

步骤2：利用临近矿区的采空区参数和观测值求取临近矿区的概率积分法模型系数；

其中，临近矿区的采空区参数包括采空区走向长D3、倾向长D1、开采厚度m、开采深度H、走向方位角

，观测值包括利用水准仪、全站仪或GPS测量的临近矿区中任意点的地表垂直下沉值W和水平移动值U；

概率积分法模型中任一点的下沉和水平移动预计公式如下：

其中，(x,y为地表任一点的坐标，erf为概率积分函数，其形式为

u为积分参数，W₀为该地质采矿条件下的最大下沉值，W₀＝m·q·cosα，α为煤层倾角，

为主要影响半径，L为倾向工作面计算长度，

D₁为工作面倾向斜长，l为走向有限开采时的计算长度，l＝D₃-s₃-s₄；

求出临近矿区的概率积分法模型的系数：

包括下沉系数q，取值范围为0.01～1；

下山拐点偏移距、上山拐点偏移距分别为s₁、s₂，取值范围均为0.05H～0.3H；

走向左、右拐点偏移距包括s₃、s₄，取值范围均为0.05H～0.3H；

走向、倾向下山、倾向上山水平移动系数分别为b、b₁、b₂，取值范围均为0.1～0.4；

走向、倾向下山、倾向上山主要影响角正切分别为tanβ、tanβ1、tanβ2，取值范围均为1～3.8；

开采影响传播角θ₀＝90°-kα，α为煤层倾角，可根据实际测量获取，其中k的取值范围为0.5～0.8，煤层倾角α，其取值范围为0～45°；

所述临近地质矿区是指与开采矿区的采煤方法和顶板管理方法相同，工作面煤矿上面的岩石力学性质、岩层分布、开采厚度和深度有70%以上的相同；

步骤3：利用开采矿区的概率积分法模型系数等于临近矿区的概率积分法模型系数，设计搜索算法的适应度函数f，其形式为：f(goaf)＝||LOS-LOS′||；

式中，goaf表示开采矿区的采空区的待测参数，即采空区的中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角，LOS为监测的开采矿区的雷达视线向形变场，从步骤1中获得，LOS′为利用概率积分法模型求得的雷达视线向形变场；

利用概率积分法模型和步骤2中求得的概率积分法模型系数计算获得开采矿区的东西、南北方向水平移动场U_E、U_N及垂直方向下沉值W，按照雷达成像几何条件转换到雷达视线向，得到计算雷达视线向形变场，即LOS′；

开采矿区的地表任一点(x′,y′)在东西、南北方向的水平移动

及下沉W′＝W(x′,y′)；

所述雷达成像几何条件由下式表征：

LOS＝Wcosθ-sinθ[U_Ncos(α_h-3π/2)+U_Esin(α_h-3π/2)]，

其中，θ为雷达卫星入射角，α_h为卫星飞行方位角，其值从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得；

则，LOS′＝W′cosθ-sinθ[U_N′cos(α_h-3π/2)+U_E′sin(α_h-3π/2)]，其中，θ为雷达卫星入射角，α_h为卫星飞行方位角，其值与从步骤1中采用InSAR技术获取地理编码后的矿区雷达视线向形变场的过程中所涉及的雷达卫星影像头文件中获得的值相同；

步骤4：利用遗传算法计算符合迭代次数的种群个体，选取适应度函数最小值对应的个体，即经过遗传算法搜索迭代得到的开采矿区采空区的参数goaf′，其参数包括开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角；

步骤A：设置种群大小为N，迭代次数为M，个体基因包括采空区中心点坐标为(X_c′,Y_c′)，走向长D₃′、倾向长D₁′、开采厚度m′、深度H′，走向方位角

随机生成开采矿区采空区参数goaf的初始种群个体，即初始个体基因，其中goaf是指开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角；

步骤B：按照步骤3根据个体基因计算得到LOS′，同时按照目标函数计算个体的适应度，按照设定的迭代终止条件，即目标函数小于目标设定阈值或者为超过迭代次数，判断迭代结果是否需要继续迭代，若不满足，则对种群个体进行选择、交叉、变异操作，得到新的种群个体，然后再重复步骤B，若满足，则输出种群中适应度函数值最小的个体，即开采矿区采空区的参数值goaf′，即通过遗传算法得出的开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角的值；

步骤5：采用模式搜索法提高经遗传算法得到的开采矿区采空区参数值的精度；

步骤A：设置格网大小、格网扩展因子及格网收缩因子，迭代终止条件即格网大小达到格网大小下限值；令goaf′作为模式搜索法中待求解参数的初始值goaf″_old，即利用模式搜索法得出的开采矿区的采空区中心点坐标、走向长、倾向长、开采厚度、开采深度及走向方位角的值，初始值包含n个参数，且令f(goaf″_old)为目标函数值A；

步骤B：将一个格网大小作为步长，根据模式搜索法生成2n个模式向量（即搜索方向），将goaf″_old在2n个方向上按照步长移动，得到2n个新的goaf″_new；

按照步骤3分别计算f(goaf″_new)，从2n个f(goaf″_new)选出值最小的f_min(goaf″_new)作为目标函数值B；比较B和A，若B大于A，则将格网按照格网扩展因子增大；否则，则将格网按照格网收缩因子缩小；

判断此时格网大小是否满足迭代终止条件，若不满足，将f_min(goaf″_new)对应的goaf″_new赋值给goaf″_old，返回步骤B；若满足，则将f_min(goaf″_new)对应的goaf″_new作为求解参数的结果，得到开采矿区采空区的参数值，即利用模式搜索法得出的开采矿区的采空区的中心点坐标、走向长、倾向长、开采厚度、开采深度、走向方位角的值，提高经遗传算法得到的开采矿区采空区参数值的精度。

3.根据权利要求2所述的基于InSAR技术的矿区开采的监测方法，其特征在于，所述步骤4中遗传算法的种群大小为20，种群生成方式为随机生成、种群选择操作的方式为轮盘赌、交叉概率为0.7、交叉方式为两点交叉、变异函数为高斯函数，迭代次数为300代。

4.根据权利要求3所述的基于InSAR技术的矿区开采的监测方法，其特征在于，所述步骤5中的模式搜索法的初始格网大小为1，扩展因子为2、收缩因子0.5、迭代终止条件即格网尺寸<10^-6。

5.根据权利要求2所述的基于InSAR技术的矿区开采的监测方法，其特征在于，所述地表下沉W⁰(x)、水平移动U⁰(x)、地表下沉W⁰(y)和水平移动U⁰(y)的预计公式如下：

1）当矿区地下开采工作面倾向充分开采，走向为有限开采时，造成的地表下沉W⁰(x)和水平移动U⁰(x)的预计公式如下：

$\left\{\begin{array}{c}{W}^0\left(x\right)=\frac{W_0}{2}\{\mathrm{erf}\left(\frac{\sqrt{\mathrm{\&pi;}}}{r}x\right)-\mathrm{erf}[\frac{\sqrt{\mathrm{\&pi;}}}{r}(x-l)\left]\right\}\\ U^0\left(x\right)=b\&CenterDot;W_0[e^{-\mathrm{\&pi;}\frac{x^2}{r^2}}-e^{-\mathrm{\&pi;}\frac{{(x-l)}^2}{r^2}}]\\ W_0=m\&CenterDot;q\&CenterDot;\cos \mathrm{\&alpha;}\end{array}\right.---\left(1\right)$

其中，x为地表点的计算横坐标值，横坐标原点，定义为煤壁向采空区或煤柱方向偏移s₃后所对应的地表点，当拐点偏移距s₃>0，则偏向采空区，反之偏向采空区相反方向，横坐标轴沿着工作面走向指向采空区一侧，erf为概率积分函数，其形式为u为积分参数，W₀为该地质采矿条件下的最大下沉值，

为主要影响半径，b为水平移动系数，α为煤层倾角，q为下沉系数，l为走向有限开采时的计算长度，其值为工作面走向长度D₃与走向左、右边界的拐点偏移距s₃、s₄之差，即：l＝D₃-s₃-s₄；

2）当矿区地下开采工作面走向充分开采，倾向有限开采时，造成的地表下沉W⁰(y)和水平移动U⁰(y)的预计公式为：

$\left\{\begin{array}{c}{W}^0\left(y\right)=\frac{W_0}{2}\{\mathrm{erf}\left(\frac{\sqrt{\mathrm{\&pi;}}}{r_1}y\right)-\mathrm{erf}[\frac{\sqrt{\mathrm{\&pi;}}}{r_2}(y-L)\left]\right\}\\ U^0\left(y\right)=W_0(b_1{e}^{-\mathrm{\&pi;}\frac{x^2}{r_1^2}}-b_2{e}^{-\mathrm{\&pi;}\frac{{(x-L)}^2}{r_2^2}})-\mathrm{ct}{\mathrm{g\&theta;}}_0{W}^0\left(y\right)\end{array}\right.---\left(2\right)$

其中，y为地表点的计算纵坐标，纵坐标原点，定义为倾向下山煤壁向采空区或煤柱方向偏移s1后，按照开采影响传播角投影到地表所对应的点，纵坐标轴沿倾向指向采空区一侧，b₁、b₂为下山、上山水平移动系数；

为下山、上山主要影响半径，

分别为下、上山的采深，s₁、s₂分别为下、上山拐点偏移距，开采影响传播角θ₀＝90°-kα，α为煤层倾角，可根据实际测量获取，其中k的取值范围为0.5～0.8，煤层倾角α，其取值范围为0～45°，L为倾向工作面计算长度，其计算表达式为： $L=(D_1-s_1-s_2)\frac{\sin ({\mathrm{\&theta;}}_0+\mathrm{\&alpha;})}{\sin {\mathrm{\&theta;}}_0},$ D₁为工作面倾向斜长；

当工作面倾向未达到充分开采时，利用公式（1）预计走向地表下沉W⁰(x)和水平移动U⁰(x)，求出的形变值需要乘上一个小于1的倾向采动程度系数

W_ym是假设走向已经达到充分开采，利用公式（2）中第一式求出W⁰(y)的最大值W_ym；

当工作面走向未达到充分开采时，倾向地表下沉和水平移动同样利用公式（2）预计，求出的形变值需乘上走向采动程度系数

W_xm是假设倾向已经达到充分开采，利用式（1）中第一式求出W⁰(x)的最大值W_xm；

求取W_xm、W_ym时，x和y的取值范围为矿区采深的正负2-3倍，利用公式（1）和公式（2）求出取值范围中走向和倾向的的最大下沉值，即为W_xm、W_ym。

Images