基于时序TM影像的淮南潘谢矿区采煤沉陷水域动态监测

王启元 肖 武,2 李素萃 张文凯 张 伟

(1.中国矿业大学(北京)土地复垦与生态重建研究所，北京市海淀区，100083；2.浙江大学公共管理学院，浙江省杭州市，310058)

摘 要利用时序TM遥感影像数据，在对比单波段阈值法、谱间关系法、归一化水体指数法和改进归一化水体指数法4种水体信息提取精度的基础上，选择改进归一化水体指数(MNDWI)法并结合GIS空间叠加技术，对潘谢矿区2005-2017年采沉陷水域进行提取并分析其时空演变特征。结果显示，2005-2017年，淮南潘谢矿区采煤沉陷水域面积增至4948.65 hm2，总增量达到4003.52 hm2，年均增长333.63 hm2，年均增长率为20.1%；分析研究区内煤炭开采作为主要贡献因素对沉陷水域面积变化影响及预测。通过线性拟合分析，预计到2020年潘谢矿区沉陷水域面积5968.12 hm2；2021-2030年沉陷水域面积年均增约为350 hm2，至2030其总量将达到9468.12 hm2。利用时序TM遥感影像对高潜水位矿区沉陷水体的时序动态变化监测，分析了其变化原因，为后期采煤沉陷区的土地复垦与治理提供了科学依据。

东部高潜水位矿区煤炭开采引起的地表沉陷积水，严重影响了当地的农业发展和居民生活。国内矿区生态环境方面的研究大多集中在生态环境修复、沉降监测以及环境影响评估等方面，对于矿区采煤沉陷积水的研究主要集中在沉陷积水动态监测以及范围确定等方面。快速且精准获取沉陷水域时空变化信息是高潜水位煤粮复合区进行环境整治与沉陷地复垦工作的前提。与传统手段相比较，遥感监测方法具有快速、精确、宏观性强和低费用等特点，在沉陷水域的长时序动态监测上具有无可比拟的优势。

国内外常见的水体信息提取方法有单波段阈值法、水体指数法、多波段谱间关系法等。Frazier等针对水体在近红外波段表现为强吸收性而植被和干土壤在该波段具有强反射性特点，提出单波段阈值法对水体进行提取。Kloiber等用TM影像的2波段加上3波段大于4波段加上5波段，提取水体信息。McFeeters提出归一化水体指数(NDWI)，该指数提取水体时植被和土壤信息能够被抑制。徐涵秋提出了改进归一化水体指数(MNDWI)来提取水体信息。彭苏萍等基于多时相TM影像，提取了淮南矿区积水沉陷面积的扩展变化信息；张瑞娅等将采区的原始地形归入地表沉陷分析，确定了采煤沉陷的积水区域范围；杨光华等利用高分遥感影像，对济宁市采煤导致的塌陷积水耕地信息用人机交互解译的方式进行了提取测算；孟磊等用TM数据分析淮南潘谢矿区1976-2006年的水体变化情况并分析了其演变的景观生态效应。利用遥感技术手段对水体信息进行提取，必须考虑到在光谱特征方面水体与其他地物的差异性，针对不同的研究区域的不同自然地理条件，在选择水体提取方法时亦有所差异。本文以淮南潘谢矿区为研究区，比较4种遥感水体信息提取方法，结合矿区塌陷范围实测资料及GIS空间叠加分析技术，实现研究区2005-2017年沉陷水域提取，揭示了潘谢矿区近十多年沉陷水域的时空演变特征极，并分析煤炭开采作为主要原因对其监测时段内的影响及趋势预测，对未来采煤沉陷水域变化趋势预测以及生态环境治理具有现实意义。

1 研究区及数据源

1.1 研究区概况

潘谢矿区位于安徽省淮南市西北部，地形平坦开阔，地处淮河冲积平原，整体地势西北高、东南低，地表高程一般为21.5～25.4 m，平均标高在23.1 m左右。所在区域属季风温暖带半湿润气候，季节性表现为夏季炎热，冬季寒冷。6～8月降雨量约占全年的40%，年平均降雨量865.5 mm，年平均蒸发量1610.44 mm，蒸发量大于降雨量，潮湿系数近似0.5。区内沟渠纵横，水利设施完善，当地潜水位埋深约为1.5 m。因此地下煤炭开采沉陷后，地表极易形成积水，如图1所示。截至2014年，淮南矿区累计原煤产量约6.7亿t，根据生产规划，到2020年，矿区累计原煤产量约11亿t；2021-2030年，原煤产量稳定在7000万t/a。

图1 淮南潘榭矿区遥感影像

1.2 数据选择及预处理

为了研究不同水体提取方法在本区域内的精度与适用性，选取研究区内丁集矿、顾桥矿和顾北矿作为实验区域对比分析4种水体提取方法中的最优法，选用2016年9月2日Landsat8 OLI影像作为数据源，然后整体提取潘谢矿区6期沉陷水域信息，见表1。所有影像来源于地理空间数据云，空间分辨率为30 m，宽幅185 km×185 km，影像包含红、绿、蓝、近红外和短波红外等9个波段。由于研究所采用的遥感影像获取时均已经过系统的辐射校正与几何校正等相关预处理工作，本文对数据的预处理工作主要包括辐射定标、大气校正以及影像裁剪等。2005-2017年6期遥感影像日期均在5月份左右，未到研究区多降水季节，降水量对于水体提取结果的影响不予考虑。

表1 研究区遥感数据基本情况

2 研究方法

基于TM影像的水体信息自动提取方法，考虑影像光谱特性的提取方法目前应用最多，其中单波段阈值法、多波段谱间关系法、水体指数法这3种方法最为广泛。本文使用4种水体信息提取方法，选择研究区内丁集、顾桥、顾北3个矿作为实验区，对实验区内所有水体进行提取，经过精度评价找到适合研究区的最优水体信息提取方法，再利用优选后的方法对研究区内2005年、2008年、2009年、2013年、2015年、2017年共6期水体进行分阶段提取，分析煤炭开采对沉陷水域变化的影响以及趋势预测。

2.1 TM影像水体信息提取方法

目前利用遥感影像进行水体信息提取的方法主要有单波段(近红外)阈值法、多波段谱间关系法和水体指数法等。单一的水体信息提取方法无法通用，需要根据研究数据情况以及区域地理条件择优选择。

2.1.1 单波段阈值法

单波段阈值法利用近红外波段水体的强吸收性以及植被和干土壤的强反射性特点，选取适当的阈值将水体信息提取出来。此方法原理简单，但水体与非水体之间存在的过渡区域容易被忽略，同时在区分细小水体和混淆在水体中的阴影上很困难，从而产生误提、漏提的现象。利用水体信息和其他地物信息在遥感图像不同波段上表现出来的不同亮度值来设定阈值，提取水体信息。依据实验，水体在TM5(近红外)波段具备很强的吸收作用，亮度值表现较低，而在该波段其他地物亮度值和水体有明显区别。使用ENVI5.2软件，对2016年9月2日Landsat8 OLI影像近红外波段灰度图像使用密度分割法确定阈值，设定水体提取阈值，利用波段运算工具，提取水体信息。

2.1.2 多波段谱间关系法

波谱间关系法可根据水体与地物波谱曲线的特征，找出其中的变化规律并设定逻辑判别式提取水体，这样能够将水体与阴影较好地区分出来，而且提取精度较高。水体具备TM2(蓝)加TM3(绿)大于TM4(红)加TM5(近红外)的特征，经过试验发现，设定一定的阈值来得到提取水体的效果更好，设定模型如下：

(1)

通过多次反复实验，选取K1及K2值，使用波段运算工具，得到水体信息提取图。

2.1.3 归一化水体指数法

水体的反射率从可见光波段到短红外波段逐级减弱，尤其在近红外和短波热红外范围内显示强烈的吸收性，因此可采用水体在可见光波段和近红外波段(NIR)的光谱差异构建水体指数，而植被在近红外波段的反射率最强，采用绿波段(GREEN)和近红外波段构建指数可以有效地抑制植被和土壤信息，NDWI指数法如下：

(2)

同样对于NDWI指数设定阈值上下限1，K2〗，可以更好地获得水体信息提取效果。

2.1.4 改进归一化水体指数法

土壤、建筑物因素的干扰在NDWI指数中被忽略，而水体绿光和近红外的波段波谱特征与土壤、建筑物几乎一致。徐秋涵等利用中红外波段(MIR)，提出MNDWI指数用于水体信息的提取，计算公式如下：

(3)

2.2 遥感影像与采煤沉陷边界空间叠加

采用水体信息提取方法中精度最高的改进归一化水体指数法分别对研究区内2005-2017年(2005年5月8日、2008年5月15日、2009年6月3日、2013年5月21日、2015年5月19日、2017年4月30日)共6期遥感影像分别进行体提取，在Arcgis软件中，使用空间叠加分析模块，将水体提取结果与淮南潘榭矿区2015年地表沉降范围叠加分析，实现对潘榭矿区采煤沉陷水域的提取。

3 结果与分析

3.1 水体信息提取方法比较及精度分析

经过前期预处理，对实验区OLI遥感影像应用各水体信息自动提取方法，得到提取结果，见表2。

表2 不同方法水体信息提取结果对比(黑色为水体)

表2从左到右依次为单波段阈值法、多波段谱间关系法、NDWI指数法和MNDWI水体指数的水体信息提取结果，黑色部分为水体信息。结合原始影像目视解译分析上述提取结果，MNDWI指数法对水体信息的提取效果最好，对水体斑块完整性及连续性提取上表现出色，并且误提现象最少；单波段阈值法对水体的提取效果次之，虽然也存在误提现象，但在水体信息连续性提取上表现出色，相比与其他两种自动提取方法效果好；水体指数法效果和谱间关系法相对较差，水体信息边界提取不连贯，存在较多误提现象。

使用ENVI5.3软件对实验区遥感原始影像目视解译，均匀选择100组样本，包含像元5562个。经统计得到水体提取混淆矩阵和Kappa统计，见表3。由表3可知，4种水体提取方法精度均达到90%以上，但MNDWI法提取研究区水体信息的用户精度、总体精度以及Kappa相关系数为各方法中最高。

表3 不同水体方法提取结果精度评价

3.2 历年水体信息提取

基于MNDWI指数法，用2005年、2008年、2009年、2013年、2015年和2017年TM/OLI影像潘谢矿区整体水域进行提取，结果如图2所示。

由图2可以看到，2005-2017年，潘谢矿区地表全域水体扩张明显，尤其在2010年之后增速加快，区内水体斑块总量不断增加，个体积水斑块多数由小变大，且斑块形状变化巨大。

图2 监测年份潘谢矿区地表水域分布

3.3 采煤沉陷水域监测

由于本文主要侧重于分析开采沉陷导致的水体变化，而矿区内仍存在大面积的自然水体，为了区分未受扰动的天然水体与开采沉陷水体，实现对研究区内采煤沉陷水域的提取监测，将淮南潘谢矿区2015年采煤沉陷边界与水域提取结果结合，利用Arcgis10.3软件进行空间叠加运算，获得2005-2017年共6期潘谢矿区沉陷区范围内的水域结果，如图3所示，并在SPSS软件中统计分析潘谢矿区监测时段的沉陷水域面积变化情况，见表4。

分析结果表明，淮南潘谢矿区沉陷水域面积由2005年的945.13 hm2增加至2017年的4948.65 hm2，累积增长量为4003.52 hm2，年均增长333.63 hm2，年均增长率为20.1%。沉陷水体总量明显呈上升趋势。其中，2005-2009年沉陷水域面积增加1459.3 hm2，增幅达154%，变化幅度相对最大；2009-2013年沉陷水域面积增加603.32 hm2，增幅为25.1%；2013-2017年沉陷水域面积增量为1940.9 hm2，增幅为64.5%，虽然增幅未至最大，但其沉陷水域面积的绝对增量在监测时段内为最大值。

表4 2005-2017年沉陷水域面积及变化值

综合图3和表4可知，2005年淮南潘谢矿区由于已采矿井数量较少，还未出现大范围的采煤沉陷积水，只有潘一矿、潘二矿、潘三矿、顾北矿、张集矿出现小范围的沉陷积水；2005-2009年，潘一矿、潘二矿、潘三矿、顾北矿、张集矿沉陷水域面积增加，其中谢集矿、张集矿扩张最为明显；2009-2013年期间，除朱集矿外，研究区内各矿沉陷水域面积均有所变化，顾桥矿水域增加明显，2013年其内部产生一大面积沉陷积水区；2013-2017年，研究区境内所有沉陷水域面积均扩有增加，且分布均匀。

图3 2005-2017年潘谢矿区沉陷水域分布

图4 累积原煤产量与沉陷水域面积相关性分析

淮南潘谢矿区2005-2017年累积原煤产量与沉陷水域面积关系如图4所示，二者相关系数r=0.94。随着累积煤炭产量的增加以及矿区沉陷水域面积不断扩大，根据淮南矿业集团生产规划，到2020年潘谢矿区累积原煤产量预计达10.85亿t；2020-2030年，原煤产量稳定在7000万t/a，累计原煤产量达17.85亿t。根据上述拟合结果预计到2020年潘谢矿区沉陷水域面积总量预计增至5968.12 hm2；2021-2030年期间，区内沉陷水域面积年均增量预计为350 hm2，到2030年沉陷水域面积将达到9468.12 hm2。

4 结论

本研究在对比分析各种水体信息方法的基础上，以淮南潘榭矿区为例，对研究区内所有水体进行了提取，结合潘谢矿区2015年采煤沉陷边界，利用GIS空间叠加，实现对研究区2005-2017年地表沉陷水体的动态监测，分析了区内煤炭开采对沉陷水域面积变化的影响并做出预测，得出以下结论：

(1)依据研究所用数据特点以及研究区地理条件，对比了单波段阈值法、谱间关系法、归一化水体指数法及改进归一化水体指数法4种水体信息提取方法的优劣，选择提取效果最好、精度最高的改进归一化水体指数法作为潘谢矿区水体信息提取方法，其总体精度达96.03%，Kappa系数为0.88。

(2)2005-2017年潘谢矿区沉陷水域的监测结果显示，2005-2017年，区内沉陷水域面积增至4948.65 hm2，总增量为4003.52 hm2，年均增长333.63 hm2，年均增长率为20.1%。其中，2005-2009年沉陷水域面积增加1459.3 hm2，增幅达154%，变化幅度相对最大，2013-2017年沉陷水域面积增量为1940.9 hm2，沉陷水域面积的绝对增量在监测时段内为最大值。

(3)通过对监测时段内的累积原煤产量与沉陷水域面积的相关性定量分析可知，随着煤炭开采量的累积增加，潘谢矿区沉陷水域面积随之不断上升，二者具有很强的相关性，相关系数为0.94。依据生产规划和拟合结果，预计到2020年，潘谢矿区沉陷水域面积总量将增至5968.12 hm2，2021-2030年，其沉陷水域面积年均增量预计为350 hm2，到2030年沉陷水域面积将达到9468.12 hm2。