基于物元可拓法的高潜水位煤矿区土地损毁程度评价

李 超，崔瑞豪，张 曦，张曦元，马晋明，王培俊

（1.永煤集团复垦公司，河南 永城 476600；2.中国矿业大学环境与测绘学院，江苏 徐州 221116；3.西北农林科技大学资源环境学院，陕西 杨凌 712100；4.中国矿业大学公共管理学院（应急管理学院），江苏 徐州 221116）

中国煤炭资源95%以上采用井工开采方式，煤炭开采造成大面积的土地下沉，并导致房屋坍塌、基础设施受损、水土流失加重、生态环境恶化，农田遭到了严重的破坏，农民的生产生活受到严重的影响[1-2]。高潜水位采煤沉陷区具有地势平坦、潜水位埋深小、可采煤层数量多、煤层厚度大和地表下沉系数大等特点，致使地表容易塌陷积水，农田减产甚至绝产，土地损毁严重[3-4]。因此，为了更好地保护生态环境，提高复垦效率，相关学者提出了“边采边复”技术[5-6]。由于“边采边复”技术是在地面沉陷前或沉陷过程中采取的复垦措施，从而促进了“末端治理”向“源头和过程控制与治理”的转变[7]，极大地提高了土地复垦率，减少了对土地的破坏。

土地损毁程度评价是复垦方案编制的重要组成部分，为土地复垦提供基础数据，其结果直接影响复垦的效果及受损农民补偿标准的制定[8-9]。国外的煤炭资源大多采用露天矿开采方式，不会产生采煤沉陷地。因此，国外对于采煤沉陷区土地损毁程度评价的研究很少[10]。随着中国对土地复垦工作重视程度不断提升，对采煤沉陷区土地损毁程度评价的研究也越来越多[11-12]。部分学者对评价指标体系及其分级等进行了探讨，常用的评价方法有模糊综合评价法、极限条件法、指数法等[13-16]，并利用GIS 手段优化评价过程[17]，实现评价结果的可视化。但如何客观有效地构建采煤沉陷区土地损毁程度评价指标体系，并完善其评价方法是重中之重。程琳琳等[18]使用可拓评价方法，建立高潜水位采煤沉陷区土地损毁评价指标的形式化模型，通过分析事物变换的规律，降低了多因素评价中的主观片面性，大大提高了评价结果真实性。该方法已广泛应用于其他领域的评价[19]，陈秋计等[20]初步验证了该方法在矿区土地复垦适宜性评价、采煤沉陷区土地损毁程度评价中的可行性，并取得了良好的效果。

在多数采煤沉陷区土地损毁程度评价中，积水情况属于必须要考虑的内容，多数土地损毁评价中对当地积水情况不清楚，对积水范围和积水时间不明确，这严重影响了土地损毁程度评价的准确性。因此，借鉴HE 等[21]、赵艳玲等[22]对华东地区沉陷区积水的识别方法，使用Landsat 8 遥感数据生成年内水频率指数（AWFI），以代表沉陷区的年度积水量，并定量描述积水情况，作为评价基础数据，使得评价结果更为准确。

目前，矿区土地损毁程度评价研究蓬勃发展，但对于东部高潜水位平原矿区的土地损毁评价研究较少。以永城市的城郊矿区和新桥矿区为例，在梳理已有研究的基础上，构建评价指标体系，划分评价单元，采用熵权法确定指标权重，针对高潜水位矿区土地损毁的特点，将极限条件法和可拓评价方法相结合，对过去数年内开采造成的沉陷区土地损毁程度进行评价，并使用无人机影像数据绘制植被长势进行评价验证[23-24]，为该区域土地复垦方案编制、受损农民补偿标准的制定等提供理论支持，为其他高潜水位采煤沉陷区土地损毁程度的评价提供借鉴。

1 研究区概况及数据来源

1.1 研究区概况

永城市属暖温带半湿润半干旱大陆性季风气候，夏季炎热多雨，冬季干燥寒冷，四季气候变化分明。永城矿区整体土地约70%属于耕地，易积水形成湖泊，耕地破坏面积大，沉陷区盐渍化趋势加重。

城郊矿区位于黄淮冲积平原的东部，地势较为平坦开阔，总体地势西北高，东南低，地面标高32～34 m，相对高差2 m 左右。据永城市气象站资料，降雨量多集中于6 月—8 月，占全年降雨量的50%；夏季多东南-南风，冬季多西北-北风，年平均风速3.4 m/s，最大风速20 m/s。井田内地表水系多为纵横交错的人工排涝沟渠，最大河流为沱河，由西北流向东南，横贯全井田。该河流为季节性河流，雨季流量增大，旱季流量剧减。新桥矿区处于华北平原与黄淮冲积平原的交接部位，地势平坦，自西向东略有倾斜。地面标高31～33 m，相对高差2 m 左右。降雨、水系等情况与城郊矿区类似。城郊矿区和新桥矿区地理区位和土地利用分类情况如图1 所示。

图1 城郊矿区与新桥矿区地理区位和土地利用分类情况Fig.1 Geographical location and land use classification of Chengjiao Mining Area and Xinqiao Mining Area

1.2 数据来源

数据资料主要包括基础矿山数据、实地调研数据以及遥感影像数据：①商丘市永城市的清华2017土地利用数据，主要获取农田土地利用位置，以确定试验田和采样点在大致位置；②商丘市永城市2020年地形图数据；③Landsat 8 TM 影像用以绘制年积水频率图；④无人机影像数据；⑤煤炭资源开发利用方案、土地复垦方案、环境监测报告等。

使用的遥感传感器是基于大疆PSDK 定制开发的MS600 PRO 无人机多光谱相机，相机镜头可以获取蓝光波段（波长450 nm、带宽27 nm）、绿光波段（波长555 nm、带宽27 nm）、红光波段（波长660 nm、带宽22 nm）、红边波段（波长710 nm、带宽10 nm）和近红外波段（波长840 nm、带宽15 nm）等波段影像。搭载MS600 PRO 的遥感平台为四旋翼无人机大疆M210-RTK，航向重叠度及旁向重叠度均设置为75%，航高90 m，航速1.9 m/s，三次飞行累计航点数81 个，航线长度约13.6 km，飞行范围约133 749 m2，预计飞行时间67 min，共计获取像片8 388 张，覆盖莲藕、芡实、水稻的主要种植范围。

空间数据包括土地利用数据、DEM（数字高程模型）、矿区分布图及Landsat 卫星影像数据等。表1展示了数据源、分辨率和用途。通过裁剪、光栅化和重采样对这些数据进行预处理，并对收集统计到的数据进行梳理检查，确保单位的规范统一；在QGIS 中实现各专题地图的空间投影坐标系一致。

表1 数据源、分辨率及用途Table 1 Data source, resolution and purpose

2 土地损毁评价

2.1 评价过程

2.1.1 评价指标的选取

根据《农用地质量分等规程》（GB/T 28407—2012）、《土地复垦方案编制规程》（TD/T 1031.3—2011）、相关研究文献和土地复垦方案案例，初步筛选了地表形变的7 个备选指标，即积水情况、下沉深度、裂缝密度、倾斜变形、曲率、水平移动和水平变形等。结合实地考察，剔除多项争议指标。研究区待评价单元主要种植水稻、芡实和莲藕等水生植物，需要从时间角度考虑积水情况，因此加入了时间因素层；东部高潜水位矿区的裂缝与西部露天煤矿的地裂缝情况不同，常出现在水下，且较为细小难以观测，因此去掉裂缝密度指标；但又因为裂缝是实际存在的，所以需要水平位移指标进行表示。最终确定的高潜水位平原型采煤沉陷区土地损毁程度评价指标体系见表2。

表2 采煤沉陷土地损毁程度评价指标体系Table 2 Evaluation index system of damage degree of coal mining subsidence land

2.1.2 指标因子获取

积水时间指标主要通过Landsat 数据生成计算获取，并在实地考察时验证准确性。使用年内水频率指数表示年积水情况，在GEE 平台上通过归一化差分植被指数（NDVI）和增强植被指数（EVI）来降低植被对水体分类的干扰，水体分类主要是使用归一化修正差分水指数（MNDWI）[25]，具体计算见式（1）～式（3）。一个像素点满足MNDWI>NDVI或MNDWI>EVI和EVI<0.1 条件时将被归类为水体，通过计算像素点有水的影像数量与全年的影像数量的值来表示积水情况。相关研究表明水体分为季节性水（0.3≤AWFI<0.75）和永久水（AWFI≥0.75），（0，0.3）范围为未积水或无效数据区。积水结果如图2 所示。

图2 研究区积水分布图Fig.2 Water distribution map of study area

式 中： ρblue、 ρgreen、 ρred、 ρnir和 ρswir分别为Landsat 8 影像的蓝色波段（450～520 nm）、绿色波段（520～600 nm）、红色波段（630～690 nm）、近红外波段（760～900 nm）和短波红外波段（1 550～1 750 nm）。

地表形变的空间指标获取主要通过MSPS 软件获取，如下沉深度、倾斜变形、水平变形和曲率。根据城郊矿区和新桥矿区实测地表移动参数，结合开采实践及上覆岩层物理力学特性，确定矿区内工作面采用概率积分法进行地表移动变形预计的参数为：下沉系数q=0.70、水平移动系数b=0.33、主要影响角正切tgβ=2.0、最大下沉角θ=90°−0.6α、拐点偏移距s=0.05H。地表移动变形预计的参数结果见表3。

表3 工作面开采地表移动变形各参数最大值Table 3 Maximum parameters of surface movement and deformation in working face mining

2.1.3 指标因子分级

1）划分评价等级。依据积水、地形变化和土壤肥力等方面差异，将高潜水位采煤沉陷区土地损毁程度评价的结果统一分为三个等级（表4），即轻度损毁、中度损毁和重度损毁，以便于该区域复垦标准及受损农民补偿标准的制定及评价结果的判断。

表4 损毁程度等级分级Table 4 Classification of damage degree

2）因子的分级标准。依据《农用地质量分等规程》（GB/T 28407—2012）、《土地复垦方案编制规程》（TD/T 1031.3—2011），结合已有研究的土地损毁程度的分级标准，根据研究区的实际情况，得到各因子的损毁程度等级划分标准，见表5。

表5 研究区损毁程度等级划分Table 5 Classification of damage degree in study area

2.1.4 确定指标权重

评价指标权重的确定是建模评价的关键环节，权重赋予合理与否直接关系到评价结果的可信度。为有效避免指标权重赋值的主观随意性，使用熵权法对各评价指标赋予权重，具体步骤见式（4）～式（6）。

熵值是一种物理计量单位，熵值越大说明数据越混乱，携带的信息越少，效用值越小，因而权重也越小。由表6 可知，积水时间熵值最小、权重最高符合一般土地损毁评价中的赋值规律。下沉值的信息熵值最大、信息效用值最小，主要是因为评价单元的下沉值较为接近。

表6 熵权值结果Table 6 The results of entropy weight

2.1.5 采用可拓理论评价损毁程度

1）数据预处理。评价指标对损毁程度的影响有正效应或负效应，即部分指标因子的数值越小损毁等级越高，部分指标数值越大损毁等级越高（如下沉深度），无法直接进行比较，因此将各评价指标进行归一化无量纲处理，定量化的指标可以根据式（7）进行归一化处理。

式中：X为指标的边界值；x¯为对应指标无量纲化后的值；Xmax和Xmin为该指标的最大值和最小值。定性指标需结合经验以及当地的实际资料来赋值，最后按照定量指标的方式进行归一化处理。

2）确定评价指标的节域及经典域。系统节域物元模型见式（8）。

式中：P为采煤沉陷区土地损毁等级评价的全部集合，即评价指标，是对应于Ci（i=1,2,···,n）的全部取值范围，也就是P的节域系统的经典域物元模型。

各指标的分级（j=1,2,3）分别表示轻度损毁、中度损毁和重度损毁，计算见式（9）。

式中，VPi为对应的指标Ci（i=1,2,···,n）的取值范围，经典域为 ⟨aPn,bPn⟩。

3）建立待评价对象的物元模型和关联函数，计算关联度。以物元表示分析得到的评价数据。式（10）为评价单元对应因子的具体数值，是待评价对象。

关联度是事物间或者因素间关联性的量度，这里描述的是各评价指标对于不同损毁程度的影响水平，计算公式见式（11）。

式中：Kji为第i个评价因子隶属于评价等级j的关联度； ρ（vi,vo ji） 为 点vi与区间vo ji的距；的模。其中， ρ（vi,vo ji）可以通过式（12）确定。

式中：y为实域（–∞，+∞）上的一点；Y=（a，b）为实域上的区间。

4）确定损毁等级。损毁等级计算见式（13）。

式中：Kj（Po）为评价单元对不同土地损毁程度等级的符合程度。依据式（14）确定土地损毁程度等级。

2.2 评价单元的划分

评价单元是具有特定土地特性和土地质量的土地评价和制图的基本单元。依据《农用地分等定级规程》，参考已有研究成果，得出常用的评价单元划分方法有图斑法、叠置法、地块法及动态网格法。首先选取叠置法将永城市土地利用现状图、DEM、开采沉陷预计的沉陷范围和矿权边界图进行叠加，其次在此基础上提取所有建筑用地、水体和农用地作为评价对象，用以表示整个采煤沉陷区地物类别，再次将其与年内水频率指数图叠加并使用缓冲区分析得出34 个评价单元，如图3（a）和图3（b）所示。水稻玉米田评价单元编号为1～10、芡实评价单元编号为11～21、莲藕池评价单元编号为22～34。最后将获取的评价单元根据位置状况相互连接绘制成图。

图3 土地利用图和年内水频率图叠加获取的评价单元图Fig.3 Evaluation unit diagram obtained by superposition of land use map and annual water frequency map

3 结果与讨论

3.1 损毁程度结果

芡实地、水稻玉米田和莲藕池各评价单元相应指标的实测值见表7，土地损毁情况如图4～图6 所示。通过上述计算得出芡实地评价单元中11 点位、15 点位、16 点位、19 点位、20 点位、21 点位为中度损毁，其他评价单元为轻度损毁。芡实地总面积为444.7 hm2，中度损毁面积为135.7 hm2，占比30.51%，即有接近三分之一的土地会有季节性积水，地形产生形变，土壤肥力下降。如图4（b）所示芡实地评价单元中中度损毁分为三块，分别是11 点位，位于西北处边界，面积最小只有9.0 hm2；15 点位和16 点位位置接近所以合并单元，其面积为48.6 hm2，位于西南边界处；19 点位～21 点位位于东部，面积最大为78.0 hm2。如图5（b）所示水稻玉米田整体属于轻度损毁状态，总面积为747.2 hm2，属于三个研究区中土地损毁情况最轻的等级，即无积水或短期积水，地形未发生明显改变，土壤肥力无影响。莲藕池总损毁面积为562.6 hm2，中度损毁面积为151.8 hm2，占比26.98%。如图6（b）所示23 点位、24 点位、27 点位、28 点位、31 点位、34 点位为莲藕池的中度损毁地区。其中，23 点位位于莲藕池的东北角，损毁面积为51.0 hm2；24 点位、27 点位和34 点位位于其西南角，总面积为45.7 hm2。

表7 研究区评价指标实测值（部分）Table 7 Measured value of evaluation index in study area（part）

图4 芡实地土地损毁总体情况和放大图Fig.4 Overall situation and enlarged map of land damage in gorgon nut land

图5 水稻玉米田土地损毁总体情况和放大图Fig.5 Overall situation and enlarged map of land damage in rice corn field

图6 莲藕池土地损毁总体情况和放大图Fig.6 Overall situation and enlarged map of land damage in lotus root pond

3.2 结果验证

研究中的土地损毁评价主要通过积水和地形的改变进行判断，因此使用土地肥力条件进行验证。土壤肥力的验证主要体现在农作物长势上，图4～图6 中无人机影像位置与损毁程度评价不同程度类型同时重叠的仅有新桥矿区的莲藕池，因此以其为例进行结果验证。将无人机多光谱影像进行预处理和指数运算，绘制代表莲藕长势的NDVI图，如图7所示。轻度损毁地区即图7（a）和图7（b）中圆点位置，放大后可明显观察出该区域大部分莲藕处于长势良好的状态。中度损毁地区即图7（a）和图7（c）中菱形位置，长势呈现明显的浅色即长势较差。该结果较符合中度损毁定义中土壤肥力下降的直观表现，表明土地损毁程度评价结果较为准确合理。

图7 莲藕长势NDVI 反演结果Fig.7 NDVI inversion results of lotus root growth

3.3 讨论

积水情况是土地损毁评价的重要影响因素，但不是唯一影响因素。对高潜水位矿区土地损毁程度评价，通常的做法是通过估计积水的范围与程度来确定损毁程度。但该经验方法既具有主观片面性与不准确性，又不能定量表示损毁范围，因此需要引入新的评价方法。因此将基于遥感影像的积水时间数据引入评价体系中，从而构建全新的高潜水位采煤沉陷区土地损毁程度评价指标体系，使用基于极限条件法与可拓物元法的评价方法，以及基于熵权法的权重确定方法，提出了采煤沉陷土地损毁指标体系，使评价结果更为客观、科学和合理。将可拓方法应用于矿区土地损毁评价，把各评价因子定量化，尽量减少人为因素的影响，且使用物元可拓模型评价出的土地损毁程度与依据积水情况的经验做法已经有了直观上的差异，克服了多因素识别评价中的主观片面性，大大提高了评价结果真实性、科学性和合理程度。

4 结 论

通过构建的高潜水位平原型采煤沉陷区土地损毁程度评价指标体系方法，对城郊矿区和新桥矿区进行土地损毁程度评价，得到以下结论。

1）研究区整体土地损毁程度较轻，其中，莲藕池中度损毁面积最大，为151.8 hm2，占莲藕池损毁总面积的26.98%；水稻玉米田损毁程度最低，整体属于轻度损毁程度。

2）评价指标体系方法的评价结果与单凭积水情况的评价方法有明显不同，引入遥感等方法后，使得研究区损毁情况更加直观、客观、科学、合理。

3）通过无人机影像绘制植被长势可以直观地区分出长势的不同，进而初步得出土壤肥力不同的区域，即中度损毁区域植被长势差于轻度损毁区域。