CN113919227A - 一种矿区生态时间累积效应点与空间累积范围识别方法 - Google Patents
一种矿区生态时间累积效应点与空间累积范围识别方法 Download PDFInfo
- Publication number
- CN113919227A CN113919227A CN202111201312.5A CN202111201312A CN113919227A CN 113919227 A CN113919227 A CN 113919227A CN 202111201312 A CN202111201312 A CN 202111201312A CN 113919227 A CN113919227 A CN 113919227A
- Authority
- CN
- China
- Prior art keywords
- data
- area
- mining
- space
- index
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000005065 mining Methods 0.000 title claims abstract description 188
- 230000000694 effects Effects 0.000 title claims abstract description 106
- 238000009825 accumulation Methods 0.000 title claims abstract description 97
- 238000000034 method Methods 0.000 title claims abstract description 67
- 238000013528 artificial neural network Methods 0.000 claims abstract description 19
- 238000011160 research Methods 0.000 claims description 70
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 59
- 239000002689 soil Substances 0.000 claims description 43
- 210000002569 neuron Anatomy 0.000 claims description 42
- 238000001556 precipitation Methods 0.000 claims description 35
- 229930002875 chlorophyll Natural products 0.000 claims description 34
- 235000019804 chlorophyll Nutrition 0.000 claims description 34
- ATNHDLDRLWWWCB-AENOIHSZSA-M chlorophyll a Chemical compound C1([C@@H](C(=O)OC)C(=O)C2=C3C)=C2N2C3=CC(C(CC)=C3C)=[N+]4C3=CC3=C(C=C)C(C)=C5N3[Mg-2]42[N+]2=C1[C@@H](CCC(=O)OC\C=C(/C)CCC[C@H](C)CCC[C@H](C)CCCC(C)C)[C@H](C)C2=C5 ATNHDLDRLWWWCB-AENOIHSZSA-M 0.000 claims description 34
- 230000001186 cumulative effect Effects 0.000 claims description 28
- 238000004364 calculation method Methods 0.000 claims description 27
- 239000010410 layer Substances 0.000 claims description 24
- 238000009304 pastoral farming Methods 0.000 claims description 19
- 230000004913 activation Effects 0.000 claims description 18
- 230000001419 dependent effect Effects 0.000 claims description 18
- 230000006870 function Effects 0.000 claims description 18
- 238000012549 training Methods 0.000 claims description 17
- 238000002310 reflectometry Methods 0.000 claims description 14
- 239000003245 coal Substances 0.000 claims description 12
- 239000012298 atmosphere Substances 0.000 claims description 9
- 230000005540 biological transmission Effects 0.000 claims description 9
- 238000009499 grossing Methods 0.000 claims description 9
- 238000010606 normalization Methods 0.000 claims description 9
- 238000011002 quantification Methods 0.000 claims description 9
- 229930002877 anthocyanin Natural products 0.000 claims description 8
- 235000010208 anthocyanin Nutrition 0.000 claims description 8
- 239000004410 anthocyanin Substances 0.000 claims description 8
- 150000004636 anthocyanins Chemical class 0.000 claims description 8
- 235000021466 carotenoid Nutrition 0.000 claims description 8
- 150000001747 carotenoids Chemical class 0.000 claims description 8
- 238000007637 random forest analysis Methods 0.000 claims description 8
- 238000005316 response function Methods 0.000 claims description 7
- 230000003595 spectral effect Effects 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 6
- 238000009795 derivation Methods 0.000 claims description 6
- 238000013139 quantization Methods 0.000 claims description 6
- 230000005855 radiation Effects 0.000 claims description 6
- 238000001228 spectrum Methods 0.000 claims description 6
- 230000001550 time effect Effects 0.000 claims description 6
- 244000144972 livestock Species 0.000 claims description 5
- 238000004519 manufacturing process Methods 0.000 claims description 5
- 238000010220 Pearson correlation analysis Methods 0.000 claims description 3
- 238000012888 cubic function Methods 0.000 claims description 3
- 238000005520 cutting process Methods 0.000 claims description 3
- 238000005259 measurement Methods 0.000 claims description 3
- 210000004205 output neuron Anatomy 0.000 claims description 3
- 239000002344 surface layer Substances 0.000 claims description 3
- 230000001360 synchronised effect Effects 0.000 claims description 3
- 230000008859 change Effects 0.000 description 7
- 238000004458 analytical method Methods 0.000 description 5
- 238000012952 Resampling Methods 0.000 description 3
- 238000000354 decomposition reaction Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 230000007613 environmental effect Effects 0.000 description 3
- 238000009827 uniform distribution Methods 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007418 data mining Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 238000005067 remediation Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 231100000027 toxicology Toxicity 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
- G06Q10/063—Operations research, analysis or management
- G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
- G06Q10/06393—Score-carding, benchmarking or key performance indicator [KPI] analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q10/00—Administration; Management
- G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
- G06Q10/063—Operations research, analysis or management
- G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations
- G06Q10/06395—Quality analysis or management
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/10—Services
- G06Q50/26—Government or public services
Landscapes
- Engineering & Computer Science (AREA)
- Business, Economics & Management (AREA)
- Human Resources & Organizations (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Economics (AREA)
- Strategic Management (AREA)
- Development Economics (AREA)
- Educational Administration (AREA)
- Tourism & Hospitality (AREA)
- Entrepreneurship & Innovation (AREA)
- General Business, Economics & Management (AREA)
- Marketing (AREA)
- Mathematical Physics (AREA)
- Operations Research (AREA)
- Mathematical Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Health & Medical Sciences (AREA)
- Health & Medical Sciences (AREA)
- General Engineering & Computer Science (AREA)
- Primary Health Care (AREA)
- Software Systems (AREA)
- Evolutionary Computation (AREA)
- Quality & Reliability (AREA)
- Game Theory and Decision Science (AREA)
- Mining & Mineral Resources (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Animal Husbandry (AREA)
- Algebra (AREA)
- Agronomy & Crop Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Marine Sciences & Fisheries (AREA)
- Databases & Information Systems (AREA)
- Medical Informatics (AREA)
- Computer Hardware Design (AREA)
Abstract
本发明公开了一种矿区生态时间累积效应点与空间累积范围识别方法，A、构建矿区生态质量指数，B、构建矿区生态扰动时空累积效应指数，C、构建驱动因子数据集与地理时空加权人工神经网络模型，进而构建矿区MESCEI时间序列数据集与MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合时间累积曲线并识别曲线拐点，通过数据拟合坐标系模型在直角坐标系中拟合空间影响曲线并识别空间影响范围。本发明建立遥感反演模型并反演得到矿区生态质量指数，构建地理时空加权人工神经网络模型，然后由此构建出MESCEI时间序列数据集与MESCEI指数空间序列数据集，实现了矿区生态的时间累积效应点与空间累积范围的识别，为矿区生态管理提供数据支持。
Description
技术领域
本发明涉及采矿领域、生态学领域、遥感及地理信息领域，尤其涉及一种矿区生态时间累积效应点与空间累积范围识别方法。
背景技术
煤炭产业在推动经济发展的同时也对生态环境造成具有生态累积效应的破坏，对矿区生态累积效应的研究是矿区生态管理的重要组成部分。生态累积效应是指当一项行动与过去、现在以及可合理预见的将来行动结合在一起时，所产生的对环境的累加影响，具有时间累积性(参见图1)和空间累积性(参见图2)(如：周其刚.环境影响评价中累积效应分析方法的探讨[J].能源环境保护,2013,27(01):60-62.)。时间累积性是指当两次扰动之间的时间间隔小于环境修复所需时间时，在时间尺度上产生的累积现象。空间累积效应指相邻扰动因素之间的空间接近度小于去除每个扰动所需的距离时，在空间尺度上产生的累积现象(如：孙金玉,卜庆伟,吴东奎,等.煤矿区生态累积效应评价研究进展[J].生态毒理学报,2019,014(005):74-82.)。由概念可知，生态累积效应存在不可自然恢复与可自然恢复的临界时间(即时间累积效应点)、空间上叠加的空间累积范围。然而，现有技术方法存在以下3个问题：(1)均使用生态参数的直接观测量开展生态累积效应分析，没有从众多因子如城镇放牧、气象等的综合影响中提取出采矿的单独影响；(2)缺少从多维度(时间维、空间维)对矿区生态累积效应进行量化建模方法；(3)对矿区生态累积效应的时间累积效应点与空间累积范围的识别问题尚未有解决方案。综上所述，目前的技术方法不能有效获取矿区生态累积效应的特征。
发明内容
针对现有技术存在的缺陷，本发明的目的在于提供一种矿区生态时间累积效应点与空间累积范围识别方法，建立遥感反演模型进行反演构建出矿区生态质量指数，构建驱动因子数据集并以目标矿区的MEQI指数、自然条件数据、人类活动数据构建地理时空加权人工神经网络模型，然后由此构建出矿区MESCEI时间序列数据集，通过数学拟合方法得到其变化曲线并进一步计算其曲线拐点，该曲线拐点即为生态累积效应的时间累积效应点，也即临界时间；建立MESCEI指数空间序列数据集，通过曲线分解的方法得到单个矿的影响范围进而识别出矿区生态累积效应的空间累积范围；实现了矿区生态累积效应的时间累积效应点与空间累积范围的识别，为矿区生态管理提供数据支持。
本发明的目的通过下述技术方案实现：
一种矿区生态时间累积效应点与空间累积范围识别方法，其方法如下：
A、构建矿区生态质量指数，方法如下：
A1、建立遥感反演模型，反演模型包括PROSAIL植被辐射传输模型、随机森林算法模型、像元二分模型、kriging模型和归一化植被指数模型，收集目标矿区包括多光谱遥感影像、土壤参数产品影像、地面实测数据在内的原始数据，通过反演模型对原始数据的生态参数进行一体化同步反演并得到参数反演数据，生态参数包括植被、土壤、大气、水，植被包括植被覆盖度、叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量，土壤包括土壤含水量，大气包括PM2.5，水包括叶绿素浓度、悬浮物浓度；
A2、将所有生态参数所对应的参数反演数据逐一进行归一化处理，归一化公式如下：
其中，Xnorm为归一化后的数据，X为原始生态反演数据，Xmin为原始生态反演数据中的最小值，Xmax为原始生态反演数据中的最大值；
A3、构建MEQI指数，将各个生态参数归一化后的数据相加得到MEQI指数，公式如下：
B、构建矿区生态扰动时空累积效应指数，方法如下：
B1、选取与目标矿区相似区域作为研究对照区，研究对照区无采矿活动且远离采矿活动区域，收集研究对照区包括多光谱遥感影像、土壤参数产品影像、地面实测数据在内的原始数据；分别采集目标矿区、研究对照区的自然条件数据与人类活动数据，自然条件数据包括降水数据、气温数据、DEM数据，人类活动数据包括放牧活动数据、城镇活动数据以及采矿活动数据；构建驱动因子数据集，驱动因子数据集包括自然条件驱动因子集、人类活动驱动因子集两类，自然条件驱动因子集包括降水驱动因子、气温驱动因子、DEM驱动因子，人类活动驱动因子集包括放牧驱动因子、城镇驱动因子以及采矿驱动因子；
B2、对目标矿区、研究对照区分别依次进行降水驱动因子、气温驱动因子、DEM驱动因子、放牧驱动因子、城镇驱动因子、采矿驱动因子量化处理；
B3、以研究对照区的MEQI指数、自然条件数据、人类活动数据构建地理时空加权人工神经网络模型，在三维空间上进行研究对照区的自然条件数据、人类活动数据的空间维度叠加及高维度拓展，高维度拓展包括时间维度拓展，并将以栅格影像格式展现的低维数据转换为高维数据立方体，最终构成生态演变大数据立方体；采用滑动立方体法对生态演变大数据立方体进行数据提取，将生态演变大数据立方体中的自然条件数据、人类活动数据、MEQI指数按照研究需求划分为自变量参数和因变量参数，自变量参数包括自然条件数据、人类活动数据，因变量参数包括MEQI指数，通过地理时空加权人工神经网络模型进行模型训练并构筑出自变量与因变量间的非线性复杂定量关系；
B4、通过训练后的地理时空加权人工神经网络模型根据目标矿区的自然条件数据、人类活动数据得到目标矿区的MEQI预测值，MEQI预测值即为目标矿区无采矿条件下的MEQI预测值；
B5、通过如下公式得到矿区MESCEI指数：MESCEI＝MEQIreal-MEQIpre；其中MEQIreal为步骤A3基于遥感反演值归一化后所得到目标矿区的MEQI指数，MEQIpre为目标矿区无采矿条件下的MEQI预测值；
C、根据步骤B5中建立的MESCEI指数、目标矿区历史的原始数据构建矿区MESCEI时间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到时间累积曲线并识别出曲线拐点，曲线拐点即为时间效应点，直角坐标系的横坐标为开采年际、纵坐标为MESCEI指数。
本发明矿区生态时间累积效应点与空间累积范围识别方法第一种优选的技术方案是：在步骤C之后还包括步骤D；
D、根据目标矿区的MEQI指数、目标矿区历史的原始数据构建MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围，直角坐标系的横坐标为到目标矿区的距离值、纵坐标为MESCEI指数。
本发明矿区生态时间累积效应点与空间累积范围识别方法第二种优选的技术方案是：在步骤C之后还包括步骤D；
D、按照步骤A、步骤B的方法得出N个目标矿区的MEQI指数，N≥2，并根据N个目标矿区的MEQI指数、N个目标矿区历史的原始数据构建MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围并识别出空间影响范围，直角坐标系的横坐标为相邻目标矿区的距离值、纵坐标为MESCEI指数均值。
本发明矿区生态时间累积效应点与空间累积范围识别方法的步骤A1优选还包括如下方法：
A11、植被参数反演：归一化植被指数计算公式如下：
植被覆盖度采用像元二分模型计算，其计算公式如下：
面向矿区小尺度的Landsat和Sentinel多源数据，包含叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量的植被参数采用PROSAIL植被辐射传输模型耦合Landsat及Sentinel系列卫星传感器光谱响应函数，结合地面实测光谱和参数数据并基于随机森林算法建立植被参数反演模型，其中冠层叶绿素含量可由叶片叶绿素含量及叶面积指数相乘计算获得；
A12、土壤参数反演：通过土壤含水量影像产品在多光谱遥感卫星数据的辅助下进行降尺度研究，获取矿区长时间序列的中小空间尺度土壤含水量产品；面向矿区场景，通过Sentinel水云模型对矿区土壤含水量反演进行优化；通过降尺度和Sentinel水云模型优化组合实现表层含水量的长时序、高分辨率反演；
A13、大气参数反演：获取地表PM2.5数据，通过kriging模型对其进行克里金插值得到连续的栅格影像数据；
A14、水参数反演：对水中的叶绿素浓度和悬浮物浓度进行反演，其公式如下：
水叶绿素浓度
其中，ρRed为红波段地表反射率，ρNIR为近红外波段地表反射率；a、b、c分别为模型的系数；
水悬浮物浓度
其中，ρRed为红波段地表反射率；ρGre为绿波段地表反射率；a、b分别为模型的系数。
本发明矿区生态时间累积效应点与空间累积范围识别方法的步骤B2优选还包括如下方法：
B21、采集研究区内自然条件数据，研究区为目标矿区或研究对照区，自然条件数据包括降水数据、气温数据、DEM数据，自然条件驱动因子集的降水驱动因子对应降水数据，气温驱动因子对应气温数据，将降水数据、气温数据与步骤A1中的植被参数按照如下公式进行皮尔逊相关性分析并得到降水驱动因子与气温驱动因子分别量化对应的皮尔逊相关系数：
采集研究区内地形地貌数据，地形地貌数据包括数字高程模型数据，从数字高程模型数据中裁剪出研究区的DEM数据并对应DEM驱动因子；
B22、获取边界数据：以Landsat影像综合提取识别城镇边界与采矿边界；
B221、按照如下公式得到放牧驱动因子量化所对应的放牧强度Xgraze：
其中，Xgraze为放牧强度，X牲畜为研究区内牲畜数量，Xarea为研究区内村落总面积；
B222、采用欧氏距离获得研究区内每个像元点到城镇边界的最短距离并结合研究区域的人口量化出城镇驱动因子，其中欧式距离的计算公式如下：
其中，n为像元点的个数，Xi为各像元点的位置，Yi为城镇的像元点位置；对城镇驱动因子的量化公式如下：
B223、采用欧氏距离获得研究区内每个像元点到采矿边界的最短距离并结合研究区域的煤炭开采量量化出采矿活动的影响，其中欧式距离的计算公式如下：
其中，n为像元点的个数，Xi为各像元点的位置，Yi为采矿边界的像元点位置；对采矿活动的量化公式如下：
其中，Xmine为采矿活动的量化结果，Xmining为每年的煤炭开采量，dist(X,Y)为栅格影像像元点到采矿边界的最短距离。
本发明矿区生态时间累积效应点与空间累积范围识别方法的步骤B3优选还包括如下方法：
B31、在生态演变大数据立方体的栅格影像上建立一个滑动窗口，滑动窗口的带宽长度为L,滑动窗口的步长为S，滑动窗口的步长S≤滑动窗口的带宽长度L，然后对滑动窗口进行高维拓展，增加一个时间维度构建三维的滑动立方体，滑动立方体的时间窗口宽度为T，滑动立方体将会对时间序列所有的整幅栅格影像进行逐像元遍历，根据建立滑动立方体的栅格影像所属的集合不同，滑动立方体单次覆盖的范围内提取的像元将被化分成单位自变量样本和单位因变量样本，当滑动立方体对整幅栅格影像遍历完成，单位样本将分别组合成为自变量参数集合与因变量参数集合；
B32、地理时空加权人工神经网络模型构筑自变量与因变量之间的关系网络，该关系网络包括输入层、隐藏层、输出层三层结构，自变量样本从输入层进行数据输入，然后在输入层传导进入隐藏层时，计算公式如下式所示：
神经元i的输出计算如下所示：
pi＝Φ(Layerj)其中，Φ为激活函数，在神经元内进行激活函数的运算，神经元中采用的激活函数为非线性的双曲正切函数，其公式如下式所示：
以对应自变量样本的样本值作为目标值ti与神经网络输出值pi进行误差计算，计算公式如下所示：
误差计算完成后，后向传播的误差信号计算公式如下所示：
本发明矿区生态时间累积效应点与空间累积范围识别方法的步骤C方法优选还包括：
C1、根据步骤B5中建立的MESCEI指数、目标矿区历史的原始数据构建矿区MESCEI时间序列数据集，矿区MESCEI时间序列数据集具有时间维的数据立方体；
C2、在MESCEI时间序列数据集中建立x×y×t滑动窗口，x与y为空间横坐标与纵坐标，t为时间维，采用此滑动窗口遍历整个MESCEI时间序列数据集实现去除尖锐噪声、平滑数据的目的；
C3、数据拟合坐标系模型在时间维度上进行MESCEI时间序列数据集的三次函数关系拟合并得到时间累积曲线，时间累积曲线为y＝ax3+bx2+c；
C4、对时间累积曲线进行二阶求导，其二阶导数由负变正的点为曲线拐点或时间效应点。
优选地，本发明步骤C2中的x×y×t滑动窗口为1×1×3滑动窗口。
优选地，本发明步骤D方法包括：
D1、根据目标矿区的MEQI指数、目标矿区历史的原始数据构建多路径的MESCEI指数空间序列数据集；
D2、在多路径的MESCEI指数空间序列数据集建立x×y滑动窗口，x与y为空间横坐标与纵坐标，采用此滑动窗口遍历整个多路径的MESCEI指数空间序列数据集实现去除尖锐噪声、平滑数据的目的；
D3、通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围，空间影响曲线为y＝a×ln(x)+b，其中x为距离矿区的距离，y为MESCEI指数，a，b分别为模型系数。
优选地，本发明步骤D方法包括：
D1、根据N个目标矿区的MEQI指数、N个目标矿区历史的原始数据构建多路径的MESCEI指数空间序列数据集；
D2、在多路径的MESCEI指数空间序列数据集建立x×y滑动窗口，x与y为空间横坐标与纵坐标，采用此滑动窗口遍历整个多路径的MESCEI指数空间序列数据集实现去除尖锐噪声、平滑数据的目的；
D3、通过数据拟合坐标系模型在直角坐标系中拟合得到单个目标矿区的空间影响曲线并识别出单个目标矿区的空间影响范围，空间影响曲线为y＝a×ln(x)+b，其中x为距离矿区的距离，y为MESCEI指数，a，b分别为模型系数；
D4、按照步骤D2、D3分别实现N个目标矿区的空间影响范围的识别并得到各个目标矿区的空间影响范围交集部分，空间影响范围交集部分即为N个目标矿区下的空间累积范围。
本发明较现有技术相比，具有以下优点及有益效果：
(1)本发明建立遥感反演模型进行反演构建出矿区生态质量指数，构建驱动因子数据集并以目标矿区的MEQI指数、自然条件数据、人类活动数据构建地理时空加权人工神经网络模型，然后由此构建出矿区MESCEI时间序列数据集，通过数学拟合方法得到其变化曲线并进一步计算其曲线拐点，该曲线拐点即为生态累积效应的时间累积效应点，也即临界时间；建立MESCEI指数空间序列数据集，通过曲线分解的方法得到单个矿的影响范围进而识别出矿区生态累积效应的空间累积范围；实现了矿区生态累积效应的时间累积效应点与空间累积范围的识别，为矿区生态管理提供数据支持。
(2)本发明构建矿区生态质量指数并通过选择与矿区具有相似的自然地理特征、人文地理特征的非矿区作为对照区域，对照区与研究区除了采矿因素外具有相似的其他特征，首次构建矿区生态扰动时空累积效应指数，得到矿区生态累积效应的特征。
(3)本发明构建矿区MESCEI时间序列数据集，通过数学拟合方法得到其变化曲线并进一步计算其曲线拐点，得到了矿区生态累积效应的时间累积效应点(即临界时间)。
(4)本发明构建MESCEI指数空间序列数据集，通曲线分解的方法得到单个矿的影响范围，识别出矿区生态累积效应的空间累积范围。
附图说明
图1为矿区生态累积效应的时间累积效应点概念示意图；
图2为矿区生态累积效应的空间累积范围概念示意图；
图3为本发明方法的技术流程图；
图4为本实施例MESCEI均值与开采年份的曲线拟合图；
图5为本实施例MESCEI指数与距矿距离的曲线拟合图；
图6为本实施例西二-胜利煤矿的MESCEI指数变化图；
图7为本实施例胜利-西二矿区生态累积效应空间累积范围图。
具体实施方式
下面结合实施例对本发明作进一步地详细说明：
实施例
如图1～图7所示，一种矿区生态时间累积效应点与空间累积范围识别方法，其方法如下：
A、构建矿区生态质量指数(英文全称：Mining Ecological Quality Index)，如图3所示，方法如下：
A1、建立遥感反演模型，反演模型包括PROSAIL植被辐射传输模型、随机森林算法模型、像元二分模型、kriging模型和归一化植被指数模型，收集目标矿区包括多光谱遥感影像、土壤参数产品影像、地面实测数据在内的原始数据，通过反演模型对原始数据的生态参数进行一体化同步反演并得到参数反演数据，生态参数包括植被、土壤、大气、水，植被包括植被覆盖度、叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量，土壤包括土壤含水量，大气包括PM2.5，水包括叶绿素浓度、悬浮物浓度；
根据本实施例的一个优选实施例，步骤A1包括如下方法：
A11、植被参数反演：归一化植被指数计算公式如下：
植被覆盖度采用像元二分模型计算，其计算公式如下：
面向矿区小尺度的Landsat和Sentinel多源数据，包含叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量的植被参数采用PROSAIL植被辐射传输模型耦合Landsat及Sentinel系列卫星传感器光谱响应函数，结合地面实测光谱和参数数据并基于随机森林算法建立植被参数反演模型，其中冠层叶绿素含量可由叶片叶绿素含量及叶面积指数相乘计算获得；
在欧洲气象卫星组织(EUMETSAT)下载Landsat-5/7/8卫星的光谱响应函数，将其转换为txt格式的文件，利用MATLAB读取光谱响应函数；然后随机生成研究区范围内可能出现的1000组植被参数，通过PROSAIL模型生成其对应的模拟光谱，结合Landsat-5/7/8的光谱响应函数，分别将得到的模拟光谱通过积分的方法得到其对应Landsat-5/7/8卫星传感器的波段上，从而形成针对不同Landsat传感器的基于PROSAIL模型的植被参数反演训练数据集。通过此训练数据集训练随机森林模型用于进行叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量、叶片等效水厚度的反演，随后在Google Earth Engine上确定合成影像的月份信息和卫星，通过训练好的随机森林模型进行参数的反演，模型的输入为绿光、红光、近红外和短波红外波段的冠层反射率，模型的输出为植被的叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量、叶片等效水厚度。其中，冠层叶绿素含量可由叶片叶绿素含量及叶面积指数相乘计算获得。PROSAIL模型的输入参数包括叶片尺度参数、冠层尺度参数、背景土壤参数以及观测几何等信息。为了尽可能广的涵盖研究区不同植被类型(灌木、草地)，研究在模拟时对大多数参数给出了一定的数值分布范围。其中，叶片叶绿素含量、干物质含量与等效水厚度的数值范围满足截短的高斯正态分布，叶面积指数和平均叶倾角给定的数值范围满足均匀分布，土壤背景光谱采用PROSEPCT-D模型内置的干湿土壤线性混合的方式获得，干土壤所占比例服从0.7-1的均匀分布。太阳天顶角和观测天顶角数值分布同样设置为均匀分布。
A12、土壤参数反演：通过土壤含水量影像产品在多光谱遥感卫星数据的辅助下进行降尺度研究，获取矿区长时间序列的中小空间尺度土壤含水量产品；面向矿区场景，通过Sentinel水云模型对矿区土壤含水量反演进行优化；通过降尺度和Sentinel水云模型优化组合实现表层含水量的长时序、高分辨率反演。
本实施例步骤A12还可以如下：第一步，将分辨率较高的土壤含水量产品数据作为训练数据，将分辨率较低的土壤含水量产品数据作为辅助数据，将辅助数据重采样至训练数据分辨率，并于训练数据的差作为标签数据。在首次降尺度过程中辅助数据为原始产品数据，训练数据通过python平台中的最近邻插值函数实现。第二步，将Landsat遥感影像各波段的空值用均值替换，并进行标准化，接着将其分别重采样至与辅助数据、训练数据相同分辨率。再将与辅助数据相同分辨率的Landsat遥感影像重采样至训练数据分辨率。这两个训练数据分辨率的Landsat遥感影像作为特征数据。第三步是RF模型训练，在python中，输入标签数据与训练数据，训练RF模型，然后将Landsat遥感影像重复第二步中的过程重采样至目标预测数据分辨率，作为第三步中RF模型的输入，得到预测数据与训练数据的差值。将差值数据与训练数据做和运算得到预测的土壤含水量数据；重复上述步骤直降至目标分辨率的土壤含水量数据。
A13、大气参数反演：获取地表PM2.5数据，通过kriging模型对其进行克里金插值得到连续的栅格影像数据；
A14、水参数反演：对水中的叶绿素浓度和悬浮物浓度进行反演，其公式如下：
水叶绿素浓度
其中ρNIR为近红外波段地表反射率，在Landsat-5/7中为波段4，在Landsat-8中为波段5；ρRed为红波段地表反射率，在Landsat-5/7中为波段3，在Landsat-8中为波段4；a、b、c分别为模型的系数；
水悬浮物浓度
其中，ρRed为红波段地表反射率；ρGre为绿波段地表反射率；a、b分别为模型的系数。
A2、将所有生态参数所对应的参数反演数据逐一进行归一化处理，归一化公式如下：
其中，Xnorm为归一化后的数据，X为原始生态反演数据，Xmin为原始生态反演数据中的最小值，Xmax为原始生态反演数据中的最大值；
A3、构建MEQI指数，将各个生态参数归一化后的数据相加得到MEQI指数，公式如下：
B、构建矿区生态扰动时空累积效应指数(英文全称：Mining EcologicalSpatiotemporal Cumulative Effect Index，简称：MESCEI)，方法如下：
B1、选取与目标矿区相似区域作为研究对照区，研究对照区无采矿活动且远离采矿活动区域，收集研究对照区包括多光谱遥感影像、土壤参数产品影像、地面实测数据在内的原始数据；分别采集目标矿区、研究对照区的自然条件数据与人类活动数据，自然条件数据包括降水数据、气温数据、DEM数据，人类活动数据包括放牧活动数据、城镇活动数据以及采矿活动数据；构建驱动因子数据集，驱动因子数据集包括自然条件驱动因子集、人类活动驱动因子集两类，自然条件驱动因子集包括降水驱动因子、气温驱动因子、DEM驱动因子，人类活动驱动因子集包括放牧驱动因子、城镇驱动因子以及采矿驱动因子；
B2、对目标矿区、研究对照区分别依次进行降水驱动因子、气温驱动因子、DEM驱动因子、放牧驱动因子、城镇驱动因子、采矿驱动因子量化处理；
根据本实施例的一个优选实施例，步骤B2包括如下方法：
B21、采集研究区内自然条件数据，研究区为目标矿区或研究对照区，自然条件数据包括降水数据、气温数据、DEM数据，自然条件驱动因子集的降水驱动因子对应降水数据，气温驱动因子对应气温数据，将降水数据、气温数据与步骤A1中的植被参数按照如下公式进行皮尔逊相关性分析并得到降水驱动因子与气温驱动因子分别量化对应的皮尔逊相关系数：
采集研究区内地形地貌数据，地形地貌数据包括数字高程模型数据，从数字高程模型数据中裁剪出研究区的DEM数据并对应DEM驱动因子。
一般来说，在降水方面，植被参数与6-8月份的累计降水量相关系数最高。在气温方面，植被参数与7-9月均温的相关性最高。因此，简单粗略运算时可以选择6-8月份的累计降水量与7-9月均温作为驱动因子。
B22、获取边界数据：以Landsat影像综合提取识别城镇边界与采矿边界；
B221、按照如下公式得到放牧驱动因子量化所对应的放牧强度Xgraze：
其中，Xgraze为放牧强度，X牲畜为研究区内牲畜数量，Xarea为研究区内村落总面积；
B222、采用欧氏距离获得研究区内每个像元点到城镇边界的最短距离并结合研究区域的人口量化出城镇驱动因子，其中欧式距离的计算公式如下：
其中，n为像元点的个数，Xi为各像元点的位置，Yi为城镇的像元点位置；对城镇驱动因子的量化公式如下：
B223、采用欧氏距离获得研究区内每个像元点到采矿边界的最短距离并结合研究区域的煤炭开采量量化出采矿活动的影响，其中欧式距离的计算公式如下：
其中，n为像元点的个数，Xi为各像元点的位置，Yi为采矿边界的像元点位置；对采矿活动的量化公式如下：
其中，Xmine为采矿活动的量化结果，Xmining为每年的煤炭开采量，dist(X,Y)为栅格影像像元点到采矿边界的最短距离。
B3、以研究对照区的MEQI指数、自然条件数据、人类活动数据构建地理时空加权人工神经网络模型，在三维空间上进行研究对照区的自然条件数据、人类活动数据的空间维度叠加及高维度拓展，高维度拓展包括时间维度拓展，并将以栅格影像格式展现的低维数据转换为高维数据立方体，最终构成生态演变大数据立方体；采用滑动立方体法对生态演变大数据立方体进行数据提取，将生态演变大数据立方体中的自然条件数据、人类活动数据、MEQI指数按照研究需求划分为自变量参数和因变量参数，自变量参数包括自然条件数据、人类活动数据，因变量参数包括MEQI指数，通过地理时空加权人工神经网络模型进行模型训练并构筑出自变量与因变量间的非线性复杂定量关系；
根据本实施例的一个优选实施例，步骤B3包括如下方法：
B31、在生态演变大数据立方体的栅格影像上建立一个滑动窗口，滑动窗口的带宽长度为L(在本实施例中，带宽长度L设为3),滑动窗口的步长为S(在本实施例中，步长S设为1)，滑动窗口的步长S≤滑动窗口的带宽长度L，然后对滑动窗口进行高维拓展，增加一个时间维度构建三维的滑动立方体，滑动立方体的时间窗口宽度为T，滑动立方体将会对时间序列所有的整幅栅格影像进行逐像元遍历，根据建立滑动立方体的栅格影像所属的集合不同，滑动立方体单次覆盖的范围内提取的像元将被化分成单位自变量样本和单位因变量样本，当滑动立方体对整幅栅格影像遍历完成，单位样本将分别组合成为自变量参数集合与因变量参数集合；
B32、地理时空加权人工神经网络模型构筑自变量与因变量之间的关系网络，该关系网络包括输入层、隐藏层、输出层三层结构，自变量样本从输入层进行数据输入，然后在输入层传导进入隐藏层时，计算公式如下式所示：
神经元i的输出计算如下所示：
pi＝Φ(Layerj)；其中，Φ为激活函数，在神经元内进行激活函数的运算，神经元中采用的激活函数为非线性的双曲正切函数，其公式如下式所示：
以对应自变量样本的样本值作为目标值ti与神经网络输出值pi进行误差计算，计算公式如下所示：
误差计算完成后，后向传播的误差信号计算公式如下所示：
B4、通过训练后的地理时空加权人工神经网络模型根据目标矿区的自然条件数据、人类活动数据得到目标矿区的MEQI预测值，MEQI预测值即为目标矿区无采矿条件下的MEQI预测值；
B5、通过如下公式得到矿区MESCEI指数：MESCEI＝MEQIreal-MEQIpre；其中MEQIreal为步骤A3基于遥感反演值归一化后所得到目标矿区的MEQI指数，MEQIpre为目标矿区无采矿条件下的MEQI预测值；
C、根据步骤B5中建立的MESCEI指数、目标矿区历史的原始数据构建矿区MESCEI时间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到时间累积曲线并识别出曲线拐点，曲线拐点即为时间效应点，直角坐标系的横坐标为开采年际、纵坐标为MESCEI指数。
根据本实施例的一个优选实施例，步骤C方法包括：
C1、根据步骤B5中建立的MESCEI指数、目标矿区历史的原始数据构建矿区MESCEI时间序列数据集，参见图4，矿区MESCEI时间序列数据集具有时间维的数据立方体；
C2、在MESCEI时间序列数据集中建立x×y×t滑动窗口(优选地，本实施例x×y×t滑动窗口为1×1×3滑动窗口)，x与y为空间横坐标与纵坐标，t为时间维，采用此滑动窗口遍历整个MESCEI时间序列数据集实现去除尖锐噪声、平滑数据的目的；本实施例使用的均值模板f(m，n，p)为：
其中，m为空间横坐标方向的系数，n为空间纵坐标方向的系数，p为时间维上的系数。
C3、数据拟合坐标系模型在时间维度上进行MESCEI时间序列数据集的三次函数关系拟合并得到时间累积曲线，时间累积曲线为y＝ax3+bx2+c；本实施例时间累积曲线为：y＝2.5×10-4×x3-6.9×10-3×x2+0.05×x+0.05。
C4、对时间累积曲线进行二阶求导，其二阶导数由负变正的点为曲线拐点或时间效应点。本实施例经计算得到的曲线拐点为8.27，即矿区生态累积效应的临界时间约在2011年。
本实施例第一种优选的技术方案是：步骤C之后还包括步骤D；
D、根据目标矿区的MEQI指数、目标矿区历史的原始数据构建MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围，直角坐标系的横坐标为相距目标矿区的距离值、纵坐标为MESCEI指数。根据本实施例的一个优选实施例，本实施例第一种技术方案优选的实施例：步骤D方法包括：
D1、根据目标矿区的MEQI指数、目标矿区历史的原始数据构建多路径的MESCEI指数空间序列数据集；
D2、在多路径的MESCEI指数空间序列数据集建立x×y滑动窗口(优选地，本实施例采用3×3滑动窗口)，x与y为空间横坐标与纵坐标，采用此滑动窗口遍历整个多路径的MESCEI指数空间序列数据集实现去除尖锐噪声、平滑数据的目的；本实施例可以采用的均值模板
为：
D3、通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围，空间影响曲线为y＝a×ln(x)+b，其中x为距离矿区的距离，y为MESCEI指数，a，b分别为模型系数。
本实施例第二种优选的技术方案是：步骤C之后还包括步骤D；
D、按照步骤A、步骤B的方法得出N个目标矿区的MEQI指数，N≥2，并根据N个目标矿区的MEQI指数、N个目标矿区历史的原始数据构建MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围并识别出空间影响范围，直角坐标系的横坐标为相邻目标矿区的距离值、纵坐标为MESCEI指数均值。在本实施例中，目标矿区为包含锡林浩特市胜利一号矿区、西二矿区的区域，选择远离采矿活动区域且同时与目标矿区具有相似基础地理环境的区域，以此作为研究对照区；并在步骤A中收集分析区域内的2003-2020年逐年的8月份Landsat-5/7/8卫星影像，GLDAS、AMSR-E、AMSR2地表土壤含水量产品数据，月累积降水数据、月平均气温数据、DEM数据，区域牲畜数量数据，研究期间的城镇人口数据，矿区煤炭产量数据。遥感影像数据、地表土壤含水量产品数据以及DEM数据可通过NASA地球科学数据平台Earthdata获取，研究区的月累积降水数据、月平均气温数据数据可通过国家气象科学数据中心(https://data.cma.cn/)获取；区域牲畜数量数据可来自与研究区的统计局，城镇人口数据可来自与区域统计年鉴，矿区煤炭产量数据可来自于当地的煤炭公司。根据本实施例的一个优选实施例，本实施例第二种技术方案优选的实施例：步骤D方法包括：
D1、根据2个目标矿区的MEQI指数、2个目标矿区历史的原始数据构建多路径的MESCEI指数空间序列数据集；
D2、在多路径的MESCEI指数空间序列数据集建立x×y滑动窗口(优选地，本实施例采用3×3滑动窗口)，x与y为空间横坐标与纵坐标，采用此滑动窗口遍历整个多路径的MESCEI指数空间序列数据集实现去除尖锐噪声、平滑数据的目的；本实施例可以采用的均值模板
为：
D3、通过数据拟合坐标系模型在直角坐标系中拟合得到单个目标矿区的空间影响曲线并识别出单个目标矿区的空间影响范围，空间影响曲线为y＝a×ln(x)+b，其中x为距离矿区的距离，y为MESCEI指数，a，b分别为模型系数；如图5、图6、图7所示，在本实施例中，选择两个样区其进行数学关系拟合分析，分别得到MESCEI指数随着据矿区距离的递增呈现对数函数关系，即：y＝-0.1038×ln(x)+1；y＝-0.17×ln(x)+1.5；对受多个矿采矿活动影响的MESCEI指数与距各矿距离之间的关系数学拟合为：y＝-0.08×ln(x1)-0.07×ln(x2)+1.17。模型可分解为：y＝-0.08×ln(x1)+0.63；y＝-0.07×ln(x2)+0.54；其中：x1、x2为到胜利矿、西二矿的距离，y为MESCEI指数。
根据统计学中小概率事件实际不可能性原理，当其距各矿的距离与MESCEI指数之间的对数关系的概率分布函数达到99％时，此点对应的距离则为各矿生态累积效应的影响范围。即：当对数函数y＝a×ln(x)+b存在x满足F(x)＝P(X≤x)＝99％条件时，x代表的距离即为矿区生态累积效应的影响范围。经计算x1＝7138.6、x2＝6075.1，即胜利矿区生态累积效应的影响范围的7138.6米，西二矿区生态累积效应的影响范围的6075.1米。
D4、按照步骤D2、D3分别实现2个目标矿区的空间影响范围的识别并得到各个目标矿区的空间影响范围交集部分，空间影响范围交集部分即为2个目标矿区下的空间累积范围(即对各矿的生态累积效应影响范围取并集，此并集范围为矿区生态累积效应的空间累积范围，参见图7)。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：其方法如下：
A、构建矿区生态质量指数，方法如下：
A1、建立遥感反演模型，反演模型包括PROSAIL植被辐射传输模型、随机森林算法模型、像元二分模型、kriging模型和归一化植被指数模型，收集目标矿区包括多光谱遥感影像、土壤参数产品影像、地面实测数据在内的原始数据，通过反演模型对原始数据的生态参数进行一体化同步反演并得到参数反演数据，生态参数包括植被、土壤、大气、水，植被包括植被覆盖度、叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量，土壤包括土壤含水量，大气包括PM2.5，水包括叶绿素浓度、悬浮物浓度；
A2、将所有生态参数所对应的参数反演数据逐一进行归一化处理，归一化公式如下：
其中，Xnorm为归一化后的数据，X为原始生态反演数据，Xmin为原始生态反演数据中的最小值，Xmax为原始生态反演数据中的最大值；
A3、构建MEQI指数，将各个生态参数归一化后的数据相加得到MEQI指数，公式如下：
B、构建矿区生态扰动时空累积效应指数，方法如下：
B1、选取与目标矿区相似区域作为研究对照区，研究对照区无采矿活动且远离采矿活动区域，收集研究对照区包括多光谱遥感影像、土壤参数产品影像、地面实测数据在内的原始数据；分别采集目标矿区、研究对照区的自然条件数据与人类活动数据，自然条件数据包括降水数据、气温数据、DEM数据，人类活动数据包括放牧活动数据、城镇活动数据以及采矿活动数据；构建驱动因子数据集，驱动因子数据集包括自然条件驱动因子集、人类活动驱动因子集两类，自然条件驱动因子集包括降水驱动因子、气温驱动因子、DEM驱动因子，人类活动驱动因子集包括放牧驱动因子、城镇驱动因子以及采矿驱动因子；
B2、对目标矿区、研究对照区分别依次进行降水驱动因子、气温驱动因子、DEM驱动因子、放牧驱动因子、城镇驱动因子、采矿驱动因子量化处理；
B3、以研究对照区的MEQI指数、自然条件数据、人类活动数据构建地理时空加权人工神经网络模型，在三维空间上进行研究对照区自然条件数据、人类活动数据的空间维度叠加及高维度拓展，高维度拓展包括时间维度拓展，并将以栅格影像格式展现的低维数据转换为高维数据立方体，最终构成生态演变大数据立方体；采用滑动立方体法对生态演变大数据立方体进行数据提取，将生态演变大数据立方体中的自然条件数据、人类活动数据、MEQI指数按照研究需求划分为自变量参数和因变量参数，自变量参数包括自然条件数据、人类活动数据，因变量参数包括MEQI指数，通过地理时空加权人工神经网络模型进行模型训练并构筑出自变量与因变量间的非线性复杂定量关系；
B4、通过训练后的地理时空加权人工神经网络模型根据目标矿区的自然条件数据、人类活动数据得到目标矿区的MEQI预测值，MEQI预测值即为目标矿区无采矿条件下的MEQI预测值；
B5、通过如下公式得到矿区MESCEI指数：MESCEI＝MEQIreal-MEQIpre；其中MEQIreal为步骤A3基于遥感反演值归一化后所得到目标矿区的MEQI指数，MEQIpre为目标矿区无采矿条件下的MEQI预测值；
C、根据步骤B5中建立的MESCEI指数、目标矿区历史的原始数据构建矿区MESCEI时间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到时间累积曲线并识别出曲线拐点，曲线拐点即为时间效应点，直角坐标系的横坐标为开采年际、纵坐标为MESCEI指数。
2.按照权利要求1所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤C之后还包括步骤D；
D、根据目标矿区的MEQI指数、目标矿区历史的原始数据构建MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围，直角坐标系的横坐标为相距目标矿区的距离值、纵坐标为MESCEI指数。
3.按照权利要求1所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤C之后还包括步骤D；
D、按照步骤A、步骤B的方法得出N个目标矿区的MEQI指数，N≥2，并根据N个目标矿区的MEQI指数、N个目标矿区历史的原始数据构建MESCEI指数空间序列数据集，通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围并识别出空间影响范围，直角坐标系的横坐标为到目标矿区的距离值、纵坐标为MESCEI指数均值。
4.按照权利要求1所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤A1包括如下方法：
A11、植被参数反演：归一化植被指数计算公式如下：
其中ρNIR为近红外波段地表反射率，ρRed为红波段地表反射率；
植被覆盖度采用像元二分模型计算，其计算公式如下：
面向矿区小尺度的Landsat和Sentinel多源数据，包含叶面积指数、叶片叶绿素含量、叶片等效水厚度、叶片类胡萝卜素含量、叶片花青素含量、叶片干物质含量、冠层叶绿素含量的植被参数采用PROSAIL植被辐射传输模型耦合Landsat及Sentinel系列卫星传感器光谱响应函数，结合地面实测光谱和参数数据并基于随机森林算法建立植被参数反演模型，其中冠层叶绿素含量可由叶片叶绿素含量及叶面积指数相乘计算获得；
A12、土壤参数反演：通过土壤含水量影像产品在多光谱遥感卫星数据的辅助下进行降尺度研究，获取矿区长时间序列的中小空间尺度土壤含水量产品；面向矿区场景，通过Sentinel水云模型对矿区土壤含水量反演进行优化；通过降尺度和Sentinel水云模型优化组合实现表层含水量的长时序、高分辨率反演；
A13、大气参数反演：获取地表PM2.5数据，通过kriging模型对其进行克里金插值得到连续的栅格影像数据；
A14、水参数反演：对水中的叶绿素浓度和悬浮物浓度进行反演，其公式如下：
水叶绿素浓度
其中，ρRed为红波段地表反射率，ρNIR为近红外波段地表反射率；a、b、c分别为模型的系数；
水悬浮物浓度
其中，ρRed为红波段地表反射率；ρGre为绿波段地表反射率；a、b分别为模型的系数。
5.按照权利要求1或2所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤B2包括如下方法：
B21、采集研究区内自然条件数据，研究区为目标矿区或研究对照区，自然条件数据包括降水数据、气温数据、DEM数据，自然条件驱动因子集的降水驱动因子对应降水数据，气温驱动因子对应气温数据，将降水数据、气温数据与步骤A1中的植被参数按照如下公式进行皮尔逊相关性分析并得到降水驱动因子与气温驱动因子分别量化对应的皮尔逊相关系数：
其中，r为皮尔逊相关系数，n为每个变量中需要进行分析的数据量，Xi为降水数据或气温数据的值，
为降水数据或气温数据的平均值，Yi为植被参数的值，
为植被参数的平均值；
采集研究区内地形地貌数据，地形地貌数据包括数字高程模型数据，从数字高程模型数据中裁剪出研究区的DEM数据并对应DEM驱动因子；
B22、获取边界数据：以Landsat影像综合提取识别城镇边界与采矿边界；
B221、按照如下公式得到放牧驱动因子量化所对应的放牧强度Xgraze：
其中，Xgraze为放牧强度，X牲畜为研究区内牲畜数量，Xarea为研究区内村落总面积；
B222、采用欧氏距离获得研究区内每个像元点到城镇边界的最短距离并结合研究区域的人口量化出城镇驱动因子，其中欧式距离的计算公式如下：
其中，n为像元点的个数，Xi为各像元点的位置，Yi为城镇的像元点位置；对城镇驱动因子的量化公式如下：
其中，Xurban为城镇活动的量化结果，Xpop为城镇的人口数，dist(X，Y)为栅格影像像元点到城镇边界的最短距离；
B223、采用欧氏距离获得研究区内每个像元点到采矿边界的最短距离并结合研究区域的煤炭开采量量化出采矿活动的影响，其中欧式距离的计算公式如下：
其中，n为像元点的个数，Xi为各像元点的位置，Yi为采矿边界的像元点位置；对采矿活动的量化公式如下：
其中，Xmine为采矿活动的量化结果，Xmining为每年的煤炭开采量，dist(X，Y)为栅格影像像元点到采矿边界的最短距离。
6.按照权利要求1或2所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤B3包括如下方法：
B31、在生态演变大数据立方体的栅格影像上建立一个滑动窗口，滑动窗口的带宽长度为L，滑动窗口的步长为S，滑动窗口的步长S≤滑动窗口的带宽长度L，然后对滑动窗口进行高维拓展，增加一个时间维度构建三维的滑动立方体，滑动立方体的时间窗口宽度为T，滑动立方体将会对时间序列所有的整幅栅格影像进行逐像元遍历，根据建立滑动立方体的栅格影像所属的集合不同，滑动立方体单次覆盖的范围内提取的像元将被化分成单位自变量样本和单位因变量样本，当滑动立方体对整幅栅格影像遍历完成，单位样本将分别组合成为自变量参数集合与因变量参数集合；
B32、地理时空加权人工神经网络模型构筑自变量与因变量之间的关系网络，该关系网络包括输入层、隐藏层、输出层三层结构，自变量样本从输入层进行数据输入，然后在输入层传导进入隐藏层时，计算公式如下式所示：
神经元i的输出计算如下所示：
pi＝Φ(Layerj)；其中，Φ为激活函数，在神经元内进行激活函数的运算，神经元中采用的激活函数为非线性的双曲正切函数，其公式如下式所示：
以对应自变量样本的样本值作为目标值ti与神经网络输出值pi进行误差计算，计算公式如下所示：
误差计算完成后，后向传播的误差信号计算公式如下所示：
7.按照权利要求1或2所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤C方法包括：
C1、根据步骤B5中建立的MESCEI指数、目标矿区历史的原始数据构建矿区MESCEI时间序列数据集，矿区MESCEI时间序列数据集具有时间维的数据立方体；
C2、在MESCEI时间序列数据集中建立x×y×t滑动窗口，x与y为空间横坐标与纵坐标，t为时间维，采用此滑动窗口遍历整个MESCEI时间序列数据集实现去除尖锐噪声、平滑数据的目的；
C3、数据拟合坐标系模型在时间维度上进行MESCEI时间序列数据集的三次函数关系拟合并得到时间累积曲线，时间累积曲线为y＝ax3+bx2+c；
C4、对时间累积曲线进行二阶求导，其二阶导数由负变正的点为曲线拐点或时间效应点。
8.按照权利要求7所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤C2中的x×y×t滑动窗口为1×1×3滑动窗口。
9.按照权利要求2所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤D方法包括：
D1、根据目标矿区的MEQI指数、目标矿区历史的原始数据构建多路径的MESCEI指数空间序列数据集；
D2、在多路径的MESCEI指数空间序列数据集建立x×y滑动窗口，x与y为空间横坐标与纵坐标，采用此滑动窗口遍历整个多路径的MESCEI指数空间序列数据集实现去除尖锐噪声、平滑数据的目的；
D3、通过数据拟合坐标系模型在直角坐标系中拟合得到空间影响曲线并识别出空间影响范围，空间影响曲线为y＝a×ln(x)+b，其中x为距离矿区的距离，y为MESCEI指数，a，b分别为模型系数。
10.按照权利要求3所述的一种矿区生态时间累积效应点与空间累积范围识别方法，其特征在于：步骤D方法包括：
D1、根据N个目标矿区的MEQI指数、N个目标矿区历史的原始数据构建多路径的MESCEI指数空间序列数据集；
D2、在多路径的MESCEI指数空间序列数据集建立x×y滑动窗口，x与y为空间横坐标与纵坐标，采用此滑动窗口遍历整个多路径的MESCEI指数空间序列数据集实现去除尖锐噪声、平滑数据的目的；
D3、通过数据拟合坐标系模型在直角坐标系中拟合得到单个目标矿区的空间影响曲线并识别出单个目标矿区的空间影响范围，空间影响曲线为y＝a×ln(x)+b，其中x为距离矿区的距离，y为MESCEI指数，a，b分别为模型系数；
D4、按照步骤D2、D3分别实现N个目标矿区的空间影响范围的识别并得到各个目标矿区的空间影响范围交集部分，空间影响范围交集部分即为N个目标矿区下的空间累积范围。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111201312.5A CN113919227B (zh) | 2021-10-15 | 2021-10-15 | 一种矿区生态时间累积效应点与空间累积范围识别方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111201312.5A CN113919227B (zh) | 2021-10-15 | 2021-10-15 | 一种矿区生态时间累积效应点与空间累积范围识别方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113919227A true CN113919227A (zh) | 2022-01-11 |
| CN113919227B CN113919227B (zh) | 2022-04-15 |
Family
ID=79240572
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111201312.5A Expired - Fee Related CN113919227B (zh) | 2021-10-15 | 2021-10-15 | 一种矿区生态时间累积效应点与空间累积范围识别方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113919227B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115526098A (zh) * | 2022-09-14 | 2022-12-27 | 国家能源投资集团有限责任公司 | 一种矿区地表植被叶面积指数遥感计算方法及电子设备 |
| CN116504327A (zh) * | 2022-09-26 | 2023-07-28 | 中国疾病预防控制中心环境与健康相关产品安全所 | 近地面o3人群暴露的时空精细化分析评估方法及系统 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113128134A (zh) * | 2021-06-17 | 2021-07-16 | 中国矿业大学(北京) | 一种矿区生态环境演变驱动因子权重量化分析方法 |
| CN113240296A (zh) * | 2021-05-19 | 2021-08-10 | 中国矿业大学 | 一种矿区生态累积效应的评价方法 |
-
2021
- 2021-10-15 CN CN202111201312.5A patent/CN113919227B/zh not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113240296A (zh) * | 2021-05-19 | 2021-08-10 | 中国矿业大学 | 一种矿区生态累积效应的评价方法 |
| CN113128134A (zh) * | 2021-06-17 | 2021-07-16 | 中国矿业大学(北京) | 一种矿区生态环境演变驱动因子权重量化分析方法 |
Non-Patent Citations (5)
| Title |
|---|
| LI JUN ETAL: "A Review of Remote Sensing for Environmental Monitoring in China", 《REMOTE SENSING 》 * |
| 张成业等: "矿区生态环境定量遥感监测研究进展与展望", 《金属矿山》 * |
| 李军等: "时空大数据支持的土地储备智能决策体系与应用研究", 《中国土地科学》 * |
| 李军等: "资源型城市长时间序列土壤含水量化分析——以锡林浩特市为例", 《测绘通报》 * |
| 邵亚琴: "基于多源动态监测数据的草原区煤电基地生态扰动与修复评价研究", 《中国优秀博硕士学位论文全文数据库(博士)工程科技Ⅰ辑》 * |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115526098A (zh) * | 2022-09-14 | 2022-12-27 | 国家能源投资集团有限责任公司 | 一种矿区地表植被叶面积指数遥感计算方法及电子设备 |
| CN116504327A (zh) * | 2022-09-26 | 2023-07-28 | 中国疾病预防控制中心环境与健康相关产品安全所 | 近地面o3人群暴露的时空精细化分析评估方法及系统 |
| CN116504327B (zh) * | 2022-09-26 | 2024-01-30 | 中国疾病预防控制中心环境与健康相关产品安全所 | 近地面o3人群暴露的时空精细化分析评估方法及系统 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113919227B (zh) | 2022-04-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113128134B (zh) | 一种矿区生态环境演变驱动因子权重量化分析方法 | |
| Danner et al. | Efficient RTM-based training of machine learning regression algorithms to quantify biophysical & biochemical traits of agricultural crops | |
| Yuan et al. | Deep learning in environmental remote sensing: Achievements and challenges | |
| Seyam et al. | Identifying the land use land cover (LULC) changes using remote sensing and GIS approach: A case study at Bhaluka in Mymensingh, Bangladesh | |
| Duveiller et al. | Crop specific green area index retrieval from MODIS data at regional scale by controlling pixel-target adequacy | |
| Hladik et al. | Salt marsh elevation and habitat mapping using hyperspectral and LIDAR data | |
| Braswell et al. | A multivariable approach for mapping sub-pixel land cover distributions using MISR and MODIS: Application in the Brazilian Amazon region | |
| CN113988626B (zh) | 一种矿区生态环境遥感综合评价指数实现方法 | |
| CN113919227B (zh) | 一种矿区生态时间累积效应点与空间累积范围识别方法 | |
| Wang et al. | Estimation of soil salt content using machine learning techniques based on remote-sensing fractional derivatives, a case study in the Ebinur Lake Wetland National Nature Reserve, Northwest China | |
| Hame et al. | Improved mapping of tropical forests with optical and SAR imagery, Part II: Above ground biomass estimation | |
| Heo et al. | Digital elevation model-based convolutional neural network modeling for searching of high solar energy regions | |
| CN114021348A (zh) | 一种精细土地利用类型的矿区植被碳汇遥感反演方法 | |
| CN113553697B (zh) | 基于长时序多源数据的煤炭开采植被扰动分析方法 | |
| CN105809148A (zh) | 基于遥感时空谱融合的作物干旱识别及风险评估方法 | |
| CN113919226B (zh) | 基于权重的采矿植被生态累积效应扰动范围识别方法 | |
| Polinova et al. | Reconstructing pre-fire vegetation condition in the wildland urban interface (WUI) using artificial neural network | |
| CN114005040A (zh) | 一种基于di的森林扰动变化遥感监测方法及装置 | |
| Liu et al. | Hyperspectral infrared sounder cloud detection using deep neural network model | |
| Zhang et al. | Winter wheat yield prediction using integrated Landsat 8 and Sentinel-2 vegetation index time-series data and machine learning algorithms | |
| CN117197668A (zh) | 基于深度学习的作物倒伏级别的预测方法及系统 | |
| Somers et al. | Endmember library approaches to resolve spectral mixing problems in remotely sensed data: Potential, challenges, and applications | |
| Al-Shammari et al. | Impact of spatial resolution on the quality of crop yield predictions for site-specific crop management | |
| Aziz et al. | An assessment of sedimentation in Terengganu River, Malaysia using satellite imagery | |
| Xu et al. | Air temperature estimation over winter wheat fields by integrating machine learning and remote sensing techniques |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220415 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |