使用ENVI时序立方体和随机森林进行鄱阳湖区单双季水稻信息提取
ENVI时序立方体
ENVI5.6.3新增的时序分析工具可以构建时序立方体。时序立方体是一个多波段栅格，其中每个波段代表同一地理区域某一时期的影像。输入多时相影像，选择不同的光谱指数，为每个时期创建一个光谱指数波段，每个波段都带有一个时间戳。创建时序立方体后，您可以像使用任何其他栅格一样使用ENVI的工具（例如监督或非监督分类）来执行时间序列分析。
ENVI5.6.3软件的试用请浏览：https://envi.geoscene.cn/envi_license
1. 数据源
采用哨兵2号Level-1C数据，下载从2022年3月到11月的鄱阳湖南部地区单景哨兵影像，共计23景。影像时间见下表。
|
20220302 |
20220312 |
20220406 |
20220411 |
|
20220421 |
20220506 |
20220516 |
20220625 |
|
20220710 |
20220804 |
20220809 |
20220819 |
|
20220824 |
20220913 |
20220918 |
20220928 |
|
20221003 |
20221013 |
20221018 |
20221023 |
|
20221102 |
20221107 |
20221112 |
|
单景哨兵2号影像覆盖鄱阳湖南部区域，下图为标准假彩色显示结果。
2. 技术流程
整体技术流程包括时序立方体构建、时序特征分析和随机森林分类。
3. 构建时序立方体
3.1. 数据打开
对于大批量数据的打开显示，可以使用“自定义显示打开影像工具”，工具安装参考：https://www.cnblogs.com/enviidl/p/16396353.html
在ENVI菜单栏选择File>Open with Custom Display，打开自定义显示加载工具。
在Open with Custom Display对话框，点击Open Files右侧添加按钮，打开选择数据对话框，在右上角输入筛选条件（可以用*通配符），按住键盘Ctrl键选择加载的数据。
在Open with Custom Display对话框，按照默认参数：加载第一个多光谱数据集，RGB选择4,3,2波段标准假彩色显示，拉伸方法按照2%线性拉伸。
3.2. 大气校正
在ENVI Modeler中，构建批处理工作流对影像进行批量QUAC快速大气校正得到地表反射率结果，批量大气校正模型如下图所示：
3.3. 构建时序立方体
在ENVI Tool工具箱，选择Spatiotemporal Analysis/Build Temporal Cube打开构建时序立方体工具，在Data Selection数据选择对话框，点击下方Select All，选择构建时序立方体的多期数据，点击OK。
注：数据应至少选择两期数据，且输入栅格必须具有有效的"Acquisition Time"采集时间元数据信息，ENVI预处理后的数据依然会保留采集时间元数据信息。
弹出Build Temporal Cube参数设置对话框，如下图所示：
在Build Temporal Cube对话框中，Overlap Operation重叠操作：
- Intersection交集，输出的时序立方体仅包括多期栅格重叠的范围（默认）
- Union联合，输出的时序立方体包括所有栅格包含的范围
Spectral Index：选择一个光谱指数指数用于生成时序立方体。此处选择归一化差值植被指数NDVI作为构建时序立方体的指数。工具会根据输入影像的中心波长自动列出可以计算的所有光谱指数。光谱指数包括但不限于植被、矿物、大气等类型：
|
明矾石指数（ALUI） |
花青素反射指数 1 (ARI1) |
花青素反射指数 2 (ARI2) |
耐大气植被指数 (ARVI) |
烧伤面积指数 (BAI) |
|
方解石指数 (CALI) |
碳酸盐指数 (CARI) |
纤维素吸收指数 (CAI) |
粘土蚀变指数 (CLAI) |
粘土矿物比例（粘土） |
|
类胡萝卜素反射指数 1 (CRI1) |
类胡萝卜素反射指数 2 (CRI2) |
差异植被指数（DVI） |
疾病水分压力指数 (DWSI) |
白云石指数 (DOLI) |
|
增强型植被指数 (EVI) |
绿帘石/绿泥石/角闪石指数 (ECAI) |
三价铁蚀变指数 (FEAI) |
亚铁 (Fe2+) 指数 (FEI) |
黑色金属矿物比例（铁） |
|
硅酸亚铁指数 (FESI) |
全球环境监测指数（GEMI） |
绿色大气阻力指数 (GARI) |
绿色叶绿素指数 (GCI) |
绿色差异植被指数 (GDVI) |
|
绿色归一化植被指数 (GNDVI) |
绿色优化土壤调整植被指数 (GOSAVI) |
绿色比率植被指数 (GRVI) |
绿土调整植被指数 (GSAVI) |
绿土调整植被指数 (GSAVI) |
|
含羟基 (OH) 变质矿物指数 1 (OHI1) |
含羟基 (OH) 变质矿物指数 1 (OHI1) |
含羟基 (OH) 变质矿物指数 3 (OHI3) |
红外植被百分比指数 (IPVI) |
氧化铁比 (Iron Oxide) |
|
高岭石指数 1 (KAI1) |
高岭石指数 2 (KAI2) |
高岭石指数 3 (KAI3) |
红土指数 (LATI) |
叶面积指数 (LAI) |
|
叶绿素指数 (LCI) |
叶水植被指数 1 (LWVI1) |
叶水植被指数 2 (LWVI2) |
菱镁矿指数 (MAGI) |
甲烷指数 (METHANE) |
|
MgOH-碳酸盐丰度指数 (MGAI) |
改良叶绿素吸收率指数 (MCARI) |
改良叶绿素吸收率指数 (MCARI) |
修正非线性指数 (MNLI) |
修正红边归一化植被指数 (MRENDVI) |
|
修正红边简单比率 (MRESR) |
修正红边简单比率 (MRESR) |
修正三角植被指数 (MTVI) |
修正三角植被指数 - 改进 (MTVI2) |
水分压力指数 (MSI) |
|
蒙脱石指数（MONI） |
白云母指数 (MUSI) |
归一化燃烧比 (NBR) |
归一化燃烧比热 1 (NBRT1) |
归一化差异累积指数 (NDBI) |
|
归一化差异红外指数 (NDII) |
标准化差异木质素指数 (NDLI) |
归一化差异氮指数 (NDNI) |
归一化差异雪指数 (NDSI) |
归一化植被指数 (NDVI) |
|
归一化差水指数 (NDWI) |
归一化多波段干旱指数 (NMDI) |
归一化多波段干旱指数 (NMDI) |
非线性指数 (NLI) |
优化土壤调整植被指数 (OSAVI) |
|
Phenitic 生物指数 (PHEI) |
光化学反射率指数 (PRI) |
Phyllic植物变迁指数 (PHAI) |
植物衰老反射指数 (PSRI) |
丙岩蚀变指数 (PRAI) |
|
富石英指数 (QRRI) |
石英/硅质岩指数 (QSRI) |
红边归一化植被指数 (RENDVI) |
红边位置指数 (REPI) |
红绿比指数 (RGRI) |
|
重整化植被指数 (RDVI) |
二氧化硅指数 (SIDI) |
简单比率 (SR) |
土壤调整植被指数 (SAVI) |
硫酸盐指数（SULI） |
|
综合绿色指数 (SGI) |
结构不敏感色素指数 (SIPI) |
结构密集色素指数 (SIPI1) |
转化叶绿素吸收反射指数 (TCARI) |
转换差值植被指数 (TDVI) |
|
三角植被指数 (TVI) |
可见大气阻力指数 (VARI) |
Vogelmann 红边指数 1 (VREI1) |
Vogelmann 红边指数 2 (VREI2) |
水带指数 (WBI) |
|
WorldView 综合指数 (WV-BI) |
WorldView 改良植物指数 (WV-VI) |
WorldView 新铁指数 (WV-II) |
WorldView 非同质特征差异 (WV-NHFD) |
WorldView 土壤指数 (WV-SI) |
|
WorldView 水指数 (WV-WI) |
|
|
|
|
Date Format：指定一个时期格式，时序立方体中的波段名称将按此日期格式命名。
Output Raster：选择时序立方体的输出路径和文件名，点击OK，生成时序立方体栅格。
4. 作物时序特征分析
在进行作物类型识别之前，首先要了解不同类型作物的生育时期和时序特征差异。
4.1. 作物生育时期
鄱阳湖区主要种植水稻，根据熟制情况分为单季稻和双季稻。鄱阳湖区双季稻生长期安排因受气候条件的制约相对固定, 大致从每年3 月下旬到10月下旬, 而单季稻生长期安排相对自由, 且全生育期略长, 一般从5月中上旬到10月上旬。近20a农业气象数据表明, 鄱阳湖区不同熟制水稻种植时间安排上略有差异, 但总体情况如下图所示。该图为鄱阳湖区年内不同熟制水稻物候历示意图[1]，图中A为早稻，B为晚稻，C为单季稻。
4.2. 时序特征
从物候历可以看出，双季稻早稻在六月中旬发育旺盛，晚稻九月中旬为发育高峰期，这期间植被指数NDVI较高；在四月和十月下旬处于育秧和收割阶段，NDVI相对较低。因此按照红通道：9月18 NDVI（或6月25日NDVI），绿通道：4月21日NDVI，蓝通道：10月23日NDVI的波段加载方案显示时序立方体。在Data Manager依次选择这三个时期的NDVI波段，点击Load Data。
加载显示的时序立方体图像如下图所示，图中粉色区域为不同的作物。
接下来，进一步观察单季稻和双季稻特征。在ENVI上方工具栏，点击Spectral Profile按钮，打开波谱曲线对话框，鼠标定位到图上粉红色区域，查看双季稻NDVI时序特征，可以看出双季稻NDVI时序有明显双峰特点。如下图所示：
之后，鼠标定位到暗绿色区域，查看NDVI时序特征，可以看到曲线呈单峰特征（单季稻）。
以上仅举例这两种典型地物的时序特征曲线。
5. 随机森林分类
识别出不同地类的时序特征曲线之后就可以进行作物类型识别，本文使用随机森林方法进行作物类型识别。
5.1. 绘制样本
随机森林分类方法首先要绘制样本。在ENVI上方工具栏，点击Region of Interest (ROI) Tool，打开ROI工具面板。添加如下几个类别：
- 单季稻
- 双季稻
- 其他
在ROI工具面板，选择Geometry>Point绘制点状ROI。
在图中寻找这几类类别并绘制ROI，如下图所示。图中左侧为双季稻样本，中间为单季稻样本，右侧为其它地类样本。
5.2. 随机森林分类
绘制好样本之后即可使用随机森林分类工具执行分类。在Toolbox工具箱，选择Machine Learning>Classification>Supercised>Random Foreast Classification随机森林分类工具。此工具对单个栅格影像执行随机森林监督分类。
在随机森林分类对话框中：
- Input Raster：选择构建的时序立方体。
- Input ROIs：选择绘制好的ROI 类别。
- Background Labels：在“背景标签”字段中，指定要用作背景的 ROI 类。这里指定Other为背景类。
- Balance Classes：选择是否平衡类别。默认为Yes，平衡类别表示指定在训练期间应将所有类视为平衡的。选择Yes有助于增加样本少的类别的提取范围。这里按照默认选Yes。
- Estimators：输入要使用的决策树的数量。估计器是算法的预测器。默认值为100。此处按照默认。
- Max Depth：指定树的最大深度。如果未指定，则扩展节点直到所有叶子都是纯的。
- Output Raster：选择分类结果的输出路径和文件名。
单击OK执行随机森林分类，得到分类结果，如下图所示：
6. 总结
本教程介绍了使用ENVI时序立方体进行作物类型识别，区分单双季水稻。此外还可以使用时序立方体进行作物健康趋势分析，揭示较大地理区域或特定田地内作物健康状况的变化。此方法也可以运用于大气、矿物等领域典型地物的特征分析。
7. 引用
[1]李鹏，姜鲁光，封志明，等. 鄱阳湖区粮食供给功能的空间格局分析[J]. 自然资源学报.
2011, 26(2): 190-200.
