ENVI学习总结（十二）——基于改进的 CASA 模型反演 NPP

远方~ · 发表于 2020-8-24 15:14

基于改进的 CASA 模型反演 NPP

数据为地理空间数据云提供的 TM 影像，其具体信息如图 1 所示。该实验中所使用的 TM 数据已进行过 QUAC 快速大气校正，来消除大气和光照等因素对地物反射的影响。Landsat 主题成像仪(TM)是 Landsat4 和 Landsat5 携带的传感器，每 16 天扫描同一地区，即 16 天覆盖全球一次。 TM 影像包含 7 个波段，波段 1-5 和波段 7 的空间分辨率为 30 米，波段 6（热红外波段）的空间分辨率为 120 米

软件安装

将文件夹“程序”下的 zhu_npp_calc.sav 文件（朱文泉教授设计的NPP模型）复制到 ENVI 安装目录的 save_add 文件夹下面；用写字板打开 ENVI 安装目录下的 envi.men 文件，在最后面加上如下内容，然后保存退出；

启动 ENVI Classic 经典版本之后，将会在 ENVI 菜单栏出现“实用函数”；

单击实用函数，选择NPP估算 ,弹出植被 NPP 计算设置窗口

一、数据准备

在软件运行之前，需要准备以下数据月平均温度栅格文件，单位为℃，由气象数据插值得到，时间范围与 NDVI 一致。月总降水量栅格文件，单位为 mm，由气象数据插值得到，时间范围与 NDVI 一致。

月太阳总辐射

栅格文件：单位为 MJ/m2，由气象数据插值得到，时间范围与 NDVI 一致。
NDVI 时间序列数据：栅格文件，由遥感数据计算得到。可以是一个时间序列，如：12 个月的 NDVI 数据。
植被类型图：栅格文件，确定各植被类型的空间分布。静态参数文件该文件配置各植被类型的 NDVI 最大值、NDVI 最小值、SR 最大值、SR 最小值以及最大光能利用率（gC/MJ）。

二、设置结果文件存放路径

设置结果存放文件夹，所有结果都会输出到该文件夹下面。默认情况下，会输出年度植被净初级生产力（npp_sum）和植被年平均覆盖率（veg_cov_mean），如果勾选了“同时输出 NPP 及植被覆盖度时间序列文件”复选框，则还会额外输出每个月的植被 NPP及植被覆盖度。

三、实验过程

1.数据来源及说明

遥感数据：来源于中国科学院资源环境科学数据中心(http://www.resdc.cn)所提供的 2001 年的 SPOT-VEGETATION NDVI 影像。该数据集经过最大值合成法 MVC（Maximum Value Composites）处理，空间分辨率为 1km，时间分辨率为逐月。
气象数据：来自中国气象科学数据共享服务网（http://cdc.cma.gov.cn/），包括 2001 年的月降水量、月平均气温和月总太阳辐射数据，共覆盖东北地区的 107 个气象站点。为保持气象和遥感数据在空间上的一致性，利用插值工具将点数据转换为空间分辨率为 1km 的面栅格数据。文中所有数据均使用以 WGS 84 为基准面的 Albers 等面积圆锥投影。
植被类型图：植被图为 2001 年中国科学院植被图编辑委员会编制的 1:100 万中国植被类型图。

2.数据处理

（1）月平均温度空间插值（克里金插值法必做）启动 ArcGIS，并加载气象（温度和降水）站点数据和东北地区边界数据，在 ArcGIS 中加载温度数据（excel 格式）。气象（温度和降水）站点数据属性表与温度数据（excel 格式）属性表建立连接。在内容列表中，右键单击图层名 NE_station.shp→打开属性表（Open Attribute Table），气象（温度和降水）站点数据属性表如图所示。

右键单击温度数据 Sheet1→打开属性表（Open），温度数据属性表如图所示

温度 excel 表中的“ID”与气象站点属性表中的站点代码“ID”相对应（实际操作中，这一字段的名称可以不同，但数据类型必须相同）。在内容列表中，右键单击想要连接的图层名 NE_station.shp→Joins and Relates（连接和关联）→Join（连接），弹出“连接数据”对话框。在连接数据对话框中，先选择用作连接依据的字段 ID；然后选择要连接到图层的表 Sheet1$；最后选择 excel 表中要用作连接依据的字段 ID，其他参数使用缺省值，单击 OK，如图所示。查看新属性表，发现增加了 2001 年逐月平均温度数据，如图所示。如果想要将连接的数据与地理要素永久保存在一起，可将数据导出到一个新要素类中，在内容列表中右键单击该图层→data→Export data，设置存储路径，单击 OK。

在执行插值处理之前，先查看菜单栏→Geoprocessing→Geoprocessing options中最后一个选项框 Results are temporary by default（创建临时结果）是否被勾选。本实验中需要保留插值结果，故不选中 Results are temporary by default 复选框。在 ArcToolBox 中，单击 Spatial Analyst Tools→Interpolation→Kriging→Batch，打开该工具的窗口。（ArcGIS 汉化版本，Spatial Analyst→插值→克里金插值→批处理）如图 15 所示，在插值窗口中，Input point features（输入点要素）中选择被用来进行插值的离散点数据 NE_station；在 Z value field（Z 值字段）中选择要参与插值的字段 Sheet1$.tem01；在 Output surface raster（输出表面栅格）中设置输出文件名称，建议按月份顺序命名（命名规则：tem2001MM.tif，tem 表示月平均温度， MM 表示月。）； Semivariogram properties（半变异函数属性）中使用缺省值 Ordinary（普通克里金法）和 Spherical（球面半变异函数）；在 Output cell size（输出像元大小）中设置输出表面的栅格大小为“1000”（和遥感影像像元大小相一致）；Search radius（搜索半径）和 Output variance of prediction raste（预测栅格数据的误差）使用缺省值。

单击右侧的“+”增加 11 行用来输入其它 11 个月的数据。在 Z value field（Z 值字段）中选择要参与插值的字段（按序选择对应的月平均温度字段）；在 Output surface raster （输出表面栅格）中设置输出文件名称。其他参数设置与第一行相同。如果某一列值相同，可直接单击该列空白处，在菜单中选择 Fill，完成一列的参数设置。在插值窗口中，单击 Environments …（环境设置）→Processing Extent（处理范围），在 Extent 的下拉菜单中选择 Same as layer 东北地区边界，单击 OK，如图所示。

2001 年 1 月东北地区平均温度插值结果如图所示。

启动 ENVI 软件，选择 File→Open，打开 12 个月的温度插值数据。在 Toolbox 工具箱中，双击 Raster Management→Layer stacking 进行波段叠加，打开 Layer Stacking Parameters 窗口。单击 Import File…按钮，弹出 Layer Stacking Input File 窗口，选中 12 个月的温度插值数据，单击 OK。输入的文件将出现在 Selected Files for Layers Stacking 列表中，如果不一致可通过 Reorder Files…按钮来调整波段顺序。数据加载进来后，会自动读出图像投影信息和像元大小。Resampling（重采样）方法使用缺省值，即 Nearest Neighbor。单击 Inclusive 和 Exclusive，选择输出文件范围。如果选择 Inclusive，输出图像的地理范围将是所有输入文件范围的并集；如果选择 Exclusive，输出图像的地理范围仅包含所有输入文件的重叠范围。此处选用缺省值 Inclusive。设置文件名及存储路径，单击 OK。波段叠加后，band1 表示 1 月的气温数据，band2 表示 2 月的气温数据。

在 Layer Manager（图层管理）中，右键单击图层名 tem2001，选择 Profiles→Spectral，可以看到某一点像元的时间序列曲线，如图所示。横坐标为月份，纵坐标为温度。

月总降水量空间插值使用克里金插值方法将月总降水量数据转换为面栅格数据，操作步骤同上。在 Layer Manager（图层管理）中，右键单击图层名 rain2001，单击 Quick Statics （Compute Statics）即可弹出统计界面。快速统计可以快速的统计图像的最大值、最小值、均值、标准差和直方图分布。如图所示，波段 12，即经插值后的 12 月份月总降水量数据最小值出现负数，这是不合理的。使用 Band Math 工具，把图像中小于 0 的值赋值为 0。

在 Toolbox 工具箱中，双击 Band Algebra→Band Math 工具，打开 Band Math 对话框

在 Enter an expression（运算表达式输入框）中输入表达式：( b1 ge 0 ) * b1 + ( b1 lt 0 ) * 0 ，关系运算符対真值（关系成立）返回值为 1，对假值（关系不成立）返回值为 0。第一个关系运算符( b1 ge 0 )是找出像元值为正或 0 的像元，乘以它们的初始值；第二个关系运算符( b1 lt 0 )是找出像元值为负的像元，乘以 0。单击 Add to List 按钮，将表达式添加到列表中。单击 OK 按钮，打开 Variables to Bands Pairings 对话框，为运算表达式中各个变量赋予图像文件或者图像波段。单击 Map Variable to Input File 按钮，为变量 b1 选择 rain2001 图像的所有波段。单击 Choose 按钮，选择文件名及路径保存结果，单击 OK 按钮，执行运算。

对 rain2001_ge0 文件进行快速统计，如图 23 所示，图像所有波段的像元最小值都大于等于 0，降水插值数据符合真实情况，可进行后续实验。

（3）月太阳总辐射空间插值

由于辐射站点较少，使用克里金方法进行插值精度不高，对月太阳总辐射数据进行插值时使用反距离加权插值方法（IDW）。其操作步骤与克里金插值方法类似，在 Output cell size（输出像元大小）中设置输出表面的栅格大小为 “1000”（和遥感影像像元大小相一致）；在 Power（幂文本框）中输入 IDW 的幂值，幂值是一个正实数，其缺省值为 2；Search radius（搜索半径）为缺省值，VARIABLE 12；Input barrier polyline features（输入障碍折线要素）为可选项。

（4）NDVI 时间序列

2001 年 1 月东北地区太阳总辐射插值结果如图 25 所示。在 ENVI 中打开 2001 年 NDVI 时间序列数据时，由于背景值为-9999，像元值范围在[0,1]（NDVI 值范围在[-1,1]，该数据已经去除由于云、积雪等影响而出现的负值），数值之间差距太大，图像呈现为黑色。在 Toolbox 工具箱中，双击 Raster Management→Edit ENVI Header 工具，在 Data Ignore Value 文本框中填入-9999，忽略背景值影响。在 Layer Manager（图层管理）中，右键单击图层名 northeast_ndvi_2001，选择 Profiles →Spectral，可以看到某一点 NDVI 时间序列曲线，如图 26 所示。横坐标为月份，纵坐标为 NDVI 值。

（5）植被类型图根据实验需求，在 2001 年中国科学院植被图编辑委员会编制的 1:100 万中国植被类型图基础上对各植被类型重新筛选、合并和编码等处理。在 ENVI 中打开经处理后的东北地区植被类型图，如图所示。植被类型主要有：针叶林、阔叶林、针阔混交林、灌丛、草地、栽培植被、沼泽、荒漠和非植被。

（6）静态参数文件生成

引用朱文泉教授研究结果，配置 9 类植被类型的 NDVImax、NDVImin、SRmax、SRmin 和 Emax（理想状态下最大光能利用率）参数。其中 NDVImax 和 SRmax 的计算需要东北地区植被类型图和 NDVI 时间序列最大值数据。

①NDVI 时间序列最大值计算方法：在 ENVI 中打开 NDVI 时间序列数据，在 Toolbox 工具箱中，双击 Band Algebra→Band Math 工具，打开 Band Math 对话框。在 Enter an expression（运算表达式输入框）中输入表达式：b1>b2>b3>b4>b5>b6>b7>b8>b9>b10>b11>b12 。单击 Add to List 按钮，将表达式添加到列表中。单击 OK 按钮，打开 Variables to Bands Pairings 对话框，为运算表达式中各个变量赋予图像文件或者图像波段。在 Variables used in expression 列表框中选择变量 b1，在 Available Bands List 中为 b1 指定一个波段。利用同样的方法分别为所有变量指定波段。单击 Choose 按钮，选择文件名及路径保存结果，单击 OK 按钮，执行运算。

b1：选择 northeast_ndvi_2001 文件的第 1 个波段
b2：选择 northeast_ndvi_2001 文件的第 2 个波段
b3：选择 northeast_ndvi_2001 文件的第 3 个波段
b4：选择 northeast_ndvi_2001 文件的第 4 个波段
b5：选择 northeast_ndvi_2001 文件的第 5 个波段
b6：选择 northeast_ndvi_2001 文件的第 6 个波段
b7：选择 northeast_ndvi_2001 文件的第 7 个波段
b8：选择 northeast_ndvi_2001 文件的第 8 个波段
b9：选择 northeast_ndvi_2001 文件的第 9 个波段
b10：选择 northeast_ndvi_2001 文件的第 10 个波段
b11：选择 northeast_ndvi_2001 文件的第 11 个波段
b12：选择 northeast_ndvi_2001 文件的第 12 个波段

② 启动 ENVI Classic 经典版本，单击实用函数→NPP 估算，在植被 NPP 计算设置窗口中，选择配置静态参数按钮。在 npp 及植被覆盖度计算静态参数设置窗口中，选择计算 NDVImax，SRmax 按钮，依次选择东北地区植被类型图和 NDVI 时间序列最大值。在输入分类精度窗口中，输入 70，单击确定。得到 NDVImax 和 SRmax 结果，如图所示。

引用朱文泉研究结果，手动输入 NDVImin、SRmin 和 Emax 数据，如表 1 所示。完成后单击导出数据按钮，设置文件名（命名规则：静态参数设置 2001）和存放路径，生成 cfg 格式的静态参数文件。

NPP 反演

启动 ENVI Classic 经典版本，单击实用函数→NPP 估算，在植被 NPP 计算设置窗口中，依次选入准备好的数据，单击完成，如图所示。

程序运行完成后，会弹出 npp 及植被覆盖度计算完成提示框（图 30）。 NPP 计算结果经过程序计算后，会生成 4 个文件，如下所示： npp_sum：年度植被净初级生产力npp_time_series：每个月的植被 NPP veg_cov_mean：植被年平均覆盖率 veg_cov_time_series：每个月的植被覆盖度

实验结果制图输出

（1）数据导入 ArcGIS

在 ArcGIS 中栅格图像进行分级显示时需要计算图像的统计特征，可以通过ArcToolBox→Data Management Tools→Raster→Raster properties→Calculate Statistics 工具计算图像统计特征。在弹出的 Calculate Statistics 窗口中，在 Input Raster Dataset 选择需要统计的栅格数据，单击 OK。

（2） NPP 分级色彩设置添加经纬网格、指北针、比例尺等成图输出结果展示：