一种基于时序遥感影像混合像元分解的耕地种植强度估算方法——以湖北省为例

刘轶青， 陶建斌， 陈 曦， 陈 濂， 卫诗琪

(地理过程分析与模拟湖北省重点实验室/华中师范大学城市与环境科学学院， 武汉 430079)

农业是人类社会生存和发展的基石，耕地是进行农业生产的基本资源和条件.随着城市化迅速发展，大量耕地被转化成建设用地[1-2].同时，随着人口持续增长、饮食结构改变和非食物用途(如生物能源)农产品消耗的增加，人类对农产品的需求不断.而在我国，市场经济的发展对农民种粮收益的冲击、农户个体追求利益最大化而增加经济作物生产，使得部分地区农作物种植频率呈现下降趋势(表现为耕地撂荒和冬闲田).如何缓解这种人地矛盾，提高现有耕地的集约化利用程度，是国际社会普遍关注的焦点问题[6].

现有耕地集约化利用的指标定义主要有种植频率和复种指数[7-12].种植频率是对种植模式的硬划分(如单季作物和双季作物)，缺乏详细的空间信息[13-15].复种指数的计算多基于统计数据，忽略了行政区内作物分布状况的空间异质性[16]，且数据存在时效性和准确性问题[17].随着遥感技术的发展，利用时间序列数据，进行大尺度的耕地集约化利用估算成为可能[18-20].但在我国由于复杂地形和复杂种植结构导致的混合像元的影响，估算结果存在不确定性.总之，现有耕地集约化利用指标及其估算方法不能满足精细作物制图的要求，如何获得高时空分辨率的数据集亟待解决.

本文基于MODIS(Moderate Resolution Imaging Spectroradiometer)时间序列数据，在时间维度构建一个特有的特征空间，基于线性光谱分解的方法进行混合像元分解来估算种植强度.本文方法为复杂地形、复杂种植结构下高时空分辨率的种植强度估算提供了一种新方法，对于耕地的集约化研究有着重要现实意义.

1 研究区域及数据

湖北省地处长江中游(108°21′42″E～116°07′50″E、29°01′53″N～33°6′47″N)，是我国主要的粮食生产基地，素有“鱼米之乡”的美誉.全省70%的耕地主要集中分布在江汉平原、鄂东沿江平原及鄂中丘陵地区；兼有水田和旱地，主要作物类型有水稻、小麦、玉米、油菜等.

图1 研究区域土地覆盖状况Fig.1 The main land-cover types in the research area

本文采用的主要数据为：2015年MOD13Q1(MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid)EVI(Enhanced Vegetation Index)数据，来源于美国国家航空航天局MODIS WEB网站(https://modis.gsfc.nasa.gov/).基于湖北省矢量行政边界，研究区所涉及图幅编号为h27v05、h27v06、h28v05、h28v06.研究数据可满足大尺度、长时间连续监测地表植被的需要.

在数据预处理上，采用MRT工具(MODIS Reprojection Tool)进行影像的拼接、重采样、重投影处理.同时，考虑到数据受到云、雪、阴影等诸多因素的影响，原始EVI曲线季节变化不明显，利用数学方法，使用TIMESAT工具重建原始数据.TIMESAT作为平滑长时间序列数据的专业滤波处理工具，集成了高斯滤波以及S-G滤波等滤波方式.本文采用的是S-G滤波算法，除达到去云、消除离异值的目的以外，还使得EVI数据在时间维度上的时序变化更加稳定.

除此之外，高分辨率遥感数据采用了Landsat 8多光谱遥感数据，空间分辨率为30 m，重返周期为16 d，数据来源于美国地质调查局(https://glovis.usgs.gov).结合农作物的物候规律，选择时相在1月中下旬至2月上旬的遥感影像，可以较好区分二者，获取较为纯净的双季农作物端元.在该时间段内，双季农作物中的第一季(主要是冬小麦和冬油菜)长势较好，而单季农作物尚未播种.在标准假彩色合成方式下，双季农作物在影像中呈现红色，单季农作物或者撂荒地呈现灰黄色.本文中，Landsat 8遥感数据主要用于端元选取的参考和样本数据的制作.相比于MODIS数据，能在较大尺度上满足精细观测的需要，提供更为丰富的地表信息.

2 研究方法

2.1 样本数据制作

本文利用Landsat 8遥感影像制作与MODIS空间分辨率一致的耕地种植频率的样本数据(图2).

利用ISODATA方法对样本区Landsat 8遥感影像进行分类，最小分类数定为10，最大分类数定为20，迭代20次.通过分类后处理，将分类结果合并成单季农作物、双季农作物和非耕地这3个类别.将样本数据与估算得到的种植强度数据在像元尺度上进行比较，以用于精度检验与评价.

图2 种植频率样本数据制作示意图Fig.2 The schematic diagram of preparing cropping frequency samples

2.2 特征空间构建

根据农作物物候规律，双季农作物的EVI曲线在一年中的第65 d和第225 d左右会各自出现一个波峰，在第145 d左右出现一个波谷.水体及建设用地在全年内EVI曲线总体数值较低，没有较大波动.天然植被在第145 d之后长势最好，EVI曲线出现一个明显的峰值.对时相EVI数据分析后发现，不同土地覆盖类型具有不同EVI时间谱曲线(图3)，农作物具有鲜明的区别于其他地物类型的季节性物候特征.这些明显的物候特征，是区分不同地物类型的基础.

本文选取主成分变换后的前两个主分量构建特征空间，特征空间呈现为三角形(图4).利用ENVI软件的N维可视化工具，建立从影像空间到特征空间的投影，从而定位、识别特征空间的端元.通过查看其时间谱曲线后，发现这3类点群分别为：双季农作物、天然植被、水体.绿色一端为天然植被，从顶点到中心，土地覆盖类型向灌丛和草地过渡；红色一端为双季农作物，趋近中心，向单季农作物过渡；蓝色一端为水体，从顶点至中心，土地覆盖类型过渡为建设用地.

图3 主要土地覆盖类型时间谱曲线Fig.3 Time series profiles of the main land-cover types

图4 利用PCA 1和PCA 2构建二维特征空间Fig.4 The two-dimensional feature space based on PCA 1 and PCA 2

2.3 端元选择

端元(End Member)代表某种纯净的地物，所选端元的质量直接影响混合像元分解的精度[21-22].根据所构建的特征空间，端元的选择应覆盖双季农作物、天然植被、水体这3类主要土地覆盖类型.考虑到不同地域同一地物类型的差异性，所选端元在研究区内尽量均匀分布且端元总数达到像元总数的3%[23].

本文利用不同农作物的物候规律，以关键物候期增强型植被指数的假彩色合成影像为底图，在时间序列 MODIS 影像上选择端元，双季农作物在影像中显示为品红色.为保证端元的纯净性，本文还利用 Landsat8遥感影像以及Google Earth影像上大块的连续的纯净地物作为参考，辅助进行端元选取.

2.4 丰度估算

本文利用全限制性线性光谱分解的方法进行丰度估算.不同于在光谱维度上进行混合像元分解，本文方法是在时间维度上进行的.通过对一年23期的MODIS时间序列遥感影像进行混合像元分解，得到双季农作物、天然植被、水体的丰度影像.线性光谱分解模型表达式(式1)及约束条件(式2、式3)如下：

(1)

(2)

(3)

式中，Rib为一年中的第b期影像上的第i像元的像素值；fki为对应于第i像元的第k个基本组分所占分量值；Ckb为第k个基本组分在一年中的第b期影像的像素值；n为像元i所包含的基本组分数目；εib为误差项.

通过解混可得到3类主要地物的丰度数据，考虑到丰度的物理意义，本文利用MATLAB软件将丰度影像的值域控制在0～1.双季农作物的丰度可代表耕地种植强度信息，0值代表撂荒地，1代表较为纯净的双季农作物像元，根据丰度影像的像素值高低可直观判断耕地种植强度大小.

图5 3种主要地物类型的端元选择(影像RGB合成方案：DOY 065， 145， 065)Fig.5 Endmember collections of the main land-cover types (RGB composite: DOY 065， 145， 065)

3 结果与分析

通过丰度估算，本文得到湖北省2015年耕地种植强度结果(图6).对关键物候期的MODIS假彩色合成影像、湖北省种植频率数据、种植强度估算结果进行目视比较，三者具有较好的空间一致性(图7).相比于MODIS原始数据，基于混合像元分解所得的种植强度数据，降低了混合像元的影响，空间分辨率得到提高.相比于湖北省种植频率数据，种植强度数据避免将耕地种植模式硬性划分为单季、双季或非耕地，所蕴含的信息更真实、丰富.

图6 湖北省2015年耕地种植强度Fig.6 Cropping intensity in Hubei Province in 2015

(a)和(b)为MODIS假彩色合成数据(影像RGB合成方案：DOY 065， 145， 065)，(c)和(d)为湖北省种植频率数据，(e)和(f)为种植强度估算结果(a) and (b) are MODIS false color composite( DOY 065，145，065)， (c) and (d) are cropping frequency in Hubei， (e) and (f) are the estimated cropping intensity图7 2015年种植强度估算结果目视比较Fig.7 Visual Comparison of cropping intensity in 2015

除目视比较之外，本文也采用定量方法进行精度检验.首先，将本文得到的种植强度估算结果与2015年中国农田熟制遥感监测数据集进行县级尺度的比较[24]，R2达到0.90.

图8 湖北省2015年县级尺度的精度检验Fig.8 Accuracy validation at county level in Hubei Province in 2015

其次，选择和样本数据同一区域的影像进行像元尺度上的验证.利用相关性分析，验证二者间的线性关系，R2达到0.84(图9).

图9 湖北省2015年像元尺度的精度检验Fig.9 Accuracy validation at pixel level in Hubei Province in 2015

4 讨论与结论

时间序列遥感影像使得耕地资源集约化利用的精确表达成为可能.本文引入“种植强度”这一概念，种植强度是亚像元尺度的农作物分布和种植频率状况的综合表征.相较于“种植频率”、“复种指数”等传统耕地集约化利用指标，种植强度避免将种植模式硬性划分为单季或双季农作物，得到更高时空分辨率的子像元水平的丰度信息，更准确地表征了农作物的种植“强度”.

利用“种植强度”这一概念，在时间维度构建特征空间，并使用线性光谱分解的方法估算耕地种植强度，在亚像元尺度对耕地集约利用程度进行表达.经检验，估算结果精度较高，证明了本文方法提取种植强度的可靠性.

本文主要探讨的是一种基于时间序列遥感数据进行混合像元分解的种植强度估算方法.不同于多数已有的在光谱维度上进行的混合像元分解，本文的混合像元分解方法是在时间维度上进行.不同的土地覆盖类型具有其特有的EVI时间谱曲线，因此，农作物具有鲜明的不同于天然植被、水体的物候特征.本文构造的特征空间，能够较好地区分农作物与其它土地覆盖类型，以及不同种植强度的农作物.基于此提取的耕地种植强度能够在亚像元级别实现对耕地集约利用程度的表达，优于传统的种植频率和复种指数，同时，可满足高时间分辨率和高空间分辨率的耕地集约利用监测要求.

农作物与天然植被的季节性物候差异是本文进行混合像元分解的理论基础.但考虑到湖北省境内农作物物候的地域性差异带来的“同物异谱”的现象[25]，以及所收集到的端元不够纯净的情况，精度在一定程度上受到影响.其次，精度验证时所利用的Landsat遥感数据，受到获取时相的限制，也无法完全反映验证区域的真实种植情况.