本发明公开了基于长时序多源数据的煤炭开采植被扰动分析方法，首先，对植被参数进行长时间尺度上高频次定量反演，根据遥感反演、统计数据获取矿区长时间尺度、连续空间的植被参数及气候气象因子、地理因子和人类活动因子数据集；然后基于采矿前的长时序多源数据，利用地理时空加权回归模型构建植被变化的理论驱动模型；最后利用上述模型预测无采矿活动条件下的植被演变过程，进而与遥感监测的采矿活动背景下的实际植被演变进行对比，分离出煤炭开采对植被的扰动量V‑MD。本发明可得到煤炭开采对植被的扰动量V‑MD，能够分离和量化煤炭开采活动对植被的影响，并揭示了不同采矿阶段的演变规律，为矿区的生态环境保护提供了理论数据支持。
1.基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：其方法如下：A、建立长时序多源数据的植被参数反演模型，遥感参数反演模型包括归一化指数模型、像元二分模型、PROSAIL植被辐射传输模型；收集包括多光谱遥感影像、地面实测数据在内的原始数据，通过遥感参数反演模型对原始数据进行反演得到植被参数，植被参数包括归一化植被指数NDVI、植被覆盖度FVC、叶面积指数LAI、植被叶绿素含量Cab；
B、构建无采矿背景下的植被响应模式：对矿区无煤炭开采活动时期的长时序逐年植被参数和驱动因子集合进行建模，驱动因子集合包括气候气象因子、地理因子和人类活动因子，具体如下：气温、降水、风速、日照对应气候气象因子，坡度、坡向、海拔对应地理因子，放牧、城镇开发、发电对应人类活动因子；通过地理时空加权回归模型建立植被参数与各驱动因子之间的定量关系并训练其定量关系；通过地理时空加权回归模型建立植被参数与各驱动因子之间的定量关系公式如下：
式中，(xi1,xi2,…,xid；yi)表示第i个观测点(ui,vi,ti),(i＝1,2,…,n)处的植被参数y和驱动因子x1,x2,…,xd的n组观测值；βk(ui,vi,ti),(k＝0,1,…,d)是第i个数据点(ui,vi,ti)处的未知回归系数；(ε1,ε2,…,εn)为独立同分布的误差项；
C、采矿植被扰动量计算：在矿区有煤炭开采活动时期所收集的原始数据中提取驱动因子集合和实际植被参数，驱动因子集合包括气候气象因子、地理因子和人类活动因子，并将该驱动因子集合输入训练后的植被参数与驱动因子之间定量关系中得到对应的植被参数预测值，然后将实际植被参数与植被参数预测值对照求差并得到植被参数差值，该植被参数差值为煤炭开采而导致的植被变化量，该植被参数差值为采矿植被扰动量，将其命名为V‑MD；通过V‑MD进行时空分布差异和演变量化特征分析与展示。
2.按照权利要求1所述的基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：步骤A中反演得到植被参数方法如下：利用归一化指数模型按照如下公式计算得到归一化植被指数NDVI：其中，ρNIR为近红外波段地表反射率，在Landsat‑5和Landsat‑7中为波段4，在Landsat‑
8中为波段5；ρRed为红波段地表反射率，在Landsat‑5和Landsat‑7中为波段3，在Landsat‑8中为波段4；
利用像元二分模型按照如下公式计算得到植被覆盖度FVC：其中，NDVI为像元的NDVI值，NDVImin为研究区内完全为裸土的像元NDVI值，NDVImax为研究区纯植被像元的NDVI值；
叶面积指数、叶片叶绿素含量计算：采用PROSAIL植被辐射传输模型耦合Landsat系列、Sentinel‑2A卫星传感器光谱响应函数，结合地面实测光谱和参数数据，基于随机森林算法建立植被参数反演模型，在Google Earth Engine平台上生产长时间序列的两种植被参数产品，两种植被参数产品分别为叶面积指数、叶片叶绿素。
3.按照权利要求1所述的基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：步骤A中还包括对反演得到的植被参数进行多源数据一致性校正与精度检验，方法如下：
对于植被参数，采用传感器配准的方法将所有反演结果都以Landsat8为基准；首先选择日期最为接近的Landsat5、Landsat7和Landsat8的一天的影像，分别计算植被参数；然后选择一个小块研究区，将该区域内的像元值提取到点，运用最小二乘原理，分别构建不同影像对应像元之间的线性模型；基于回归线性模型实现不同传感器之间的配准，最后都以Landsat8为基准；多源结果的一致性校正完成之后，再与地面实测数据进行检验。
4.按照权利要求1所述的基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：
根据加权最小二乘方法，第i个观测点的回归系数估计值 为：第i个观测点因变量的拟合值 为：
式中，X为驱动因子矩阵 Xi是X矩阵的第i行；y为植被参数矩阵 Wi表示空间核函数矩阵 常用的核函数包括高斯核函数和近高斯核函数；
调整型高斯核函数下，观测点j对观测点i影响的权重值：调整型近高斯核函数下，观测点j对观测点i影响的权重值：式中:dij表示点(uj,vj,tj)到点(ui,vi,ti)的时空距离；dN表示点(ui,vi,ti)到最近邻第N个点的时空距离；用时空椭球坐标系来计算时空距离：其中，λ和μ为时间和空间距离的比例调整系数。
5.按照权利要求4所述的基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：步骤C中将该驱动因子集合输入模型中预测得到对应的植被参数预测值V预测。
6.按照权利要求5所述的基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：步骤C中采矿植被扰动量V‑MD计算方法如下：V‑MD＝V实际‑V预测，其中，V实际为矿区有煤炭开采活动时期通过遥感反演得到的实际植被参数，V预测为模型预测得到的植被参数预测值。
基于长时序多源数据的煤炭开采植被扰动分析方法
技术领域
[0001] 本发明涉及生态学领域、煤炭领域、遥感及地理信息领域，尤其涉及一种基于长时序多源数据的煤炭开采植被扰动分析方法。
背景技术
[0002] 现有关于煤炭开采对植被的扰动分析主要分为两大类，一是在时间尺度上，研究植被受煤炭开采扰动的变化规律；二是在空间尺度上，通过拟合植被参数随离矿距离的变
化曲线来界定植被受煤炭开采扰动的范围。但现有研究均假定直接监测到的植被状态变化
代表了采矿的影响，但事实上，植被状态变化是自然条件、采矿及其它人类活动等综合作用
的结果，因此当前研究没有将采矿影响和其他因素的影响分离开，导致无法揭示采矿活动
对植被的单独影响规律。
发明内容
[0003] 针对现有技术存在的不足之处，本发明的目的在于提供一种基于长时序多源数据的煤炭开采植被扰动分析方法，基于地理时空加权回归与多源大数据相结合，以局部多元
回归代替整体回归，且充分考虑了各因子的复杂性和时空异质性，能够剔除气候因子、地理
条件、人类活动等多因素对植被的耦合影响，分离并量化煤炭开采对植被的单独影响，即可
得到煤炭开采对植被的扰动量V‑MD，为矿区的生态环境保护提供了理论数据支持。
[0004] 本发明的目的通过下述技术方案实现：
[0005] 一种基于长时序多源数据的煤炭开采植被扰动分析方法，基于长时序多源数据的煤炭开采植被扰动分析方法，其特征在于：其方法如下：
[0006] A、建立长时序多源数据的植被参数反演模型，遥感参数反演模型包括归一化指数模型、像元二分模型、PROSAIL植被辐射传输模型；收集包括多光谱遥感影像、地面实测数据
在内的原始数据，通过遥感参数反演模型对原始数据进行反演得到植被参数，植被参数包
括归一化植被指数NDVI、植被覆盖度FVC、叶面积指数LAI、植被叶绿素含量Cab；
[0007] B、构建无采矿背景下的植被响应模式：对矿区无煤炭开采活动时期的长时序逐年植被参数和驱动因子集合进行建模，驱动因子集合包括气候气象因子、地理因子和人类活
动因子，具体如下：气温、降水、风速、日照对应气候气象因子，坡度、坡向、海拔对应地理因
子，放牧、城镇开发、发电对应人类活动因子；通过地理时空加权回归模型建立植被参数与
各驱动因子之间的定量关系并训练其定量关系；
[0008] C、采矿植被扰动量计算：在矿区有煤炭开采活动时期所收集的原始数据中提取驱动因子集合和实际植被参数，驱动因子集合包括气候气象因子、地理因子和人类活动因子，
并将该驱动因子集合输入训练后的植被参数与驱动因子之间定量关系中得到对应的植被
参数预测值，然后将实际植被参数与植被参数预测值对照求差并得到植被参数差值，该植
被参数差值为煤炭开采而导致的植被变化量，该植被参数差值为采矿植被扰动量，将其命
名为V‑MD；通过V‑MD进行时空分布差异和演变量化特征分析与展示。
[0009] 本发明优选的技术方案如下：本发明步骤A中反演得到植被参数方法如下：
[0010] 利用归一化指数模型按照如下公式计算得到归一化植被指数NDVI：
[0011]
[0012] 其中，ρNIR为近红外波段地表反射率，在Landsat‑5和Landsat‑7中为波段4，在Landsat‑8中为波段5；ρRed为红波段地表反射率，在Landsat‑5和Landsat‑7中为波段3，在
Landsat‑8中为波段4；
[0013] 利用像元二分模型按照如下公式计算得到植被覆盖度FVC：
[0014]
[0015] 其中，NDVI为像元的NDVI值，NDVImin为研究区内完全为裸土的像元NDVI值，NDVImax为研究区纯植被像元的NDVI值；
[0016] 叶面积指数、叶片叶绿素含量计算：采用PROSAIL植被辐射传输模型耦合Landsat系列、Sentinel‑2A卫星传感器光谱响应函数，结合地面实测光谱和参数数据，基于随机森
林算法建立植被参数反演模型，在Google Earth Engine平台上生产长时间序列的两种植
被参数产品，两种植被参数产品分别为叶面积指数、叶片叶绿素。
[0017] 本发明优选的技术方案如下：本发明步骤A中还包括对反演得到的植被参数进行多源数据一致性校正与精度检验，方法如下：
[0018] 对于植被参数，采用传感器配准的方法将所有反演结果都以Landsat8为基准；首先选择日期最为接近的Landsat5、Landsat7和Landsat8的一天的影像，分别计算植被参数；
然后选择一个小块研究区，将该区域内的像元值提取到点，运用最小二乘原理，分别构建不
同影像对应像元之间的线性模型；基于回归线性模型实现不同传感器之间的配准，最后都
以Landsat8为基准；多源结果的一致性校正完成之后，再与地面实测数据进行检验。
[0019] 本发明优选的技术方案如下：本发明步骤B中地理时空加权回归模型建立植被参数与驱动因子之间的定量关系模型公式如下：
[0020]
[0021] 式中，(xi1,xi2,…,xid；yi)表示第i个观测点(ui,vi,ti),(i＝1,2,…,n)处的植被参数y和驱动因子x1,x2,…,xd的n组观测值；βk(ui,vi,ti),(k＝0,1,…,d)是第i个数据点
(ui,vi,ti)处的未知回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定均值为0，方
2
差为σ；
[0022] 根据加权最小二乘方法，第i个观测点的回归系数估计值 为：
[0023]
[0024] 第i个观测点因变量的拟合值 为：
[0025]
[0026] 式中，X为驱动因子矩阵 Xi是X矩阵的第i行；y为植被参数矩阵 Wi表示空间核函数矩阵 常用的核函数包括高
斯核函数和近高斯核函数；
[0027] 调整型高斯核函数下，观测点j对观测点i影响的权重值：
[0028]
[0029] 调整型近高斯核函数下，观测点j对观测点i影响的权重值：
[0030]
[0031] 式中:dij表示点(uj,vj,tj)到点(ui,vi,ti)的时空距离；dN表示点(ui,vi,ti)到最近邻第N个点的时空距离；用时空椭球坐标系来计算时空距离：
其中，λ和μ为时间和空间距离的比例调整系数。
[0032] 本发明优选的技术方案如下：本发明步骤C中将该驱动因子集合输入模型中预测得到对应的植被参数预测值V预测。
[0033] 本发明优选的技术方案如下：本发明步骤C中采矿植被扰动量V‑MD计算方法如下：
[0034] V‑MD＝V实际‑V预测，其中，V实际为矿区有煤炭开采活动时期通过遥感反演得到的实际植被参数，V预测为模型预测得到的植被参数预测值。
[0035] 本发明较现有技术相比，具有以下优点及有益效果：
[0036] (1)本发明提出了分离采煤干扰两阶段处理方法，基于地理时空加权回归与多源大数据相结合，以局部多元回归代替整体回归，且充分考虑了各因子的复杂性和时空异质
性，能够剔除气候因子、地理条件、人类活动等多因素对植被的耦合影响，分离并量化煤炭
开采对植被的单独影响，即可得到煤炭开采对植被的扰动量V‑MD，为矿区的生态环境保护
提供了理论数据支持。
[0037] (2)本发明能够揭示不分离干扰因素所无法呈现出的规律，在一定程度上解决了从矿区植被变化的多因素耦合影响中如何单独分离并量化煤炭开采影响的难题。
[0038] (3)本发明能够分离和量化煤炭开采活动对植被的影响，剖析了采矿影响的时空分布差异，并揭示了不同采矿阶段的演变规律，为矿区的生态环境保护提供了理论数据支
持。
附图说明
[0039] 图1为本发明的流程原理图。
具体实施方式
[0040] 下面结合实施例对本发明作进一步地详细说明：
[0041] 实施例一
[0042] 如图1所示，一种基于长时序多源数据的煤炭开采植被扰动分析方法，其方法如下：
[0043] A、建立长时序多源数据的植被参数反演模型，遥感参数反演模型包括归一化指数模型、像元二分模型、PROSAIL植被辐射传输模型；收集包括多光谱遥感影像、地面实测数据
在内的原始数据，通过遥感参数反演模型对原始数据进行反演得到植被参数，植被参数包
括归一化植被指数NDVI、植被覆盖度FVC、叶面积指数LAI、植被叶绿素含量Cab；
[0044] 步骤A中反演得到植被参数方法如下：
[0045] 利用归一化指数模型按照如下公式计算得到归一化植被指数NDVI：
[0046]
[0047] 其中，ρNIR为近红外波段地表反射率，在Landsat‑5和Landsat‑7中为波段4，在Landsat‑8中为波段5；ρRed为红波段地表反射率，在Landsat‑5和Landsat‑7中为波段3，在
Landsat‑8中为波段4；
[0048] 利用像元二分模型按照如下公式计算得到植被覆盖度FVC：
[0049]
[0050] 其中，NDVI为像元的NDVI值，NDVImin为研究区内完全为裸土的像元NDVI值，NDVImax为研究区纯植被像元的NDVI值；
[0051] 叶面积指数、叶片叶绿素含量计算：采用PROSAIL植被辐射传输模型耦合Landsat系列、Sentinel‑2A卫星传感器光谱响应函数，结合地面实测光谱和参数数据，基于随机森
林算法建立植被参数反演模型，在Google Earth Engine平台上生产长时间序列的两种植
被参数产品，两种植被参数产品分别为叶面积指数、叶片叶绿素。
[0052] B、构建无采矿背景下的植被响应模式：对矿区无煤炭开采活动时期的长时序逐年植被参数和驱动因子集合进行建模，驱动因子集合包括气候气象因子、地理因子和人类活
动因子，具体如下：气温、降水、风速、日照对应气候气象因子，坡度、坡向、海拔对应地理因
子，放牧、城镇开发、发电对应人类活动因子；通过地理时空加权回归模型建立植被参数与
各驱动因子之间的定量关系并训练其定量关系；
[0053] 步骤B中地理时空加权回归模型建立植被参数与驱动因子之间的定量关系模型公式如下：
[0054]
[0055] 式中，(xi1,xi2,…,xid；yi)表示第i个观测点(ui,vi,ti),(i＝1,2,…,n)处的植被参数y和驱动因子x1,x2,…,xd的n组观测值；βk(ui,vi,ti),(k＝0,1,…,d)是第i个数据点
(ui,vi,ti)处的未知回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定均值为0，方
2
差为σ；
[0056] 根据加权最小二乘方法，第i个观测点的回归系数估计值 为：
[0057]
[0058] 第i个观测点因变量的拟合值 为：
[0059]
[0060] 式中，X为驱动因子矩阵 Xi是X矩阵的第i行；y为植被参数矩阵 Wi表示空间核函数矩阵 常用的核函数包括高
斯核函数和近高斯核函数；
[0061] 调整型高斯核函数下，观测点j对观测点i影响的权重值：
[0062]
[0063] 调整型近高斯核函数下，观测点j对观测点i影响的权重值：
[0064]
[0065] 式中:dij表示点(uj,vj,tj)到点(ui,vi,ti)的时空距离；dN表示点(ui,vi,ti)到最近邻第N个点的时空距离；用时空椭球坐标系来计算时空距离：
其中，λ和μ为时间和空间距离的比例调整系数。地理时空加权回归模型输入模型数据的长
时序气候气象因子、地理因子和人类活动因子、植被参数如表1所示：
[0066]
[0067]
[0068] 表1输入模型数据
[0069] C、采矿植被扰动量计算：在矿区有煤炭开采活动时期所收集的原始数据中提取驱动因子集合和实际植被参数，驱动因子集合包括气候气象因子、地理因子和人类活动因子，
并将该驱动因子集合输入训练后的植被参数与驱动因子之间定量关系中得到对应的植被
参数预测值，本实施例将该驱动因子集合输入模型
中预测得到对应的植被参数预测值V预测。
[0070] 然后将实际植被参数与植被参数预测值对照求差并得到植被参数差值，该植被参数差值为煤炭开采而导致的植被变化量，该植被参数差值为采矿植被扰动量，将其命名为
V‑MD；通过V‑MD进行时空分布差异和演变量化特征分析与展示。
[0071] 本实施例步骤C中采矿植被扰动量V‑MD计算方法如下：
[0072] V‑MD＝V实际‑V预测，其中，V实际为矿区有煤炭开采活动时期通过遥感反演得到的实际植被参数，V预测为模型预测得到的植被参数预测值。
[0073] 优选地，本实施例步骤A中还包括对反演得到的植被参数进行多源数据一致性校正与精度检验，方法如下：
[0074] 对于植被参数，采用传感器配准的方法将所有反演结果都以Landsat8为基准；首先选择日期最为接近的Landsat5、Landsat7和Landsat8的一天的影像，分别计算植被参数；
然后选择一个小块研究区，将该区域内的像元值提取到点，运用最小二乘原理，分别构建不
同影像对应像元之间的线性模型；基于回归线性模型实现不同传感器之间的配准，最后都
以Landsat8为基准；多源结果的一致性校正完成之后，再与地面实测数据进行检验。
[0075] 实施例二
[0076] 如图1所示，一种基于长时序多源数据的煤炭开采植被扰动分析方法，其方法如下：
[0077] A、建立长时序多源数据的植被参数反演模型，遥感参数反演模型包括归一化指数模型、像元二分模型、PROSAIL植被辐射传输模型；收集分析区域内包括多光谱遥感影像(包
括Landsat‑5、Landsat‑7、Landsat‑8卫星影像产品，Sentinel‑2A影像产品)、地面实测数据
在内的原始数据，通过遥感参数反演模型对原始数据进行反演得到植被参数，植被参数包
括归一化植被指数NDVI、植被覆盖度FVC、叶面积指数LAI、植被叶绿素含量Cab；
[0078] A1、在Google Earth Engine平台，加载多光谱遥感影像(包括Landsat‑5、Landsat‑7、Landsat‑8卫星影像产品，Sentinel‑2A影像产品)，采用归一化植被指数模型
选取公式中对应影像波段，进行计算归一化植被指数；其中ρNIR为近红
外波段地表反射率，在Landsat‑5和Landsat‑7中为波段4，在Landsat‑8中为波段5；ρRed为红
波段地表反射率，在Landsat‑5和Landsat‑7中为波段3，在Landsat‑8中为波段4，最后获取
分析区域内的归一化植被指数栅格影像数据。
[0079] A2、在ArcGIS平台，加载步骤A1获取的NDVI数据，再利用分区域统计脚本提取NDVImax与NDVImin，然后采用栅格计算器构建像元二分模型 其
中，NDVI为像元的NDVI值，NDVImax为研究区内完全为裸土的像元NDVI值，NDVImin为研究区纯
植被像元的NDVI值，最后输出获取分析区域内的植被覆盖度栅格影像数据。
[0080] B、以步骤A中反演得到的植被覆盖度FVC参数作为因变量，构建无采矿背景下气候气象因子、地理因子和人类活动因子等多种驱动因子与植被覆盖度FVC的定量关系。
[0081] B1、通过国家气象数据网站和研究区统计年鉴获得如下表2所示数据。
[0082]
[0083]
[0084] 表2实例二中输入模型自变量
[0085] B2、对步骤B1中获取的众多驱动因子做多重共线性检验，并与步骤A得到的植被覆盖度FVC做多元逐步线性回归，筛选得到除了煤炭开采以外与植被覆盖度相关性最高的驱
动因子。
[0086] B3、将步骤B2得到的驱动因子作为自变量输入时空地理加权回归模型，训练得到植被覆盖度FVC与其各类驱动因子的关系模型：
[0087]
[0088] 式中，(xi1,xi2,…,xid；yi)表示第i个观测点(ui,vi,ti),(i＝1,2,…,n)处的植被参数y和驱动因子x1,x2,…,xd的n组观测值；βk(ui,vi,ti),(k＝0,1,…,d)是第i个数据点
(ui,vi,ti)处的未知回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定均值为0，方
2
差为σ；
[0089] 根据加权最小二乘方法，第i个观测点的回归系数估计值 为：
[0090]
[0091] 第i个观测点因变量的拟合值 为：
[0092]
[0093] 式中，X为驱动因子矩阵 Xi是X矩阵的第i行；y为植被参数矩阵 Wi表示空间核函数矩阵 常用的核函
数包括高斯(Gauss)核函数和近高斯(Bisquare)核函数。
[0094] 调整型高斯核函数下，观测点j对观测点i影响的权重值：
[0095]
[0096] 调整型近高斯核函数下，观测点j对观测点i影响的权重值：
[0097]
[0098] 式中:dij表示点(uj,vj,tj)到点(ui,vi,ti)的时空距离；dN表示点(ui,vi,ti)到最近邻第N个点的时空距离；用时空椭球坐标系来计算时空距离：
其中，λ和μ为时间和空间距离的比例调整系数。
[0099] C、剥离采煤对植被的干扰，根据步骤B训练得到的定量关系，输入矿区有采矿活动背景下的气候气象因子、地理因子和人类活动因子数据，可以预测出没有煤矿开采情况下
矿区逐年的植被覆盖度FVC。预测的植被覆盖度FVC与遥感监测的植被覆盖度实际FVC可认
为是剔除其他影响后，由煤炭开采导致的植被覆盖度变化量，将其定义为采矿植被扰动量
V‑MD(Vegetation–Mining Disturbance)。
[0100] C1、采矿植被扰动量V‑MD的计算方法为：V‑MD＝V实际‑V预测，其中，V实际为矿区有煤炭开采活动时期通过遥感反演得到的实际植被参数，V预测为模型预测得到的植被参数预测值。
按照上述方法，可以得到叶面积指数LAI、植被叶绿素含量Cab分别的采矿植被扰动量。
[0101] 以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。
