TIMESATv3.3 EVI物候提取

本文是在兰大付阳师兄的TIMESAT软件操作教程的基础上，结合自己实际操作经历，以及同学和网友分享的经验，撰写的个人学习记录。由于刚刚开始学习植被物候的遥感提取，对于许多机制还没有深刻的理解，难免存在错漏，希望能和知友们多多交流，感谢！

1. 准备数据

1.1 下载数据，转换格式

首先在NASA EARTHDATA官网上下载数据。以MOD13Q1植被指数数据集为例，首先在MRT软件中批量拼接、裁剪，将hdf文件批量转换为TIMESAT可以打开的dat或img格式（关于MRT软件的使用方法请看虎鲸妹妹：MRT批量提取hdf文件特定波段并拼接影像）；

MOD13Q1产品时间分辨率16d，时间跨度3年（TIMESAT要求至少3年的数据才能得到中间的1年物候），每年23期数据，空间分辨率250m；像元值为16位整形（务必确认数据类型，可以在arcgis里右键属性pixel depth找到）；

1.2 制作数据路径list

创建一个txt文档，第一行为图像的总期数，列出所有图像的路径（在Excel里就能快速实现）：

2. 数据处理

运行Matlab，在命令行输入大写的TIMESAT，回车，调出TIMESAT；

2.1 输入图像预览 【TSM_Imageview】【数据预览不是必要步骤，初次学习建议尝试预览，以后可以直接跳过】

打开刚才的路径文档 list.txt；

选中其中一张图像，填写正确的行列号（行列不要填反）；选择正确的像元值类型（选错会导致无法加载）；

如果不清楚像元值类型，可以在下载数据的产品说明页面查找，一般来说MODIS的植被指数都是16位整形（对原始数据进行诸如乘以0.0001系数的运算，会导致数据深度变成32或者64位）；

输入行列数，点击draw，逐幅查看图像是否存在异常（数据本身存在的问题，比如缺失某一期，或拼接时出现很大缝隙等，会影响到后续物候信息提取）。

2.2 设置滤波和提取方式【TSM_GUI 】

点击load data，选中一小块区域（不建议全局查看，会卡住或者报错），查看曲线拟合情况，自行选择一种合适的拟合方法：

目前很多物候提取的文章使用的都是S-G方法，这是一种高阶多项式拟合的滑动窗口平滑法，相对能保留一定的细节；D-L法也比较常用，与S-G方法相比更不容易被异常值/极值干扰；

如果采用季节性振幅法，需要进行生长季开始/结束的阈值设定：

TIMESAT默认使用的是第一个Seasonal amplitude 季节振幅法，也就是动态阈值法，需要结合自己的研究区实际，拟定一个合适的动态阈值，可以参考前人研究确定，TIMESAT软件的作者建议可以用0.2（20%振幅）试试，不合适可以再改：

2.3 添加质量文件和土地利用文件（选填）保存设定

不要忘了将处理的行列号范围设置成全图！（行列号的开头是1不是0）保存为.set格式文件；

如果有像元质量文件，可以添加进来，同样是以一个list.txt的形式，制作方法和上述EVI制作方法类似，此处不赘述；

注意，如果像元质量文件中，质量差的像元比重太高，可能会导致预览界面报错，建议多点几下plot the next series按钮，看看报错的概率，如果只是偶尔几个像元报错，那么质量文件应该还是可以用的；

像元质量文件的权重可以结合TIMSAT3.3用户手册和遥感产品的说明书/用户手册决定，例如MOD13Q1的质量文件说明：

TIMESAT用户手册中关于质量文件和权重的建议：有云的像元权重设置成0.1，混合的设置成0.5，质量很好的设置成1.0：

如果有土地利用类型文件，也可以添加进来：

在 Data processing 中打开刚才保存的.set格式文件【如果nodata像元过多，data process这步可能会报错，可能需要退回数据准备的第一步，适当调整输入的图像范围。比如，我的研究区位于沿海，nodata值比较多，于是想到使用一个规则的矩形矢量掩膜裁剪经MRT融合后的数据，使得图幅包含研究区，并且Nodata值较少】

运行，耐心等待，得到以下结果，表示拟合文件创建成功；

2.4 后处理post-processing部分

注意：以下5个按钮都属于预览结果，不是必要步骤。第一次使用TIMESAT可以操作一下，熟练了以后可以直接跳到【TSF_seas2img】这一步。

2.4.1【TSM_fileinfo】的功能：打开上一步生成.tpa .tts文件查看信息，例如年份数，期数，行列数等，看看数据是否正常；

2.4.2【TSM_printseasons】的功能：指定处理范围内一个小区域，输出一小部分拟合结果，供用户预览。可以输一个比较小的范围，逐个像元按回车查看一下拟合结果：

2.4.3【TSM_viewfits】的功能：查看拟合曲线，与上一步相似，结合2.3步骤里生成的拟合文件，预览拟合曲线：

接下来是TSF系列的按钮：

2.4.4【TSF_fit2time】的功能：创建ASCII码格式的拟合时间序列。打开2.3步骤生成的fit.tts格式的数据，弹出命令窗，输入一个小范围，命名一个.txt格式的文件；

运行成功后，在路径（Matlab当前工作路径）下找到生成的txt文件，改成合适的后缀就可以在TSM_imageview 里查看【这步我有一点小问题，还没能成功显示图像，但好像这一步不能显示的话，关系也不大，可以继续后续操作，以最终结果能否在ENVI里成功显示为准】

2.4.5【TSF_fit2img】的功能：查看拟合图像，还是打开之前的fit.tts文件，根据cmd窗口里建议的值输入就可以；

我的数据图像总共有69期，可以选择其中的第35期查看（因为35期大概对应第二年的峰值），像元值为拟合后的EVI值：

这一步起到的是预览图像的作用；

2.4.1~2.4.5 这些按钮只起到预览功能，嫌麻烦可以直接跳过。

2.4.6 TSF_seas2img的功能：输出物候参数图像【重点】

按提示选择需要输出的参数（比如1=生长季始期，3=生长季长度），提取所需的时间段（比如我选择第20~50期，保证能完整覆盖第二年），保存为16位整形（考虑到一般情况下日期用整数表达，但有的文献也用小数点，也可以保存为32位浮点）

【请记住这一步自己选择了什么类型的数据，因为这一步输出的是不含各种地理信息和像元深度等信息的二进制文件，因此在ENVI里打开的时候，要手动选择数据类型，否则不能显示！】

根据3.3使用手册里的说法，对于只有一个生长季的植被来说，我们应该打开season1图查看，在TSM_imageview中打开查看如下：

右边图例中的数字代表了期数，需要转换成天数。由于3年的数据只能提取出中间一年的物候，第一年的物候是无效的，需要减去第一年的期数，再乘以每一期的天数得到具体的天数。举个栗子，上图中某个像元的像元值为30，则它代表的天数为[30-（23+1）]*16天。天数换算在arcgis里利用栅格计算器（raster calculator）就可以完成。

由于得到的season1文件是二进制文件，无法直接在arcgis中打开，利用ENVI可以将其转为arcgis能够打开的tif或者dat格式（如果转出来的tif格式打不开，可以转成dat试试）；

以我的ENVI 5.3x64位为例：

open as——binary——选中刚才生成的物候结果二进制文件（如begin_s1）

需要我们手动给这个二进制文件添加信息，点击header info——input header info from——open——找一个区域（行列数）完全一致的tif或者dat数据放进去，为season1提供地理信息（我用的就是evilist里面的第一个.dat文件）。

datatype这一行要注意，我使用的原始数据是16位整形，但在TIMESAT输出是选择的是32位浮点型，因此直接使用原始数据的头文件是会报错的，需要手动修改。行列数目也要写对。File Type最好选unknown，选tif和dat有可能会报错。另外文件路径里不要有中文！

如果第一次设置错误，导致第二次直接无法在ENVI中打开season1文件，是由于第一次设置的时候自动生成了.hdr头文件（相当于它自动保存了你刚才错误的操作，导致你第二次打开的时候ENVI自动读取了错误的设置）。请在文件目录里找到自动生成的同名.hdr文件删掉，再重复上图操作。

然后在ENVI里 File——save as选中该图层，将图片输出为tif或者dat格式，就可以在arcgis里打开了：

如果你有很多地理信息等参数完全相同的season1文件，需要重复上述的添加信息的步骤，那么可以不必在ENVI里一个个打开输出，只需要复制之前生成的.hdr头文件，把名字改成一样的就可以了（批量修改文件名可以通过Python或者Matlab等实现，数量不多的话可以手工复制改一下）。

有同学问我为什么提取后的物候期在ArcGIS中打开会出现负值，典型的负值如-1和-2是用户在TSF_seas2img这一步里按照系统提示自己定义的缺省值（本教程里有这个设置的截图），所以-1和-2是正常的，只需要用ArcGIS的按属性值提取（extract by attributes）去掉就可以。按属性值提取时，建议参考前人研究和自己研究区的实际物候情况，设置一个物候期开始/结束的正常范围，排除掉明显不正常的像元。

【2021.7.20第12次更新】