多因子影响下洱海流域NDVI时空变化规律探究

郑舒元，李亚强，董亚坤，海 燕，王建雄

1.云南农业大学 水利学院，云南 昆明 650201；

2.云南省高校农业遥感与精准农业工程研究中心，云南 昆明 650201

0 引 言

洱海作为云南省九大高原湖泊之一，其流域为生态系统的重要组成部分。洱海流域生态环境变化对环境整体有着重要影响[1]，而植被作为流域生态系统的重要组成部分，生长状态及覆盖度变化能够切实反映生态系统的变化。植被变化与光照、温度、降雨量、土壤状况等自然环境要素密切相关[2-4]，与不同自然环境要素间并不是简单的单一对应关系。不同要素对植被变化影响程度不同，在不同程度的共同作用下，植被生长状态与植被覆盖会发生不同变化。

归一化植被指数（NDVI）能够反映植被生长状态及植被覆盖，是反映转化生长状态的优良指标，原理简单易操作，是目前应用最广泛的植被指数。国内外许多作者利用不同因子对NDVI进行拟合，探究不同因子作用下NDVI的变化规律。任荣仪等利用MODIS数据，对贵州省NDVI与温度、降水数据进行分析[5]，发现贵州省2000-2018年NDVI受温度影响强于降水；邓玉娇等对广东省NDVI与气温、降水、日照等气象因子进行分析，发现广东省2000-2018年月平均NDVI与气温、降水、日照时数相关性显著[6]；李元春等对甘南和川西北地区草地植被NDVI变化与人类活动、气象因素进行相关性分析，发现2000-2018年温度和降水对研究区草地NDVI主要呈正向驱动，人类活动在2000-2015年对草地NDVI呈负向影响，2006-2017年呈正向影响[7]；曹俊涛等对三亚市NDVI与降水量、日照时数进行相关性分析，发现2010-2019年三亚市NDVI与降水量整体呈正相关，与日照时数呈负相关[8]。然而以上研究多侧重于NDVI与少数因子之间的分析，多种因素对NDVI的影响尚未深入进行探究。因此，本文利用2010-2015年Landsat TM/OLI影像与洱海流域月均温（℃）、月降雨量（mm）、月风速（m/s）、月蒸发量（mm）、月海平面气压（hPa）、月总云量（tcc）、月日照强度（J/m2.d）、月紫外强度（J/m2.d）、月短波辐射（W/m2）、月长波辐射（W/m2）、月空气湿度及月40 cm土壤湿度（kg/m2）数据，采用主成分分析法（Principal Component Analysis，PCA）对各因子与NDVI的相关性进行分析，并分别采用BP神经网络与偏最小二乘对NDVI进行反演，取拟合效果较好的模型对洱海流域NDVI变化进行预测。

1 研究区概况

研究区域位于云南省中部偏西洱海流域，地处北纬25°25′～26°10′，东经99°32′～100°27′。流域面积2607.75 km2，湖面高程196 6m时湖面面积为252 km2，一般湖水面积约246 km2，蓄水量约29.5×108m3，呈狭长形。流域内河流呈现聚集状汇入中心湖泊，境内有大小河溪共117条，主要入湖沟渠125条，辖大理市、洱源县16个乡镇，具有供水、灌溉、发电、调节气候、渔业航运及旅游七大功能，是云南省九大高原湖泊之一。流域内干湿季分明，5-10月为雨季，11月至次年4月为旱季，多年年均温度及降雨量分别为15℃和 1048 mm。

2 数据来源及研究方法

2.1 数据来源

遥感影像来自Google_earth_engine，时间跨度为2010年1月-2015年12月的Landsat TM/OLI数据，时间分辨率为1个月，空间分辨率为30 m，并进行辐射定标、大气校正、去云处理、NDVI合成。2010-2011年影像为Landsat5 TM数据，2012年数据为Landsat 7 ETM+，2013-2015年采用Landsat 8 OLI日间成像数据。Landsat 7数据自2003年6月以来，由于扫描线矫正器（SCL）故障导致Landsat 7影像数据出现空白区域，下载的影像缺失区域呈条带状，利用ENVI的Landsat gapfill插件，使用三角插值方法对2012年的影像数据进行填充修复。洱海流域月均温、月降雨量、月风速、月蒸发量、月海平面气压、月总云量、月日照强度、月紫外强度、月短波辐射、月长波辐射、月空气湿度、月40 cm土壤湿度数据来自美国国家气候数据中心（NCDC）的公开FTP服务器（ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-lite/）与中国气象数据网（http://data.cma.cn/）。

2.2 研究方法

2.2.1 一元线性回归法

一元线性回归法在植被长时序列变化趋势的研究中较为常用，时间x与NDVI数据y的线性数学模型为：

式中，k为变化趋势斜率；b为未知常数ε为随机误差。利用NDVI值(xi，yi) (i=1，2，，11)求解未知参数k的公式为：

式中，k为NDVI值长时间序列变化趋势，k＞0表示植被呈上升趋势，k＜0则表示植被为减少趋势；n为长时间序列数。

2.2.2 主成分分析法

主成分分析（Principal Com-ponent Analysis，PCA）是一种将原来众多具有一定相关性的变量，通过线性组合转换成一组互不相关的综合变量，用几个综合变量代替原来众多变量的一种多元统计分析方法[9-10]，数学模型如下：

式中，α11，……，αpn为X的协方差阵特征值对应的特征向量；ZX1，……，ZXP为原变量经过标准化处理后的值。

2.2.3 BP神经网络

BP（back propagation）神经网络是1986年由Rumelhart和McClelland为首的科学家提出的概念，是一种按照误差逆向传播算法训练的多层前馈神经网络，是应用最广泛的神经网络。在输入层与输出层之间增加若干层神经元，这些神经元称为隐单元，它们与外界没有直接的联系，但其状态的改变，能影响输入与输出之间的关系，每一层可以有若干个节点。

BP神经网络的学习过程分为两个部分，即信号的正向传播和误差的反向传播。正向传播是信号由输入层传输到隐含层，经隐含层处理单元计算后传向输出层；误差反向传播是指当实际输出信号与期望输出信号不符时，误差信号就会沿着误差减小的方向由输出层经隐含层向输入层传输，将误差分摊给各层所有单元，进而修正输入层与隐含层、隐含层与输出层之间的权值[11-15]。

在正向传播过程中计算所有节点的误差平方和，得：

式中，m为输出层的节点数；dj为输出层节点j期望输出；yj为输出层节点j的实际输出。如果E小于规定值，则进行下一个样本训练，否则需进行误差逆向传播，并调整权值的大小。

如此反复训练，使网络输出的误差在达到预先设定的学习次数之前减少到可接受的程度。这样的过程也称作BP神经网络的训练（图1）。

图1 BP神经网络的学习过程Fig.1 Learning process of BP neural network

基于BP神经网络的洱海流域NDVI变化，利用研究时间段内经过标准化处理的12个影响因子数据作为输入参数，NDVI值作为仿真输出。隐含层节点对神经网络预测效率有很大影响，若隐含层节点太少，则容错性较差，无法得出正确成果，造成欠拟合；隐含层节点数过多则会增加学习时间，影响学习效率，出现过拟合现象。隐含层节点数通常根据经验公式计算，常用公式为：

式中，m为输入层节点数；n为输出层节点数；s为隐含层节点数；k为1～10间的整数。在模型中m=12，n=1，经过3～5的计算，隐含层节点数为6，经过3～6计算，隐含层节点数为4～13，通过实际计算对比，确认隐含层节点数为4的时候收敛速度最快，因此建立12-4-1的BP神经网络。

2.2.4 偏最小二乘法

偏最小二乘回归分析是一种通过使误差平方和最小化，找到数据的最佳匹配函数的回归模型，这种模型提供一种多对多线性回归建模的方法。其原理为：设有q个因变量{y1，y2，……，yq}与p个自变量{x1，x2，……，xp}，为研究因变量和自变量的统计关系，观测n个样本点，由此构成自变量与因变量的数据集。偏最小二乘回归分别从数据集中提取成分t1和u1，t1和u1应尽可能好地代表数据表X和Y，同时自变量的成分t1对因变量的成分u1又有很强的解释能力。在第一个成分t1和u1被提取后，偏最小二乘回归分别实施X对t1的回归以及Y对u1的回归，若回归方程已经达到满意的精度，则算法中止；否则，将利用X被t1解释后的残余信息以及Y被u1解释后的残余信息进行第二轮的成分提取，经过一系列重复学习，直到精度达到一个较为满意的值，最终X共提取了m个成分t1，t2，……，tm。偏最小二乘回归将利用yk（k=1，2，……，q）对t1，t2，……，tm的回归构成因变量yk对自变量x1，x2，……，xp的回归方程[16-22]。建模整体可归纳为以下4个步骤。

1）数据标准化处理。为了加速权重参数的收敛，消除特征之间的差异性，样本点的集合重心与坐标原点重合，经过标准化的X数据矩阵记为E0=（E01，……，E0p）n×p经过标准化的Y的数据矩阵记为F0=（F01，……，F0q）n×q。

2）从E0中提取第一成分t1。t1为x1，x2，……，xp的线性组合，必须尽可能代表解释变量矩阵E0，使t1方差最大化，计算公式如下：

式中，t1为从E0中提取的第一个主成分；w1为E0的第一个主轴，且向量的模为1；r(xj,y)为解释变量与被解释变量之间的相关系数。

3）进行E0对t1的回归，得到两变量的回归系数与方程如下：

式中，P1为E0对第一个成分t1的回归系数；E1为回归方差的残差矩阵。而后用残差矩阵E1重复步骤2）、3），以此来提取r个主成分t1，t2，，tm，直到不满足交叉有效性检验标准为止。交叉有效性是新增成分tr，代表对模型预测效果是否有显著改善，计算公式为：

式中，yi为变量y的第i个样本点的真实值；为采用全部样本点对k个成分建立回归方程后第i个样本点的预测值；k(-i)为选取k个成分时将第i个样本点带入排除第i个样本点回归方程的预测值；为交叉有效性，当≥0.0975时，新增成分tm对预测模型有显著改善。

4）建立F0与t1，t2，……，tm的回归方程，与的标准化变量回归方程分别为：

2.2.5 验证方法

BP神经网络与偏最小二乘拟合成果的良好程度衡量指标有不同方法，通常采用计算累加差的方式对比累加差大小或计算拟合优度R2，但无法同时对误差的敛散程度以及拟合优劣进行反映。为此，本文提出了一种基于比值的反演效果分析法（Analysis of prediction results based on ratio，ABR），即将拟合值与真值视为一组正交向量。如图2所示，在坐标系里以y轴为真值，x轴为拟合值，真值与拟合值的比值k能够表示拟合度。若拟合结果等于真值，则；若拟合结果大于真值，则；若拟合结果小于真值，则。图中，直线y0表示拟合结果等于真实值，y1表示拟合结果小于真实值，y2表示拟合结果大于真实值，即k1＞k0＞k2。

图2 不同拟合成果与真实值对比Fig.2 Comparison of different fitting results and real values

若φ越大，则说明拟合数据离散程度越大，拟合效果偏差较大；若φ越小，则说明拟合效果偏差较小，拟合成果越好。

利用ABR分析法解算φ，对其准确率进行计算，采用平均相对误差绝对值δ与R2对评价结果进行验证。δ在评价精度误差分析时是一个重要指标，能够计算预测结果的相对误差率，同时得到各点的相对误差绝对值，计算公式如下：

式中，n为数据个数；y为对应的因变量，将预测得到的NDVI值与原始数据带入进行解算即可得到其平均相对误差绝对值。R2用来解释yi变量的方差比例大小，用于衡量预测值对于真实值的拟合好坏程度，R2取值区间为(-∞，1]，值越接近1表示拟合程度越好，R2计算公式如下：

式中，yi'为第i个月的NDVI预测值；yi为第i个月NDVI真值；为i个月NDVI真值的均值。

3 结果与分析

3.1 洱海流域NDVI变化分析

3.1.1 洱海流域NDVI变化规律分析

洱海流域NDVI月际变化能够直观反映出1年内当地植被覆盖状况与生长状况，2010-2015年洱海流域NDVI值变化如图3所示。

图3 2010-2015年洱海流域NDVI变化趋势Fig.3 Variation trend of NDVI in Erhai lake basin from 2010 to 2015

由目视解译可知，洱海流域2010-2015年NDVI月际变化呈现一定规律，总体趋势为先降后升。为探究空间变化规律及影响空间变化规律的因子，将洱海流域2010-2015年NDVI数据利用一元线性回归法计算，结果见表1和图4。

表1 一元线性回归法分析植被变化结果Tab.1 Analysis results of vegetation change by single variable linear regression method

图4 一元线性回归法植被空间变化分布Fig.4 Spatial change distribution of vegetation by linear regression method

根据合成结果，2010-2015年洱海流域绝大部分植被呈现明显改善趋势，洱海西岸、洱海南部和东南部城镇区域、洱海流域西北部植被整体呈现退化趋势。通过与原始影像及行政区划图进行对比，推断植被覆盖变化是因为洱海流域大理市、大理镇、银桥镇、喜洲镇、上关镇对洱海沿岸土地开发利用及作物种植造成的，表明人类活动对植被覆盖有着直接影响。

通过对2010-2015年每个月的NDVI均值进行拟合可得到如图5所示的曲线，发现1-2月趋于平稳，3-4月有下降趋势，5-7月处于较低水平，8月开始上升，9-12月维持在高水平，呈现一定规律性波动。

图5 洱海流域NDVI月合成拟合曲线Fig.5 Monthly synthetic fitting curve of NDVI in Erhai lake basin

洱海流域NDVI月变化趋势可用傅里叶级数进行拟合，变化趋势与正弦函数相似，拟合函数为：

式中，x为月份；y为当月NDVI值；R2=0.8931，拟合效果良好。为验证上述NDVI月变化规律与实际变化情况的一致性，在此用2010-2015年洱海流域6年72个月的NDVI数据再次进行拟合（图6）。

图6 洱海流域NDVI逐月份拟合曲线Fig.6 Monthly fitting curve of NDVI in Erhai lake basin

如图6所示，利用傅里叶级数对2010-2015年NDVI数据进行拟合，所得拟合函数为：

式中，x为月份；y为当月NDVI值；R2=0.5266＞0.5，所得拟合曲线能够代表其变化规律。由以上结果分析可得，洱海流域NDVI值从12月至次年5月呈现下降趋势，6-11月呈现上升，拟合曲线谷值出现在5月左右，峰值出现在11月左右。

3.1.2 洱海流域NDVI影响因子分析

影响洱海流域NDVI的因素多且复杂，但有大小之分，确定不同因子的影响大小对变化分析与变化规律反演有重要意义。由于选用12种不同的影响因子，可使用PCA法进行降维分析，得出的主成分系数能够解释不同因子的影响大小，通过对12个因子6年的逐月数据降维分析，所得主成分两个，主成分表达式为：

对两组主成分系数进行可视化，分析不同因子对洱海流域NDVI值的影响大小（图7）。可见气温、海平面气压、紫外强度与短波辐射为影响洱海流域NDVI变化的主导因子，其中气温、紫外强度与短波辐射为正向影响因子，海平面气压为负向影响因子。

图7 2010-2015年洱海流域NDVI变化影响主导因素Fig.7 Main factors influencing NDVI change in Erhai lake basin from 2010 to 2015

3.2 洱海流域NDVI变化预测

3.2.1 基于BP神经网络的洱海流域NDVI变化预测

利用MATLAB2016b进行网络仿真，取(-1，1)之间的随机数初始化网络权值w，阈值θ，输入学习速度γ以及期望误差ε，输入标准化数据α，数据矩阵T，其中训练样本占总样本70%，验证样本占总样本15%，测试样本占总样本15%。通过预测模型进行多次训练，直到训练结果与实际值有较高符合度，再对预测模型输入测试样品数据可得到预测结果。

如图8所示，训练网络预测值均处于较高水平，预测结果是较为准确的，网络性能也处于良好水平，BP神经网络拟合成果与原始数据变化规律呈现高度相似，可认为拟合结果符合要求，可对洱海流域NDVI变化做出预测。通过将测试样品带入预测模型，并与期望值相比较的可视化成果如图8所示。

图8 BP神经网络预测结果与原始数据比对图Fig.8 Comparison of BP neural network prediction results and original data

3.2.2 基于偏最小二乘的洱海流域NDVI变化预测

根据PLS建模流程，对原始数据进行标准化，利用2010-2015年自变量数据与因变量数据构建数据矩阵，提取主成分并计算回归系数从而构建回归方程，通过以上流程完成整个OLS预测模型的构建。结果显示提取两个主成分即可满足回归需要，回归方程由一组回归系数与一个常数项构成，通过计算所得回归方程如下：

利用构建好的偏最小二乘模型，将2010-2015年12个影响因子逐月数据带入预测模型，对2010-2015年NDVI值逐月预测检验，将所得结果可视化（图9）。

图9 PLS预测结果与原始数据比对图Fig.9 Comparison between PLS prediction results and original data

由目视解译可得PLS预测成果与原始数据变化规律高度相似，拟合效果较为平滑，但存在部分预测值与原始数据差值过大的问题。

3.2.3 BP神经网络预测效果与PLS预测效果对比

数据解算完毕，利用ABR分析法对两种不同方法的预测效果进行评价。BP神经网络预测结果评估如图10所示，偏最小二乘预测结果评估如图11所示。

图10 BP神经网络预测结果评估Fig.10 BP neural network prediction result evaluation

图11 偏最小二乘预测结果评估Fig.11 Evaluation of partial least squares prediction results

通过对图10、图11目视解译，可直观得出BP神经网络拟合结果偏差量整体小于偏最小二乘。根据3个描述性统计量计算结果，对两种方法预测结果的评估见表2，2010-2015年洱海流域NDVI值BP神经网络预测效果整体优于偏最小二乘预测效果，使用BP神经网络模型能够更准确地预测洱海流域NDVI值的变化。

表2 BP神经网络预测与偏最小二乘预测效果对比Tab.2 Comparison of BP neural network prediction and partial least squares prediction

4 结 论

以2010-2015年的Landsat TM/OLI为研究数据，对月变化趋势进行研究，运用主成分分析法分析不同因子对洱海流域NDVI值变化的影响大小，再利用BP神经网络与偏最小二乘回归法对洱海流域的NDVI值进行预测，结论如下。

1）洱海流域2010-2015年NDVI值变化呈现傅里叶级数波状变化，6-10月植被长势最好，11月NDVI值达到最大。

2）2010-2015年洱海流域绝大部分植被呈现明显改善趋势，洱海西岸、洱海南部和东南部城镇区域、洱海流域西北部植被整体呈现退化趋势，植被严重退化土地面积为126.23 km2，轻微退化面积为133.35 km2，分别占总流域面积的4.81%与5.09%。

3）气温、海平面气压、紫外强度与短波辐射对洱海流域NDVI变化影响较大，其中气温、紫外强度与短波辐射对NDVI值起到较大的正向推动作用，海平面气压对NDVI值起到负增长的作用。

4）实验结果表明，BP神经网络预测模型比偏最小二乘预测模型平均相对误差高出2.44%，R2高出0.3123，准确率∅高出0.1275，BP神经网络对洱海流域2010-2015年NDVI值变化预测效果更优。