CN110766299B - 一种基于遥感数据的流域植被变化分析方法 - Google Patents

Abstract

本发明为一种基于遥感数据的流域植被变化分析方法，包括以下步骤：第一步骤，选用NASA发布的植被遥感数据，提取研究流域相应数据；第二步骤，植被覆盖度计算，选用像元二分模型法；第三步骤，基于Savitzky‑Golay算法的数据重构；第四步骤，通过将NDVI数据拟合到双逻辑斯蒂函数来获取地表物候；第五步骤，对获取的物候特征值进行Mann‑Kendall秩次相关趋势检验。本发明能够有效地对流域植被变化进行定量化的分析评价，简单可靠。

Description

一种基于遥感数据的流域植被变化分析方法

技术领域

本发明属于一种流域植被变化分析方法，具体涉及一种基于遥感数据的流域植被变化分析方法。

背景技术

河流的水沙情势的深入分析意义重大，它可以关乎下游宽河段治理的方向，关乎水资源的配置和利用策略，关乎水沙调控工程的布局和整体治理方略的确定。以黄河为例，21世纪以来黄河水沙情势发生史无前例的剧烈变化，而对黄河水沙变化机理的深入剖析需要对流域内各个影响因子的全面评价，其中植被是不可忽略的重要部分。森林、草地与耕地等地表类型都可算作为植被，在环境中发挥了重要的作用，如固土、防沙、缓流、截流等。但当下黄河水沙研究中植被的部分多囿于水土保持功能，没有充分利用植被在自然界中的指示能力，植被包含的丰富信息没有被完全挖掘。而且退耕还林举措虽然减沙效果显著，但却对黄河径流是否有消极影响，近年来有不少争议。有观点认为退耕还林政策的继续推进将对社会和自然环境造成更大的伤害，确保当地的粮食与水资源的供应应当置于首要位置，与不断改善植被覆盖情况相比，保持与气候条件、水文情况和土壤侵蚀水平相当的植被覆盖度则更有利于当地的可持续发展。在此背景下，急需一种更加全面完善的流域植被变化分析方法。

发明内容

为了解决上述流域植被分析不全面的问题，本发明提出一种基于遥感数据的流域植被变化分析方法。

本发明的技术方案如下：

一种基于遥感数据的流域植被变化分析方法，其特征在于包括以下步骤：

第一步骤，选用NASA发布的植被遥感数据，提取研究流域相应数据；

第二步骤，植被覆盖度计算，选用像元二分模型法，像元二分模型基本方程如下：

NDVI＝FVC·NDVIveg+(1-FVC)NDVIsoil，

植被覆盖度基本计算公式如下：

FVC＝(NDVI-NDVImin)/(NDVImax-NDVImin)

其中，FVC代表植被覆盖度；NDVIveg代表完全被植被覆盖地区的NDVI值；NDVIsoil为裸土或是没有植被覆盖地区的NDVI值；NDVImax和NDVImin分别为区域内的NDVI最大值和最小值；

第三步骤，基于Savitzky-Golay算法的数据重构，

Savitzky-Golay滤波一般方程如下：

其中NDVI是NDVI初始数据，NDVI^*是经过重构处理的数据序列，W_i是在一次局部拟合平滑过程中滤波器内第i个初始数据的权重值，下标j是数据在初始数据序列的位置，N表示一个滑动窗口内处理的数据数量，m是滑动窗口宽度的一半大小，N＝2m+1；

第四步骤，通过将NDVI数据拟合到双逻辑斯蒂函数来获取地表物候，拟合方程如下：

其中f(t)是特定日期t的NDVI值；v₁是全年的背景NDVI水平；v₂是全年NDVI的振幅。m₁、n₁和m₂、n₂分别是一对参数，两个参数m和n用于确定NDVI增大阶段和减小阶段的总体斜率和基本相位；

第五步骤，对获取的物候特征值进行Mann-Kendall秩次相关趋势检验，

在公式中，

var(S)＝n(n-1)(2n+5)/18

式中：1<j<i<n，β所代表的是斜率，β为正表示上升，β为负表示下降，值的大小表示趋势显著与否；

零假设H₀：β＝0；当|Z|>Z_1-α/2，拒绝H₀假设；Z>Z_1-α/2时，序列有显著上升趋势；Z<-Z_1-α/2时，序列有显著下降趋势；式中：Z_1-α/2为标准正态方差，α为显著性检验水平。

优选的，第一步骤中选用NASA发布的第3代NDVI数据集GIMMS NDVI 3g的植被遥感数据。

优选的，第二步骤中，

优选的，第三步骤中，NDVI数据序列重构的步骤如下：

(1)用Savitzky-Golay滤波方法获得NDVI0序列的长期变化趋势拟合线，需要确定拟合多项式的次数b与拟合窗口的大小m，这里的两个主要参数根据已有的滤波效果指数研究取b＝4，m＝3，将对NDVI⁰进行第一次滤波后的序列记为NDVI^tr；

(2)评估每一个原始数据的质量即可靠度，根据数据点与包络线上的对应值的差值来计算权重W_i，根据偏大的数值监测误差更小、可信度更高的假设，原始数据值越大，权重W_i更大。W_i的计算公式为：

式中，

d_max是d_i的最大值；

(3)生成新的NDVI序列，仍依据遥感中偏大的数值可信度更高的基础假设，选择重构后数据序列与原始NDVI数据序列中相同位置对应的两个数据点中较大的一个作为新序列的值

(4)用同样的方法再次拟合上一步得到的

序列，拟合后的序列结果标记为

k＝1代表第一次拟合；

(5)对重构后的数据序列的质量和滤波的效果进行评价，滤波效果指数F_k能够体现重构后的数据序列接近上包络线的程度，F_k计算方法如下：

式中，

表示第k次滤波后得到的数据序列中第i个NDVI值，

是未经处理序列中相同位置的NDVI初始值，W_i是(2)中计算得到的相对应位置的数据权重值；

循环第(3)-(5)步，得到新的NDVI序列

第6步：退出循环的条件：F_k-1≥F_k≤F_k+1。

优选的，第四步骤中，该模型共有6个参数，全部参数需要率定，该非线性最小二乘优化过程选用Levenberg-Marquardt算法求解，计算拟合后曲线的斜率，斜率最大值处对应的时间为生长期开始时间Start-of-season，SOS，斜率最小处对应的时间为生长期结束时间End-of-season，EOS，两者之差为生长期长度Growing Season Length，GSL。

优选的，采用分区物候分析和/或逐项元物候分析。

进一步优选的，所述分区物候分析采用五个分区。

本发明的有益技术效果：

本发明的方法，第一步是确立植被遥感数据，优选采用NASA发布的第3代NDVI数据集GIMMS NDVI 3g，该数据集汇编自美国国家海洋和大气管理局(NOAA)卫星上的高级超高分辨率辐射计(AVHRR)传感器采集的NDVI图像。5种不同的卫星覆盖了34年的数据：NOAA-7，9，11，14和16。NDVI图像来自AVHRR通道1和2图像，分别对应红色(0.58～0.68μm)和红外波长(0.73～1.1μm)。时间跨度为1982—2015年，时间分辨率为15d，每月2景，上半月和下半月各一景，空间分辨率为0.083°×0.083°。

图像是通过最大值合成(Maximum Value Compositing，MVC)技术获得的，该方法最大限度地减少了大气气溶胶和云的影响，且处理过程校正了天顶角、校准损失，轨道漂移，火山爆发等影响因素，误差减小，精度更高。该数据集虽然没有其他信道信息作为参考，仅为单一NDVI信息，但已经根据其前身Pathfinder AVHRR Land(PAL)数据集提出了若干改进。一方面改进数据处理方式，包括导航，传感器校准和平流层气溶胶的大气校正。另一方面主要是通过经验模式分解(EMD)技术来修正NOAA的轨道漂移。GIMMS NDVI 3g数据集为目前可获得的时间跨度最长的NDVI序列数据集，适用于大范围植被覆盖变化的长期监测。本文研究区黄河流域约7.95×10⁵km²，占到了中国国土面积的8.3％，且长时间跨度的遥感数据能更加充分地体现植被时空演化变化。

图像中每个像元的NDVI值可以被认为是植被覆盖组分的NDVI值和没有植被覆盖的裸土组分的NDVI值加权平均的结果。理论上，NDVIsoil应该是非常接近零的相对稳定值，但实际上，由于大气干扰、云层污染和地表湿度温度差异等影响因子，NDVIsoil的值会有所波动。同理，NDVIveg也不能通过一个定值简单确定。

在这个植被覆盖度计算模型中，NDVIsoil和NDVIveg的可以通过上述公式计算，计算参数有FVCmax、FVCmin、NDVImax及NDVImin。其中取FVCmax＝1和FVCmin＝0。

植被覆盖度计算公式可简化如下：

FVC＝(NDVI-NDVImin)/(NDVImax-NDVImin) (2.3)

区域内的NDVI最大值和最小值分别记为FVCmax和FVCmin。然而，因为遥感中噪声的干扰难以完全剔除，且即使在植被覆盖度均达到了100％的情况下，植物种类不同也引起NDVI值的差异，所以直接取NDVI序列中的最大值和最小值作为参数的误差较大。通常，NDVImax和NDVImin取一定置信区间上的极值。在本发明中，置信水平取95％，即NDVImin是取一年内所有像元NDVI累计频率为5％的值，NDVImax则是累计频率为95％的NDVI值。

相对于单逻辑斯蒂函数和其他拟合方法，双逻辑斯蒂函数有以下优越性：

(1)双逻辑斯蒂函数拟合不需要将一年的NDVI时间序列划分为两个阶段，可以实现年尺度上的变化整体拟合，平滑且连续性好。

(2)双逻辑斯蒂函数拟合的六个参数有明确的实际意义，在运用非线性拟合算法确定参数时易于确定初始值范围。

(3)双逻辑斯蒂函数能够很好地刻画NDVI时序曲线中的平台阶段，对应于植被冬季的休眠期或夏季持续的光合作用高峰期。这些平台阶段的拟合需要很多组谐波才能够拟合完善，这对于逐个像元的黄河流域全局分析来说过于复杂费时。

(4)黄河流域内存在一定范围的农业区，由于农业区中存在轮种的情况，两个种植周期会导致像元NDVI年变化曲线会出现两个波峰，单逻辑斯蒂函数无法实现双峰NDVI模式的拟合，而双逻辑斯蒂函数能够精准识别出第一个波峰的上升阶段和第二个波峰的下降阶段，实现两个波峰的相连。

Mann-Kendall非参数统计方法是数据趋势检测中一个非常有用且已经得到广泛应用的方法，其主要特点有：(1)不针对特定参数；(2)不需要严格假定变量分布，能够分析没有明确分布、含随机性的数据序列；(3)极端值对结果的干扰小，对数据降噪的要求减弱，并且允许数据序列有空白值；(4)适用于微小值的趋势分析，这是因为分析过程是相对于数量级而不是数字本身的大小。

Mann-Kendall检验是世界气象组织推荐的一种水文气象数据分析方法。该方法能够有效地确定某一个自然过程有无显著的变化趋势。本文研究的植被物候期没有明确的分布规律，因此Mann-Kendall秩次相关检验具有显著的优越性。

本发明双逻辑斯蒂函数对NDVI时序曲线的拟合效果较好，能够准确刻画NDVI变化的平台期，结合基于Savitzky-Golay滤波器的数据重构方法和最大斜率法能够得到较为稳定可靠的物候期识别结果，适用面广，且不需要很多物候识别方法要求的经验结合。

附图说明

图1为本发明实施例1黄河流域植被覆盖度变化；

图2为本发明实施例1黄河流域五分区示意图；

图3为本发明的实施例1根据图2的五分区NDVI 1982～2016年时间变化曲线；

图4为本发明的实施例1生长期总时长作时间变化曲线的线性回归分析；

图5为本发明的实施例2生长期开始时间(SOS)对比；

图6为本发明的实施例2生长期结束时间(EOS)对比；

图7为本发明的实施例2生长期长度(GSL)对比；

图8为本发明的实施例2生长期长度(GSL)变化趋势；

图9为本发明的实施例2生长期长度(GSL)变化速率

具体实施方式

为了更清楚理解本发明内容，将结合附图1-9和具体实施方式详细说明。

实施例1

本实施例涉及地处96°～119°E，32°～42°N的黄河流域，面积广大，流域面积约7.95×105km²，包含4.2×104km²的内流区。黄河流域地跨我国九省(自治区)，流域内地势总体特征为为西高东低，根据海拔可近似为三个台阶，第一级台阶主要是海拔在3000m以上的青藏高原地区，该区域是黄河的主要来水区，水土流失程度很轻；中部海拔在1000～2000m之间黄土高原是黄河流域台阶的第二级，该区域大面积有严重的水土流失情况，是黄河的主要来沙区。黄河闻名于世的高含沙量主要是源于其流域内大面积的黄土地区，达到了640,000km²，，其中212,000km²的严重侵蚀区域提供了黄河中大约90％的泥沙。第三级阶梯是由黄河下游冲积平原和鲁中丘陵组成，该区域的产流产沙量在黄河水沙总量中占比很小。研究区内各地气候差距很大，干旱、半干旱、半湿润及湿润四大类气候类型皆有，半湿润区与半干旱区的分界线大体与400mm降雨等值线一致，不同的地貌和气候条件造就了黄河流域丰富多样的植被类型，高山草甸、灌木林、森林及农耕地等类型都能在黄河流域中找到，各地差异极大。

黄河约为其流域1.07亿人提供淡水，约占中国总人数的8.7％，是我国西北和华北地区主要水源。黄河水文形势与当地自然环境与社会经济紧密相连，环境与发展之间的矛盾使黄河流域一直备受关注。

1、植被覆盖度变化

选择NASA发布的第3代NDVI数据集GIMMS NDVI 3g的植被遥感数据，提取黄河流域NDVI数据，依据像元二分模型估算植被覆盖度，计算式如下：

FVC＝(NDVI-NDVImin)/(NDVImax-NDVImin)

其中，FVC代表植被覆盖度；NDVIveg代表完全被植被覆盖地区的NDVI值；NDVIsoil为裸土或是没有植被覆盖地区的NDVI值；NDVImax和NDVImin分别为区域内的NDVI最大值和最小值。

计算得到34年每年24景的植被覆盖度结果，图1选取了1982年、1999年、2015年七月下半月的计算结果。选择了1982、1999、2015三组图像是由于1982年是能够获得的最早数据，2015年组是较新的完整一年数据，中间插入1999年是不仅是因为这是这段时间的中点位置，更是因为退耕还林政策在当年实施，黄土高原植被在此人为干预下快速改变。

根据前述的步骤进行数据的处理。由于是通过像元二分模型进行计算，NDVIveg和NDVIsoil这两个直接关系覆盖度计算结果精度的参数是直接通过统计获取，分别取当年NDVI数据累计频率为95％与5％的值，没有与实测数据进行对比校准，覆盖度估算结果与实测结果可能存在一定的偏差，但由于所有处理方法与标准是相同的，各地覆盖度计算结果为NDVI数据的线性拉伸，得到的各地覆盖度对比情况是可靠的。图像中植被变化尤为突出的区域是位于黄土高原的主要产沙区，1982年到1999年该区域植被覆盖度有所提高，在1999年退耕还林政策施行之后，可看到2015年黑色面积相对于1999年大幅度拓展，植被覆盖度增长速度加快。各地径流试验小区开展的研究已经证明，无论植被类型如何，当区域植被覆盖度达60～70％时，水土流失的程度是非常轻微的，而植被覆盖度在10％～20％之间时，植被基本没有水土保持的作用。根据图1呈现的计算结果，在1982年，黄河流域主要产沙区植被覆盖度达60％的区域不足30％，而在2015年，黄土高原中黄河主要产沙区内70％左右的土地植被覆盖度达到了60％以上，可见该地区水土保持能力相较于1982年大幅度提高，这应当是进入21世纪以来黄河含沙量骤降的一个关键性原因。

2、分区物候分析

本发明旨在从物候角度提取黄河流域植被特征值，并以此来耦合黄河水沙变化。为此，提出以植被生长期长度作为一个代表性指标。植被生长期由植物种类、地理位置和气候条件等多方面因素决定，与其他植被指数相比，具有综合性的优势，生长期不仅能够在一年内植物生长状况，还对气温、降水等条件因素敏感，能够反映气候变化，同时能够体现人类活动的影响，如耕作方式和土地类型的改变等。在研究河流产沙中，植被生长期也能够直接反映植被起作用的时长。

通过选择一景中研究区域内NDVI最大值作为代表值来实现由面到点的转换，由区域内NDVI最大值随时间的变化序列来提取代表生长期。在每一年同一时期下，区域内NDVI最大值直接由植物量、植物种类和雨热条件的影响，选取最大值的稳定性与可靠性较高。在相同情况下，NDVI最大的值植被生长状况最好，不仅仅在物种上有代表性，对气候条件的反应也最佳。尽管忽略了NDVI在空间尺度上的变异性，但可以通过区域的合理划分在一定程度上减小这种空间差异。

由于黄河流域面积广大，各地特征不同，对黄河流域的植被研究最终是要结合到黄河水沙研究中，根据现今黄河水沙变异性可将黄河流域划分为五个区域，如图2所示。

提取每个区域每景NDVI最大值组成代表NDVI时序曲线(见图3)，再通过上述的数据重构方法与物候提取方法计算该区域的生长期，以此作为该区域植被物候层面特征值。

先对各个区域代表NDVI时序曲线进行基于Savitzky-Golay滤波器的数据重构，Savitzky-Golay滤波一般方程如下：

接着拟合到双逻辑斯蒂函数来获取地表物候，拟合方程如下：

拟合需要率定6个参数，该非线性最小二乘优化过程选用普遍使用的Levenberg-Marquardt算法求解。计算拟合后曲线的斜率，斜率最大值处对应的时间为生长期开始时间(Start-of-season，SOS)，斜率最小处对应的时间为生长期结束时间(End-of-season，EOS)，两者之差生长期长度(Growing Season Length，GSL)。

观察图3可得，五个区域NDVI时序曲线均表现了非常好的周期性与稳定性。在人类活动剧烈的黄河下游区域#2，时序曲线出现了显著的双峰型，能够体现黄河下游冲积平原夏玉米、冬小麦套作的情况。

对各个区域的NDVI最大值时间变化序列应用前文介绍的数据重构方法与物候识别方法，得到五个区域1982～2016年每年的生长期开始时间、生长期结束时间和生长期总时长。表1罗列了五个区域34年生长期长度计算结果，仅#2区域在2006年出现了无法识别的情况，其余均得到了有效值。在#2区域在2006年初春出现了与夏季相当的明显波峰干扰，双逻辑斯蒂函数无法有效拟合造成了识别错误，该波峰在34年间仅出现一次，气候改变不会出现如此剧烈的波动，可能是该区域的人类活动所致，也有可能仅仅是原始数据的误差造成，该误差出现频率很低，在接受范围以内。

对每个区域的生长期总时长作时间变化曲线，为了确定各个区域生长期变化是否有显著趋势，采用了线性回归分析与Mann-Kendall检验进行趋势分析。

Mann-Kendall秩次相关趋势检验计算式如下：

在公式中，

var(S)＝n(n-1)(2n+5)/18

式中：1<j<i<n，β所代表的是斜率，β为正表示上升，β为负表示下降，值的大小表示趋势显著与否；

零假设H₀：β＝0；当|Z|>Z_1-α/2，拒绝H₀假设；Z>Z_1-α/2时，序列有显著上升趋势；Z<-Z_1-α/2时，序列有显著下降趋势；式中：Z_1-α/2为标准正态方差，α为显著性检验水平。

检验结果证明各个区域生长期在一定范围内波动且在几个区域中具有显著趋势性(表1)。

如图4的线性回归分析中，五个区域生长期的总体趋势均为延长，其中1，2，5在0.95的置信区间下表现了显著性，P值分别为1.55×10^-5，3.22×10^-3及0.0488。在Mann-Kenndal检验中，五个区域总体趋势均为延长，1，2，5三区域同样在0.95的置信度下表现了显著趋势。

表1分区物候识别结果——生长期长度

五个区域的延长趋势在一定程度上印证了近几十年来黄河流域气温总体上升的趋势。#5区域为黄河源地区，人为因素影响小，显著的生长期延长趋势集中体现了气候条件的影响。#1、#2区域表现了极为显著的特征，同时延长的速率最高，这两个区域的生长期延长中人类活动的影响为主要因素。1#区域主要为黄土高原区，含大面积的黄河流域产沙区，近几十年来大规模水土保持举措在此实施，尤其是1999年来的退耕还林和封山禁牧政策使该区域的植被覆盖度以极快的速度增长。植被类型的改变能在物候期中有突出体现。与20世纪后期相比，如今黄土高原的林草植被结构已经有了极大的改变，过去黄土高原除了裸地外，植被主要是人工种植的乔木和灌木植物，有物种单一且地被物密度小的特点。在退耕还林政策下，近十几年黄土高原的新增植被主要来自自然修复形成的草木植物或是灌木植被，郁闭度高、物种多样、植被紧贴地表，近两年的实地调查已经可以看到草灌下的枯落物和苔藓。#2区域为黄河下游平原区，以灌区为主，是《国家主体功能区规划》中划定的全国重要粮食生产区，冬小麦是当地主要的粮食作物，且通常都与夏玉米进行套作，该地区近30年粮食产量连续增长，生长期直接受耕种方式和耕种作物的影响。#2区域生长期的延长说明了一年内有效耕种时间的延长，播种时间提前，收割时间延迟，农田的利用效率提高，体现了耕作技术的进步。最后，上述结果有力地说明了该物候指数提取方法的可行性与可靠性。

实施例2

逐像元物候分析

由于实施例1划分的五个区域内存在气候条件和地貌条件差异较大的情况，以点概面必定会造成很多信息的丢失。为了探究黄河流域空间域上的分布，实施前述的第三步骤-第五步骤对黄河流域每个像元进行了物候期分析，像元总数11698个。每年的物候识别结果包括生长期开始时间(Start-of-season，SOS)、结束时间(End-of-season，EOS)、总时长(Growing Season Length，GSL)，将三者分别导出为tiff图像文件以观察各地差异，得到1982～2015年共34组图像文件。接着对每个像元的物候期变化进行Mann-Kendall检验，判断三种特征值的数值变化。对于生长期长度，将时长有显著增长趋势(数值显著增加)的像元记为1，显著缩短趋势(数值显著减小)的像元记为-1，无明显趋势的像元标记为0，生长期开始时间与生长期结束时间相同，在当年发生的时间显著推迟(数值显著增加)的像元记为1，发生的时间显著提前(数值显著减小)的像元记为-1，无明显趋势的像元记为0，导出趋势分布图。与此同时记录每个像元检验过程中评价变化速率的β，导出变化速率图。

生长期分布总体上西部短、东部长与农业区生长期长的现象，这与黄河流域地势西高东低及从西到东气候逐渐湿润的特点一致。黄河源区地处高原，气候干旱寒冷，以草地为主，生长期开始晚、结束早，总时间相对于东部地区较短。根据黄河流域水资源评价成果，黄土高原1956～2000年系列年降水量为300～800mm，水面蒸发量为700～1200mm，由西北向东南，年降水量递增，水面蒸发量递减。。在黄土高原物候识别结果中，灌区以外的区域生长期长度呈现了与降水量一致的东南长西北低的现象，相同的趋势在一定程度上说明了植被生长期与降水的密切关系。从34年物候期识别结果中依旧选择了1982、1999、2015三组图像进行对比(图5，图6，图7)，可以直观地观察到在1999年以后，深色区域(较长生长期)在黄土高原中部产沙区及黄河下游地区的延展与颜色加深，说明了这两个区域植被条件的大幅度改善。这种变化与植被覆盖度的改变(图1)呼应，证明了该物候识别方法的可行性。

提取三个时间节点的计算结果进行对比能够大致判断出一些比较突出的整体变化，而对于单个像元难以判断趋势，且生长期对气候条件的变化较为敏感，这体现为较大波动幅度。为了对生长期变化进行更为精细的分析，对每一个像元进行了趋势分析。

图8为M-K检验分析后的GSL趋势分布图，黑色点表示在95％置信水平下，SOS数值有显著减小趋势，即生长期开始时间提前，白色点则反之，表示SOS数值显著增大，即生长期开始时间推迟，灰色表示趋势不显著。图9为GSL在M-K分析中得到的β值分布图，表征黄河流域生长期开始时间的年变化速率。生长期持续时间(GSL)变化趋势分布与结束时间(EOS)非常接近，由于GSL＝EOS-SOS，而SOS的变化幅度普遍要小于EOS，因此GSL的趋势与EOS趋同。与前文1982、1999及2015三时间节点GSL图像对比得到黄河主要产沙区及黄河下游生长期延长的结论一致，95％以上的生长期显著延长点(黑色点)分布在这些区域。

生长期持续时间(GSL)的变化速率分布同生长期结束时间(EOS)速率分布基本一致。对比黄河流域降雨及蒸发等值线图，生长期缩短速率大于0.5的区域95％落在降雨量小于400mm的干旱半干旱区，生长期延长速率大于0.5的区域有95％以上落在重点治理的黄土丘陵沟壑区及黄河下游地区，黄河源区颜色普遍较浅，变化幅度小于0.5。物候期延长的最大幅度达到了5.5，缩短的最大幅度为-2.3左右，绝对值不到最大延长速率的二分之一。由于叠加了生长期开始时间变化，生长期总时间变化速率的范围大于结束时间的变化范围。

近年来我国西北干旱区NDVI和植被覆盖度均呈下降的趋势，而我国东部湿润地区NDVI和植被覆盖度则主要为上升趋势。结合在东部季风区温度变化主导、在西北干旱半干旱区植被的变化降雨因素主导的论断，以及森林年变化主要受温度影响，草原植被与荒漠等干旱半干旱地区的植被类型主要受降水因素控制的观点，可以推测黄河流域植被物候变化在气候角度的原因。很多研究都证明了黄河流域近40年来气温的上升趋势与降雨的减少趋势。增温对植物生长的影响具有两面性：积极作用是更高的温度能够促进植被的光合作用，提高水分利用率，从而对植被的生长有益。消极作用是随着温度上升，植物的水分消耗也增加，干旱的情况下植物缺水更为严重，从而植被生长受阻。升温对植物的影响两面性可能是黄河流域生长期变化呈现两个分区的气候方面一个原因，黄土高原中的湿润半湿润地区，植物生长所需的水分较为充足，因此气温上升促进植物光合作用，提高水分利用率，进而有利于该地区植物的生长，生长期延长；而黄河源和黄土高原西北部地处干旱半干旱区，增温使植物的蒸腾量增加，本就缺水的地区更为干旱，从而阻碍了当地植物的生长。与此同时，黄河流域的降水也呈现了减少的趋势，干旱半干旱区的植被对降水情况更为敏感，降水减少使植物生长环境更加恶劣，而对较为湿润的地区，降雨减少造成的生长负面影响要小于升温引起的正面影响，从而生长期呈延长趋势。这能够从气候角度解释黄河流域生长期变化趋势呈现两个区域且分界线接近400mm等降水量线的现象。

综上，由以上实施例可以看出，双逻辑斯蒂函数对NDVI时序曲线的拟合效果较好，能够准确刻画NDVI变化的平台期，结合基于Savitzky-Golay滤波器的数据重构方法和最大斜率法能够得到较为稳定可靠的物候期识别结果，适用面广，且不需要很多物候识别方法要求的经验结合。提取流域内物候特征值时采用区域内NDVI最大值实现空间域上由点到面的压缩是一种简单可行的方法。实施例2计算了黄河流域植被生长期开始时间(SOS)、结束时间(EOS)及总时长(GSL)三个指标并进行了趋势检验，结合植被覆盖度变化进行了分区物候分析与逐像元物候分析。结果证明本发明能够有效地对流域植被变化进行定量化的分析评价，简单可靠。

Claims (2)

1.一种基于遥感数据的流域植被变化分析方法，其特征在于，针对地处96°～119°E，32°～42°N的黄河流域，包含4.2×10⁴ km²的内流区，包括以下步骤：

第一步骤，选用NASA发布的第3代NDVI数据集GIMMS NDVI 3g的植被遥感数据，提取研究流域相应数据；

第二步骤，植被覆盖度计算，选用像元二分模型，像元二分模型基本方程如下：

，

植被覆盖度计算公式如下：

其中，FVC代表植被覆盖度；NDVIveg代表完全被植被覆盖地区的NDVI值；NDVIsoil为裸土或是没有植被覆盖地区的NDVI值；NDVImax和NDVImin分别为区域内的NDVI最大值和最小值；

计算得到每年24景的植被覆盖度结果，至少选取1982年、1999年和2015年七月下半月的计算结果，NDVIveg和NDVIsoil这两个参数是直接通过统计获取，分别取当年NDVI数据累计频率为95％与5％的值；

第三步骤，从物候角度提取黄河流域植被特征值，并以此来耦合黄河水沙变化，提出以植被生长期长度作为一个代表性指标；通过选择一景中研究区域内NDVI最大值作为代表值来实现由面到点的转换，由区域内NDVI最大值随时间的变化序列来提取代表生长期；

根据黄河水沙变异性将黄河流域划分为五个区域，提取每个区域每景NDVI最大值组成代表NDVI时序曲线，再通过数据重构方法与物候提取方法计算该区域的生长期，以此作为该区域植被物候层面特征值；

对各个区域代表NDVI时序曲线基于Savitzky-Golay滤波方法进行数据重构，

Savitzky-Golay滤波方法一般方程如下：

其中，

是经过重构处理的数据序列，

是在一次局部拟合平滑过程中滤波器内第i个初始数据的权重，下标j是数据在初始数据序列的位置，N表示一个滑动窗口内处理的数据数量，m是拟合窗口的大小，

；

其中，NDVI数据序列重构的步骤如下：

(1)用Savitzky-Golay滤波方法获得

序列的长期变化趋势拟合线，确定拟合多项式的次数b与拟合窗口的大小m，这里的两个参数根据已有的滤波效果指数研究取b＝4，m＝3，将对

序列进行第一次滤波后的数据序列记为

；

(2)评估每一个原始数据的质量即可靠度，根据数据点与包络线上的对应值的差值来计算权重

，

的计算公式为：

式中，

，

是

中的最大值；

(3)生成新的NDVI序列，选择第一次滤波后的数据序列与原始NDVI数据序列中相同位置对应的两个数据点中较大的一个作为新序列的值：

(4)用同样的方法再次拟合上一步得到的

序列，拟合后的序列标记为

，k=1代表第一次拟合；

(5)对拟合后的序列的质量和滤波的效果进行评价，滤波效果指数

能够体现拟合后的序列接近上包络线的程度，

计算方法如下：

式中，

表示第k次滤波后得到的数据序列中第i个NDVI值，

是未经处理序列中相同位置的NDVI初始值，

是计算得到的相对应位置的权重；

循环第(3) - (5)步，得到新的NDVI序列

：

；

(6)退出循环的条件：

且

；

第四步骤，通过将NDVI数据拟合到双逻辑斯蒂函数来获取地表物候特征值，拟合方程如下：

其中

是特定日期t的NDVI值；

是全年的背景NDVI水平；

是全年NDVI的振幅，

、

和

、

分别是一对参数，用于确定NDVI增大阶段和减小阶段的总体斜率和基本相位，该拟合方程共有6个参数，全部参数需要率定，非线性最小二乘优化过程选用Levenberg-Marquardt算法求解，计算拟合后曲线的斜率，斜率最大值处对应的时间为生长期开始时间SOS，斜率最小处对应的时间为生长期结束时间EOS，两者之差为生长期长度GSL；

第五步骤，对获取的物候特征值进行Mann-Kendall秩次相关趋势检验，

在公式中，

式中：1<i<j<c，β所代表的是斜率，β为正表示上升，β为负表示下降，值的大小表示趋势显著与否；

零假设H₀：β＝0；当|Z|>Z_1-α/2，拒绝H₀假设；Z>Z_1-α/2时，序列有显著上升趋势；Z<-Z_1-α/2时，序列有显著下降趋势；式中：Z_1-α/2为标准正态方差，α为显著性检验水平；

或根据黄河水沙变异性进行逐项元物候分析，像元总数11698个，每年的物候识别结果包括生长期开始时间SOS、结束时间EOS、生长期长度GSL，将三者分别导出为tiff图像文件，接着对每个像元的物候期变化进行Mann-Kendall秩次相关趋势检验，判断三种特征值的数值变化，对于生长期长度，将时长有显著增长趋势的像元记为1，显著缩短趋势的像元记为-1，无明显趋势的像元标记为0，生长期开始时间与生长期结束时间相同，在当年发生的时间显著推迟的像元记为1，发生的时间显著提前的像元记为-1，无明显趋势的像元记为0，导出趋势分布图，与此同时记录每个像元检验过程中评价变化速率的

，导出变化速率图，对每一个像元进行趋势分析。

2.根据权利要求1所述的方法，其特征在于第二步骤中，

植被覆盖度计算公式如下：

。

Images