CN109901240A - 一种探测露头残煤自燃区域的方法 - Google Patents

Abstract

本发明构建了一种探测露头残煤自燃区域的方法，该方法首先获取研究区的地质图和遥感影像。然后，对研究区的地质图进行矢量化，从中提取煤层分布区；对多期遥感影像进行预处理，再利用决策树方法对预处理后影像进行分类，利用分类后变化检测方法提取露头煤和裸地增加的区域；利用温度反演算法对预处理后遥感影像进行温度反演，并对地表温度进行分级，从中提取地表高温异常区；最后将煤层分布区、露头煤和裸地增加区、高温异常区三者进行空间叠置分析，提取三者重叠区域作为露头残煤自燃区域。该方法可以快速、准确地提取煤自燃区域，对煤自燃区域生态修复、残煤回收提供科学依据和技术支撑。

Description

一种探测露头残煤自燃区域的方法

技术领域

本发明涉及矿山生态环境监测的技术领域，尤其是涉及一种探测露头残煤自燃区域的方法。

背景技术

露头残煤主要是采掘揭露、煤层天然露头产生的，暴露在空气中，引起煤自燃现象。露头残煤自燃不但造成了大量煤炭资源的浪费，而且对生态环境造成了巨大的破坏，对矿区生活健康和生产安全造成了威胁。对煤自燃区域生态环境进行监测及修复，需将煤自燃区域提取出来。

遥感技术具有覆盖面广、成本低、获取便捷便利快捷的特点，因此被广泛应用于煤自燃方面的探测。如，《一种基于地表裂隙和温度的煤火勘探方法》，发明公开专利号为CN102176070B；《--一种煤田自燃的遥感监测方法》，发明公开专利号为CN104598874A；《一种基于卫星遥感的煤矿火灾监测方法》，发明公开专利号为CN102074093A；《热红外与雷达遥感联合探测煤火的方法》，发明公开专利号为CN108983309A。这些发明均是基于煤自燃的温度异常原理，适合大型的天然煤田火区。相对煤田火区，露头残煤自燃区规模小、分布零散，且伴随采掘扰动，地表有显著的地表覆盖变化特征，这给传统的温度异常探测带来了不利影响。对于露头残煤自燃区域的精准探测，目前仍缺乏有效方法。

发明内容

为满足上述需求，本发明提供了一种探测露头残煤自燃区域的方法，为解决上述问题提供了一种可行性方案，为进行矿区生态环境修复，保证矿区的生活健康和生产安全提供了有力的技术支撑。

本发明采用的技术方案是一种探测露头残煤自燃区域的方法，包括以下步骤：

步骤1)获取研究区的煤田地质图，对煤田地质图进行矢量化，从矢量化的煤田地质图中选取为煤层的区域，该区域为研究区的煤层分布区，记为A；

步骤2)获取研究区的多期遥感影像，对影像进行预处理，得到预处理后影像；利用决策树方法对预处理后影像进行分类，得到分类后影像；再利用分类后变化检测方法提取露头煤和裸地增加的区域，记为B；其过程主要包括步骤2-1)至步骤2-3)；

步骤2-1)对遥感影像进行预处理，包括几何校正、辐射校正、图像融合及图像增强；

步骤2-2)利用决策树方法对预处理后影像进行分类；所述决策树方法的分类指标包括高程、坡度、归一化植被指数及光谱特征；利用决策树方法将预处理后影像分为林地、草地、耕地、建筑用地、裸地和露头煤六种地表覆盖类型；

步骤2-3)利用分类后变化检测方法提取露头煤和裸地增加的区域；所述分类后变化检测方法的具体步骤是：首先选取两个相邻时期的分类后影像，记为f和g，然后分别从f和g中提取出地表覆盖类型为露头煤和裸地的区域，记为f₁和g₁，其他区域为f₂和g₂；其次，叠置f₁和g₁，从g₁中扣除g₁与f₁重合区域，得到两个相邻时期露头煤和裸地增加的区域，记为B；

步骤3)利用温度反演算法对预处理后遥感影像进行温度反演，得到地表温度影像；再利用自然间隔点分割法将地表温度影像划分为低、较低、中、较高、高五个级别；提取级别为高的区域作为高温异常区，记为C，其过程主要包括步骤3-1)至步骤3-4)；

步骤3-1)将预处理影像的像元DN值转转换为辐射亮度值，转换公式为：

式中L_λ为辐射亮度值，单位为W·m^-2·sr^-1·μm^-1，DN为像元灰度值，L_max，L_min分别为相应波段灰度值分别为255和1时所对应的光谱辐射率，单位为W·m^-2·sr^-1·μm^-1；

步骤3-2)将辐射亮度值转换为亮度温度值，转换公式为：

T_rad＝K₂/ln(K₁/L_λ+1) (公式2)

式中T_rad为亮度温度值，K₁为常数，单位为W·m^-2·sr^-1·μm^-1；K²为常数，单位为K；

步骤3-3)求取地表比辐射率；自燃区域，即地表覆盖类型为裸地、露头煤类型的区域，地表比辐射率等于常数0.970；非自燃区域，即地表覆盖类型为林地、草地、耕地、建筑用地的区域，地表比辐射率的计算公式为：

式中ε_b为地表覆盖类型为建筑用地的地表比辐射率；ε_n为自然表面，即地表覆盖为林地、草地、耕地的地表比辐射率；P_V为植被覆盖度，计算公式为：

P_V＝[(NDVI-NDVI_min)/(NDVI_max-NDVI_min)]² (公式4)

式中NDVI为归一化植被指数，NDVI_max是NDVI的最大值，NDVI_min是NDVI的最小值NDVI；当像元的NDVI大于NDVI_max时，P_V取值为1；当NDVI小于NDVI_min时，P_V取值为0；

步骤3-4)求取地表温度，计算公式为：

式中，T_s为地表温度，λ为热红外波段的中心波长， c为光速，h为普朗克常数；σ为波尔兹曼常数；ε为地表比辐射率；

然后，利用自然间隔点分割法确定地表温度的自然间断点，利用自然间断点将地表温度分为低、较低、中、较高、高五个级别；提取地表温度级别为高的区域作为高温异常区，记为C；

步骤4)最后将煤层分布区A、露头煤和裸地增加区B、高温异常区C三者进行空间叠置分析，提取A、B、C三者重叠的区域作为露头残煤自燃区域。

本发明的技术效果在于：

(1)利用露头残煤自燃区域地质构造、地表覆盖变化、热异常特点，实现对露头残煤自燃区域的精准探测，排除了地面采掘扰动、露头残煤规模小、露头残煤分布零散等不利因素，保证了露头残煤自燃区域提取结果的可靠性。

(2)利用遥感影像实现露头残煤自燃区域的精准探测，不但成本较低，耗费时间较短，而且覆盖范围较广，可节约大量人力、物力。

(3)本发明最终提供露头残煤燃烧中心的地图，可以为煤自燃的环境监测、露头煤的开发和管理、生态精准修复提供直接支撑。

附图说明

图1为本发明的技术流程图；

图2为实施例的2002年和2006年的露头煤和裸地增加区图；

图3为实施例的2006年和2010年的露头煤和裸地增加区图；

图4为实施例的2002年份地表温度分级结果图；

图5为实施例的2006年份地表温度分级结果图；

图6为实施例的2010年份地表温度分级结果图；

图7为实施例的将2002-2006和2006-2010裸地和露头煤增加区域，2002年、2006年、2010年高温异常区叠加在一起的图。

具体实施方式

本发明中研究区为有露头残煤燃烧区域，遥感影像数据来源较广泛，多光谱及高光谱数据均可，例如：Landsat、QuickBird、sentinel、Worldview、SPOT等数据。在本实施例中，以Landsat影像为例，并结合图1进行详细描述，具体步骤如下：步骤1：选定车道沟、炭窑坪、石炭沟、寺儿沟四个采矿区作为研究区，四个采矿区的相对位置如图2和图3所示。获取研究区的煤田地质图，利用ArcGIS软件对煤田地质图进行矢量化，从矢量化的煤田地质图中选取为煤层的区域，并重采样至与遥感影像同一尺度，该区域为研究区的煤层分布区，记为A。

步骤2：获取研究区的多期遥感影像，对影像进行预处理，得到预处理后影像。利用决策树方法对预处理后影像进行分类，得到分类后影像。再利用分类后变化检测方法提取露头煤和裸地增加的区域，记为B，其过程综合使用ENVI与ArcGIS软件完成，主要包括步骤2-1至步骤2-3。

步骤2-1，利用ENVI软件对遥感影像进行预处理，包括几何校正、辐射校正、图像融合及图像增强。

步骤2-2，利用决策树方法对预处理后影像进行分类。所述决策树方法的分类指标包括高程、坡度、归一化植被指数及光谱特征。利用决策树方法将预处理后影像分为林地、草地、耕地、建筑用地、裸地和露头煤六种地表覆盖类型。决策树方法的分类过程是：设定高程阈值H1，归一化植被指数阈值NDVI1、NDVI2、NDVI3，坡度阈值s1，第五波段光谱阈值t1。对所有遥感影像的像元进行判断：当坡度小于H1、且归一化植被指数小于NDVI1时，该像元为耕地；当坡度小于H1、且归一化植被指数大于NDVI1时，该像元为建筑用地；当坡度大于H1、归一化植被指数小于t2、且第五波段光谱值小于t1时，该像元为露头煤；当坡度大于H1、归一化植被指数小于t2、且第五波段光谱值大于t1时，该像元为裸地；当坡度大于H1、归一化植被指数大于t2、且归一化植被指数大于t3时，该像元为林地；当坡度大于H1、归一化植被指数大于t2、归一化植被指数小于t3时、且坡度小于s1时，该像元为耕地；当坡度大于H1、归一化植被指数大于t2、归一化植被指数小于t3时、且坡度大于s1时，该像元为草地。

步骤2-3，利用分类后变化检测方法提取露头煤和裸地增加的区域。所述分类后变化检测方法的具体步骤是：首先选取两个相邻时期的分类后影像，记为f和g。然后利用ArcGIS软件分别从f和g中提取出地表覆盖类型为露头煤和裸地的区域，记为f₁和g₁，其他区域为f₂和g₂。其次，叠置f₁和g₁，从g₁中扣除g₁与f₁重合区域，得到两个相邻时期露头煤和裸地增加的区域，记为B。本实施例中考虑了两组相邻时期，分为为2002和2006年、2006和2010年。图2和图3中分别给出了2002年、2006年以及2006年、2010年的矿区内露头煤和裸地区域。2002年与2006年的露头煤和裸地叠加(如图2所示)、2006年与2010年的露头煤和裸地叠加(如图3所示)，叠加区域就是相邻时期露头煤和裸地增加的区域。

步骤3：利用温度反演算法对预处理后遥感影像进行温度反演，得到地表温度影像。再利用自然间隔点分割法将地表温度影像划分为低、较低、中、较高、高五个级别。提取级别为高的地表温度影像作为高温异常区，记为C，其过程主要包括步骤3-1至步骤3-4。

步骤3-1，利用经过几何校正的影像像元的DN值计算辐射亮度值。计算公式为：

式中L_λ为辐射亮度值，单位为W·m^-2·sr^-1·μm^-1，DN为像元灰度值，L_max，L_min分别为相应波段灰度值分别为255和1时所对应的光谱辐射率，单位为W·m^-2·sr^-1·μm^-1。针对Landsat影像，对于TM第6波段：L_max＝15.303W·m^-2·sr^-1·μm^-1，L_min＝1.238W·m^-2·sr^-1·μm^-1；对于ETM+第6波段低增益波段L_max＝17.040W·m^-2·sr^-1·μm^-1，L_min＝0，高增益波段L_max＝12.650W·m^-2·sr^-1·μm^-1，L_min＝3.200W·m^-2·sr^-1·μm^-1。将其带入公式1中，计算TM/ETM+第6波段辐射亮度值的公式为：

对于TM6，公式为：

L_λ＝0.05537DN+1.1826 (公式6)

对于ETM+6，公式为：

L_λL＝0.067DN-0.067

L_λH＝0.0372DN+3.16 (公式7)

式中L_λh，L_λL分别为低增益和高增益波段的辐射亮度。

步骤3-2，利用辐射亮度值计算亮度温度值，计算公式如下：

T_rad＝K₂/ln(K₁/L_λ+1) (公式2)

式中L_λ为上面求取的辐射亮度值，T_rad为量度温度值，K₁为常数，单位为W·m^-2·sr^-1·μm^-1。K₂为常数，单位为K。对于TM数据，K₁＝607.76W·m^-2·sr^-1·μm^-1，K₂＝1260.56K；对于ETM+数据K₁＝666.09W·m^-2·sr^-1·μm^-1，K₂＝1282.71K。

步骤3-3，计算地表比辐射率。自燃区域，即地表覆盖类型为裸地、露头煤类型的区域，地表比辐射率等于常数0.970。非自燃区域，即地表覆盖类型为林地、草地、耕地、建筑用地的区域，地表比辐射率的计算公式为：

式中ε_b为地表覆盖类型为建筑用地的地表比辐射率；ε_n为自然表面，即地表覆盖为林地、草地、耕地的地表比辐射率；P_V为植被覆盖度，计算公式为：

P_V＝[(NDVI-NDVI_min)/(NDVI_max-NDVI_min)]² (公式4)

式中NDVI为归一化植被指数，NDVI_max是NDVI的最大值，NDVI_min是NDVI的最小值NDVI。当像元的NDVI大于NDVI_max时，P_V取值为1；当NDVI小于NDVI_min时，P_V取值为0。

步骤3-4，求取地表温度，计算公式为：

式中，T_s为地表温度，λ为热红外波段的中心波长， c为光速，h为普朗克常数；σ为波尔兹曼常数；ε为地表比辐射率。

然后，利用自然间隔点分割法确定地表温度的自然间断点，利用自然间断点将地表温度分为低、较低、中、较高、高五个级别。提取地表温度级别为高的区域作为高温异常区，记为C。图4、图5和图6为实施例的研究区域2002、2006、2010三个年份地表温度分级结果，图中1、2、3、4、5分别代表低、较低、中、较高、高五个级别。

步骤4：利用ArcGIS软件对煤层分布区A、露头煤和裸地增加区B、高温异常区C三者进行空间叠置分析，提取A、B、C三者重叠区域作为露头残煤自燃区域。图7中将2002-2006裸地和露头煤增加区域，2006-2010裸地和露头煤增加区域、2002年高温异常区、2006年高温异常区、2010年高温异常区叠加在一起，重叠区域即是露头残煤自燃区域，在露头残煤自燃区域标记了燃烧中心，如图7所示。该图件可以为露头残煤的开发、环境影响评价、土地复垦提供支撑。

Claims (3)

1.一种探测露头残煤自燃区域的方法，其特征在于，包括以下步骤：

步骤1)获取研究区的煤田地质图，对煤田地质图进行矢量化，从矢量化的煤田地质图中选取为煤层的区域，该区域为研究区的煤层分布区，记为A；

步骤2)获取研究区的多期遥感影像，对影像进行预处理，得到预处理后影像；利用决策树方法对预处理后影像进行分类，得到分类后影像；再利用分类后变化检测方法提取露头煤和裸地增加的区域，记为B；

步骤3)利用温度反演算法对预处理后遥感影像进行温度反演，得到地表温度影像；再利用自然间隔点分割法将地表温度影像划分为低、较低、中、较高、高五个级别；提取级别为高的区域作为高温异常区，记为C；

步骤4)最后将煤层分布区A、露头煤和裸地增加区B、高温异常区C三者进行空间叠置分析，提取A、B、C三者重叠的区域作为露头残煤自燃区域。

2.根据权利要求1所述的一种探测露头残煤自燃区域的方法，其特征在于，步骤2)中提取露头煤和裸地增加的区域B的方法包括以下步骤：

步骤2-1)对遥感影像进行预处理，包括几何校正、辐射校正、图像融合及图像增强；

步骤2-2)利用决策树方法对预处理后影像进行分类；所述决策树方法的分类指标包括高程、坡度、归一化植被指数及光谱特征；利用决策树方法将预处理后影像分为林地、草地、耕地、建筑用地、裸地和露头煤六种地表覆盖类型；

步骤2-3)利用分类后变化检测方法提取露头煤和裸地增加的区域；所述分类后变化检测方法的具体步骤是：首先选取两个相邻时期的分类后影像，记为f和g，然后分别从f和g中提取出地表覆盖类型为露头煤和裸地的区域，记为f₁和g₁，其他区域为f₂和g₂；其次，叠置f₁和g₁，从g₁中扣除g₁与f₁重合区域，得到两个相邻时期露头煤和裸地增加的区域，记为B。

3.根据权利要求1所述的一种探测露头残煤自燃区域的方法，其特征在于，步骤3)中提取高温异常区C的方法包括以下步骤：

步骤3-1)将预处理影像的像元DN值转转换为辐射亮度值，转换公式为：

式中L_λ为辐射亮度值，单位为W·m^-2·sr^-1·μm^-1，DN为像元灰度值，L_max，L_min分别为相应波段灰度值分别为255和1时所对应的光谱辐射率，单位为W·m^-2·sr^-1·μm^-1；

步骤3-2)将辐射亮度值转换为亮度温度值，转换公式为：

T_rad＝K₂/ln(K₁/L_λ+1) (公式2)

式中T_rad为亮度温度值，K₁为常数，单位为W·m^-2·sr^-1·μm^-1；K₂为常数，单位为K；

步骤3-3)求取地表比辐射率；自燃区域，即地表覆盖类型为裸地、露头煤类型的区域，地表比辐射率等于常数0.970；非自燃区域，即地表覆盖类型为林地、草地、耕地、建筑用地的区域，地表比辐射率的计算公式为：

式中ε_b为地表覆盖类型为建筑用地的地表比辐射率；ε_n为自然表面，即地表覆盖为林地、草地、耕地的地表比辐射率；P_V为植被覆盖度，计算公式为：

P_V＝[(NDVI-NDVI_min)/(NDVI_max-NDVI_min)]² (公式4)

式中NDVI为归一化植被指数，NDVI_max是NDVI的最大值，NDVI_min是NDVI的最小值NDVI；当像元的NDVI大于NDVI_max时，P_V取值为1；当NDVI小于NDVI_min时，P_V取值为0；

步骤3-4)求取地表温度，计算公式为：

式中，T_s为地表温度，λ为热红外波段的中心波长， c为光速，h为普朗克常数；σ为波尔兹曼常数；ε为地表比辐射率；

然后，利用自然间隔点分割法确定地表温度的自然间断点，利用自然间断点将地表温度分为低、较低、中、较高、高五个级别；提取地表温度级别为高的区域作为高温异常区，记为C。