本发明涉及生态地质环境保护领域，尤其是涉及一种基于煤炭资源开发的生态地质环境类型划分方法。
背景技术：
煤炭资源是一种重要的自然资源，同时也是许多工业诸如钢铁、水泥、化工等能源及材料的基础来源，占据中国一次性能源消费结构的比重超过70％。随着中国东部地区煤炭资源逐渐枯竭，煤炭工业发展的重点迅速向中国西部地区转移。在未来10年，中国西部五省包括山西、陕西、内蒙、宁夏、新疆的煤炭产量将超过中国煤炭总产量的70％。然而，中国西部地区多年平均降雨量稀少，蒸发量巨大，属于干旱-半干旱地区，水资源贫乏，生态环境脆弱。且近几年来，随着该区煤炭资源大规模开采，尤其埋深较浅的首采煤层的开采，带来了一系列的矿山环境地质问题，地下水位下降格外凸显，产生了井泉干涸，地表径流量减少，河流两岸的流域生态严重破坏等直接后果，致使生态地质环境质量一直下降。尤其近些来，生态地质环境问题引起了人民和国家高度重视，为此，在2014年国家“973”计划项目指南中将“我国西部生态脆弱区煤炭科学规模开发与水资源保护”列入了“能源科学领域”资助研究方向之一。
生态地质环境是研究地质环境与生态的关系，包括各种地质体、地质过程、环境变化、生物效应和生物活动(主要是人类活动)对地质环境的影响。在中国西部干旱半干旱生态脆弱区，大规模煤炭开采活动将对潜水含水层水资源的赋存产生重大影响，由于采煤引起地表开裂和沉陷，可能发生严重的漏水，导致潜水水位大幅下降。潜水水位的下降将进一步影响地表植被，因为植物将无法吸收潜水含水层的水分。结果，如果潜水水位持续下降，生态地质环境可能恶化。因此，地表径流和松散砂层潜水是连接地质环境与生态环境的桥梁，具有重要的生态功能。然而，不同类型的生态地质环境对采煤活动的敏感性也是不同的。生态地质环境较差的地区对煤矿开采活动的敏感性较差，而生态地质环境较好的地区对即使小规模的采矿活动都具有很强的敏感性。因此，根据不同的生态地质环境特征分类的生态地质环境类型是必要的。这样的分类可以为保护宝贵的潜水含水层水资源，维护原本脆弱的生态地质环境，以及为矿区规划及开采方式选择等工作提取基础依据，对实现干旱-半干旱区生态环境保护开采具有重要意义。
现阶段国内外的技术多集中于采煤活动完成后造成的对原有的地表地质环境或生态环境破坏后的监测、评价及修复措施，缺少在实施采煤活动前，对将开采区域内地表不同的地质环境与生态环境相结合考虑，进行不同生态地质环境类型的区分，以根据不同的生态地质环境条件，为具体的煤炭资源开采活动提供，以达到既实现对煤炭资源的开发，又可以尽可能的减少对地表生态地质环境的破坏，并为之后的矿区地表生态地质环境的修复治理奠定必要的基础，实现煤炭资源开发与生态地质环境保护的协调发展。
影响生态地质环境的因素翻多复杂，并且各个因素之间相互关联、相互影响，各个因素对生态地质环境的整体影响作用不同，加之影响生态地质环境的各个因素大都具有数据模糊性和评价标准模糊性等特点。因此，利用模糊数学的理论与方法，以arcgis和matlab为计算平台，采用构建划分生态地质环境类型的层次结构模型，对西部干旱半干旱生态脆弱地区的生态地质环境划分不同的类型。
对划分结果产生的权重系数的计算有：客观法和主观法。客观法主要有以下几种：熵权法、主成分分析法、均方差法等；主观法主要有以下几种：直接打分法、专家打分法、层次分析法、环比评分法、对比排序法等。然而，生态地质环境类型划分中涉及的相关因素指标多为没有确切数值只有模糊性，一般无法满足客观法的计算要求。对于这一类问题，由于层次分析法作为一种把定性分析和定量分析相结合的系统分析方法可以把复杂的问题层次化，将定性条件定量化，多广泛采用主观法中的层次分析法，然而传统的层次分析法中要求进行一致性检验，但一致性检验困难，不允许决策者出现较大的不一致的情况，其实从行为决策分析的角度来分析，应该允许决策者出现较大的不一致的情况。模糊德尔菲层次分析法是综合层次分析法、模糊评价原理以及德尔菲群体决策法，是一种能够让决策者充分参与权重确定和分析的决策方法，形成一个交互式的权重向量决策分析过程，最终确定出决策者满意的群体决策权重向量.这种决策交互的过程可以在层次结构的任意单准则下进行，并且这种方法允许决策者做出不合理的判断，判断矩阵不需要一致性检查。因此，对于划分指标需要进行更加精确的定量化评价，才能为生态地质环境的合理划分提供更加精确的科学依据。
聚类的定义是由everitt在1974年提出的，他指出聚类就是以某种方式把数据分成指定数目的簇类，最终使得同一簇类中的元素非相似性尽可能的小，不同簇类之间的元素的非相似性尽可能的大。工程中许多问题的解决都用到了聚类方法，像统计学、图像处理等等，这些聚类算法主要有基于模型的聚类算法、划分式聚类算法、层次聚类算法等等，每种算法都有各自的特点，工程实际问题的多样化和复杂化决定了没有一种算法可以解决所有问题，但随着计算机技术的发展，数字计算和程序实现已经越来越容易操作，所以基于目标函数的聚类方法得到了更加深远的发展和推广，模糊聚类就属于该类算法，它是以k均值聚类为基础，引入了模糊理论，而在模糊c均值聚类算法中加入各属性的权重，就形成了属性加权模糊c均值聚类算法，该方法更为科学准确。
技术实现要素：
鉴于上述的分析，本发明旨在提供一种基于煤炭资源开发的生态地质环境类型划分方法，为保护宝贵的潜水含水层水资源，维护原本脆弱的生态地质环境，以及为矿区规划及开采方式选择等工作提取基础依据，对实现干旱-半干旱区生态环境保护开采具有重要意义。
本发明的目的主要是通过以下技术方案实现的：
一种基于煤炭资源开发的生态地质环境类型划分方法，包括以下步骤：
步骤一、获取区域生态、水文、地质资料；
步骤二、建立生态地质环境类型划分的层次结构模型；
步骤三、根据步骤一获取的资料和步骤二中建立的层次结构模型，选取影响生态地质环境的相关因素作为划分指标，并获取待划分区域内生态地质环境类型划分层次结构模型中所有参与类型划分的划分指标所对应的生态、水文及地质数据；
步骤四、将步骤三中所获取划分指标的相关数据转换成浮点型数据；
步骤五、利用归一化函数对步骤四中所述浮点型数据进行无量纲化处理；
步骤六、采用模糊德尔菲层次分析法分析计算出各个划分指标的权重系数；
步骤七、将步骤五中无量纲化数据和步骤六中所述权重系数结合，利用加权模糊c均值聚类方法对影响因素进行叠加聚类计算；
步骤八、根据步骤七中的聚类计算结果和各划分指标的生态、水文及地质特征进行分析判别，确定不同生态地质环境类型，得到生态地质环境类型分区图。
进一步地，步骤二所述层次结构模型包括目标层和指标层，所述的目标层为生态地质环境类型划分的总目标，所述指标层为所有参与类型划分的指标。
进一步地，步骤五用于无量纲化处理的归一化函数为：
式中，fi为每个划分指标中第i个无量纲化处理后的数据，a和b分别为归一化范围的下限和上限，每个划分指标中都有n个数据，xi是每个划分指标中第i个无量纲化前的原始数据，max(xi)和min(xi)为各划分指标原始数据的最大值和最小值。
进行无量纲化处理能够去除量纲对后面步骤中聚类计算的影响。
进一步地，所述归一化范围的下限a＝0，所述归一化范围的上限b＝1。
进一步地，步骤六具体为：利用模糊德尔菲层次分析法，通过向有关生态、水文、地质方面的专家咨询，并结合t.l.saaty1-9标度法对各划分指标进行相对生态地质环境整体重要性评分，建立群体的模糊判断矩阵，确定群体模糊权重向量，最后单准则权重分析计算出各个划分指标的权重系数。
进一步地，步骤六具体包括以下步骤：
步骤6.1、设有m个要判断的划分指标以及n个相关领域的咨询专家，通过德尔菲专家调查法，相关领域咨询专家在某个准则下对指标层中的划分指标相对目标层的相对重要性程度的打分，第k个专家对第i个划分指标fi以及第j个划分指标fj两个划分指标之间的相对重要程度判断bij·k，其中i＝1,2,……m，j＝1,2,……m，k＝1,2……n，确定第k个专家的两两比较判断矩阵b(k)＝[bij·k]；
其中，bij·k＝pi·k/pj·k，pi·k为第k个专家对第i个划分指标相对于目标层重要性的打分值；pj·k为第k个专家对第j个划分指标相对于目标层重要性的打分值；
步骤6.2、构建用三角模糊数表示全部相关领域咨询专家的群体两两模糊判断矩阵c：
c＝[αij,βij,γij]＝[b1b2…bm]
式中，所述判断矩阵由αij，βij，γij三个计算元素组成，其中i＝1……m，j＝1……m，αij≤βij≤γij，αij,βij,γij∈[1/9,1]∪[1,9]，所述计算元素αij，βij和γij由下式确定：
αij＝min(bij·k),k＝1,2,...,n，
γij＝max(bij·k),k＝1,2,...,n，
其中，k＝1,2……n，n为相关领域咨询专家的总数，min(bij·k)为全部相关领域咨询专家打分结果的最小值，geomean(bij·k)为全部相关领域咨询专家打分结果的几何平均数，max(bij·k)为全部相关领域咨询专家打分结果的最大值；
由此构造的全部相关领域咨询专家的群体两两模糊判断矩阵：
步骤6.3、对于所有划分指标中的任意一个划分指标fi，计算群体模糊权重向量过程中涉及的过程计算向量ri：
则确定任意一个划分指标fi群体模糊权重向量为：
式中，符号
关于三角模糊数运算关系说明：
设a＝[a1,a2,a3]和b＝[b1,b2,b3]为两个正三角模糊数，根据三角模糊数理论：
其中a1,a2,a3和b1,b2,b3分别为两组任意实数。
步骤6.4、对于任意一个划分指标fi的群体模糊权重向量为：
其中，
则任意一个划分指标fi指标的权重系数wi经归一化处理后为：
进一步地，步骤七包括以下步骤：
步骤7.1、给定包含n个d维向量数据的待聚类样本集合x，x＝{x1,x2,x3,……xn}，将样本集合分成c个簇类gi(i＝1,…,c)，i为第i个簇类，从样本数据中随机选取c个数据点作为初始的聚类中心，xk＝{xk1,xk2,xk3,…,xkd}t∈rd(k＝1,…c)，xkj为数据点xk的第j维属性的赋值，给定加权指数m、目标函数迭代终止阈值ε和迭代终止最大次数l的值；
步骤7.2、计算各个样本内的数据点与聚类中心的加权欧式距离dw-ij；
步骤7.3、计算每个样本内的数据相对于每个簇类的隶属度；
步骤7.4、计算新的聚类中心矩阵p；
步骤7.5、重复步骤7.2、7.3和7.4，对于各个样本指标内每个数据点，当第t次的迭代计算出新的聚类中心矩阵p(t)与第t+1次迭代计算出新的聚类中心矩阵p(t+1)的差值小于给定的迭代终止阈值ε，即||p(t+1)-p(t)||≤ε，或迭代次数达到给定的最大次数l时，停止计算。
进一步地，步骤7.1中，加权指数m＝2；迭代终止阈值ε取值0.001到0.01。
进一步地，步骤7.2包括以下步骤：
步骤7.2.1、包含n个样本数据点xk(k＝1,…,n)的样本集合x＝{x1,x2,x3,…,xn}，分成c个簇类gi(i＝1,…,c)，从各样本数据点xk(k＝1,…,n)中任意选出c个数据点作为每个簇类的初始簇中心，xk＝{xk1,xk2,xk3,…,xkd}t∈rd(k＝1,…c)，其中xkj为数据点xk的第j维属性的赋值，分别计算每个样本内各数据点到初始簇类中心ci(i＝1,…c)的距离，计算各样本内数据点到初始簇类中心的误差平方和；
步骤7.2.2、对每个样本内数据点和初始簇类中心的欧式距离dki＝||xk-ci||乘以在步骤6.4中计算得出的权重系数wi加以修正，则：
欧式距离
加权欧式距离dw-ij＝d||xj-ci||w＝[(xj-ci)tw2(xj-ci)]1/2，
其中，权重向量w由步骤6.4所述权重系数wi组成，所述权重向量w＝[w1,w2,…,wi]t,(i＝1……d)，所述权重向量中权重系数wi需满足下式：
wi≥0,i＝{1,2,…,d}且
进一步地，步骤7.3包括以下步骤：
步骤7.3.1、新的评价聚类性能的误差平方和准则函数，即新的加权目标函数为：
其中，
步骤7.3.2、利用拉格朗日乘子法求解，构造出的新的拉格朗日函数为：
式中，u为模糊加权划分矩阵，p为新的聚类中心矩阵，uij是第j个数据点对簇类gi的聚类隶属度，ci是相应的模糊向量集的聚类中心，λj为n个约束式的拉格朗日乘子；
结合约束条件
步骤7.3.3、一个数据点对某一簇类的归属是按照隶属度最大原则来确定的，所述数据点归属于隶属度最大的簇类，表达式为：
本发明有益效果如下：
本发明基于煤炭资源开发的生态地质环境类型划分方法，该方法是将我国西北煤炭资源丰富、生态地质环境却十分脆弱的干旱半干旱地区划分不同的生态地质环境类型，并绘制出生态地质环境类型分区图。为保护宝贵的潜水含水层水资源，维护原本脆弱的生态地质环境，以及为矿区规划及开采方式选择等工作提取基础依据，对实现干旱-半干旱区生态环境保护开采具有重要意义。
本发明能够根据现有的生态水文地质资料，快速有效地划分出不同的生态地质环境类型，确定不同类型的生态地质环境的生态地质特征及其对煤炭资源开采活动的敏感性，从而为保护类似干旱半干旱地区珍惜的潜水资源，维护好脆弱的生态环境同时选择恰当的采煤方法实现煤炭资源的开发利用提供科学依据，对西北生态环境脆弱区保水采煤具有重要意义。
本发明对将开采区域内地表不同的地质环境与生态环境相结合考虑，进行不同生态地质环境类型的区分，以根据不同的生态地质环境条件，为具体的煤炭资源开采活动提供，以达到既能够实现对煤炭资源的开发，又可以尽可能的减少对地表生态地质环境的破坏，并为之后的矿区地表生态地质环境的修复治理奠定必要的基础，实现煤炭资源开发与生态地质环境保护的协调发展。
本发明中，上述各技术方案之间还可以相互组合，以实现更多的优选组合方案。本发明的其他特征和优点将在随后的说明书中阐述，并且，部分优点可从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过说明书、权利要求书以及附图中所特别指出的内容中来实现和获得。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图1为本发明方法实施流程图；
图2为待划分区域生态地质环境类型划分的层次结构模型；
图3-为生态地质环境类型中植被指数专题图；
图4为生态地质环境类型中地表高程专题图；
图5为生态地质环境类型中地形坡度专题图；
图6为生态地质环境类型中地表岩性专题图；
图7为生态地质环境类型中地貌类型专题图；
图8为生态地质环境类型中水系河网影响程度专题图；
图9为生态地质环境类型中植被指数归一化专题图；
图10为生态地质环境类型中地表高程归一化专题图；
图11为生态地质环境类型中地形坡度归一化专题图；
图12为生态地质环境类型中地表岩性归一化专题图；
图13为生态地质环境类型中地貌归一化专题图；
图14为生态地质环境类型中水系河网影响程度归一化专题图；
图15为生态地质环境类型区划图。
具体实施方式
下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。
以下参照附图1，并举具体实施例对本发明作进一步的说明。
如图1所示，为一种基于煤炭资源开发的生态地质环境类型划分方法，包括如下步骤：
1、收集区域水文、地质、水文地质资料；
2、建立生态地质环境类型划分的层次结构模型，包括目标层和指标层，所述的目标层为生态地质环境类型划分的总目标，所有的参与类型划分的指标作为指标层；
3、根据步骤1获取的资料及步骤2建立的层次结构模型，选取影响生态地质环境的相关因素作为划分指标，并获取待划分区域内生态地质环境类型划分层次结构模型中所有参与类型划分的指标所对应的生态、水文及地质数据；
4、将步骤3所获取划分指标的相关数据在arcgis中处理成matlab软件能够读取的浮点型.flt数据；
5、在matlab中利用归一化函数将步骤4获得的划分指标的浮点型数据进行无量纲化处理，去除量纲对后面步骤中聚类计算的影响；
6、利用模糊德尔菲层次分析法，通过向有关生态、水文、地质方面的专家咨询，并结合t.l.saaty1-9标度法对各划分指标进行相对生态地质环境整体重要性评分，建立群体的模糊判断矩阵，确定群体模糊权重向量，最后单准则权重分析计算出各个划分指标的权重系数；
7、利用加权模糊c均值聚类方法，将步骤5获得的各个划分指标无量纲化数据与步骤6确定的各个划分指标相对生态地质环境整体重要性的权重系数相结合在matlab中进行聚类计算，输出不同的聚类计算结果，并以文本文件(.txt)的形式存储；
8、将步骤7中计算得出的以文本文件(.txt)存储的聚类结果在arcgis中软件中打开，结合在步骤7中计算得出的各因素聚类中心值，并根据各划分指标的生态、水文及地质特征进行分析判别，确定不同生态地质环境类型，得到生态地质环境类型分区图。
本实施例步骤1具体为：通过遥感图像提取植被指数(ndvi)，所选用影像为landsat8卫星遥感数据，根据研究区范围，选用两幅数据经图像镶嵌而成，卫星过境采集数据时，研究区天气晴朗，天空没有覆盖大面积的云层，因而两幅图全图云量较低，成像质量高，图像清晰，分辨率均为30米。基于30m的数字高程模型数据，利用arcgis10.5空间分析功能，提取研究区的高程和坡度，通过实地踏勘，及多年地质资料的积累，整理所需的生态、水文、地质资料。
本实施例步骤2中，生态地质环境类型划分作为目标层，植被归一化指数(f1)、地表高程(f2)、地形坡度(f3)、地表岩性(f4)、地貌类型(f5)、水系河网(f6)作为划分指标，形成待划分区域的生态地质环境的层次结构模型，如图2所示。
接着步骤2提取所述的6个划分指标对应的生态、水文、地质数据，继续执行步骤3。
步骤3中，将待划分区域的生态、水文、地质数据导入到arcgis中，建立各指标单因素图层，如图3-图8。
步骤4中，在arcgis10.5中将评价因子中shp格式的数据转换成grid格式的栅格数据，进而转换成matlab可识别的.flt浮点型数据，转化结果包含两个文件，一个为hdr扩展名的头文件，包含了栅格左下角的x，y坐标、栅格大小、栅格的行数和列数等信息，另一个为flt扩展名的浮点数据。
对待划分生态地质环境类型区域的各单项指标数据图层进行栅格化处理，将待评价区域划分为n个基础评价单元，n＝682*903＝615846个基础单元。
步骤5中，在matlab中，使用read_agaschdr函数读取待划分区域各指标，用normalization函数对因子做归一化去量纲处理，各划分指标归一化处理后如图9-图14。
normalization函数：
式中，fi为每个划分指标中第i个无量纲化处理后的数据，a和b分别为归一化范围的下限和上限，xi是每个划分指标中第i个无量纲化前的原始数据，max(xi)和min(xi)为各划分指标原始数据的最大值和最小值。
步骤6包括以下步骤：
(601)t.l.saaty1-9标度法对各划分指标进行相对生态地质环境整体重要性评分：
(602)建立两两比较判断矩阵
(603)构建群体的模糊判断矩阵
(604)确定群体模糊权重向量
w1＝[0.0630.0980.158]w2＝[0.0570.0920.177]
w3＝[0.0900.1430.228]w4＝[0.1770.2850.437]
w5＝[0.1800.2780.419]w6＝[0.0620.1040.173]
(605)各划分指标权重系数
步骤7中，对聚类函数custom_fcm进行改进，在计算欧氏距离过程中加入属性权重wi，设置聚类参数，对上述归一化因子进行聚类分析。在matlab处理之后，用fprintf函数对结果进行后处理，先将文件读入时获取的栅格左下角的x，y坐标和栅格行列数等信息参数重新写入头文件，再输出计算的栅格数值，进而将计算结果转为ascii数据，利用arcgis软件读取ascii文件，将其转换成栅格文件，输出生态地质环境类型区划图，如图15。
本发明涉及了一种基于煤炭资源开发的生态地质环境类型划分方法，该方法是将我国西北煤炭资源丰富、生态地质环境却十分脆弱的干旱半干旱地区划分不同的生态地质环境类型，并绘制出生态地质环境类型分区图。本发明方法首先是在对区域生态、水文、地质等相关资料调查的基础之上，将影响生态地质环境的诸多因素收集整理，并利用归一化函数将各个因素无量纲化；其次，利用模糊德尔菲层次分析法确定每个因素对生态地质环境影响的权重系数；再次，以matlab为计算平台，利用加权模糊c均值聚类法对各个影响因素进行叠加聚类计算，得出三种不同的聚类结果；最后，利用arcgis将聚类结果进行图像处理，通过每个因素的聚类中心值分析判断确定不同的生态地质环境类型。本发明能够根据现有的生态水文地质资料，快速有效地划分出不同的生态地质环境类型，确定不同类型的生态地质环境的生态地质特征及其对煤炭资源开采活动的敏感性，从而为保护类似干旱半干旱地区珍惜的潜水资源，维护好脆弱的生态环境同时选择恰当的采煤方法实现煤炭资源的开发利用提供科学依据，对西北生态环境脆弱区保水采煤具有重要意义。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。