植被覆盖率envi_增强型植被指数EVI

建议大家问问题时，详细描述问题，最好是写上我怎么做了，但是存在…导致失败了。最后，多读几遍检查句子是否存在语病或歧义。

1. 植被覆盖度（Fractional Vegetation Cover, FVC）遥感反演方法

植被覆盖度定义：
一般是指植被（包括叶、茎、枝）在地面的垂直投影面积占统计区总面积的百分比。

遥感反演方法常用的包括：像元二分法、回归模型法、机器学习算法等

1.1 像元二分法：

FVC = (VI – VIsoil)/ ( VIveg – VIsoil)

式中：VI为植被指数，通常选用归一化差分植被指数，VIsoil选定的裸土端元VI值，VIveg为植被端元VI。

什么意思呢？我们来举一个简单的例子：

我们在沙漠地区选中一点，（图为MOD13Q1，NDVI值由小到大依次为：橙-黄-绿渐变），该点NDVI值为0.059，那么该值可以作为VIsoil的值，即选定的裸土端元VI值。

同样，在植被茂盛的地区

我们可以选取0.911200，该值作为VIveg的值，即选定的植被端元VI值。

现在，我们来看在专业软件（ENVI或者ArcGIS）像元二分法的操作步骤：

（1）数据源获取，所需数据：NDVI、EVI等植被指数，本文以归一化差分植被指数NDVI为例。

本文使用的数据为MODIS的MOD13Q1植被产品，产品具体介绍省略。该产品提供250m的NDVI数据，可以从LPDAAC网址中下载，也可以在GEE平台中调用。

为了方便，我在GEE平台上下载，代码如下（传不上，需要的话私信，一起交流学习）

得到的是以下结果，采用最大值合成的方法，也就是把6~8月的所有MOD13Q1 NDVI影像以最大值的方式合成为单一影像。

影像打开结果：

NDVI正常的值范围在0~1之间，但不免会有一些负值（云雪雨等，少部分土壤也有可能为负值），但这些数量非常少，可以忽略不计：

注意：

1）MOD13Q1需要进行拼接、裁剪等预处理。原始数据像元值为整型，需要乘比例系数0.0001，即NDVI值范围处于[-1,1]之间。

2）在某些软件中，加载MODIS影像会出现一片白/黑的现象，原因可能是存在背景值，或者是MOD13Q1数据统计值异常（存在着一个很大或很小的值）。

需要去除异常值的情况：

此时需要通过波段计算去除异常值：

b1 为植被指数波段，

lt （ less than）-1 指的是小于-1的值；*0表示赋值为0；

（b1 ge（greather equal ）-1 and b1 le（less equal）1）

指 b1 大于等于-1且小于等于1的值；

*b1指赋值为b1

现在，大家的影像应当是正常的影像，我们可以开始用像元二分法进行FVC反演。本文最开始，在遥感影像上随机选的两点，实际上是一种纯植被/裸土端元值选取的方法，可以直接使用这两个值。同时，还有另一种方法，即置信区间法选取 纯植被/裸土端元端元值。我们通过ENVI软件，右击影像，点击“快速统计”。该操作能让你快速了解影像像元值分布（最大值，最小值，均值等）。

置信区间法，是将影像像元值由大到小或由小到大排列。5%的置信区间，则表示排列后第5和第95分位数，在“快速统计”中为第五列。纯裸土端元像元值为：0.01344013（置信区间5%），0.00278434（置信区间2%）。

设置好纯裸土端元值，我们再看看纯植被端元值。纯植被端元像元值为：0.63147597（置信区间5%，找到第95分位数），（置信区间2%，找到98分位数）。

总结：端元值选取问题，及各自优缺点

做法1：基于置信区间法选取端元值

做法2：目视解译高空间分辨率遥感影像，手动选择土壤或植被像元，并记录其端元值（理论上可以）

做法1：选取5%置信区间，裸土端元的NDVI值应该在-0.1~0.2之间，纯植被NDVI值应该在0.8以上，裸土端元的低估会造成FVC的高估[1]，表中纯植被端元选择为0.63会造成FVC估算不准确，这和QTP植被生长季密切相关。因此，应该结合地域、季节，慎重选择端元值。

[1] Montandon L M, Small E E. The impact of soil reflectance on the quantification of the green vegetation fraction from NDVI – ScienceDirect[J]. Remote Sensing of Environment, 2008, 112( 4):1835-1845.

土壤端元：0.039

植被端元： 0.858

最后，使用波段计算计算植被覆盖度，ENVI软件中Band Math，

键入公式

输入输出路径，即可计算出区域尺度的植被覆盖度：

最后，根据FVC的定义，应当把计算出的FVC值归为0到1之间（再次利用Band Math）

后面的部分读者想知道的话请评论留言，我一定尽快更新，有问题欢迎私聊（2023.10.14）

剩下两种方法（经验回归模型和机器学习）均需要实测FVC数据。

1.2 经验回归模型法：植被指数、光谱波段线性回归。

经验回归模型法指的是FVC和某一遥感影像的波段或植被指数之间呈现一定的关系。

FVC反演值=a×NDVI+b（最简单的一次线性函数，也可以为指数、二次函数等等）

在最小二乘条件下，能得到散点唯一的条拟合直线。此时，我们可以对一张遥感影像使用Band Math：

B1为植被指数数据或者反射率数据，而0.1和0.9实际上是拟合直线的斜率以及截距。通过计算我们可以得到区域尺度的植被覆盖度数据。

1.3 机器学习方法：随机森林、支持向量机、多元自适应回归样条、高斯过程回归、神经网络、卷积神经网络等等

机器学习方法也很简单，我们需要收集实测FVC数据（表中为GLASS），以及辅助数据（列2：列7）。

在MATLAB中把实测FVC数据作为目标数据，把辅助数据作为驱动数据，使用70%（不唯一）的随机样本来训练模型，余下30%的数据用来测试模型。

最后，我们训练得到模型，并测试模型的准确性（通常为决定系数，以及均方根误差）。

在编程软件中，将遥感影像叠加，顺序需要和训练模型时的顺序一致，例如（BLUE、RED、NIR、NDVI、SAVI、EVI）。叠加后实际上是7个波段的影像集合，一个像元拥有七个波段的信息。在编程软件中，遥感影像可以看成是一个矩阵，X，Y：（行和列，也可以为坐标），”Z”信息为像元携带的值（x1~x7）。模型f（x1~x7）即可计算出区域尺度的FVC。（写的匆忙，后期再改）

1.3 机器学习方法（2024 03 11）MATLAB

思路：

（1）构建数据集（包括测试集和训练集）

例如：M~O列为驱动数据；P列为目标数据

（2）构建反演模型（以一个简单的随机森林回归模型为例）
反演模型在MATLAB里面有很多，很多CSDN博主也有很多模型的代码，只需运行出来

clear                   % 清空变量
clc                     % 清空命令行

TRAIN = xlsread('G:\\TR.xlsx');%%训练集目录
TEST= xlsread('G:\\TR.xlsx');  %%测试集目录
P_train = TRAIN(:, 2: 4)';     %%训练集读取xlsx全部的行，2~4列 （驱动数据）
T_train = TRAIN(:,1)';         %%训练集读取xlsx全部的行，1列 （目标数据）
M = size(P_train, 2);

P_test = TEST(: ,2: 4)';
T_test = TEST(: , 2)';
N = size(P_test, 2);
 %%  转置以适应模型
p_train = P_train'; p_test = P_test';
t_train = T_train'; t_test = T_test';

%%  训练模型
trees = 500;                                     % 决策树数目
leaf  = 50;                                       % 最小叶子数
OOBPrediction = 'on';                             % 打开误差图
OOBPredictorImportance = 'on';                    % 计算特征重要性
Method = 'regression';                            % 分类还是回归
net = TreeBagger(trees, p_train, t_train, 'OOBPredictorImportance', OOBPredictorImportance,...
      'Method', Method, 'OOBPrediction', OOBPrediction, 'minleaf', leaf);
importance = net.OOBPermutedPredictorDeltaError;  % 重要性

%%  模拟
t_sim1 = predict(net, p_train);
t_sim2 = predict(net, p_test );

%%  均方根误差
error1 = sqrt(sum((t_sim1' - t_train).^2) ./ M);
error2 = sqrt(sum((t_sim2' - t_test ).^2) ./ N);

%%  绘图
figure
plot(1: M, t_train, 'r-*', 1: M, t_sim1, 'b-o', 'LineWidth', 1)
legend('真实值','预测值')
xlabel('预测样本')
ylabel('预测结果')
string = {'训练集预测结果对比'; ['RMSE=' num2str(error1)]};
title(string)
xlim([1, M])
grid

figure
plot(1: N, t_test, 'r-*', 1: N, t_sim2, 'b-o', 'LineWidth', 1)
legend('真实值','预测值')
xlabel('预测样本')
ylabel('预测结果')
string = {'测试集预测结果对比';['RMSE=' num2str(error2)]};
title(string)
xlim([1, N])
grid

%%  绘制误差曲线
figure
plot(1 : trees, oobError(net), 'b-', 'LineWidth', 1)
legend('误差曲线')
xlabel('决策树数目')
ylabel('误差')
xlim([1, trees])
grid

%%  绘制特征重要性
figure
bar(importance)
legend('重要性')
xlabel('特征')
ylabel('重要性')

%%  相关指标计算
%  R2
R1 = 1 - norm(t_train - t_sim1')^2 / norm(T_train - mean(t_train))^2;
R2 = 1 - norm(t_test -  t_sim2')^2 / norm(T_test -  mean(t_test ))^2;

disp(['训练集数据的R2为：', num2str(R1)])
disp(['测试集数据的R2为：', num2str(R2)])

%  MAE
mae1 = sum(abs(t_sim1' - t_train)) ./ M;
mae2 = sum(abs(t_sim2' - t_test )) ./ N;

disp(['训练集数据的MAE为：', num2str(mae1)])
disp(['测试集数据的MAE为：', num2str(mae2)])

%  MBE
mbe1 = sum(t_sim1' - t_train) ./ M ;
mbe2 = sum(t_sim2' - t_test ) ./ N ;

disp(['训练集数据的MBE为：', num2str(mbe1)])
disp(['测试集数据的MBE为：', num2str(mbe2)])

完毕后，用save model保存模型即可。

（3）制作空间分布图。

先加载模型

放入波段数据，在这里为了便于分析，我将一幅影像的波段重命名为：Blue（2001.tif）、Green（2002.tif）、Red（2003.tif）、NIR（2004.tif），将其依次读取并做成一个有空间位置信息的大矩阵。

Blue（2001.tif）、Green（2002.tif）、Red（2003.tif）、NIR（2004.tif）为4个驱动数据，对应一个目标数据（反演结果，也是FVC）。

我感觉可以这么想：4张影像（例如Blue、Green、Red、NIR）叠加。影像大小200*200

那么现在是200行*200列*4

将其平铺为40000（行）*4（列），

取出第1行（data1）放入随机森林回归模型predict（net，data）

随机森林回归模型（Blue、Green、Red、NIR）——>预测FVC值，这个值就是第一行第一列预测的像元值。逐个预测，依次填回影像。最后可以得到FVC机器学习反演结果，带经纬度信息。

2. 植被覆盖度产品数据

GLASS、GEOV1、GEOV2、GEOV3、MuSyQ FVC使用及处理教程

GEOV3全球数据很大，但可以定制区域，详细过程目前还没时间写。

2.1 MuSyQ FVC产品数据

来源：国家地球系统科学数据中心首页申请后通过邮箱FTP获取

参考文献：

Zhao J , Li J , Liu Q , et al., 2020. Estimating fractional vegetation cover from leaf area index and clumping index based on the gap probability theory. International Journal of Applied Earth Observation and Geoinformation 90:102112.

数据时空分辨率为4d/500m，时间范围：2001~2019，空间范围：全球

HDF5(H5) 可以通过MATLAB、Python打开，并转化为tif

2.2 GEOV3、GEOV2、GEOV1

来源：哥白尼全球土地服务中心(https://land.copernicus.eu/global/products/fcover)

参考文献：

Fuster, B.; Sánchez-Zapero, J.; Camacho, F.; García-Santos, V.; Verger, A.; Lacaze, R.; Weiss, M.; Baret, F.; Smets, B. Quality Assessment of PROBA-V LAI, FAPAR and FCOVER Collection 300 m Products of Copernicus Global Land Service. Remote Sens. 2020, 12, doi:10.3390/rs12061017.

Baret, F.; Weiss, M.; Lacaze, R.; Camacho, F.; Makhmara, H.; Pacholcyzk, P.; Smets, B. GEOV1: LAI and FAPAR Essential Climate Variables and FCOVER Global Time Series Capitalizing over Existing Products. Part1: Principles of Development and Production. Remote Sens. Environ. 2013, 137, 299–309, doi:10.1016/j.rse.2012.12.027.

Verger, A.; Baret, F.; Weiss, M. Biophysical Variables from VEGETATION-P Data. MultiTemp 2013 7th Int. Work. Anal. Multi-temporal Remote Sens. Images. IEEE, 2013, 10–13.

3. 植被覆盖度动态变化分析

3.1 趋势分析

所需数据：研究区时间序列影像，下面补充一些常用的影像处理知识。

打开建模工具

插入循环工具，循环栅格（Rasters）

将工具拖入模型内

放入栅格：选择的一定是循环文件

掩膜文件：可以是矢量，也可以是栅格（研究区/感兴趣区域）

%name%表示输出文件的名字与放入文件的名字一模一样。

点击开始运行即可

得到了23张研究区的时间序列影像。

3.1.1 一元线性回归

3.1.2 Sen+MK

需要的数据：研究区的时间序列影像，所有影像行列必须一致，坐标系统，空间分辨率等基本信息完全一致。

放入MATLAB里，注意修改参数信息，（代码C站里有很多）。

Sen趋势计算结果：

最后再重分类：

先看到MODIS_Sen.tif的像元值分布

栅格值小于0的（斜率小于0）的栅格赋值为1；栅格值大于0的（斜率大于0）的栅格赋值为2.

重分类结果：1表示该区域FVC有减小趋势，2表示该区域FVC有增加趋势。

3.2 相关分析，偏相关分析以及复相关分析

例如：探究A区域植被覆盖度变化与温度、降水之间的关系。

相关分析：放入2000~2023年FVC时间序列影像及2000~2023年平均气温等数据，计算出的结果为相关系数。

详细可以参考：（2024.3.14）

生态参数反演（植被覆盖度FVC）手把手教会

3.3 基于地理探测器，探究影影响FVC动态变化的因子

3.3.1 创建渔网

（1）选择栅格

DEM数据完整度高，栅格像元基本无缺失值，故选取DEM数据创建渔网。

DEM数据采用的是Copernicus DEM。

在该网站有其详细介绍：https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_DEM_GLO30

The Copernicus DEM is a Digital Surface Model (DSM) which represents the surface of the Earth including buildings, infrastructure and vegetation. This DEM is derived from an edited DSM named WorldDEM&trade, i.e. flattening of water bodies and consistent flow of rivers has been included. Editing of shore- and coastlines, special features such as airports and implausible terrain structures has also been applied.

The WorldDEM product is based on the radar satellite data acquired during the TanDEM-X Mission, which is funded by a Public Private Partnership between the German State, represented by the German Aerospace Centre (DLR) and Airbus Defence and Space. More details are available in the dataset documentation.

Earth Engine asset has been ingested from the DGED files.

DEM数据一般需要分块下载，GEE平台可以完成对这些块块的拼接与裁剪。这里给出GIS中如何拼接（镶嵌）。

例如：当前有两块区域。