基于权重的采矿植被生态累积效应扰动范围识别方法与流程

1.本发明涉及采矿遥感数据识别处理领域，尤其涉及一种基于权重的采矿植被生态累积效应扰动范围识别方法。


背景技术：

2.矿产资源开采通过对岩层的开挖搬运对矿区植被造成强烈影响，具有生态累积效应，破坏了当地的自然生态系统，这种影响存在扰动的空间地理范围，即采矿植被生态累积效应扰动范围(以下简称采矿对植被的扰动范围)。植被作为生态系统的生产者，在生态系统中处于关键地位，识别采矿活动对植被的扰动范围具有重要的现实意义。当前对矿区扰动范围的识别方法是基于植被指数(或者植被参数)的直接观测，并利用傅里叶分析(例如：闫凤飞.基于fourier分析的煤炭开采对植被的影响范围研究[d].中国地质大学(北京),2020)、函数型主成分分析(例如：袁涛,倪璇,周伟.煤炭开发区植被扰动时空效应及影响范围界定：以宁东矿区为例[j].地学前缘,2021,28(04):110-117)、缓冲区分析(例如：倪璇.煤炭开发扰动区植被时空效应及影响范围识别[d].中国地质大学(北京),2020)、趋势线拟合(例如：李晶,韩颖,杨震,苗辉,殷守强.基于温度植被干旱指数的兖州煤田煤炭开采影响边界遥感提取[j].农业工程学报,2018,34(19):258-265)等方法分别对植被指数进行分析。现有技术普遍存在如下技术缺陷：第一，植被指数的变化是由气候气象、地形地貌、人类活动三类数据共同作用的结果，其中气候气象包括降水和气温，地形地貌包括数字高程模型，人类活动包括放牧、城镇、采矿活动等。而现有有关矿区扰动范围的研究(傅里叶分析、函数型主成分分析、缓冲区分析、趋势线拟合，还如：yang y,erskine p d,lechner a m,et al.detecting the dynamics of vegetation disturbance and recovery in surface mining area via landsat imagery and landtrendr algorithm[j].journal of cleaner production,2018,178(mar.20):353-362)并没有考虑诸如气温、地形、放牧等其他因素对植被的影响，只是单纯的对植被指数进行分析，将多因素耦合的结果作为采矿因素影响的结果，导致采矿扰动范围识别存在较大误差。第二，现有研究是在研究区域的某几个方向选取样区，对每个样区的植被指数进行分析，拟合趋势线找到趋于平稳的点作为采矿对植被扰动的阈值，参见图2；而实际情况中，趋势线趋于平稳的点无法确定一定是采矿对植被扰动的阈值，目前对阈值确定的方法仍缺乏一定的合理性。第三，现有研究是将各个方向上样区的阈值点位置连接成一个闭合的曲线作为矿区的扰动范围，然而相邻样区阈值点位置的连接区域属于研究的盲区，阈值点与阈值点之间具有空间不连续性；并且现有研究对于阈值点与阈值点之间的连接方式也没有合理的说明；因此阈值点与阈值点连接构成的矿区扰动范围存在一定的误差，如图3所示。第四，矿区开采范围会随着时间的变化而改变，现有方法仅考虑了研究年份一年的影响范围，没有考虑矿区采矿活动对植被影响的时间异质性。
[0003]
综上所述，现有研究只是单纯的对植被指数进行分析，没有考虑诸如气温、地形、放牧等其他因素对植被的影响，将多因素耦合的结果作为采矿因素影响的结果；目前常见
方法都不能准确获取采矿活动对植被的扰动范围。


技术实现要素：

[0004]
针对现有技术存在的不足之处，本发明的目的在于提供一种基于权重的采矿植被生态累积效应扰动范围识别方法，利用高维数据下的因子定权法量化出气候气象、地形地貌、人类活动三类驱动因子的权重，避免了多因素耦合造成的扰动范围识别误差，利用显著性检验法能够最终识别采矿对植被的扰动范围，为挖掘矿区采矿活动对生态环境的影响机制、保护矿区生态环境等提供数据支持。
[0005]
本发明的目的通过下述技术方案实现：
[0006]
一种基于权重的采矿植被生态累积效应扰动范围识别方法，其方法如下：
[0007]
a、收集包括landsat系列卫星影像产品、sentinel-2a影像产品在内的研究区原始数据，优选地，本发明研究区原始数据为栅格化影像。landsat系列卫星影像产品对应为landsat系列遥感影像，sentinel-2a影像产品对应为sentinel-2a遥感影像；
[0008]
a1、植被参数反演：首先采用如下公式计算归一化植被指数ndvi：其中，ρ
nir
为近红外波段地表反射率，在landsat-5/7中为波段4，在landsat-8中为波段5；ρ
red
为红波段地表反射率，在landsat-5/7中为波段3，在landsat-8中为波段4；接着采用像元二分模型来计算植被覆盖度，其计算公式如下：其中，fvc为像元的植被覆盖度，ndvi为像元的ndvi值，ndvi
min
为研究区内完全为裸土的像元ndvi值，ndvi
max
为研究区纯植被像元的ndvi值；植被覆盖度为反演后的植被参数；
[0009]
b、构建驱动因子数据集，驱动因子数据集包括气候气象驱动因子集、地形地貌驱动因子集、人类活动驱动因子集三类，气候气象驱动因子集包括降水驱动因子与气温驱动因子，地形地貌驱动因子集包括dem驱动因子，人类活动驱动因子集包括放牧驱动因子、城镇驱动因子以及采矿驱动因子；
[0010]
b1、采集研究区内气候气象数据，气候气象数据包括降水数据与气温数据，气候气象驱动因子集的降水驱动因子对应气候气象数据的降水数据，气温驱动因子对应气候气象数据的气温数据，将气候气象数据与步骤a1中的植被参数按照如下公式进行皮尔逊相关性分析并得到降水驱动因子与气温驱动因子分别量化对应的皮尔逊相关系数：其中，r为皮尔逊相关系数，n为每个变量中需要进行分析的数据量，xi为降水数据或气温数据的值，为降水数据或气温数据的平均值，yi为植被参数的值，为植被参数的平均值；
[0011]
b2、采集研究区内地形地貌数据，地形地貌数据包括数字高程模型数据，从数字高程模型数据中裁剪出研究区域的dem数据并对应dem驱动因子；
[0012]
b3、获取地理行政边界数据并辅以landsat影像综合提取识别行政边界与采矿边界；
[0013]
b31、按照如下公式得到放牧驱动因子量化所对应的放牧强度x
graze
：其中，x
graze
为放牧强度，x
牲畜
为研究区域内牲畜数量，x
area
为研究区域内村落总面积；
[0014]
b32、采用欧氏距离获得研究区域内每个像元点到城镇边界的最短距离并结合研究区域的人口量化出城镇驱动因子，其中欧式距离的计算公式如下：其中，n为像元点的个数，xi为各像元点的位置，yi为城镇的像元点位置；对城镇驱动因子的量化公式如下：其中，x
urban
为城镇活动的量化结果，x
pop
为城镇的人口数，dist(x,y)为栅格影像像元点到城镇边界的最短距离；
[0015]
b33、假设m0为没有采矿的年份，m1为开始采矿的年份，m2为结束采矿的年份，则对m
1-m2时期采矿驱动因子所对应的采矿活动数据进行量化；采用欧氏距离获得研究区域内每个像元点到采矿边界的最短距离并结合研究区域每年的煤炭开采量量化出采矿活动的影响，其中欧式距离的计算公式如下：其中，n为像元点的个数，xi为各像元点的位置，yi为采矿边界的像元点位置；对采矿活动的量化公式如下：其中，x
mine
为采矿活动的量化结果，x
mining
为每年的煤炭开采量，dist(x,y)为栅格影像像元点到采矿边界的最短距离；
[0016]
b4、数据归一化处理：将所有驱动因子的量化结果按照如下公式进行归一化处理，其中，x
norm
为归一化后的数据，x为各驱动因子量化后的数据，x
min
为各驱动因子量化后数据中的最小值，x
max
为各驱动因子量化后数据中的最大值；
[0017]
优选地，步骤b1的气候气象数据来源包括中国气象局在内的气候气象数据；步骤b2中数字高程模型数据来源包括地理空间数据云；步骤b3中地理行政边界数据来源包括国家地理信息公共服务平台。
[0018]
c、在三维空间上进行m
1-m2时期内驱动因子数据集的空间维度叠加及高维度拓展，m
1-m2时期内驱动因子数据集所对应的数据为m
1-m2时期研究区原始数据，高维度拓展包括时间维度拓展，并将以栅格影像格式展现的低维数据转换为高维数据立方体，最终构成生态演变大数据立方体；
[0019]
d、采用滑动立方体法对生态演变大数据立方体进行数据提取，将生态演变大数据立方体中的驱动因子按照研究需求划分为自变量参数集合和因变量参数集合两个集合，自变量参数集合包括气候气象、地形地貌、人类活动数据，因变量参数集合包括植被参数；构建地理时空加权人工神经网络模型，通过地理时空加权人工神经网络模型进行模型训练，
在每个时空节点上都构筑出自变量与因变量间的非线性复杂定量关系；
[0020]
e、构建高维数据下的驱动因子定权法量化各驱动因子的权重，并首先对m
1-m2时期下每个空间位置的驱动因子计算其在高维神经网络传播中的放大率，然后对所有量化因子权重进行归一化，并构筑生态演变驱动因子的权重立方体；
[0021]
f、按照步骤b的方法得到m
0-m1时期下的驱动因子数据集及除采矿驱动因子外的驱动因子量化数据，m
0-m1时期为无开采活动时期，m
0-m1时期内驱动因子数据集所对应的数据为m
0-m1时期研究区原始数据；按照b33方法将m
1-m2时期所对应的采矿活动数据代入m
0-m1时期中进行m
0-m1时期的采矿驱动因子所对应采矿活动数据的假设性量化，并得到m
0-m1时期采矿驱动因子虚拟量化数据；根据m
0-m1时期采矿驱动因子虚拟量化数据、m
0-m1时期下除采矿驱动因子外的驱动因子量化数据按照步骤e的驱动因子定权法得到m
0-m1时期采矿驱动因子虚拟权重；
[0022]
g、对m
0-m1时期与m
1-m2时期的采矿驱动因子权重进行像元显著性检验判断栅格影像像元是否受采矿活动扰动影响显著；由此得到m
1-m2时期单个像元采矿对植被扰动影响的显著度情况，然后对研究区域栅格影像的所有像元进行遍历，获得研究区域内受采矿扰动影响显著的区域，进而确定采矿对植被扰动的影响范围。
[0023]
为了更好地实现本发明，本发明步骤d还包括如下方法d11：d11、cv交叉验证法：采用k倍折叠交叉验证法将生态演变大数据立方体的数据集分成k份，循环抽取1份作为验证数据集，其他k-1份作为训练集，进行k次循环，取k次训练的平均mse作为损失函数，如下式所示，以损失函数值最小的滑动窗口带宽值作为最优带宽值，滑动窗口的带宽长度l取最优带宽值；其中，yi与分别代表模型的预测值与标签。
[0024]
优选地，本发明步骤g中像元显著性检验方法如下：g1、设定原假设与备择假设的两种假设模型，具体如下：原假设h0为分析像元不受采矿的扰动影响：μ≤μ0；备择假设h1为分析像元受采矿的扰动影响：μ＞μ0。其中，μ表示m
1-m2时期采矿驱动因子权重的真值，该真值为k次实验结果的均值逼近，μ0表示原假设m
0-m1时期采矿驱动因子虚拟权重的均值；
[0025]
g2、按照如下公式计算检验统计量：其中，z为检验统计量，为k次实验后m
1-m2时期采矿驱动因子权重样本的均值，μ0为原假设m
0-m1时期采矿驱动因子虚拟权重的均值，σ为m
0-m1时期采矿驱动因子虚拟权重数据的总体标准差，n为m
1-m2时期采矿驱动因子权重数据的样本量；总体标准差σ的公式如下：其中，xi为单个像元k次实验的m
1-m2时期采矿驱动因子权重值，为k次实验后m
1-m2时期采矿驱动因子权重样本的均值，n为m
1-m2时期采矿驱动因子权重数据的样本量；
[0026]
g3、确定拒绝域的形式：当检验统计量取到拒绝原假设h0区域中的值时，则该区域为拒绝域，拒绝域的边界点即为临界点，所以z≥z
α
为拒绝域，z＝z
α
为临界点；其拒绝域公式
如下：其中，z为检验统计量，z
0.05
为α等于0.05时检验统计量的临界值，对于右边检验来说查找标准正态分布表可知该值为1.96，为k次实验后m
1-m2时期采矿驱动因子权重样本的均值，μ0为原假设m
0-m1时期采矿驱动因子虚拟权重的均值，σ为m
0-m1时期采矿驱动因子虚拟权重数据的总体标准差，n为m
1-m2时期采矿驱动因子权重数据的样本量；若存在z≥z
0.05
的情况即z的值落在了拒绝域中，则在显著性水平α＝0.05下拒绝了原假设h0，即认为该像元受采矿的扰动影响显著。
[0027]
优选地，本发明步骤g中采矿对植被扰动的影响范围确定如下：
[0028]
g4、对整个研究区域逐像元执行采矿活动对植被扰动的影响显著度判断，并获得每个像元是否属于受采矿扰动影响显著的区域，将这些受采矿扰动影响显著的像元选择出来；g5、确定整个研究区域受采矿扰动影响显著的最外围像元位置并连接成闭合的曲线，该曲线为采矿对植被扰动的影响范围。
[0029]
优选地，本发明步骤c具体方法如下：c1、在统一空间位置下，将研究区m
1-m2时期不同驱动因子的栅格影像数据和矢量数据在三维空间的z轴方位上，进行数据的空间维度叠加与高维度拓展，完成由低维空间数据至高维数据立方体的转换，其中栅格影像数据和矢量数据可视为在空间直角系统的x轴与y轴上拓展的二维平面，最终构成此空间位置下的生态演变大数据立方体；c2、在生态演变大数据立方体的三维空间中，每一层代表一个驱动因子，驱动因子包括降水驱动因子、气温驱动因子、dem驱动因子、放牧驱动因子、城镇驱动因子以及采矿驱动因子，单层驱动因子的厚度代表该驱动因子的时间序列数据，查询生态演变大数据立方体上的某一点时，可以获取到该驱动因子的时间序列数据变化趋势。
[0030]
优选地，本发明步骤d中滑动立方体法如下：d1、在生态演变大数据立方体的栅格影像上建立一个滑动窗口，滑动窗口的带宽长度为l，滑动窗口的步长为s，滑动窗口的步长s≤滑动窗口的带宽长度l，然后对滑动窗口进行高维拓展，增加一个时间维度构建三维的滑动立方体，滑动立方体的时间窗口宽度为t，滑动立方体将会对时间序列所有的整幅栅格影像进行逐像元遍历，根据建立滑动立方体的栅格影像所属的集合不同，滑动立方体单次覆盖的范围内提取的像元将被化分成单位自变量样本和单位因变量样本，当滑动立方体对整幅栅格影像遍历完成，单位样本将分别组合成为自变量参数集合与因变量参数集合。优选地，本发明步骤d中地理时空加权人工神经网络模型方法如下：d2、地理时空加权人工神经网络模型构筑自变量与因变量之间的关系网络，该关系网络包括输入层、隐藏层、输出层三层结构，自变量样本从输入层进行数据输入，然后在输入层传导进入隐藏层时，计算公式如下式所示：其中，w
ij
为神经元i和j之间的连接权值，pi为神经元i的输出，sj为与神经元j有外向连接的神经元集合；神经元i的输出计算如下所示：pi＝φ(layerj)；其中，φ为激活函数，在神经元内进行激活函数的运算，神经元中采用的激活函数为非线性的双曲正切函数，其公式如下式所示：其中f(x)为神经
元激活后的传递值，x为神经元激活前的参数值；
[0031]
以对应自变量样本的样本值作为目标值ti与神经网络输出值pi进行误差计算，计算公式如下所示：其中ri是目标值，pi是输出神经元i的输出，n是目标值的数目，di为时空权重值；其中时空权重值计算如下所示：其中u0，v0，t0为滑动立方体范围内中心像元的三维坐标值，l为滑动立方体的带宽长度，t为滑动立方体的时间窗口宽度；
[0032]
误差计算完成后，后向传播的误差信号计算公式如下所示：
[0033]
其中pj是神经元j的输出，rj是神经元j的目标值，w
jk
是神经元j和k之间的连接权重，δk是神经元k的误差信号，layerj是神经元j的网络输入，并且φ’是激活函数的导数。
[0034]
优选地，本发明步骤e中驱动因子定权法如下：
[0035]
e1、从训练时所用自变量参数集合随机选取其中1个自变量因子，给其的训练数据添加一个偏置增量δx构建一个待定权自变量因子x δx，与其他自变量因子作为已经训练好的地理时空加权人工神经网络模型的输入层，通过模型计算得到对应空间位置下的因变量因子y δy；
[0036]
e2、根据步骤e1所计算得到的待定权因子x δx，与对应因变量因子y δy，计算网络传播放大率w(也即待定驱动因子的权重值)，其计算公式为：
[0037]
e3、进行n次循环，重复步骤d1，直至所有自变量因子权重都已计算，将所有自变量因子的权重进行求和，最终得到一张因子权重之和的栅格影像，然后进行权重的归一化计算，如下式所示：其中wi为自变量因子i的权重值，g(wi)为自变量因子i归一化后的权重。
[0038]
优选地，本发明步骤a1中的ρ
nir
在landsat-5/7中为波段4，在landsat-8中为波段5，ρ
red
在landsat-5/7中为波段3，在landsat-8中为波段4；步骤a1得到的landsat-5/7所对应归一化植被指数ndvi采用最小二乘拟合方法统一拟合校正至landsat-8，接着再计算植被覆盖度。
[0039]
本发明较现有技术相比，具有以下优点及有益效果：
[0040]
(1)本发明利用高维数据下的因子定权法量化出气候气象、地形地貌、人类活动三类驱动因子的权重，其中气候气象包括气温和降水，地形地貌包括数字高程模型，人类活动包括放牧、城镇和采矿活动数据，然后剥离出采矿驱动因子的影响权重进行分析，分析采矿
因子的影响权重，这种权重是对采矿影响的分离，剔除了气温、放牧、放牧活动等其他因素对植被的影响，避免了多因素耦合造成的扰动范围识别误差，旨在识别采矿对植被的扰动范围，为挖掘矿区采矿活动对生态环境的影响机制、保护矿区生态环境等提供数据支持。
[0041]
(2)本发明通过对无开采活动时期的采矿驱动因子权重做出假设，由于无开采活动时期并没有采矿数据，所以要将有开采活动时期的采矿数据代入到无开采活动时期中，量化出无开采活动时期的采矿驱动因子虚拟权重数据，然后对有开采活动时期的采矿驱动因子权重进行显著性检验分析，来获取采矿对植被扰动的影响范围。
[0042]
(3)本发明是对整个研究区域栅格影像中的像元进行分析，具有空间连续性，减少误差的产生，而且考虑到多因素耦合的情况，单独提取分析采矿驱动因子的影响，结果更具有合理性。
[0043]
(4)本发明得到的驱动因子影响权重是考虑到像元级的分析，栅格影像上的每个像元都有各驱动因子的权重值，通过显著性检验方法得到的采矿对植被扰动的影响范围具有空间连续性，有效避免了常规方法在分析样区阈值点与阈值点之间的空间不连续造成的研究盲区，并且避免了常规方法需要人为设置阈值点与阈值点间连接方式而引入的误差。
[0044]
(5)本发明利用显著性检验方法得到的采矿对植被扰动的影响范围更具有准确性，在显著性水平α下，求出原假设(分析像元不受采矿的扰动影响)为真时得到样本的概率，若α《0.05，则分析像元受采矿的扰动影响显著，影响显著的像元集合就是采矿对植被扰动明显的区域范围；本发明方法避免了常规方法(傅里叶分析\缓冲区分析\函数型主成分分析\趋势线拟合)在分析样区拟合趋势时人为确定阈值引入的误差。
附图说明
[0045]
图1为本发明的流程原理图；
[0046]
图2为现有研究中阈值点确定的不合理性示意图；
[0047]
图3为现有研究中阈值点与阈值点之间的研究盲区示意图；
[0048]
图4是本发明中矿区生态演变大数据立方体的示意图；
[0049]
图5为实施例中地理时空加权人工神经网络模型的可视化表达示意图；
[0050]
图6为本发明无开采活动时期的驱动因子集合示意图；
[0051]
图7为本发明无开采活动时期的采矿权重结果示意图；
[0052]
图8为本发明m
0-m1时期采矿驱动因子权重的正态分布以及拒绝域示意图；
[0053]
图9为本发明采矿对植被扰动的影响范围示意图；
[0054]
图10为本发明与现有研究在剔除多因素耦合方面的区别效果示意图；
[0055]
图11为本发明与现有研究在阈值点确定方面的区别效果示意图；
[0056]
图12为本发明与现有研究在空间连续性方面的区别效果示意图；
[0057]
图13为实施例中部分年份放牧驱动因子量化结果示意图；
[0058]
图14为实施例中部分年份城镇驱动因子量化结果示意图；
[0059]
图15为实施例中部分年份采矿驱动因子量化结果示意图；
[0060]
图16为实施例中部分年份基于权重采矿对植被扰动的影响范围图。
具体实施方式
[0061]
下面结合实施例对本发明作进一步地详细说明：
[0062]
实施例
[0063]
如图1～图16所示，一种基于权重的采矿植被生态累积效应扰动范围识别方法，其方法如下：
[0064]
a、收集包括landsat系列卫星影像产品、sentinel-2a影像产品在内的研究区原始数据，优选，本实施例研究区原始数据为栅格化影像。landsat系列卫星影像产品对应为landsat系列遥感影像，sentinel-2a影像产品对应为sentinel-2a遥感影像。本实施例确定研究矿区为锡林浩特市胜利一号矿区，根据各矿业公司煤炭开采量数据，由于2004-2020年为采矿活动时期，则选择无开采活动的时间段1990-2003年，有开采活动的时间段2004-2020进行研究。采集研究矿区的1990-2020年逐年的landsat-5、landsat-7、landsat-8卫星影像产品以及sentinel-2a影像产品，利用收集到的数据对植被参数进行反演。
[0065]
a1、植被参数反演：在google earth engine(简称gee)平台，加载landsat-5、landsat-7、landsat-8卫星影像产品；首先采用如下公式计算归一化植被指数ndvi：其中，ρ
nir
为近红外波段地表反射率，在landsat-5/7中为波段4，在landsat-8中为波段5。ρ
red
为红波段地表反射率，在landsat-5/7中为波段3，在landsat-8中为波段4。接着采用像元二分模型来计算植被覆盖度，其计算公式如下：其中，fvc为像元的植被覆盖度，ndvi为像元的ndvi值，ndvi
min
为研究区内完全为裸土的像元ndvi值，ndvi
max
为研究区纯植被像元的ndvi值；植被覆盖度为反演后的植被参数。
[0066]
根据本实施例的一个优选实施例，本实施例步骤a1中的ρ
nir
在landsat-5/7中为波段4，在landsat-8中为波段5，ρ
red
在landsat-5/7中为波段3，在landsat-8中为波段4。步骤a1得到的landsat-5/7所对应归一化植被指数ndvi采用最小二乘拟合方法统一拟合校正至landsat-8，接着再计算植被覆盖度。本实施例最小二乘拟合方法原理如下：采样截取研究区包含多种地表类型的小块区域，选择landsat-7/8相近日期影像反演得到的植被参数结果，随机选取约500个随机点，运用最小二乘原理构建两种影像得到植被参数的数学拟合模型；应用此模型将步骤a1中由landast-7得到的ndvi结果校正到landsat-8得到的结果上。
[0067]
b、构建驱动因子数据集，驱动因子数据集包括气候气象驱动因子集、地形地貌驱动因子集、人类活动驱动因子集三类，气候气象驱动因子集包括降水驱动因子与气温驱动因子，地形地貌驱动因子集包括dem驱动因子，人类活动驱动因子集包括放牧驱动因子、城镇驱动因子以及采矿驱动因子。
[0068]
b1、采集研究区内气候气象数据，气候气象数据包括降水数据与气温数据，气候气象驱动因子集的降水驱动因子对应气候气象数据的降水数据，气温驱动因子对应气候气象数据的气温数据。本实施例可以通过中国气象局、google earth engine(简称gee)平台、地
理空间数据云、锡林浩特市的统计年鉴以及各矿业公司煤炭开采量年统计数据获取驱动因子并进行量化。将气候气象数据与步骤a1中的植被参数按照如下公式进行皮尔逊相关性分析并得到降水驱动因子与气温驱动因子分别量化对应的皮尔逊相关系数：
[0069]
其中，r为皮尔逊相关系数，n为每个变量中需要进行分析的数据量，xi为降水数据或气温数据的值，为降水数据或气温数据的平均值，yi为植被参数的值，为植被参数的平均值；本实施例可以按照皮尔逊相关性分析的公式可以分别进行降水驱动因子与气温驱动因子的量化。
[0070]
b2、采集研究区内地形地貌数据，地形地貌数据包括数字高程模型数据，从数字高程模型数据中裁剪出研究区域的dem数据并对应dem驱动因子。本实施例在arcgis平台上，对数字高程模型数据(digital elevation model，dem)即从地理空间数据云中获取的aster gdem数据集利用裁剪工具，裁剪出研究区域内的dem数据并对应作为驱动因子。
[0071]
b3、获取地理行政边界数据并辅以landsat影像综合提取识别行政边界与采矿边界。
[0072]
b31、按照如下公式得到放牧驱动因子量化所对应的放牧强度x
graze
：其中，x
graze
为放牧强度，x
牲畜
为研究区域内牲畜数量，x
area
为研究区域内村落总面积；本实施例部分年份放牧强度的量化结果如图13所示。
[0073]
本实施例在arcgis平台上，对人类活动数据即放牧、城镇以及采矿活动数据进行量化，首先加载由国家地理信息公共服务平台(天地图)中下载的乡镇级行政边界以及由landsat系列影像目视解译识别出1990-2020年的城镇以及2004-2020年的采矿边界，并通过锡林浩特市的统计年鉴与当地的煤炭公司获取城镇人口以及煤炭开采量数据。对放牧活动的量化如下：从锡林浩特市统计局中获取锡林浩特市各村落每年的牛、马、羊数据，并且依据牲畜之间的数量转化关系，将牛和马的数量全部转换为羊的数量，转换公式为x
羊
＝5*x
牛/马
，其中x
羊
指的是羊的数量，x
牛/马
指的是牛和马的数量，x
牲畜
＝x
羊
x
牛/马
。
[0074]
b32、采用欧氏距离获得研究区域内每个像元点到城镇边界的最短距离并结合研究区域的人口量化出城镇驱动因子，其中欧式距离的计算公式如下：
[0075]
其中，n为像元点的个数，xi为各像元点的位置，yi为城镇的像元点位置。获得研究区内每个像元点到城镇边界的最短距离然后结合锡林浩特市人口量化出城镇活动的影响，对城镇驱动因子的量化公式如下：
[0076]
其中，x
urban
为城镇活动的量化结果，x
pop
为城镇的人口数，dist(x,y)为栅格影像像元点到城镇边界的最短距离，本实施例部分年份城镇活动量化结果如图14所示。
[0077]
b33、假设m0为没有采矿的年份，m1为开始采矿的年份，m2为结束采矿的年份，则对m1-m2时期采矿驱动因子所对应的采矿活动数据进行量化(本实施例研究区域2004-2020年才有煤炭开采活动，则对2004-2020年的采矿活动数据进行量化)。采用欧氏距离获得研究区域内每个像元点到采矿边界的最短距离并结合研究区域每年的煤炭开采量量化出采矿活动的影响，其中欧式距离的计算公式如下：其中，n为像元点的个数，xi为各像元点的位置，yi为采矿边界的像元点位置。对采矿活动的量化公式如下：其中，x
mine
为采矿活动的量化结果，x
mining
为每年的煤炭开采量，dist(x,y)为栅格影像像元点到采矿边界的最短距离；本实施例部分年份采矿活动量化结果如图15所示。
[0078]
b4、数据归一化处理：将所有驱动因子的量化结果按照如下公式进行归一化处理，其中，x
norm
为归一化后的数据，x为各驱动因子量化后的数据，x
min
为各驱动因子量化后数据中的最小值，x
max
为各驱动因子量化后数据中的最大值。
[0079]
优选地，本实施例步骤b1的气候气象数据来源包括中国气象局在内的气候气象数据。本实施例步骤b2中数字高程模型数据来源包括地理空间数据云。本实施例步骤b3中地理行政边界数据来源包括国家地理信息公共服务平台。
[0080]
c、在三维空间上进行m
1-m2时期内驱动因子数据集的空间维度叠加及高维度拓展，m
1-m2时期内驱动因子数据集所对应的数据为m
1-m2时期研究区原始数据，高维度拓展包括时间维度拓展，并将以栅格影像格式展现的低维数据转换为高维数据立方体，最终构成生态演变大数据立方体。
[0081]
根据本实施例的一个优选实施例，本实施例步骤c具体方法如下：
[0082]
c1、本实施例步骤a和b所提供的以栅格影像作为基本存储单元的参数定量遥感植被反演数据，构成矿区长时间尺度、连续空间的生态因子及气候气象、地形地貌、人类活动三类驱动因子数据集。在统一空间位置下，将研究区m
1-m2时期(本实施例为2004-2020年)不同驱动因子的栅格影像数据和矢量数据在三维空间的z轴(垂直于二维平面的高程轴)方位上，进行数据的空间维度叠加与高维度拓展，完成由低维空间数据至高维数据立方体的转换，其中栅格影像数据和矢量数据可视为在空间直角系统的x轴与y轴上拓展的二维平面，最终构成此空间位置下的生态演变大数据立方体。
[0083]
c2、在生态演变大数据立方体的三维空间中，每一层代表一个驱动因子(在本实例中为植被覆盖度)，驱动因子包括降水驱动因子、气温驱动因子、dem驱动因子、放牧驱动因子、城镇驱动因子以及采矿驱动因子，单层驱动因子的厚度代表该驱动因子的时间序列数据(如从2004年到2020年，每一年该生态参数因子的栅格影像数据)，查询生态演变大数据立方体上的某一点时，可以获取到该驱动因子的时间序列数据变化趋势。
[0084]
d、采用滑动立方体法对生态演变大数据立方体进行数据提取，将生态演变大数据立方体中的驱动因子按照研究需求划分为自变量参数集合和因变量参数集合两个集合，自变量参数集合包括气候气象、地形地貌、人类活动数据，因变量参数集合包括植被参数。构建地理时空加权人工神经网络模型，通过地理时空加权人工神经网络模型进行模型训练，
在每个时空节点上都构筑出自变量与因变量间的非线性复杂定量关系。
[0085]
根据本实施例的一个优选实施例，本实施例步骤d还包括如下方法d11：
[0086]
d11、cv交叉验证法(cross validation，cv)：采用k倍折叠交叉验证法将生态演变大数据立方体的数据集分成k份(此处k值可变，按照经验模型一般定义为10)，循环抽取1份作为验证数据集，其他k-1份作为训练集，进行k次循环，取k次训练的平均mse作为损失函数，如下式所示，以损失函数值最小的滑动窗口带宽值作为最优带宽值，滑动窗口的带宽长度l取最优带宽值。
[0087]
其中，yi与分别代表模型的预测值与标签。
[0088]
本实施例经过迭代训练，反向传播修正连接权重之后，得到了一个描述由若干自变量与单个因变量关系的地理时空加权神经网络模型gtwnet。
[0089]
根据本实施例的一个优选实施例，本实施例步骤d中滑动立方体法如下：
[0090]
d1、在生态演变大数据立方体的栅格影像上建立一个滑动窗口(本实施例可以基于python环境)，滑动窗口的带宽长度为l(本实例中的滑动窗口的带宽长度l可以为定值，通过指定带宽长度或者临近栅格像元数目的方式，获取相应大小的滑动窗口，也可以为自适应值，通过输入带宽范围采用交叉验证法(cv)来决定最佳带宽)，滑动窗口的步长为s(本实例中的滑动窗口的带宽长度s为定值1)，滑动窗口的步长s≤滑动窗口的带宽长度l，然后对滑动窗口进行高维拓展，增加一个时间维度构建三维的滑动立方体，滑动立方体的时间窗口宽度为t(本实例中的时间窗口宽度为一个定值t，取决于栅格影像时间序列的时间长度)，滑动立方体将会对时间序列所有的整幅栅格影像进行逐像元遍历，根据建立滑动立方体的栅格影像所属的集合不同，滑动立方体单次覆盖的范围内提取的像元将被化分成单位自变量样本和单位因变量样本(例如栅格影像属于自变量参量集合，则在此栅格影像上建立的滑动立方体单次覆盖的范围提取的像元将成为自变量样本)，当滑动立方体对整幅栅格影像遍历完成，单位样本将分别组合成为自变量参数集合与因变量参数集合。
[0091]
根据本实施例的一个优选实施例，本实施例步骤d中地理时空加权人工神经网络模型方法如下：
[0092]
d2、地理时空加权人工神经网络模型构筑自变量与因变量之间的关系网络(本实施例可基于pytorch搭建三层人工神经网络结构-输入层、隐藏层、输出层三层结构并构筑自变量与因变量之间的关系)，该关系网络包括输入层、隐藏层、输出层三层结构(如图5所示)，自变量样本从输入层进行数据输入，然后在输入层传导进入隐藏层时，计算公式如下式所示：
[0093]
其中，w
ij
为神经元i和j之间的连接权值，pi为神经元i的输出，sj为与神经元j有外向连接的神经元集合。
[0094]
神经元i的输出计算如下所示：
[0095]
pi＝φ(layerj)；其中，φ为激活函数，在神经元内进行激活函数的运算，神经元中采用的激活函数为非线性的双曲正切函数，其公式如下式所示：
[0096]
其中f(x)为神经元激活后的传递值，x为神经元激活前的参数值。以对应自变量样本的样本值作为目标值ti与神经网络输出值pi进行误差计算，计算公式如下所示：其中ri是目标值，pi是输出神经元i的输出，n是目标值的数目，di为时空权重值。其中时空权重值计算如下所示：
[0097]
其中u0，v0，t0为滑动立方体范围内中心像元的三维坐标值，l为滑动立方体的带宽长度，t为滑动立方体的时间窗口宽度。
[0098]
误差计算完成后，后向传播的误差信号计算公式如下所示：
[0099]
其中pj是神经元j的输出，rj是神经元j的目标值，w
jk
是神经元j和k之间的连接权重，δk是神经元k的误差信号，layerj是神经元j的网络输入，并且φ’是激活函数的导数。
[0100]
e、构建高维数据下的驱动因子定权法量化各驱动因子的权重，并首先对m
1-m2时期下每个空间位置的驱动因子计算其在高维神经网络传播中的放大率，然后对所有量化因子权重进行归一化，并构筑生态演变驱动因子的权重立方体。
[0101]
根据本实施例的一个优选实施例，本实施例步骤e中驱动因子定权法如下：
[0102]
e1、从训练时所用自变量参数集合(自变量数目n根据研究需求决定，n＞＝1)随机选取其中1个自变量因子，给其的训练数据添加一个偏置增量δx(在发明中，即给该因子栅格影像的每一个像元都要加上一个偏置增量
△
x)构建一个待定权自变量因子x δx，与其他自变量因子作为已经训练好的地理时空加权人工神经网络模型的输入层，通过模型计算得到对应空间位置下的因变量因子y δy。
[0103]
e2、根据步骤e1所计算得到的待定权因子x δx，与对应因变量因子y δy，计算网络传播放大率w(即待定驱动因子的权重值)，其计算公式为：其中，w为待定驱动因子的权重值，
△
y为因变量因子的偏置，
△
x为自变量因子的偏置。
[0104]
e3、进行n次循环，重复步骤d1，直至所有自变量因子权重都已计算，将所有自变量因子的权重进行求和，最终得到一张因子权重之和的栅格影像，然后进行权重的归一化计算，如下式所示：其中wi为自变量因子i的权重值，g(wi)为自变量因子i归一化后的权重。
[0105]
本实施例按照步骤c的流程，将植被参数(即植被覆盖度)遥感反演数据更换为驱动因子归一化权重数据，其它流程与步骤c类似，此处不再赘述，最终得到以驱动因子归一化权重数据为主体的生态演变驱动因子权重立方体ew-cub。本实施例步骤e可获取有开采
活动时(即2004-2020年)各驱动因子的权重归一化结果，由于人工神经网络的普遍特性其结果具有微小的随机性，为了避免单次实验的偶然性，本实施例在同一栅格像元中进行10次gtwann-w方法实验，取10次采矿驱动因子权重归一化结果的平均值为最终的结果值
[0106]
f、按照步骤b的方法得到m
0-m1时期下的驱动因子数据集及除采矿驱动因子外的驱动因子量化数据，m
0-m1时期为无开采活动时期，m
0-m1时期内驱动因子数据集所对应的数据为m
0-m1时期研究区原始数据。按照b33方法将m
1-m2时期所对应的采矿活动数据代入m
0-m1时期中进行m
0-m1时期的采矿驱动因子所对应采矿活动数据的假设性量化，并得到m
0-m1时期采矿驱动因子虚拟量化数据。本实施例根据步骤b可以量化1990-2003年(m
0-m1时期)的气温、降水、高程、放牧和城镇活动驱动因子数据；由于1990-2003年(m
0-m1时期)并没有采矿活动，所以要将2004-2020年(m
1-m2时期)有开采活动时期的煤炭开采量数据代入到1990-2003年中进行假设性采矿活动数据的量化。根据m
0-m1时期采矿驱动因子虚拟量化数据、m
0-m1时期下除采矿驱动因子外的驱动因子量化数据按照步骤e的驱动因子定权法得到m
0-m1时期采矿驱动因子虚拟权重。本实施例输入采矿数据就能计算出采矿驱动因子的权重，所以对于1990-2003年无开采活动时期的采矿驱动因子也能计算出权重，这种权重就是引入的噪声(即模拟采矿环境下无开采活动时期的采矿权重)。具体方法如下：
[0107]
f1、量化无开采活动时期(1990-2003年)的驱动因子数据：基于arcgis平台，根据步骤b量化1990-2003年的驱动因子数据，构建无开采活动时期的驱动因子集合。由于1990-2003年并没有煤炭开采活动，所以在同一采矿范围中将2004-2020年的煤炭开采量代入到1990-2003年中，对1990-2003年的采矿驱动因子进行假设性采矿活动数据的量化。
[0108]
首先假设1990-2003年存在煤炭开采活动，将2004-2020年的采矿边界以及每年的煤炭开采量代入到1990-2003年中，使用欧氏距离获得研究区域内每个像元点到采矿边界的最短距离然后结合研究区域每年的煤炭开采量量化出采矿活动的影响。其中欧式距离的计算公式如下：其中，n为像元点的个数，xi为各像元点的位置，yi为采矿边界的像元点位置。
[0109]
对采矿活动的量化公式如下：其中，x
mine
为采矿活动的量化结果，x
mining
为每年的煤炭开采量，dist(x,y)为栅格影像像元点到采矿边界的最短距离。
[0110]
f2、计算无开采活动时期(1990-2003年)的采矿权重(噪声)：得到1990-2003年的驱动因子量化数据后，作为输入数据代入到步骤e高维数据下的因子定权法(gtwann-w)中计算1990-2003年各驱动因子权重，获得假设存在煤炭开采条件下的采矿驱动因子虚拟影响权重即噪声数据。
[0111]
g、对m
0-m1时期与m
1-m2时期的采矿驱动因子权重进行像元显著性检验判断栅格影像像元是否受采矿活动扰动影响显著。由此得到m
1-m2时期单个像元采矿对植被扰动影响的显著度情况，然后对研究区域栅格影像的所有像元进行遍历，获得研究区域内受采矿扰动影响显著的区域，进而确定采矿对植被扰动的影响范围。
[0112]
根据本实施例的一个优选实施例，本实施例步骤g中像元显著性检验方法如下：
[0113]
g1、设定原假设与备择假设的两种假设模型，通常情况下，实验结果α》0.05表示差异性不显著；α《0.05表示差异性显著；其中，α表示用于确定假设检验结果的参数，称为显著性水平；对本实施例来说，若α《0.05，则分析像元受采矿的扰动影响显著。具体如下：原假设h0为分析像元不受采矿的扰动影响：μ≤μ0；备择假设h1为分析像元受采矿的扰动影响：μ＞μ0。其中，μ表示m
1-m2时期采矿驱动因子权重的真值，该真值为k次实验结果的均值逼近，μ0表示原假设m
0-m1时期采矿驱动因子虚拟权重的均值；对于本实施例，其中，μ表示2004-2020年采矿驱动因子权重的真值，该真值可以用多次实验结果的均值逼近，本发明定为10次实验，μ0表示原假设1990-2003年采矿驱动因子虚拟权重的均值。
[0114]
g2、给定显著性水平α以及确定检验统计量：经过发明人反复的实验与研究，发现1990-2003年采矿驱动因子虚拟权重即噪声数据集呈现正态分布，本实施例中给定显著性水平α为0.05。首先利用公式计算出2004-2020年的检验统计量，其中，z为检验的统计量，为10次实验后2004-2020年采矿驱动因子权重样本的均值，μ0为原假设1990-2003年采矿驱动因子虚拟权重的均值，σ为1990-2003年采矿驱动因子虚拟权重数据的总体标准差，n为2004-2020年采矿驱动因子权重数据的样本量。总体标准差的公式为其中xi为单个像元10次实验的2004-2020年采矿驱动因子权重值，为10次实验后2004-2020年采矿驱动因子权重样本的均值，n为2004-2020年采矿驱动因子权重数据的样本量。
[0115]
给定显著性水平α以及确定检验统计量：经过发明人反复的实验与研究，发现1990-2003年采矿驱动因子虚拟权重即噪声数据集呈现正态分布，本实施例中给定显著性水平α为0.05。首先利用公式计算出2004-2020年的检验统计量，其中，z为检验的统计量，为10次实验后2004-2020年采矿驱动因子权重样本的均值，μ0为原假设1990-2003年采矿驱动因子虚拟权重的均值，σ为1990-2003年采矿驱动因子虚拟权重数据的总体标准差，n为2004-2020年采矿驱动因子权重数据的样本量。总体标准差的公式为其中xi为单个像元10次实验的2004-2020年采矿驱动因子权重值，为10次实验后2004-2020年采矿驱动因子权重样本的均值，n为2004-2020年采矿驱动因子权重数据的样本量。
[0116]
根据本实施例的一个优选实施例，本实施例步骤g中采矿对植被扰动的影响范围确定如下(本实施例遍历整个研究区域栅格影像的像元，确定采矿对植被扰动的影响范围：根据步骤f可以得到2004-2020年单个像元采矿对植被扰动影响的显著度情况，在此基础上对研究区栅格影像的所有像元进行遍历，获得研究区域受采矿扰动影响显著的区域，确定
采矿对植被扰动的影响范围)：
[0117]
g4、对整个研究区域逐像元执行采矿活动对植被扰动的影响显著度判断，并获得每个像元是否属于受采矿扰动影响显著的区域，将这些受采矿扰动影响显著的像元选择出来。
[0118]
g5、确定整个研究区域受采矿扰动影响显著的最外围像元位置并连接成闭合的曲线，该曲线为采矿对植被扰动的影响范围；本实施例部分年份采矿影响范围如图16所示。
[0119]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。