基于空间自相关区域植被净初级生产力时空分异探测方法与流程

本发明涉及植被生产力探测领域，特别是涉及一种空间自相关区域植被净初级生产力时空分异探测，尤其涉及一种利用卫星遥感数据快速、高效、低廉、准确探测大区域尺度植被净初级生产力时空分异特征的方法。

背景技术：

气候和人类活动对植被生产力水平区域分布有很大影响。准确探测植被净初级生产力空间分布对评估陆地生态系统的结构和功能具有重要意义，净初级生产力(netprimaryproductivity)简称为npp，传统的npp时空分布情况监测主要依赖于单一像元属性值的计算和比较，忽略了各部分之间的空间相关性，研究结果表示的是多数点的变化情况，这导致在不同的研究尺度下对同一地区的变化描述存在差异，甚至相反。在当前全球气候变化频繁与人类活动加剧的背景下，急需发明一种快速、高效、低廉、准确的标准方法来探测和表达区域植被净初级生产力时空分异特征。

技术实现要素：

本发明的目的是提供一种空间自相关区域植被净初级生产力时空分异探测方法，以解决所述现有技术存在的问题。

为实现所述目的，本发明提供了如下方案：本发明提供一种空间自相关区域植被净初级生产力时空分异探测方法，包括以下步骤：

s1、获取研究区域研究年份的mod17a3数据

下载研究区域研究时段对应的modis陆地产品mod17a3数据，然后对mod17a3数据进行格式转换、地图投影转换，即hdf格式转换为tiff格式，正弦曲线地图投影转换为wgs84/albers地图投影；

s2、提取研究区域研究时段的植被npp图

将步骤s1中的mod17a3处理后的数据，用于研究区边界的矢量图对处理后的数据进行掩膜提取，空间分辨率1km×1km，数据格式为tiff格式，投影为wgs84/albersequalareaconic；

s3、计算研究区域研究时段的逐年植被npp变化量

依据步骤s2中提取的研究区域研究时段植被npp图，利用后一年的植被npp减去前一年的植被npp，生成研究区域研究时段的逐年植被npp差值图；

s4、计算研究区域研究时段植被npp的全局莫兰指数i和高/低值聚类指数g(d)及相关统计量

根据步骤s2中提取的研究时段的植被npp，计算全局莫兰指数i，公式如下：

其中，xi、xj分别为点i和点j的观测值；是相关观测值的均值；wij代表空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；n为观测值的数量。再计算每个全局莫兰指数值对应的期望e、p值和z得分；

计算高/低值聚类指数g(d)，公式如下：

其中，wij(d)为空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；xi、xj分别为点i和点j的观测值。再计算每个高/低值聚类指数对应的期望e、p值和z得分；

s5、判断研究区域数据样本的整体聚类特征

根据步骤s4计算的研究区域研究时段的全局莫兰指数i，高/低值聚类指数g(d)及相关统计量，检验数据能否拒绝零假设，并判断是否存在聚类以及聚类类型：若全局莫兰指数i为正则表示数值相近的npp之间呈聚集趋势，若全局莫兰指数i为负则表示呈离散趋势；若高/低值聚类指数g(d)大于期望值则表示研究区域的聚类类型为高值聚类，若高/低值聚类指数g(d)小于期望值则表示研究区域的聚类类型为低值聚类；

s6、计算研究区域研究时段植被npp的热点分析指数g*(d)

根据步骤s2中提取的研究区域研究时段的植被npp，计算研究区域研究时段的热点分析指数g*(d)，公式如下：

其中，wij(d)为空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；xj为点j处的npp观测值。再计算热点分析指数g*(d)对应的期望e、p值和z得分；

s7、判断研究区域研究时段植被npp的热点分布特征

根据步骤s6得到的研究区域研究时段植被npp热点分析指数g*(d)估计值，判断研究区域研究时段植被npp的热点分布特征：如果z(gi*)为正且在显著区间，表明点i周围的值属高值空间集聚(热点区)；反之，如果z(gi*)为负且在显著区间，则表明点i周围的值属低值空间集聚(冷点区)；

s8、计算研究区域研究时段逐年植被npp变化量的局部莫兰指数ii

根据步骤s3中计算的研究区域研究时段的逐年植被npp变化量，计算局部莫兰指数ii，公式如下：

其中，xi为区域i的观测值；为相关观测值的均值；wij为空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；n为观测值的数量；

s9、判断研究区域研究时段逐年植被npp变化量的空间聚类特征

根据步骤s8中得到的研究区域研究时段逐年植被npp变化量局部莫兰指数ii计算结果，判断研究区域研究时段逐年植被npp变化量的空间聚类特征，包括hh—高值围绕高值、hl—高值围绕低值、lh—低值围绕高值、ll—低值围绕低值四种类型，其中hh和ll这两种为正相关模式，hl和lh为负相关模式。

本发明公开了以下技术效果：

本发明基于空间自相关指数对空间信息的表达能力，构建了探测植被npp时空分异的方法，弥补了传统方法用多数点的趋势概括总体来描述植被净初级生产力变化特征，从而导致在不同的研究尺度下对同一地区的变化描述存在差异的不足。本发明能够准确表达植被npp的整体聚类特征，并探测出局部空间特征，以便将大尺度区域的研究讨论的对象转换为具有空间相关性的几部分，实现了基于卫星遥感数据的大区域范围植被净初级生产力时空分异和演变特征的快速、高效、低廉、准确探测与表达。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明方法的流程示意图；

图2为2001年内蒙古呼伦贝尔市与乌兰察布市植被npp热点分析中冷点和热点占比示意图

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

为使本发明的所述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

参照图1-2本发明提供一种基于空间自相关区域植被净初级生产力时空分异探测方法，包括如下步骤：

s1、获取内蒙古2000-2014年的mod17a3数据

下载内蒙古2000-2014年modis陆地产品mod17a3数据，使用mrt软件将mod17a3数据进行格式转换、地图投影转换，即hdf格式转换为tiff格式，正弦曲线地图投影转换为wgs84/albers地图投影；

s2、提取内蒙古2000-2014年植被npp图

依据所述步骤s1中获取的mod17a3陆地产品数据，用内蒙古边界矢量图对处理后的mod17a3数据进行掩膜提取，空间分辨率1km×1km，数据格式为tiff格式，投影为wgs84/albersequalareaconic；

s3、计算内蒙古地区逐年植被npp变化量

依据所述步骤s2提取的2000-2014年的植被npp图，利用后一年的植被npp减去前一年的植被npp，生成内蒙古2000-2014年的逐年植被差值图；

s4、计算内蒙古2000-2014年植被净初级生产力(npp)的全局莫兰指数i和高/低值聚类指数g(d)及相关统计量

根据步骤s2中提取的内蒙古2000-2014年植被npp，计算全局莫兰指数i，公式如下：

其中，xi、xj分别为点i和点j的观测值；是相关观测值的均值；wij代表空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；n为观测值的数量。再计算每个全局莫兰指数对应的期望e、p值和z得分；

计算高/低值聚类指数g(d)，公式如下：

其中，wij(d)为空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；xi、xj分别为点i和点j的观测值。再计算每个高/低值聚类指数对应的期望e、p值和z得分；

s5、判断内蒙古植被净初级生产力的整体聚类特征

根据步骤s4计算的内蒙古2000-2014年的全局莫兰指数i，高/低值聚类指数g(d)及相关统计量，检验数据能否拒绝零假设，并判断是否存在聚类以及聚类类型：若全局莫兰指数i为正则表示数值相近的npp之间呈聚集趋势，若全局莫兰指数i为负则表示呈离散趋势；若高/低值聚类指数g(d)大于期望值则表示内蒙古地区的聚类类型为高值聚类，若高/低值聚类指数g(d)小于期望值则表示内蒙古地区的聚类类型为低值聚类；

s6、计算内蒙古2000-2014年植被npp的热点分析指数g*(d)

根据步骤s2中提取的内蒙古2000-2014年植被npp，计算热点分析指数g*(d)，公式如下：

其中，xj为点j处的npp观测值；wij(d)为空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0。再计算热点分析指数g*(d)对应的期望e、p值和z得分；

s7、判断内蒙古2000-2014年植被npp的热点分布特征

依据所述步骤s6计算的内蒙古2000-2014年植被npp的热点分析指数g*(d)，判断内蒙古植被npp的热点分布特征：如果z(gi*)为正且在显著区间，表明点i周围的值属高值空间集聚(热点区)；反之，如果z(gi*)为负且在显著区间，则表明点i周围的值属低值空间集聚(冷点区)；

s8、计算内蒙古2000-2014年逐年植被npp变化量的局部莫兰指数ii

根据步骤s3中计算的内蒙古2000-2014年逐年植被npp变化量，计算局部莫兰指数ii，公式如下：

其中，xi为区域i的观测值；为相关观测值的均值；wij为空间权重矩阵，若点i和点j空间相邻，矩阵第i行第j列对应的元素为1，不相邻则为0；n为观测值的数量；

s9、判断内蒙古2000-2014年的逐年植被npp变化空间聚类特征

依据所述步骤s8得到的npp内蒙古2000-2014年逐年植被npp变化量局部莫兰指数ii计算结果，判断内蒙古逐年植被净初级生产力npp变化量的空间聚类特征。

在本发明的描述中，需要理解的是，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。