CN113836740A - 一种高潜水位采煤沉陷积水区历史空间信息计算方法 - Google Patents

Abstract

本发明的空间信息计算方法将遥感方法与沉陷预计结合，一方面利用遥感方法提取研究区域遥感影像的积水边界，另一方面利用工作面预计参数、地质采矿条件，通过概率积分法预计地表下沉盆地地形。将遥感影像提取的积水边界与预计下沉地形通过坐标进行匹配，确定积水区水面高程，进一步计算沉陷积水区水深，水下地形，积水面积，积水体积，沉陷区库容各项空间信息，绘制水下地形图。相对比传统的沉陷区积水空间信息获取方法，本发明方便快捷，结果精度高，可计算过往积水空间信息数据，可以为采煤沉陷区治理提供可靠的数据支持，提高采煤沉陷区水土资源的利用价值，改善采煤沉陷区生态环境，缓解了采煤沉陷区人地矛盾。

Description

一种高潜水位采煤沉陷积水区历史空间信息计算方法

技术领域

本发明涉及矿区观测技术领域，涉及一种高潜水位采煤沉陷积水区历史空间信息计算方法。

背景技术

长期大规模、高强度的煤炭资源开采，造成土地沉陷、耕地损毁、房屋开裂、基础设施损坏、生态破坏等一系列问题。同时，地下开采引起大规模的地表沉陷后，受大气降水、地下水渗流等因素影响，地表形成大面积的沉陷积水区。中国东部地区矿区具有地下潜水位高、松散层厚、煤层厚的特点，沉陷区积水问题尤为严重。目前针对此类因为采煤形成的积水区，常见的治理模式包括：充填复垦恢复土地、鱼塘改造发展水产养殖、下沉盆地改造为水库、因地制宜修建特色公园等。在生态环境治理时，缺少积水区不同时期的空间数据会造成充填复垦无法确定充填量、改造鱼塘无法确定养殖范围、改造水库无法确定储蓄水量。因此沉陷积水区的三维空间信息在沉陷区治理以及开采沉陷危害评估时起了决定性作用，所以获取沉陷积水区空间信息是矿区观测工作中的重要一环。

现有方法均存在一些问题：

(1)现场实测的问题

现场实测多为利用测量船通过声纳、测深仪、GPS等装置，观测积水区空间信息。现场实测的过程费时费力、花费成本较高，且只能实测积水区实时的空间信息，不能获得积水区以往的历史空间信息数据，无法研究积水区动态变化情况。

(2)单纯采用沉陷预计方法计算积水区空间信息

传统的积水区空间信息的预计方法，是基于地下潜水位高程，结合下沉预计结果，来判断水面在下沉盆地中的位置，从而计算积水区空间信息。但是由于地下潜水位常伴随降水量、蒸发量改变而动态变化，同时潜水位观测站数量少、分布散，所以沉陷区详细潜水位难以确定，只能取整体区域的潜水位高度，导致预测结果精度低。传统方法虽可以通过历史潜水位高度计算过去的沉陷积水区空间信息，但准确度较低，不能完全满足于工程需要。

(3)遥感方法的问题

遥感技术在采煤沉陷区的应用多为通过遥感影像提取积水范围，而难以获取沉陷积水区的水下地形、水深、积水体积等三维空间信息。

发明内容

(一)解决的技术问题

由于高潜水位矿区采煤沉陷积水区空间信息难以获取，而新兴的遥感技术也难以准确获取积水区三维空间信息。针对上述问题，本发明提供一种基于遥感影像与沉陷预计方法融合的、用于高潜水位采煤沉陷积水区历史空间信息的计算方法。

(二)技术方案

一种高潜水位采煤沉陷积水区历史空间信息计算方法，包括如下步骤：

S1，搜集沉陷积水区所在矿区地质采矿条件，包括采煤沉陷区的地层组成、地质构造、水文地质资料、矿区沉陷预计参数和矿区煤炭开采资料；

S2，利用Knothe时间函数结合采煤沉陷预计方法预计不同时刻地表的下沉情况，再通过原始地形计算沉陷后地表高程，并根据沉陷盆地高程结果，生成下沉盆地高程等值线；

S3，获取积水区需要计算空间信息同时期的遥感影像，对遥感影像进行辐射定标、大气校正和图像裁剪预处理；再利用水体指数法进行积水边界提取，得到具有坐标信息的沉陷积水区的积水边界；

S4，将S2计算所得沉陷区地形与S3得到的遥感影像提取的沉陷积水区边界根据坐标信息进行匹配，得到遥感方法提取的积水边界在预计的下沉地形高程等值线中的位置；

S5，判断是否可以直接确定积水边界处高程，即判断提取的积水区边界各个位置在下沉地形等值线中高程是否相等；

S6，若积水边界在下沉地形等值线中的各处高程相等，则可直接确定水面高程；若受地形不规则等因素的影响，提取的积水边界各处高程在下沉地形等值线中不相等，无法确定水面高程，此时重新采用一种边界离散插值法确定积水边界处高程；

S7，确定水面高程后，结合S2中计算得到的沉陷区地形，计算积水区空间信息，包括水深，积水面积，积水体积，库容，水下地形。

优选的，S1中采矿条件搜集的参数包括平均采深H、采厚m、煤层倾角α、预计工作面的走向长L、倾向长l、开采速度v；

其中矿区沉陷预计参数包括下沉系数q、水平移动系数b、主要影响半径r、开采影响传播角θ₀、拐点偏移距s。

优选的，S2中先计算地表不同时刻下沉情况，再利用原始地形高程减去地表下沉计算沉陷后地表高程，计算地表不同时刻下沉情况采用公式：

式中W(x,y,T)为T时刻地表点(x,y)的下沉值，t为工作面开采时间，等于T时刻减去工作面开始开采的时刻；e为自然常数，c是时间比例系数，m为煤层采厚，q为下沉系数，α为煤层倾角，r为主要影响半径，r＝H/tanβ,H为采深，tanβ是主要影响角正切，D是地下开采范围，

是二重积分变量；

计算沉陷后地表高程采用公式：

H(x,y,T)＝H₀(x,y,T)-W(x,y,T)；

式中H₀(x,y,T)为地表点(x,y)原始高程，W(x,y,T)为T时刻地表点(x,y)下沉值，H(x,y,T)是T时地表点(x,y)的高程；

其中，下沉盆地高程等值线的单位为m。

优选的，S3中遥感方法提取沉陷积水区的积水边界，先选取遥感影像，经过辐射定标、大气校正和图像裁剪预处理，采用改进的归一化差异水体指数法即MNDWI，获取沉陷积水区的积水边界，公式为：

式中ρ_Green为绿光波段的反射值、ρ_MIR为短波红外波段的反射值。

优选的，S6中，边界离散插值法确定积水边界处高程是将遥感方法提取的沉陷区积水边界随机离散为n个点，利用预计的下沉地形结果利用克里金插值法对离散点的高程进行插值，得到各离散点的高程，计算各积水边界离散点高程的平均值作为积水边界处高程，即水面高程；

其中，计算各积水边界离散点高程的平均值采用公式：

式中

为各插值点处的高程值，

为T时刻计算后得到的积水边界高程值。

优选的，S7中计算沉陷积水区内任意点的水深，得到水下地形，采用公式：

式中h(x_i,y_i,T)为T时地表任意点(x_i,y_i)位置水深，

为T时刻水面高程，H(x_i,y_i,T)是T时地表任意点(x_i,y_i)高程。

优选的，S7中根据S3中得到的沉陷盆地高程等值线，计算各高程等值线截面面积，采用公式：

确定水面高程后，计算水面面积，采用公式：

上述式中

为T时刻水面面积，

为积水区内地表坐标。

优选的，S7中利用S3得到的沉陷区高程地形，根据T时刻的水面高程选取沉陷积水区域，将积水区域划分为一定间距的方格区域，从地形图得到每个方格的水深，对划分的每个方格按四角水深取平均值；将沉陷积水区分成N个方格区域，利用方格区域的面积乘以水深，即为该区域的积水体积，方格宽度越小，计算精度越高；将每个方格区域的积水体积空间积分得到积水区的体积，计算公式为：

式中

为积水体积，

为T时刻剖分的各个方格区域面积，

为T时刻各子区域的水深。

优选的，S7中利用沉陷区高程地形，将沉陷区域划分为一定间距的方格区域，从地形图得到每个方格的高程，方格按四角高程取平均值；将沉陷区分成M个方格区域，利用方格区域的面积乘以方格高程到目标高程的差值，即为该方格区域的蓄水体积，目标高程即蓄水的水位；

通过对每个方格区域的蓄水体积空间积分得到沉陷区的蓄水体积，即库容，计算公式为：

式中V_库为沉陷区库容，H(x_i,y_i,T)是T时地表任意点(x_i,y_i)高程，

为T时刻下沉盆地剖分的方格区域的面积，H_e是目标高程、即蓄水水位。

(三)有益效果

本发明提供了一种高潜水位采煤沉陷积水区历史空间信息计算方法，具有以下优点：

(1)不需要现场实测，方便快捷，成本低；

(2)相比传统计算方法，提高了计算结果的准确性；

(3)可以计算出历史不同时期的沉陷积水区空间信息；

(4)为高潜水位采煤沉陷区积水区的土地复垦规划和评估等研究提供科学、准确的数据支持与指导。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的，保护一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图；

图2为实施例中沉陷盆地动态下沉预计图；

图3为实施例中动态下沉预计叠加原始地形计算下沉盆地地形图；

图4为实施例中遥感影像提取的积水边界图；

图5为实施例中积水边界与下沉地形匹配叠加图；

图6为实施例中积水边界随机离散图；

图7为实施例中重新确定的水面高程在下沉盆地地形中的位置图；

图8为实施例中利用计算所得空间信息结果绘制的地形图。

具体实施方式

下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本发明的描述中，需要说明的是，如出现术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等，其所指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的设备或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，如出现术语“第一”、“第二”、“第三”，其仅用于描述目的，而不能理解为指示或暗示相对重要性。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，如出现术语“安装”、“相连”、“连接”，应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

实施例

以某矿工作面地表沉陷积水区为例计算2020年8月中旬积水区空间信息。

S1，获取示例区工作面信息：工作面1走向长1124m，倾向长为175m，平均采深为400m、工作面平均采厚m＝4.1m，煤层倾斜角a＝8°，开采时间为2013年10月—2015年2月；工作面2走向长1116m，倾向长127m，平均采深为400m、工作面平均采厚m＝3m，煤层倾斜角a＝11°，开采时间为2016年10月—2017年8月；工作面3走向长426m，倾向长为90m，平均采深为350m、工作面平均采厚m＝3.7m，煤层倾斜角a＝10°，开采时间为2015年3月—2015年6月；工作面4走向长276m，倾向长为132m，平均采深为420m、工作面平均采厚m＝3.7m，煤层倾斜角a＝10°，开采时间为2015年8月—2015年11月。获取矿区的预计参数：下沉系数q＝0.98、主要影响角正切tanβ＝2.16、水平移动系数b＝0.31，拐点偏移距为10m，影响传播角86.3°。

S2，根据所获取数据信息对工作面进行基于Knothe时间函数的开采沉陷动态预计，获得沉陷盆地动态下沉预计值，如图2。工作面原始地表平均高程为27.42m，叠加原始地形计算下沉盆地地形，得到高程等值线图，如图3。

S3，获取沉陷积水区2020年8月21日landsat8遥感影像，经过预处理后利用水体指数法方法提取积水边界如图4所示。

S4，将提取的积水边界于下沉地形匹配叠加，如图5。

S5、S6，将S2中得到的高程等值线与S3中得到的积水边界经过坐标匹配比对，可以看到积水边界高程不相等，因此重新确定积水边界处高程。

先将积水边界随机离散为150个离散点，如图6，并利用下沉地形结果利用克里金插值法进行插值，确定各个离散点的高程；

再利用公式计算水面高程；结果为25.64m，在下沉盆地地形中的位置如图7；

S7，确定水面高程后，利用公式计算积水区空间信息。其中最大水深为2.474m，积水区面积为464439m²，积水体积为529078m³，沉陷区库容为1750871m³；

利用计算所得空间信息结果绘制2020年8月中旬示例区地形图，如图8。

在本说明书的描述中，参考术语“一个实施例”、“示例”、“具体示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上公开的本发明优选实施例只是用于帮助阐述本发明。优选实施例并没有详尽叙述所有的细节，也不限制该发明仅为的具体实施方式。显然，根据本说明书的内容，可作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地理解和利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。

Claims (9)

1.一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，包括如下步骤：

S1，搜集沉陷积水区所在矿区地质采矿条件，包括采煤沉陷区的地层组成、地质构造、水文地质资料、矿区沉陷预计参数和矿区煤炭开采资料；

S2，利用Knothe时间函数结合采煤沉陷预计方法预计不同时刻地表的下沉情况，再通过原始地形计算沉陷后地表高程，并根据沉陷盆地高程结果，生成下沉盆地高程等值线；

S3，获取积水区需要计算空间信息同时期的遥感影像，对遥感影像进行辐射定标、大气校正和图像裁剪预处理；再利用水体指数法进行积水边界提取，得到具有坐标信息的沉陷积水区的积水边界；

S4，将S2计算所得沉陷区地形与S3得到的遥感影像提取的沉陷积水区边界根据坐标信息进行匹配，得到遥感方法提取的积水边界在预计的下沉地形高程等值线中的位置；

S5，判断是否可以直接确定积水边界处高程，即判断提取的积水区边界各个位置在下沉地形等值线中高程是否相等；

S6，若积水边界在下沉地形等值线中的各处高程相等，则可直接确定水面高程；若受地形不规则等因素的影响，提取的积水边界各处高程在下沉地形等值线中不相等，无法确定水面高程，此时重新采用一种边界离散插值法确定积水边界处高程作为水面高程；

S7，确定水面高程后，结合S2中计算得到的沉陷区地形，计算积水区空间信息，包括水深，积水面积，积水体积，库容，水下地形。

2.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S1中采矿条件搜集的参数包括平均采深H、采厚m、煤层倾角α、预计工作面的走向长L、倾向长l、开采速度v；

其中矿区沉陷预计参数包括下沉系数q、水平移动系数b、主要影响半径r、开采影响传播角θ₀、拐点偏移距s。

3.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S2中先计算地表不同时刻下沉情况，再利用原始地形高程减去地表下沉计算沉陷后地表高程，计算地表不同时刻下沉情况采用公式：

式中W(x,y,T)为T时刻地表点(x,y)的下沉值，t为工作面开采时间，等于T时刻减去工作面开始开采的时刻；e为自然常数，c是时间比例系数，m为煤层采厚，q为下沉系数，α为煤层倾角，r为主要影响半径，r＝H/tanβ,H为采深，tanβ是主要影响角正切，D是地下开采范围，

是二重积分变量；

计算沉陷后地表高程采用公式：

H(x,y,T)＝H₀(x,y,T)-W(x,y,T)；

式中H₀(x,y,T)为地表点(x,y)原始高程，W(x,y,T)为T时刻地表点(x,y)下沉值，H(x,y,T)是T时地表点(x,y)的高程；

其中，下沉盆地高程等值线的单位为m。

4.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S3中遥感方法提取沉陷积水区的积水边界，先选取遥感影像，经过辐射定标、大气校正和图像裁剪预处理，采用改进的归一化差异水体指数法即MNDWI，获取沉陷积水区的积水边界，公式为：

式中ρ_Green为绿光波段的反射值、ρ_MIR为短波红外波段的反射值。

5.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S6中，边界离散插值法确定积水边界处高程是将遥感方法提取的沉陷区积水边界随机离散为n个点，利用预计的下沉地形结果利用克里金插值法对离散点的高程进行插值，得到各离散点的高程，计算各积水边界离散点高程的平均值作为积水边界处高程，即水面高程；

其中，计算各积水边界离散点高程的平均值采用公式：

式中

为各插值点处的高程值，

为T时刻计算后得到的积水边界高程值。

6.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S7中计算沉陷积水区内任意点的水深，得到水下地形，采用公式：

式中h(x_i,y_i,T)为T时地表任意点(x_i,y_i)位置水深，

为T时刻水面高程，H(x_i,y_i,T)是T时地表任意点(x_i,y_i)高程。

7.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S7中根据S3中得到的沉陷盆地高程等值线，计算各高程等值线截面面积，采用公式：

确定水面高程后，计算水面面积，采用公式：

上述式中

为T时刻水面面积，

为积水区内地表坐标。

8.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S7中利用S3得到的沉陷区高程地形，根据T时刻的水面高程选取沉陷积水区域，将积水区域划分为一定间距的方格区域，从地形图得到每个方格的水深，对划分的每个方格按四角水深取平均值；将沉陷积水区分成N个方格区域，利用方格区域的面积乘以水深，即为该区域的积水体积，方格宽度越小，计算精度越高；将每个方格区域的积水体积空间积分得到积水区的体积，计算公式为：

式中

为积水体积，

为T时刻剖分的各个方格区域面积，

为T时刻各子区域的水深。

9.根据权利要求1所述的一种高潜水位采煤沉陷积水区历史空间信息计算方法，其特征在于，S7中利用沉陷区高程地形，将沉陷区域划分为一定间距的方格区域，从地形图得到每个方格的高程，方格按四角高程取平均值；将沉陷区分成M个方格区域，利用方格区域的面积乘以方格高程到目标高程的差值，即为该方格区域的蓄水体积，目标高程即蓄水的水位；

通过对每个方格区域的蓄水体积空间积分得到沉陷区的蓄水体积，即库容，计算公式为：

式中V_库为沉陷区库容，H(x_i,y_i,T)是T时地表任意点(x_i,y_i)高程，

为T时刻下沉盆地剖分的方格区域的面积，H_e是目标高程、即蓄水水位。

Images