CN106226764A - 一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法 - Google Patents

Abstract

本发明公开了一种基于D‑InSAR的煤矿开采地沉陷区域的测定方法，该方法以D‑InSAR技术为基础，对影像进行一系列地解译分析，得到地表移动形变图、下沉等值线图、地表三维形变图，求取地表移动盆地的角量参数和概率积分法预计参数。与传统测量方法相比，本发明的基于D‑InSAR的煤矿开采地沉陷区域的测定方法得到的距离向移动形变精度更高，效率高、成本低、监测周期短等优势，可在一定程度上取代传统测量方法。

Description

一种基于D-InSAR的煤矿开采地沉陷区域的测定方法

技术领域

本发明涉及的是沉陷区域测定方法的技术领域，尤其涉及的是一种基于D-InSAR的煤矿开采地沉陷区域的测定方法。

背景技术

目前，概率积分法预计参数获取主要有两种方法：①通过实测地表移动资料反演预计参数；②在没有实测资料可借鉴的情况下，参照临近矿区或规程上的预计参数经验值。概率积分法预计参数值选取的准确与否直接关系到矿区地表移动和变形预计的精度。

合成孔径雷达差分干涉测量(Differential Interferometric SAR,D-InSAR)是利用合成孔径雷达(Synthetic Aperture Radar,SAR)获取的相位信息提取地表形变信息的一门技术。D-InSAR具有全天候、全天时的工作特点和快速、高精度、大区域对地测量的技术优势,可以探测到地震、冰川漂移、地面沉降等自然灾害引起的地表微小形变,精度可达厘米级至毫米级。

发明内容

本发明的目的在于克服现有技术的不足，提供了一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，以提供一种精度高、效率高、成本低、监测周期短的煤矿开采地沉陷区域的测量方法，该方法能够测量煤矿开采地沉陷区的角量参数和概率积分法预计参数。

本发明是通过以下技术方案实现的：

一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，包括以下步骤：

(1)采集待研究的煤矿开采地沉陷区域地表形变发生前后的两幅SAR影像，将该两幅SAR影像先进行干涉处理，然后利用外部DEM数据模拟该研究区域的地形相位并剔除模拟地形相位，最后进行相位噪声处理，获得差分干涉图；

(2)对差分干涉图先进行相位解缠，然后将差分干涉图中的相位形变转换成地表形变，最后利用SRTM3DEM数据对差分干涉图进行地理编码，得到在大地坐标系中的最终形变图；

(3)将最终形变图与采煤工作面一同叠加入ArcGIS软件中进行形变分析，得到地表形变等值线图、采煤工作面沉陷盆地走向线剖面图、沿采煤工作面走向线的下沉曲线、沿采煤工作面倾向线的下沉曲线、三维地表形变图以及地表形变等高线图；

(4)地表形变等值线图、采煤工作面沉陷盆地走向线剖面图、沿采煤工作面走向线的下沉曲线、沿采煤工作面倾向线的下沉曲线、三维地表形变图以及地表形变等高线图计算研究区域地表形变的参数，所述地表形变的参数包括角量参数和概率积分法预计参数，但不限于此，其中，所述角量参数包括最外边界和边界角，所述概率积分法预计参数包括主要影响角、下沉系数、拐点偏距和水平移动系数。

所述角量参数的计算方法为：

A、最外边界：地表形变等值线图中的最外层等值线即为待测区域的地表移动盆地的最外边界；

B、边界角：边界角：在采煤工作面走向线的下沉曲线上，利用倾斜、曲率和水平变形值的公式，分别求出地表移动盆地的最外边界上倾斜变形值为3mm/m、曲率变形值为0.2mm/m²、水平变形值为2mm/m的点，取走向线与最外边界的交点为最外边一点，计算最外边一点至开采边界的连线与水平线在煤柱一侧的夹角，即为边界角。

进一步地，上述步骤B中，倾斜、曲率和水平变形值的计算方法为：

$i_{n-n+1}=\frac{w_{n+1}-w_n}{l_{n-n+1}},mm/m$

式中，i_n-n+1为n号点至n+1号点的倾斜变形值，l_n-n+1为n号点至n+1号点的水平距离；w_n+1,w_n分别为n+1号点和n号点的下沉量；

$K_{n+1~n~n-1}=\frac{i_{n+1~n}-i_{n~n-1}}{\frac{l_{n+1~n}+l_{n~n-1}}{2}}=\frac{2(i_{n+1~n}-i_{n~n-1})}{l_{n+1~n}+l_{n~n-1}},mm/m^2$

式中，K_{n+1～n～n-1}表示n+1、n、n-1号点之间的曲率变形值，i_n+1～n、i_n～n-1分别表示n+1号点至n号点和n号点至n-1号点的倾斜变形值；l_n+1～n、l_n～n-1分别表示n+1号点至n号点和n号点至n-1号点的水平距离；

${\mathrm{\ϵ}}_{n+1-n}=\frac{{\left(l_{n+1~n}\right)}_m-{\left(l_{n+1~n}\right)}_0}{{\left(l_{n+1~n}\right)}_0},mm/m$

式中，ε_n+1-n表示n+1号点至n号点的水平变形值，(l_n+1～n)₀、(l_n+1～n)_m分别表示n+1号点至n号点在首次观测时和m次观测时的水平距离。

所述概率积分法预计参数的计算方法为：

a、主要影响角：

步骤一、主要影响半径r为由沿采煤工作面走向线的下沉曲线的拐点到最大下沉点的水平距离，主要影响半径的值随采样深度的变化而变化，主要影响半径的确定，可由地表形变等值线图和采煤工作面沉陷盆地走向线剖面图中求取；

步骤二：计算主要影响角的正切值tanβ＝H/r，其中，H表示采样深度；

b、下沉系数：

$q=W_{cm}/(M\·cos\∂)$

式中，q为下沉系数，W_cm是沉陷盆地的最大下沉值，可根据沿采煤工作面走向线或沿采煤工作面倾向线的下沉曲线的最大下沉点得到，M是煤层的厚度，取煤层的法线采厚，是煤层倾角；

c、拐点偏距：将沿采煤工作面走向线的下沉曲线的拐点投影到煤层上得到计算边界，量取实际开采边界与计算边界之间沿煤层的距离，即为拐点偏距；

d、水平移动系数：

$b=\frac{U_{cm}}{W_{cm}}$

式中，b表示水平移动系数，U_cm表示最大水平移动值，U_cm是通过SARscape中的MAI功能获取的方位向的形变，即水平位移量。

本发明相比现有技术具有以下优点：本发明提供了一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，该方法空间分辨率高、检测快速，能周期性提供连续空间覆盖的地面沉降信息，观测精度达到厘米甚至毫米级，可监测地表微小形变。该方法还可反演某一矿区的概率积分法预计参数，且具有周期短、效率高、精度高、可获得全盆地沉降信息的优点。

附图说明

图1为地表形变等值线图；

图2为采煤工作面沉陷盆地走向线剖面图；

图3为沿采煤工作面走向线的下沉曲线；

图4为沿采煤工作面倾向线下沉曲线；

图5为三维地表形变图；

图6为地表形变等高线图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例1

本实施例提供的一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，包括以下步骤：

一、影像解算

二轨法是利用研究区域地表形变发生前后的两幅SAR影像生成干涉图，然后利用外部DEM(数字高程模型)数据模拟该区域的地形相位，并从干涉图中剔除模拟的地形相位得到研究区域的地表形变相位信息，本发明利用二轨法获得地表形变相位信息，具体包括以下步骤：

(1)采集研究区域地表形变发生前后的两幅SAR影像，将该两幅SAR影像先进行干涉处理，然后利用外部DEM数据模拟该研究区域的地形相位并剔除模拟地形相位，最后进行相位噪声处理，获得差分干涉图，本实施例中，SAR影像采用的是RS2-3m超精细影像，DEM数据采用美国NASA(美国国家航空航天局)发布的SRTM-3(航天飞机雷达制图任务)DEM数据；

(2)对差分干涉图先进行相位解缠，然后将差分干涉图中的相位形变转换成地表形变，最后利用SRTM3DEM数据对差分干涉图进行地理编码，得到在大地坐标系中的最终形变图；

(3)将最终形变图与采煤工作面一同叠加入ArcGIS软件中进行形变分析，得到如图1-6所示的地表形变等值线图、采煤工作面沉陷盆地走向线剖面图、沿采煤工作面走向线的下沉曲线、沿采煤工作面倾向线的下沉曲线、三维地表形变图以及地表形变等高线图，通过三维地表形变图可以明显地看出由于采煤引起的沉降漏斗。

二、求解角量参数

边界角是地表移动盆地最外边界至采空区边界的连线与水平线在煤柱一侧的夹角，一般地表移动盆地的最外边界是由下沉为10mm的点所圈定的边界，本发明中，可通过地表形变等值线图求出边界角，具体包括以下步骤：

A：求地表移动盆地的最外边界：

地表移动盆地的最外边界即下沉10mm的点所圈定的边界，本发明中，地表形变等值线图中的最外层等值线即为待测区域的地表移动盆地的最外边界；

B：求边界角：

目前，我国采用的临界形变值包括：倾斜、曲率和水平变形分别为i＝3mm/m、K＝0.2mm/m²、ε＝2mm/m的临界形变值，根据临界值，在采煤工作面走向线的下沉曲线上，利用倾斜i、曲率K和水平ε变形值的公式，分别求出地表移动盆地的最外边界上i＝3mm/m、K＝0.2mm/m²、ε＝2mm/m的点，取走向线与最外边界的交点为最外边一点，计算最外边一点至开采边界的连线与水平线在煤柱一侧的夹角，即为边界角；

倾斜i、曲率K和水平ε变形值的公式为：

$i_{n-n+1}=\frac{w_{n+1}-w_n}{l_{n-n+1}},mm/m$

式中，i_n-n+1为n号点至n+1号点的倾斜变形值，l_n-n+1为n号点至n+1号点的水平距离；w_n+1,w_n分别为n+1号点和n号点的下沉量；

$K_{n+1~n~n-1}=\frac{i_{n+1~n}-i_{n~n-1}}{\frac{l_{n+1~n}+l_{n~n-1}}{2}}=\frac{2(i_{n+1~n}-i_{n~n-1})}{l_{n+1~n}+l_{n~n-1}},mm/m^2$

式中，K_{n+1～n～n-1}表示n+1、n、n-1号点之间的曲率变形值，i_n+1～n、i_n～n-1分别表示n+1号点至n号点和n号点至n-1号点的倾斜变形值；l_n+1～n、l_n～n-1分别表示n+1号点至n号点和n号点至n-1号点的水平距离；

${\mathrm{\ϵ}}_{n+1-n}=\frac{{\left(l_{n+1~n}\right)}_m-{\left(l_{n+1~n}\right)}_0}{{\left(l_{n+1~n}\right)}_0},mm/m$

式中，ε_n+1-n表示n+1号点至n号点的水平变形值，(l_n+1～n)₀、(l_n+1～n)_m分别表示n+1号点至n号点在首次观测时和m次观测时的水平距离。

三、求预计参数：

a、求主要影响角正切值：

步骤一、主要影响半径r为由沿采煤工作面走向线的下沉曲线的拐点到最大下沉点的水平距离，主要影响半径的值随采样深度的变化而变化，主要影响半径的确定，可由图1的地表形变等值线图和图2的采煤工作面沉陷盆地走向线剖面图中求取，例如，从图1的地表形变等值线图及图2的采煤工作面沉陷盆地走向线剖面图中可判断出2013.11.13-12.07沉陷盆地的主要影响半径为312m；

步骤二：计算主要影响角的正切值tanβ＝H/r，其中，H表示采样深度，r为主要影响半径。

b：计算下沉系数q：

$q=W_{cm}/(M\·cos\∂)$

式中，W_cm是沉陷盆地的最大下沉值，可根据沿采煤工作面走向线的下沉曲线的最大下沉点得到，M是煤层的厚度，取煤层的法线采厚，是煤层倾角。

c：求拐点偏距S：

将走向线下沉曲线(如图3)的拐点投影到煤层上得到计算边界，量取实际开采边界与计算边界之间沿煤层的距离，即为拐点偏距S。

确定实测下沉曲线的拐点：拐点处的下沉值为W_cm/2。

d：求水平移动系数b：

水平移动系数b反映地表最大水平移动值U_cm与最大下沉值W_cm的比例关系，即：

$b=\frac{U_{cm}}{W_{cm}}$

式中，最大水平移动值U_cm是通过SARscape中的MAI功能获取的方位向的形变，即水平位移量，W_cm表示沉陷盆地的最大下沉值。

四、精度评定：

SAR是侧视雷达系统，成像几何属于斜距投影类型。由于合成孔径雷达影像在距离向和方位向具有完全不同的几何特征，故将其成像几何特征分离开来。根据成像几何特征，距离向形变主要体现地表的下沉，方位向形变表征在地表水平移动变形，水平移动变形更为微小、复杂。SAR影像的距离分辨率和方位分辨率由于成像卫星、波段、模式的不同存在较大差异。

(1)根据地表形变等值线图，地表移动盆地中下沉10mm的点是可以准确找到，因为影像采用的是3m分辨率的超精细影像，所以下沉10mm的边界线的确定误差在3m以内，相较与传统测量30m的监测点而言，精度得到极大提高；

(2)由高精度的移动盆地的边界点，所求的边界角、移动角等角值参数精度也必定会提高；

(3)下表1为本发明的基于D-InSAR技术获得的相关参数与传统测量获得的相关参数的精度对比分析：

表1两种观测方式获取的相关参数的精度分析

从表1中可以看出，本发明在煤矿开采地沉陷区域的测定中求得的参数的精度高于传统测量方法，说明本发明所述方法可在一定程度上取代传统测量方法。

Claims (5)

1.一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，其特征在于，包括以下步骤：

(1)采集待研究的煤矿开采地沉陷区域地表形变发生前后的两幅SAR影像，将该两幅SAR影像先进行干涉处理，然后利用外部DEM数据模拟该研究区域的地形相位并剔除模拟地形相位，最后进行相位噪声处理，获得差分干涉图；

(2)对差分干涉图先进行相位解缠，然后将差分干涉图中的相位形变转换成地表形变，最后利用SRTM3DEM数据对差分干涉图进行地理编码，得到在大地坐标系中的最终形变图；

(3)将最终形变图与采煤工作面一同叠加入ArcGIS软件中进行形变分析，得到地表形变等值线图、采煤工作面沉陷盆地走向线剖面图、沿采煤工作面走向线的下沉曲线、沿采煤工作面倾向线的下沉曲线、三维地表形变图以及地表形变等高线图；

(4)地表形变等值线图、采煤工作面沉陷盆地走向线剖面图、沿采煤工作面走向线的下沉曲线、沿采煤工作面倾向线的下沉曲线、三维地表形变图以及地表形变等高线图计算研究区域地表形变的参数。

2.根据权利要求1所述的一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，其特征在于，所述步骤(4)中，地表形变的参数包括地表移动盆地的角量参数和概率积分法预计参数，所述角量参数包括最外边界和边界角，所述概率积分法预计参数包括主要影响角、下沉系数、拐点偏距和水平移动系数。

3.根据权利要求2所述的一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，其特征在于，所述步骤(4)中，角量参数的计算方法为：

A、最外边界：地表形变等值线图中的最外层等值线即为待测区域的地表移动盆地的最外边界；

B、边界角：在采煤工作面走向线的下沉曲线上，利用倾斜、曲率和水平变形值的公式，分别求出地表移动盆地的最外边界上倾斜变形值为3mm/m、曲率变形值为0.2mm/m²、水平变形值为2mm/m的点，取走向线与最外边界的交点为最外边一点，计算最外边一点至开采边界的连线与水平线在煤柱一侧的夹角，即为边界角。

4.根据权利要求3所述的一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，其特征在于，所述步骤B中，倾斜、曲率和水平变形值的计算方法为：

$i_{n-n+1}=\frac{w_{n+1}-w_n}{l_{n-n+1}},mm/m$

式中，i_n-n+1为n号点至n+1号点的倾斜变形值，l_n-n+1为n号点至n+1号点的水平距离；w_n+1,w_n分别为n+1号点和n号点的下沉量；

$K_{n+1~n~n-1}=\frac{i_{n+1~n}-i_{n~n-1}}{\frac{l_{n+1~n}+l_{n~n-1}}{2}}=\frac{2(i_{n+1~n}-i_{n~n-1})}{l_{n+1~n}+1_{n~n-1}},mm/m^2$

式中，K_{n+1～n～n-1}表示n+1、n、n-1号点之间的曲率变形值，i_n+1～n、i_n～n-1分别表示n+1号点至n号点和n号点至n-1号点的倾斜变形值；l_n+1～n、l_n～n-1分别表示n+1号点至n号点和n号点至n-1号点的水平距离；

${\mathrm{\ϵ}}_{n+1-n}=\frac{{\left(l_{n+1~n}\right)}_m-{\left(l_{n+1~n}\right)}_0}{{\left(l_{n+1~n}\right)}_0},mm/m$

式中，ε_n+1-n表示n+1号点至n号点的水平变形值，(l_n+1～n)₀、(l_n+1～n)_m分别表示n+1号点至n号点在首次观测时和m次观测时的水平距离。

5.根据权利要求2所述的一种基于D-InSAR的煤矿开采地沉陷区域的测定方法，其特征在于，所述步骤(4)中，概率积分法预计参数的计算方法为：

a、主要影响角：

步骤一、在叠加进ArcGIS中的D-InSAR影像图的地表形变等值线图和采煤工作面沉陷盆地走向线剖面图中获取主要影响半径；

步骤二：计算主要影响角的正切值tanβ＝H/r，其中，H表示采样深度，r为主要影响半径；

b、下沉系数：

$q=W_{cm}/(M\·\cos \∂)$

式中，q为下沉系数，W_cm是沉陷盆地的最大下沉值，可根据沿采煤工作面走向线的下沉曲线的最大下沉点得到，M是煤层的厚度，取煤层的法线采厚，是煤层倾角；

c、拐点偏距：将沿采煤工作面走向线的下沉曲线的拐点投影到煤层上得到计算边界，量取实际开采边界与计算边界之间沿煤层的距离，即为拐点偏距；

d、水平移动系数：

$b=\frac{U_{cm}}{W_{cm}}$

式中，b表示水平移动系数，U_cm表示最大水平移动值，U_cm是通过SARscape中的MAI功能获取的方位向的形变，即水平位移量。