CN115344815A - 考虑植被空间变化的天然径流变化归因分析方法和系统 - Google Patents

Abstract

本发明公开了一种考虑植被空间变化的天然径流变化归因分析方法和系统，通过构建模型从网格尺度区分归一化植被指数NDVI受气候变化和人类活动影响下的地理空间变化，在此基础上，分解下垫面参数对气候变化和人类活动的影响，分别计算气候变化和人类活动引起的归一化植被指数的贡献率，从而可以更加准确地计算出气候变化对研究区域的天然径流变化的贡献率、人类活动对研究区域的天然径流变化的贡献率。对天然径流变化进行准确归因。可更科学地得出在任意空间尺度上径流受气候变化和人类活动的贡献率。本发明提高了天然径流变化归因的准确性，同时完善了气候变化和人类活动对流域河川径流变化影响的定量归因方法体系。

Description

考虑植被空间变化的天然径流变化归因分析方法和系统

技术领域

本发明涉及流域径流变化归因识别技术，尤其是一种考虑植被空间变化的天然径流变化归因分析方法和系统。

背景技术

近年来，变化环境下水文循环与水资源演化机制的研究逐渐成为热点，受气候变化和人类活动的影响，不同流域的径流也呈现复杂多变的水文情势，定量归因气候变化和人类活动对径流变化影响有十分重要的意义。

Budyko假设方法是统计分解变化的有效方法，在气候变化和人类活动引起的河流径流中得到了广泛应用。在径流变化归因研究中，径流变化通常归因于气象要素和下垫面参数的变化，其中气象要素主要包括降水、蒸散发等；下垫面参数主要由土壤质地、地形因素和土地利用所决定，其中土壤质地和地形的变化可忽略不计，而土地利用为影响下垫面参数最主要的因素。传统基于Budyko假设的径流归因研究认为气候变化对径流变化的影响主要来源于降水、蒸散发等气象要素，人类活动对径流变化的影响主要由土地利用变化造成的，即土地利用变化对下垫面参数改变引起的。

但经申请人研究发现，土地利用的变化既受气候变化的影响，也受人类活动的影响，如地表植被覆盖作为众多流域土地利用的表现形式，植被覆盖变化同时受诸如降水、蒸散发等气候要素以及植树造林、水土保持等人类活动的影响。因此，应认识到下垫面参数的变化是气候变化和人类活动综合影响下的结果，而现状研究中将下垫面参数变化完全归因为人类活动的认知是错误的。

上述研究导致归因结果中由气候变化引起的径流量变化量未包括气候变化通过影响下垫面参数引起的径流变化量，而人类活动引起的径流量变化过多地考虑了下垫面参数变化的影响，这可能引起传统基于Budyko假设的径流变化归因结果的失准，即气候变化引起的径流变化量比真实数据偏小，而人类活动引起的径流量变化量偏大。另外，下垫面参数在特定时空尺度下是一个定值，当前还没有将下垫面参数变化量做进一步归因的研究，评估下垫面参数变化在特定时空尺度下受气候变化和人类活动的影响也是一个难点，这也对基于Budyko假设的径流变化归因方法的改进带来挑战。

发明内容

发明目的：一方面，提供一种考虑植被空间变化的天然径流变化归因分析方法，以解决现有技术存在的上述问题，例如通过进一步归因下垫面参数的变化量，对气候变化和人类活动对天然径流变化的贡献率公式进行修正，改进了原有Budyko假设在径流归因研究中的不足，以提高径流变化归因结果的准确性。另一方面，提供一种基于上述方法的分析系统。

技术方案：根据本申请的一个方面，考虑植被空间变化的天然径流变化归因分析方法，包括如下步骤：

步骤S1、采集研究区域预定时间范围内的长时间序列数据，并将预定时间范围至少分为连续的第一时期和第二时期；分别计算第一时期的归一化植被指数平均值NDVI₁和第二时期的归一化植被指数平均值NDVI₂；并计算其差值∆NDVI_a，其中∆NDVI_a=NDVI₂-NDVI₁；

步骤S2、构建归一化植被指数与气象要素的地理回归关系，在约定下垫面参数恒定的情况下，分别计算第一时期归一化植被指数平均值NDVI_1s和第二时期的归一化植被指数平均值NDVI_2s；并计算其差值∆NDVI_C，即为气候变化引起的归一化植被指数平均值变化量；其中，∆NDVI_C =NDVI_2s-NDVI_1s；接着计算人类活动引起的归一化植被指数平均值变化量∆NDVI_H；∆NDVI_H =∆NDVI_a-∆NDVI_C；然后分别计算气候变化和人类活动引起的归一化植被指数的贡献率η_C和η_H；η_C=∆NDVI_C/ NDVI_a；η_H=∆NDVI_H/ NDVI_a；

步骤S3、基于天然径流对变化量的弹性系数构建天然径流的全微分方程；并通过全微分方程计算第二时期与第一时期的天然径流量平均值之差dR，以及降雨量变化、潜在蒸散发量变化、下垫面参数变化导致天然径流量的改变量∆R_p=ε_p（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆R_L=ε_L（R/L）×∆L；

步骤S4、基于上述计算结果计算气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH：

η_Rc=(ΔR_P+ΔR_ET0)/ΔR×100%+η_C×ΔR_L/ΔR×100%；

η_RH=ΔR_L/ΔR×100%-η_C×ΔR_L/ΔR×100%=η_H×ΔR_L/ΔR×100%；

其中，ΔR=ΔR_P+ΔR_ET0+η_C×ΔR_L+η_H×ΔR_L；

η_C×ΔR_L/ΔR×100%为修正项，即气候变化通过影响下垫面参数引起径流变化的贡献率；

η_H×ΔR_L/ΔR×100%为人类活动通过影响下垫面参数引起的径流变化的贡献率。

根据本申请的一个方面，所述步骤S1进一步为：

步骤S11、划定研究区域，并获取研究区域的地理数据；

步骤S12、采集研究区域预定时间范围内的长时间序列数据，所述长时间序列数据为中分辨率成像光谱仪数据。

根据本申请的一个方面，所述步骤S1中将预定时间范围至少分为连续的第一时期和第二时期的过程进一步为：

步骤S1a、判断研究区域是否存在影响下垫面参数的事件，若存在，将事件发生的时间作为第一时期和第二时期的分界点；若不存在，进入步骤S1b、

步骤S1b、构建预定长度的时间区间，并在预定时间范围内滑动，形成至少两个间隔大于阈值的时期，记为第一检测时期和第二检测时期；在预定时间范围内选择至少一个时期作为标准时期；分别计算第一检测时期、第二检测时期和标准时期的归一化植被指数平均值；

步骤S1c、分别计算第二检测时期与标准时期的归一化植被指数平均值的差值，记为第一变化值，第一检测时期与标准归一化植被指数平均值的差值，记为第二变化值，

计算第一变化值与第二变化值差值，并与阈值进行比较，若大于阈值，将第一时期和第二时期的起点时间作为下垫面参数改变的待检时刻，放入待检时刻集合。

根据本申请的一个方面，所述步骤S2进一步为：

步骤S21、构建地理加权回归模型，采集研究区域内第一时期的降雨量、平均气温和空间位置参数，采用局部加权最小二乘法求解所述地理加权回归模型，获得第一时期的归一化植被指数并计算平均值NDVI_1s；

步骤S22、采集第二时期的降雨量、平均气温和空间位置参数，设定归一化植被指数与气象要素的地理回归关系在第一时期和第二时期相同；构建地理加权回归模型，并采用局部加权最小二乘法求解所述地理加权回归模型，获得第二时期的归一化植被指数并计算平均值NDVI_2s；

步骤S23、计算第二时期和第一时期的归一化植被指数并计算平均值的差值∆NDVI_C，∆NDVI_C =NDVI_2s-NDVI_1s；即为即为气候变化引起的归一化植被指数平均值变化量；

步骤S24、计算人类活动引起的归一化植被指数平均值变化量∆NDVI_H；∆NDVI_H =∆NDVI_a-∆NDVI_C；

步骤S25、分别计算气候变化和人类活动引起的归一化植被指数的贡献率η_C和η_H；η_C=∆NDVI_C/ NDVI_a；η_H=∆NDVI_H/ NDVI_a。

根据本申请的一个方面，所述步骤S3进一步为：

步骤S31、构建研究区域的水平衡量方程R=P-ET-∆S；

其中，R为径流量，P为降水量，ET为实际蒸散发量，∆S为自然流域下的陆地水储量变化量；

步骤S32、构建干旱指数的函数，F（φ）=ET/P，其中φ为干旱指数，

其中，

，n为下垫面参数；

步骤S33、根据天然径流对变量的弹性系数的定义，构建天然径流的微分方程；

；

ε_P、ε_ET0、ε_L分别为天然径流对降水、潜在蒸散发和下垫面参数的弹性系数；若弹性系数为正值，表示径流随着该变量的增加而增加，若弹性系数为负值，表示径流随着该变量的增加而减小；

步骤S33、用天然径流第二时期和第一时期的多年平均径流量之差dR表示第一时期到第二时期年径流的变化dR=R₂-R₁，R₂和R₁分别为第二时期和第一时期的年径流量；

步骤S34、计算降雨量变化、潜在蒸散发量变化、下垫面参数变化导致天然径流量的改变量∆R_p=ε_p（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆R_L=ε_L（R/L）×∆L。

根据本申请的一个方面，还包括步骤S5：

步骤S51、获取研究区地图及各处的气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH；

步骤S52、分别采用气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH渲染研究区域地图；

步骤S53、查找研究区域地图上的分界线，并将气候变化对研究区域的天然径流变化的贡献率η_Rc超出阈值、人类活动对研究区域的天然径流变化的贡献率η_RH超出阈值的区域分别采用颜色标注在不同的图层上。

根据本申请的一个方面，所述步骤S53进一步为：

步骤S53a、读取渲染后的研究区域地图，针对每一颜色通道，计算相邻两个像素之间的像素差，并计算该像素差与所述相邻两个像素之间最大像素值之间的比例，即为相差比i；

步骤S53b、针对每一颜色通道，计算相差比大于阈值的相邻像素；

步骤S53c、连接相邻像素，并判断相邻像素是否构成连续的分界线；若构成，查找气候变化对研究区域的天然径流变化的贡献率η_Rc超出阈值的一侧或人类活动对研究区域的天然径流变化的贡献率η_RH超出阈值的一侧，并采用颜色分别标注在不同的图层。

根据本申请的一个方面，所述步骤53c还包括：

若不构成，依次判断相邻像素组成线段的相邻端点之间的欧式距离是否小于阈值，若小于，则连接两个端点，并作为分界线。

根据本申请的一个方面，还包括步骤S54：

判断分界线是否构成封闭区域，或分界线与研究区域边缘线是否构成封闭区域；若构成，采用颜色标注。

在本申请的另一实施例中，还提供一种考虑植被空间变化的天然径流变化归因分析系统，包括：

至少一个处理器；以及

与至少一个所述处理器通信连接的存储器；其中，

所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述任一项实施例所述的考虑植被空间变化的天然径流变化归因分析方法。

有益效果：通过构建地理加权回归模型（GWR）从网格尺度区分NDVI受气候变化和人类活动影响下的地理空间变化，以归因下垫面参数对气候变化和人类活动的影响，在此基础上，对气候变化和人类活动对天然径流变化的贡献率公式进行修正，改进了传统Budyko假设在径流归因研究中的不足，同时可更准确地得出在任意空间尺度上径流受气候变化和人类活动的贡献率。

附图说明

图1是本发明的整体实现过程示意图。

图2是本发明步骤S1的流程图。

图3是本发明将预定时间范围至少分为连续的第一时期和第二时期的过程示意图。

图4是本发明步骤S2的流程图。

图5是本发明步骤S3的流程图。

图6是本发明步骤S5的流程图。

图7是本发明步骤S53的流程图。

具体实施方式

结合图1描述本发明的主要技术过程和技术细节。为了解决现有技术存在的问题，申请人对现有技术进行了归集、分析和研究。在现有的研究中，多数文献将下垫面参数的变化归因于人类活动。在一些文献中，仅仅通过降水量弹性系数、潜在蒸发量弹性系数和水热耦合参数的弹性系数进行归因分析，这种分析不能够精确反应物理过程，准确性不够高，在这些分析中，也并未考虑采用归一化植被指数（NDVI）来分解影响因素。在计算时，主要通过一段时间内的平均值来表示，所以为了表述方便，NDVI用来表示归一化植被指数平均值。

根据本申请的一个方面，考虑植被空间变化的天然径流变化归因分析方法，包括如下步骤：

步骤S1、采集研究区域预定时间范围内的长时间序列数据，并将预定时间范围至少分为连续的第一时期和第二时期；分别计算第一时期的归一化植被指数平均值NDVI₁和第二时期的归一化植被指数平均值NDVI₂；并计算其差值∆NDVI_a，其中∆NDVI_a=NDVI₂-NDVI₁。

在这一步骤中，首先根据采集的数据，判断是否存在人工改造下垫面的工程项目，例如水土保持和植树造林等工程项目，如果存在的话，分别确定上述工程项目的时间范围，尤其是工程项目的时间起点。例如在某个研究区域，可以把植被变化分成两个阶段来考虑，在第一时期，称之为自然生长期，在第二时期，称之为生态改造期。在两个时期的分界点，即为植被变化的转折时刻，自该时刻开始，研究区域开始了大规模植被工程。植被同时受到自然生长和人类社会改造的影响。因此，在研究区区域，前后两个阶段平均归一化植被指数平均值NDVI可以被计算出来，同时，可以得到两者的差值。两者的差值，即为研究区域不同阶段中气候变化导致的差值与人类活动变化导致的差值之和。气候变化，主要受降水、温度等要素的影响，人类活动变化主要是受到大规模的人类活动，例如水土保持、植树造林等植被工程。换句话说，归一化植被指数平均值NDVI的变化主要有两个主要原因，一个是研究区域的气候变化，比如研究区域，在不同的时间段上，气候发生了较大的变化（例如区域气候异常等等），导致NDVI数值发生了变化。另一个是人类活动对研究区域的影响，比如进行了植树造林等等工程。

步骤S2、构建归一化植被指数与气象要素的地理回归关系，在约定下垫面参数恒定的情况下，分别计算第一时期归一化植被指数平均值NDVI_1s和第二时期的归一化植被指数平均值NDVI_2s；并计算其差值∆NDVI_C，即为气候变化引起的归一化植被指数平均值变化量；其中，∆NDVI_C =NDVI_2s-NDVI_1s；接着计算人类活动引起的归一化植被指数平均值变化量∆NDVI_H；∆NDVI_H =∆NDVI_a-∆NDVI_C；然后分别计算气候变化和人类活动引起的归一化植被指数的贡献率η_C和η_H；η_C=∆NDVI_C/NDVI_a；η_H=∆NDVI_H/ NDVI_a。

在这一步骤中，由于现有资料表明，归一化植被指数的平均值NDVI与降水、平均气温等气象因素有着密切的关系，可以采用地理回归模型，例如表征地理要素非平稳空间相关性的地理加权回归模型来描述归一化植被指数与降水和平均气温之间的地理回归关系。在这一模型中，已不存在人类活动的影响，即不考虑水土保持和退耕还林等人类活动影响，NDVI与降水、平均气温等气象要素的地理回归关系在某一时间段内保持不变，从而在第一时期和第二时期内，NDVI与不同气象要素的回归关系可以采用同样的回归模型来计算。

因此，通过上述回归模型，基于采集的研究区域的相关数据，可以获得各个时期内归一化植被指数，并求得平均值。基于上述假设和模型，可以计算气候变化造成的归一化植被指数平均值的变化量。通过步骤S1获得的总变化量，就可以算出人类活动引起的归一化植被指数变化量。通过上述变化量，可以算出各自的贡献率。也就是气候变化和人类活动对植被变化影响的贡献率。

根据气候变化和人类活动对植被变化影响的贡献率，可以对研究区域的NDVI归因分析。比如在某一区域，每升高100mm将导致NDVI增加值超过0.1。在一些区域气温的增加，导致NDVI数值的降低，可以是由于气温增加会促进光和作用，同时促进NDVI的增加。在另外一些区域，气温的增加导致NDVI的降低，可能是由于气温已经超过作物的最适宜的温度，气温的增加导致营养物质消耗过快，从而导致NDVI降低。

步骤S3、基于天然径流对不同要素变化量的弹性系数构建天然径流的全微分方程；并通过全微分方程计算第二时期与第一时期的天然径流量平均值之差dR，以及降雨量变化、潜在蒸散发量变化、下垫面参数变化导致天然径流量的改变量∆R_p=ε_p（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆R_L=ε_L（R/L）×∆L。

弹性系数定义为各变量变化导致流域径流量的变化程度，即降水量、潜在蒸散发量或下垫面参数等变量增加1%导致流域年径流量相对于多年平均值的变化百分比，径流对变量x的弹性系数一般表示为ε_x=（dR/dx）×（x/R）。因此，径流对降水量、实际蒸散发量量和下垫面参数的变化量可以用如下全微分表示：

dR/R=ε_P（dP/P）+ε_ET0（dET0/ET0）+ε_L（dL/L）；利用水热耦合平衡方程，可以得到各个变量的弹性系数，对于一个特定的研究区域，在一段时间内，可以认为各个变量的弹性系数为常数。具体不再详述。若弹性系数为正值，表示径流随着该变量的增加而增加，若弹性系数为负值，表示径流随着该变量的增加而减小。

基于上述数据，从自然生长期到生态改造期的年径流的变化可以用径流前后两时段的多年平均径流量之差dR。换句话说，可以求得第一时期与第二时期内的多年平均径流量之差。基于类似的思路，可以求得多年平均降水量、潜在蒸散发量和下垫面参数的变化量。

根据水量平衡方程R=P-E=P-f（P，ET），实际蒸散发是关于降水量、潜在蒸散发量的函数，径流变异点前后由于下垫面参数的不同导致函数表达式的不同，突变点前后径流差可表示为∆R_p=ε_p（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆R_L=ε_L（R/L）×∆L。

综上，即计算出各变量导致的径流改变量。从而可以计算出气候变化和人类活动对研究区域天然径流变化的贡献率。

步骤S4、基于上述计算结果计算气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH。

η_Rc=(ΔR_P+ΔR_ET0)/ΔR×100%+η_C×ΔR_L/ΔR×100；

η_RH=ΔR_L/ΔR×100%-η_C×ΔR_L/ΔR×100%=η_H×ΔR_L/ΔR×100%；

其中，ΔR=ΔR_P+ΔR_ET0+η_C×ΔR_L+η_H×ΔR_L；

η_C×ΔR_L/ΔR×100%为修正项，即气候变化通过影响下垫面参数引起径流变化的贡献率；

η_H×ΔR_L/ΔR×100%为人类活动通过影响下垫面参数引起的径流变化的贡献率。

在这一步骤中，基于上述步骤的计算结果，可以获得最终气候变化对研究区域的天然径流变化的贡献率η_Rc和人类活动对研究区域的天然径流变化的贡献率η_RH。

如图2所示，根据本申请的一个方面，所述步骤S1进一步为：

步骤S11、划定研究区域，并获取研究区域的地理数据。

步骤S12、采集研究区域预定时间范围内的长时间序列数据，所述长时间序列数据为中分辨率成像光谱仪数据。

在该实施例中，首先根据需求划定研究区域，并获取研究区域的地理数据，并采集该研究区域内的基础数据，比如可以采用中分辨率成像光谱仪数据。当然，可以根据实际情况采集其他数据集。采用不同的数据集，其结果可能存在偏差，但趋势应该是相同的。在其他实施例中，可以采用不同的数据集进行计算，并进行比较，给出最能够准确表征径流变化的数据集。

在研究中发现，有些研究区域的下垫面参数变化时刻是清楚的，比如某些研究区域，从某一年份开始进行植树造林计划，但是也有些区域的变化时刻是不清楚的，例如有些区域，由于降水等因素导致树木生长速度加快，森林覆盖率增加，下垫面的参数发生了明显的变化。换句话说，在不清楚什么时刻发生变化时，需要有一种可信的方法来判断研究区域是否有需要分成几个时期。在这种情况下，提供如下解决方案。

如图3所示，根据本申请的一个方面，所述步骤S1中将预定时间范围至少分为连续的第一时期和第二时期的过程进一步为：

步骤S1a、判断研究区域是否存在影响下垫面参数的事件，若存在，将事件发生的时间作为第一时期和第二时期的分界点；若不存在，进入步骤S1b。

步骤S1b、构建预定长度的时间区间，并在预定时间范围内滑动，形成至少两个间隔大于阈值的时期，记为第一检测时期和第二检测时期；在预定时间范围内选择至少一个时期作为标准时期；分别计算第一检测时期、第二检测时期和标准时期的归一化植被指数平均值。

步骤S1c、分别计算第二检测时期与标准时期的归一化植被指数平均值的差值，记为第一变化值，第一检测时期与标准归一化植被指数平均值的差值，记为第二变化值。

计算第一变化值与第二变化值差值，并与阈值进行比较，若大于阈值，将第一时期和第二时期的起点时间作为下垫面参数改变的待检时刻，放入待检时刻集合。

在这一方案中，通过建立一个滑动时间窗口的方法，从时间区间上选取两个间隔的时期，分别计算两个时期内的归一化植被指数的平均值，然后根据平均值与标准时期的平均值的变化量的差值是否大于阈值，来判断是否发生了变化，如果发生变化，则采用同样的方法，查找上述变化量的发生时刻。因此，通过上述方法，可以针对一个不确定是否存在分界点的时间区域进行检测和计算，解决了研究区域中，人工无法确定下垫面参数变化时刻的问题。

比如，当前具有某个研究区域内1972年至2020年的数据，但是并不确定该区域的植被变化存在几个阶段，每个阶段的时间分界点在哪一年。即使知道某年开始进行了植树造林运动，也没法确定在这之前下垫面已经发生了变化。在这种情况下，如何将研究阶段进行划分，是首先需要解决的问题。因此可以将1972~1976年作为标准时期，然后以10年的时间窗口在1972~2020上滑动。比如滑动后形成1982~1991和2002~2011两个时间窗口。就可以分别计算两个时期内的NDVI平均值与标准时期内的NDVI平均值的差值，作为第一变化值和第二变化值，如果第一变化值和第二变化值之间的差值大于阈值，则说明第二时期内的下垫面参数发生了变化，且变化时刻在第一时期的最终时刻至第二时期的最终时刻之间。将时间窗口向后滑动，比如滑动至1995~2004，再次判断是否发生变化，以及变化量是否超过阈值，则可以判断下垫面参数发生变化的时刻。

随后调整前一个时间窗口，向后滑动，比如滑动至1977~1986，判断是否发生变化以及发生变化的量是否超过阈值。在此基础上，判断是否还有下垫面发生变化的分界点。如果有，则根据分界点的数量分成若干个时期，分别为第一时期、第二时期、第三时期……。然后采用实施例1所述的方法来计算各个时期的最终气候变化对研究区域的天然径流变化的贡献率η_Rc和人类活动对研究区域的天然径流变化的贡献率η_RH。计算方法与实施例1中的方法相似，在此不再详述。

在进一步的实施例中，还可以采用如下方法先大致判断下垫面参数的时间变化范围。具体如下：

步骤S100、针对研究区域的基本情况，在每一年度中选取至少N个典型时间点，其中N≥8；

步骤S101、读取研究区域内每一年度各个典型时间点的基础数据图数据，并基于每一典型时间点构建跨年度的基础数据集，并按照时间顺序划分为若干子集；

步骤S102、计算每一基础数据集中各年度的基础数据图的子集数据图的相似度；针对对应每一典型时间点的基础数据集，构建权重矩阵；

步骤S103、计算各个子集内每一年度各个典型时间点总的相似度，并判断是否超出阈值；如果超出阈值，则标记为发生变化。

其中，所述基础数据图数据包括SPOT NDVI、MODIS NDVI和GIMMS NDVI。

在一个具体的实施例中，比如某一研究区域有1982~2016年的数据，针对每一年，按照季度（即四个典型时间点，即第一时间点至第四时间点）提取至少2张或3张图片，以反应当季的植被情况。还可以从当季的若干张图中选择典型图片或者对当季的图片进行加权取值，从而获得典型图片。基于上述方法，选取每一年度各个典型时间点的基础数据图，比如1982~2016年中，第一典型时间点的数据图，构建每一时间点的基础数据集；

按照时间顺序，从基础数据集中找出至少2个子集，比如第一典型时间点的1982~1996和2002~2016的子集。其余时间点的数据集处理方式与此相同。

然后计算两个子集中图片的整体相似度。比如第一典型时间点第一个子集和第二个子集的图像的相似度。如此反复，计算出所有时间点的图像的相似度。

然后基于构建的不同时间点的权重矩阵，计算出总的相似度。比如，第一典型时间点至第四典型时间点分别给予权重参数，即可得到整体的相似度。

如果整体的相似度超出阈值，比如低于某个数值，则说明植被可能发生变化。需要进行后续的计算过程。

如图4所示，根据本申请的一个方面，所述步骤S2进一步为：

步骤S21、构建地理加权回归模型，采集研究区域内第一时期的降雨量、平均气温和空间位置参数，采用局部加权最小二乘法求解所述地理加权回归模型，获得第一时期的归一化植被指数并计算平均值NDVI_1s；

步骤S22、采集第二时期的降雨量、平均气温和空间位置参数，设定归一化植被指数与气象要素的地理回归关系在第一时期和第二时期相同；构建地理加权回归模型，并采用局部加权最小二乘法求解所述地理加权回归模型，获得第二时期的归一化植被指数并计算平均值NDVI_2s；

步骤S23、计算第二时期和第一时期的归一化植被指数并计算平均值的差值∆NDVI_C，∆NDVI_C =NDVI_2s-NDVI_1s；即为即为气候变化引起的归一化植被指数平均值变化量；

步骤S24、计算人类活动引起的归一化植被指数平均值变化量∆NDVI_H；∆NDVI_H =∆NDVI_a-∆NDVI_C；

步骤S25、分别计算气候变化和人类活动引起的归一化植被指数的贡献率η_C和η_H；η_C=∆NDVI_C/ NDVI_a；η_H=∆NDVI_H/ NDVI_a。

需要说明的是，本实施例的主要思路在上述实施例中，已经进行了描述。本实施例对其不同点进行描述。在本实施例中，在步骤S21和步骤S22中，通过地理加权回归模型对研究区域进行建模，然后再进行计算。在求解的过程中，需要对回归系数进行计算，此时，可以通过最小二乘法进行求解。

如图5所示，根据本申请的一个方面，所述步骤S3进一步为：

步骤S31、构建研究区域的水平衡量方程R=P-ET∆S；

其中，R为径流量，P为降水量，ET为实际蒸散发量，∆S为自然流域下的陆地水储量变化量；

步骤S32、构建干旱指数的函数，F（φ）=ET/P，其中φ为干旱指数，

其中，

，n为下垫面参数；

步骤S33、根据天然径流对变量的弹性系数的定义，构建天然径流的微分方程；

；

ε_P、ε_ET0、ε_L分别为天然径流对降水、潜在蒸散发和下垫面的弹性系数；若弹性系数为正值，表示径流随着该变量的增加而增加，若弹性系数为负值，表示径流随着该变量的增加而减小；

步骤S33、用天然径流第二时期和第一时期的多年平均径流量之差dR表示第一时期到第二时期年径流的变化dR=R2-R1，R2和R1分别为第二时期和第一时期的年径流量；

步骤S34、计算降雨量变化、潜在蒸散发量变化、下垫面参数变化导致天然径流量的改变量∆R_p=ε_p（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆RL=ε_L（R/L）×∆L。

这一实施例的内容，在实施例1中进行了表述，需要说明的是，干旱指数函数的计算中，是与上述实施例不同的。

如图6所示，根据本申请的一个方面，还包括步骤S5：

步骤S51、获取研究区域地图及各处的气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH；

步骤S52、分别采用气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH渲染研究区域地图；

步骤S53、查找研究区域地图上的分界线，并将气候变化对研究区域的天然径流变化的贡献率η_Rc超出阈值、人类活动对研究区域的天然径流变化的贡献率η_RH超出阈值的区域分别采用颜色标注在不同的图层上。

在计算完成后，需要将计算结果渲染于图片上，从而判断研究区域中哪些地方受到人类活动影响，以及受到多大的影响。为此，通过构建与贡献率对应的颜色区间，形成映射关系，然后渲染在研究区域地图中。需要注意的是，在该实施例中，将两个影响因素的贡献率进行分层渲染和叠加渲染，形成多个图层或颜色通道。能够便于使用者快速对各个因素的影响区域。从而做出更为准确的分析和判断。

另外，对于缺少数据的区域，采用插值法进行填充，插值法可以选用现有的各种方法，在此不再详述。

在研究过程中发现，在区域划分的时候，会出现色块不一致，或者重叠等问题，导致分界线不清晰，为此，给出如下技术方案。

如图7所示，根据本申请的一个方面，所述步骤S53进一步为：

步骤S53a、读取渲染后的研究区域地图，针对每一颜色通道，计算相邻两个像素之间的像素差，并计算该像素差与所述相邻两个像素之间最大像素值之间的比例，即为相差比i；

步骤S53b、针对每一颜色通道，计算相差比大于阈值的相邻像素；

步骤S53c、连接相邻像素，并判断相邻像素是否构成连续的分界线；若构成，查找气候变化对研究区域的天然径流变化的贡献率η_Rc超出阈值的一侧或人类活动对研究区域的天然径流变化的贡献率η_RH超出阈值的一侧，并采用颜色分别标注在不同的图层。

所述步骤53c还包括：

若不构成，依次判断相邻像素组成线段的相邻端点之间的欧式距离是否小于阈值，若小于，则连接两个端点，并作为分界线。

步骤S54：

判断分界线是否构成封闭区域，或分界线与研究区域边缘线是否构成封闭区域；若构成，采用颜色标注。

在这一实施例中，通过计算每一颜色通道中相邻像素之间的像素差，也就是像素梯度，判断相邻像素之间是否为分界线。从而能够快速查找出分界线，从而对图像进行准确分割。解决了目前图像渲染过程中存在的难题。

在本申请的另一实施例中，还提供一种考虑植被空间变化的天然径流变化归因分析系统，包括：

至少一个处理器；以及

与至少一个所述处理器通信连接的存储器；其中，

所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述任一项实施例所述的考虑植被空间变化的天然径流变化归因分析方法。

所述考虑植被空间变化的天然径流变化归因分析系统包括计算机设备和配置在计算机设备中的软件指令。

所述计算机设备包括通过系统总线相互通信连接存储器、处理器、网络接口。其中，本技术领域技术人员可以理解，这里的计算机设备是一种能够按照事先设定或存储的指令，自动进行数值计算和信息处理的设备，其硬件包括但不限于微处理器、专用集成电路、可编程门阵列、数字处理器、嵌入式设备等。

所述计算机设备可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述计算机设备可以与用户通过键盘、鼠标、遥控器、触摸板或声控设备等方式进入人机交互。

所述存储器至少包括一种类型的可读存储介质，所述可读存储介质包括内存、硬盘、随机访问存储器、只读存储器、磁盘、光盘等。在一些实施例中，所述存储器也可以是所述计算机设备的外部存储设备，例如该计算机设备上配备智能存储卡、SD卡等。当然，所述存储器还可以既包括所述计算机的内部存储单元也包括其外部存储设备。本实施例中，所述存储器通常用于存储安装于所述计算机设备的操作系统和各类应用软件，例如运行上述方法的计算机可读指令等。

所述处理器在一些实施例中可以是中央处理器、控制器、微控制器或其他数据处理芯片。该处理器通常用于控制所述计算机设备的总体操作。本实施例中，所述处理器用于运行所述存储器存储的计算机可读指令或者处理数据，例如运行上述方法的计算机可读指令。

总之，归一化植被指数（NDVI）变化可大致反映出下垫面参数的变化。将以NDVI表征的下垫面参数的变化分解为受气候变化和人类活动影响下的变化量，可更为准确对研究流域天然径流的变化进行归因。本实施例首先通过GWR方法对像元尺度下NDVI变化进行归因，依次分解下垫面参数受气候变化和人类活动的变化量，然后结合Budyko方法对天然径流变化进行归因。归一化植被指数（NDVI）作为公认的反映大尺度地表植被覆盖和生长状况的有效指标，通常用来表征耕地、林地、草地植被活动的强弱。

本发明针对当前天然径流变化归因研究中将下垫面参数变化完全归因于人类活动的错误认知，应对原有基于Budyko假设的径流变化归因方法进行改进，以提出更为合理有效的归因方法。一方面，在原有气候变化对研究区域天然径流变化的贡献率公式中加入一个修正项，该修正项为气候变化通过影响下垫面参数引起的径流变化贡献率，同时，在原有人类活动引起径流变化贡献率的公式中应减去该修正项，即得到人类活动通过影响下垫面参数引起的径流变化贡献率；另一方面，为准确求解该修正项，应对下垫面参数变化量进一步归因，结合研究区域土地利用特点，通过一定技术措施评估下垫面参数变化在特定时空尺度下分别受气候变化和人类活动的影响。通过改进原有Budyko假设在径流归因研究中的不足，以提高径流变化归因结果的准确性，是目前需要解决的问题。

以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。

Claims (10)

1.考虑植被空间变化的天然径流变化归因分析方法，其特征在于，包括如下步骤：

步骤S1、采集研究区域预定时间范围内的长时间序列数据，并将预定时间范围至少分为连续的第一时期和第二时期；分别计算第一时期的归一化植被指数平均值NDVI₁和第二时期的归一化植被指数平均值NDVI₂；并计算其差值∆NDVI_a，其中∆NDVI_a=NDVI₂- NDVI₁；

步骤S2、构建归一化植被指数与气象要素的地理回归关系，在约定下垫面参数恒定的情况下，分别计算第一时期归一化植被指数平均值NDVI_1s和第二时期的归一化植被指数平均值NDVI_2s；并计算其差值∆NDVI_C，即为气候变化引起的归一化植被指数平均值变化量；其中，∆NDVI_C =NDVI_2s-NDVI_1s；接着计算人类活动引起的归一化植被指数平均值变化量∆NDVI_H；∆NDVI_H =∆NDVI_a-∆NDVI_C；然后分别计算气候变化和人类活动引起的归一化植被指数的贡献率η_C和η_H；η_C=∆NDVI_C/ NDVI_a；η_H=∆NDVI_H/ NDVI_a；

步骤S3、基于天然径流对不同因素变化量的弹性系数构建天然径流的全微分方程；并通过全微分方程计算第二时期与第一时期的天然径流量平均值之差dR，以及降雨量变化、潜在蒸散发变化、下垫面参数变化导致天然径流量的改变量∆R_p=ε_P（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆R_L=ε_L（R/L）×∆L；

步骤S4、基于上述计算结果计算气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH：

η_Rc=(ΔR_P+ΔR_ET0)/ΔR×100%+η_C×ΔR_L/ΔR×100%；

η_RH=ΔR_L/ΔR×100%-η_C×ΔR_L/ΔR×100%=η_H×ΔR_L/ΔR×100%；

其中，ΔR=ΔR_P+ΔR_ET0+η_C×ΔR_L+η_H×ΔR_L；

η_C×ΔR_L/ΔR×100%为修正项，即气候变化通过影响下垫面参数引起径流变化的贡献率；

η_H×ΔR_L/ΔR×100%为人类活动通过影响下垫面参数引起的径流变化的贡献率。

2.根据权利要求1所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，所述步骤S1进一步为：

步骤S11、划定研究区域，并获取研究区域的地理数据；

步骤S12、采集研究区域预定时间范围内的长时间序列数据，所述长时间序列数据为中分辨率成像光谱仪数据。

3.根据权利要求2所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，所述步骤S1中将预定时间范围至少分为连续的第一时期和第二时期的过程进一步为：

步骤S1a、判断研究区域是否存在影响下垫面参数的事件，若存在，将事件发生的时间作为第一时期和第二时期的分界点；若不存在，进入步骤S1b、

步骤S1b、构建预定长度的时间区间，并在预定时间范围内滑动，形成至少两个间隔大于阈值的时期，记为第一检测时期和第二检测时期；在预定时间范围内选择至少一个时期作为标准时期；分别计算第一检测时期、第二检测时期和标准时期的归一化植被指数平均值；

步骤S1c、分别计算第二检测时期与标准时期的归一化植被指数平均值的差值，记为第一变化值，第一检测时期与标准归一化植被指数平均值的差值，记为第二变化值，

计算第一变化值与第二变化值差值，并与阈值进行比较，若大于阈值，将第一时期和第二时期的起点时间作为下垫面参数改变的待检时刻，放入待检时刻集合。

4.根据权利要求1所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，所述步骤S2进一步为：

步骤S21、构建地理加权回归模型，采集研究区域内第一时期的降雨量、平均气温和空间位置参数，采用局部加权最小二乘法求解所述地理加权回归模型，获得第一时期的归一化植被指数并计算平均值NDVI_1s；

步骤S22、采集第二时期的降雨量、平均气温和空间位置参数，设定归一化植被指数与气象要素的地理回归关系在第一时期和第二时期相同；构建地理加权回归模型，并采用局部加权最小二乘法求解所述地理加权回归模型，获得第二时期的归一化植被指数并计算平均值NDVI_2s；

步骤S23、计算第二时期和第一时期的归一化植被指数并计算平均值的差值∆NDVI_C，∆NDVI_C =NDVI_2s-NDVI_1s；即为气候变化引起的归一化植被指数平均值变化量；

步骤S24、计算人类活动引起的归一化植被指数平均值变化量∆NDVI_H；∆NDVI_H =∆NDVI_a-∆NDVI_C；

步骤S25、分别计算气候变化和人类活动引起的归一化植被指数的贡献率η_C和η_H；η_C=∆NDVI_C/ NDVI_a；η_H=∆NDVI_H/ NDVI_a。

5.根据权利要求4所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，所述步骤S3进一步为：

步骤S31、构建研究区域的水平衡量方程R=P-ET-∆S；

其中，R为径流量，P为降水量，ET为实际蒸散发量，∆S为自然流域下的陆地水储量变化量；

步骤S32、构建干旱指数的函数，F（φ）=ET/P，其中φ为干旱指数，

其中，

，n为下垫面参数；

步骤S33、根据天然径流对变量的弹性系数的定义，构建天然径流的微分方程；

；

ε_P、ε_ET0、ε_L分别为天然径流对降水、潜在蒸散发和下垫面参数的弹性系数；若弹性系数为正值，表示径流随着该变量的增加而增加，若弹性系数为负值，表示径流随着该变量的增加而减小；

步骤S33、用天然径流第二时期和第一时期的多年平均径流量之差dR表示第一时期到第二时期年径流的变化dR=R₂-R₁，R₂和R₁分别为第二时期和第一时期的年径流量；

步骤S34、计算降雨量变化、潜在蒸散发量变化、下垫面参数变化导致天然径流量的改变量∆R_p=ε_p（R/P）×∆P、∆R_ET0=ε_ET0（R/ET0）×∆ET0、∆R_L=ε_L（R/L）×∆L。

6.根据权利要求5所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，还包括步骤S5：

步骤S51、获取研究区地图及各处的气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH；

步骤S52、分别采用气候变化对研究区域的天然径流变化的贡献率η_Rc、人类活动对研究区域的天然径流变化的贡献率η_RH渲染研究区域地图；

步骤S53、查找研究区域地图上的分界线，并将气候变化对研究区域的天然径流变化的贡献率η_Rc超出阈值、人类活动对研究区域的天然径流变化的贡献率η_RH超出阈值的区域分别采用颜色标注在不同的图层上。

7.根据权利要求6所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，所述步骤S53进一步为：

步骤S53a、读取渲染后的研究区域地图，针对每一颜色通道，计算相邻两个像素之间的像素差，并计算该像素差与所述相邻两个像素之间最大像素值之间的比例，即为相差比i；

步骤S53b、针对每一颜色通道，计算相差比大于阈值的相邻像素；

步骤S53c、连接相邻像素，并判断相邻像素是否构成连续的分界线；若构成，查找气候变化对研究区域的天然径流变化的贡献率η_Rc超出阈值的一侧或人类活动对研究区域的天然径流变化的贡献率η_RH超出阈值的一侧，并采用颜色分别标注在不同的图层。

8.根据权利要求7所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，所述步骤53c还包括：

若不构成，依次判断相邻像素组成线段的相邻端点之间的欧式距离是否小于阈值，若小于，则连接两个端点，并作为分界线。

9.根据权利要求7所述的考虑植被空间变化的天然径流变化归因分析方法，其特征在于，还包括步骤S54：

判断分界线是否构成封闭区域，或分界线与研究区域边缘线是否构成封闭区域；若构成，采用颜色标注。

10.一种考虑植被空间变化的天然径流变化归因分析系统，其特征在于，包括：

至少一个处理器；以及

与至少一个所述处理器通信连接的存储器；其中，

所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现权利要求1~9任一项所述的考虑植被空间变化的天然径流变化归因分析方法。

Images