基于FastICA算法和MODIS数据的水稻面积提取

耿利宁，景元书，杨沈斌，浩宇

(南京信息工程大学1．江苏省农业气象重点实验室;2．应用气象学院，江苏南京210044)

0 引言

水稻是世界主要粮食作物之一，全世界近一半人口以水稻为食。随着大多数水稻种植国家城市化进程加速，水稻种植面积逐年减少，已经关系到世界的和平与稳定。为此，采用当前先进的遥感技术进行长期、大区域的水稻生产监测成为保障水稻安全生产的重要手段，而提出动态、可靠、有效的水稻面积遥感监测方法是其中的重要内容。

近十多年来，学者们提出了不同水稻面积的遥感监测方法，面对大区域的水稻动态监测，光学数据仍是水稻面积提取研究和应用的主要数据源，其中MODIS数据发挥重要作用。MODIS的光谱观测通道大大增强了对地表类型的识别能力和对不同生态系统的观测能力(邓学良等，2009;陈爱军等，2012)，在大区域水稻遥感监测中同样得到广泛应用(Boschetti et al．，2009;Sakamoto et al．，2009;Xiao et al．，2005)。而根据 Lay et al．(1983)提出的相似性指数，顾晓鹤等(2008)、杨沈斌等(2012)提出了通过建立参考水稻EVI(Enhanced Vegetation Index，增强植被指数)生长曲线和计算相似性指数的方法，从多时相MODIS数据中提取水稻面积。

然而，MODIS数据的空间分辨率较低，最高为250 m。较低的空间分辨率给大区域水稻遥感监测带来的严重的混合像元问题，尤其是在稻田分散且面积较小的种植区。针对混合像元问题，目前常用的方法有线性混合模型分解、模糊分类模型分解以及盲源分解。其中，独立成分分析(Independent Component Analysis，ICA)属于盲源分解方法，而FastICA算法则是一种快速的ICA算法。它作为一种新兴的信号处理方法，已在越来越广泛的领域获得成功应用 (胡茂桂和王劲峰，2010;李根等，2014)。而应用于农业遥感方面，该算法不仅能有效分解出不同作物类型，还能通过计算得出不同作物的丰度(Wang and Chang，2006;Ozdogan，2010)，已有研究表明，地物的类型越复杂，计算得到的丰度结果越不准确，同时针对同一作物而言，不同的田间管理水平也会导致不同的丰度结果。

基于已有利用MODIS数据提取水稻面积方法的研究，本文拟在使用相似性指数计算获得潜在水稻像元基础上，采用FastICA算法，对潜在水稻像元的归一化差分植被指数INDV(Normalized Difference Vegetation Index)时相曲线进行分解(Xiao et al．，2005)，计算水稻丰度，获取研究区水稻种植面积和分布。同时，结合统计数据和地面观测样方资料，对水稻面积提取方法进行有效性验证。

1 研究区域和数据

1.1 研究区域概况

江苏、安徽、江西3省隶属长江中下游水稻生态区。其中，江西省属于江南丘陵平原双季稻种植区，江苏、安徽的南部地区属于长江中下游平原双单稻种植区，而江苏和安徽北部地区属于华北单季稻种植区。但随着近十多年来的城市化和劳动力转移，江苏和安徽大部分地区仅种植一季中稻和晚稻。

1.2 研究数据

1.2.1 MODIS数据

从 UGGS EROS 数据中心(http://glovis．usgs．gov/)下载了覆盖整个研究区域的2010年46景8 d合成地表反射率数据(MOD09A1)和1景(2009年)土地利用类型(MCD12Q1)MODIS产品数据。MOD09A1数据经过大气校正、气溶胶校正以及卷云处理，空间分辨率为500 m，包含7个波段，数据还包括两个表征质量的波段QA(Quality Assurance)。MCD12Q1数据包含16个类别的土地覆盖类型，本文选用其中的Croplands以及Cropland/Natural Vegetation Mosaic类型作为耕地像元(图1)。

图1 研究区耕地分布与GPS样点Fig．1 Distribution of cultivated land and samples of GPS in the study area

1.2.2 辅助数据

通过国际科学数据服务平台(http://datamirror．csdb．cn/index．jsp)获取了 6 景覆盖研究区域的Landsat TM数据(2005年1景、2006年1景、2007年2景和2010年2景)，空间分辨率为30 m;通过中国资源卫星应用中心获取了2010年3景覆盖研究区域的HJ1A-CCD1影像，空间分辨率为30 m。利用获取的研究区水稻和同期作物(玉米、大豆和棉花)GPS样方资料，结合Google Earth和上述影像，从研究区内额外建立了150个水稻样方数据。样方数据中江苏60个、安徽50个和江西40个，分布如图1(江西30个早稻样方中有20个样方为双季早晚稻)所示。

另外，获取了2010年江苏省、安徽省、江西省统计年鉴数据，提取了三省水稻种植面积数据。从各省近10 a农业气象观测站获取了各观测点水稻生育期资料。表1给出了各省观测到的主要水稻生长期最早至最晚日期的基本情况。

2 研究方法

2.1 水稻面积提取方案

首先在耕地类型像元基础上提取满足(ILSW+0.05)≥INDV(ILSW为地表水分指数，Land Surface Water Index)的水田像元(Xiao et al．，2005);然后，利用样方的INDV时相曲线，建立符合研究区各水稻生长特征的水稻参考生长曲线，对提取的参考曲线进行前后时相滑动，并分次计算每个耕地像元INDV时相曲线与参考曲线的相似性指数，选取多次结果中最低值得到一景相似性指数影像，并结合样方的指数范围确定阈值，筛选符合阈值范围的像元作为水稻像元(顾晓鹤等，2008;杨沈斌等，2012)。然后，获取水稻像元的INDV数据，再利用FastICA算法分解INDV数据，获取每个像元中含水稻的百分比即水稻丰度，以提取水稻面积。最后，利用统计和样方资料对结果进行验证。

2.2 植被指数的计算与滤波重构

ILSW和INDV的计算可参考 Huete et al．(1997)以及Xiao et al．(2006)。由于研究区水稻生长季多云雾雨天气，需要对MODIS数据进行去云滤波和噪声抑制处理。为此，参考杨沈斌等(2012)的数据预处理方法，分别对ILSW时序曲线采用LMF(基于局部最大滤波算法)、INDV时序曲线采用TIMESAT程序中的Savitzky-Golay算法进行滤波和噪声抑制。两种算法都是在给定时相窗口内进行的滤波算法，利用窗口范围内数值的平均值或者线性组合替换缺失或者不良数据。原始数据与滤波重构后的数据对比如图2所示。由图2可见，INDV与ILSW数据都得到了较好的平滑并保留了原始的变化特性。

表1 各省主要水稻生长期的基本情况Table 1 Basic situation of the rice growth period in different parts of the study area

图2 INDV和ILSW时相曲线滤波前后的对比Fig．2 Comparison of INDVand ILSWbetween before and after filtering

2.3 水田像元的获取

已有研究认为(Xiao et al．，2005)，判别水田像元可通过判断该像元的ILSW与INDV之间是否存在(ILSW+T)≥INDV的关系来确定，其中T值通常取值0.05。因此，针对耕地像元，获取移栽期时相内满足ILSW+0.05≥INDV的像元为水田像元。

2.4 水稻像元提取

2.4.1 水稻参考生长曲线的建立

符合水田像元提取特征的地物还可能包含湿地或周期性滩地等。为此，需建立符合研究区各水稻生长特征的水稻参考生长曲线。

根据样方数据中的纯像元获取各地区的INDV时相曲线，求多个纯像元的平均值，得到能代表该区域的水稻参考INDV时相曲线，如图3所示。其中，江苏地区用一条参考曲线代表，安徽地区皖北皖南各用一条参考曲线代表，而江西地区用三条参考曲线代表早、中、晚三类水稻。通过与表1的对比可以发现，选取的参考INDV时相曲线符合真实的水稻物候期。

2.4.2 相似性指数的计算

相似性指数SNDVI的原理和计算可参考杨沈斌等(2012)一文。相似性指数值越小表明与参考INDV曲线越相似，是水稻像元的可能性越大。而本文研究区域较大，使得一些水稻像元的时相曲线开始和结束时间与参考INDV时相曲线存在差异，相似性指数值过大而被当作非水稻。为了避免这一不同种植日期带来的干扰，对参考INDV时相曲线进行了前后一个时相的滑动(图略)，多次进行相似性指数计算。

2.4.3 确定水稻像元

图3 水稻参考INDV时相曲线 a．江苏省;b．安徽省北部;c．安徽省南部;d．江西省早稻;e．江西省中稻;f．江西省晚稻Fig．3 Reference INDVtemporal profiles of paddy rice a．single rice of Jiangsu province;b．single rice of northern Anhui province;c．single rice of southern Anhui Province;d．early rice of Jiangxi province;e．middle rice of Jiangxi province;f．late rice of Jiangxi province

叠加多次计算得到的相似性指数影像，选取多波段数据中各个像元的最低值，即与参考INDV时相曲线最相近的结果，综合成一景以相似性指数为特征波段的影像。获取此景影像中样方的相似性数值，确定符合样方的相似性指数范围，其中，江苏省指数范围为0.2～2.1;安徽为0.1～2.0;江西早稻为0.08～1.88，中稻为0.2～1.5，晚稻为0.1～3.0。利用样方获取的阈值范围，筛选出影像中符合范围的像元，确定为水稻像元。

2.5 基于FastICA算法的水稻像元分解

在进行FastICA算法之前，需要对数据进行预处理来简化ICA问题。预处理主要是中心化和白化。设定输入的数据X是观测得到的INDV数据，中心化即为零平均，白化则是使其各行数据互不相关，且各自的方差等于1，使得X的协方差矩阵等于单位矩阵，本文预处理通过MATLAB中的zscore函数完成。

通常FastICA算法的输入数据是以像元个数作为信号个数，像元的时序数作为样本个数，而本文的应用中，不同作物的生长时间序列可能并不符合ICA分析中提及的统计独立，因此对原始数据进行转置，这是一种特别的 FastICA算法(Ozdogan，2010)。使用MATLAB中的FastICA工具箱来完成计算，计算完成后会分别得到混合矩阵A、解混矩阵W以及估计出的独立成分CI，其中A中第i列代表第i个独立成分的时间序列，CI中第i行代表该独立成分信号在整个图像中的强弱分布。一般而言，结果中独立成分的排列顺序符合独立性从强至弱的规律(Wang and Chang，2006)。

FastICA算法在独立成分个数很多的情况下，无法彻底分解所有独立成分，丰度计算的误差也就很大(Ozdogan，2010)。而本文在水稻像元基础上进行的FastICA算法等于间接的对原始数据进行了降维，减少了误差。由于FastICA算法运行输出的结果代表的是每个独立信号占总信号能量的比重，不能直接用于表示丰度。为此，参考 Wang and Chang(2006)一文中的丰度计算公式计算了水稻丰度。

3 结果与分析

3.1 研究区水稻像元

根据水稻相似性指数和样方数据提取出各省水稻像元的分布(图略)，其中颜色越深代表水稻相似性指数值越小，暗示该像元水稻所占比重越大。将各省的结果进行叠加，选取所有种植时期水稻的像元，形成研究区水稻像元分布影像(图4)。

图4显示，江苏水稻主要种植分布在该省中北部地区;安徽水稻主要为沿淮河、长江两岸分布;江西水稻分布也符合沿水沿江分布的特征，但多集中在鄱阳湖周围以及赣江、信江、抚河等流域。早晚双季稻以鄱阳湖为中心的南昌、宜春和吉安东部、抚州北部、上饶西部;而中稻分布较离散，主要集中种植在九江东部以及宜春、吉安、抚州、南昌地区。与Xiao et al．(2005)绘制的对应区域水稻分布结果相比，本文提取的水稻像元分布与其吻合较好。

图4 水稻像元分布Fig．4 Distribution of paddy rice pixel

3.2 研究区水稻像元分解和水稻丰度

为减少非水稻种植期间不同的耕地使用情况带来的干扰，FastICA算法仅对水稻生长季内的INDV时相曲线进行信号分解。

引入降维过程中使用的表征独立成分个数的特征量 NVD(Virtual Dimensionality)(Chang，2003)。针对江苏地区，选取NVD=8时，独立成分个数占总变量的95.12%，分解得出的INDV时相曲线中，符合水稻生长特性曲线的有两条(图5a)，从图中可以看出，两条曲线代表的生长季符合水稻生长季，INDV1较INDV2略长，二者峰值相近，在时间横轴上有所偏移，INDV1曲线在9月出现明显下凹。同样的，安徽地区取NVD=7，占总变量94.14%，共分解出4条曲线符合水稻生长特性的曲线(图5b)，4条曲线同样符合水稻生长季，峰值相近，其中INDV3在9月出现明显下凹且较其余3条曲线生长季长。江西地区早稻取NVD=6，占总变量90.26%，中稻取NVD=8，占总变量92.5%，晚稻取NVD=6，占总变量92.77%，早、中、晚稻各分解出一条INDV曲线(图5c)，3条曲线峰值相近，在时间横轴上，早稻结束后晚稻的曲线开始增长，中稻的曲线位于早晚稻之间，符合江西地区早、中、晚稻的生长季，其中代表晚稻的曲线9月出现明显下凹。

使用FastICA算法分解得到各地区不同的INDV时相曲线后，获取符合水稻生长特性曲线对应的CI矩阵中的行数据，即曲线代表的水稻成分在不同像元的信号强度值，进行丰度的计算，获取每个像元对应的丰度值，绘制代表该曲线的丰度图，三省的水稻丰度如图6所示。

图5 用FastICA提取的INDV时相曲线 a．江苏省;b．安徽省;c．江西省Fig．5 INDVtemporal profiles by FastICA a．Jiangsu province;b．Anhui province;c．Jiangxi province

图 6 江苏省(a，b;对应图 5a 中 INDV1、INDV2)、安徽省(c，d，e，f;对应图 5b 中 INDV1、INDV2、INDV3、INDV4)、江西省(g，h，i;对应图 5c中 INDV早、INDV中、INDV晚)的水稻丰度图Fig．6 The abundance fraction maps of paddy rice in(a，b)Jiangsu province(corresponding to INDV1and INDV2in Fig．5a，respectively)，(c，d，e，f)Anhui province(corresponding to INDV1，INDV2，INDV3and INDV4in Fig．5b，respectively)，and(g，h，i)Jiangxi province(corresponding to INDVEarly，INDVMiddleand INDVLatein Fig．5c，respectively)

由图6可见，江苏地区水稻高丰度区域主要集中在苏中江淮平原，洪泽湖东南沿岸，南通地区以及苏南少部分地区，苏北地区水稻丰度较小。安徽地区水稻高丰度地区集中在淮河流域和长江流域沿岸，皖北地区以及皖南山区水稻丰度很小。江西地区早晚双季稻的高丰度区域明显，主要集中在以鄱阳湖为中心的南昌、宜春和吉安地区;而中稻高丰度区域分布离散，主要集中种植在九江东部以及宜春、吉安、抚州、南昌地区。通过云量统计数据发现，3条出现下凹的曲线代表的区域(图6a、6e、6i)9月有大量云层覆盖，符合实际情况。

3.3 面积提取和验证

在获得水稻丰度图后，通过叠加计算各个丰度图中每个像元的丰度值，进行水稻面积的提取，不同曲线之间采用多次丰度计算结果的相加。对计算结果，分别结合2010年统计年鉴资料以及30个覆盖3省的实测样方资料进行了对比验证，与统计年鉴的对比结果如表2，与样方对比结果如图7。

表2 提取的水稻面积与统计的水稻面积的比较Table 2 Comparison between retrieved rice areas and statistical rice areas

图7 样方面积SS与FastICA提取面积SF的比较(单位:102hm2)Fig．7 Comparision of rice area between samples data(SS)and FastICA results(SF)(units:102hm2)

由表2可见，提取的江苏一季稻面积约为192.9×104hm2，与统计面积相对误差为13.6%;安徽地区水稻提取面积为144.2×104hm2，误差百分比为12.1%;江西地区出提取面积分别为早稻61.7×104hm2，中稻 24 ×104hm2，晚稻 69.7×104hm2，误差过大，平均48.5%左右。由图7可见，FastICA提取的面积SF与样方面积SS平均相对误差小于15%，散点与趋势线的确定系数为0.79。

4 讨论

4.1 相似性算法的缺陷

相似性指数是绝对值叠加的结果，容易受到个别离差偏大点的影响。同时，混合像元的存在使得包含较为混杂地物的像元表现出来的曲线与参考曲线差异很大，即使含有水稻也可能被剔除。相似性算法在今后的研究中还需完善与改进，以获得更准确的提取结果。

4.2 NVD的确定

NVD是降维过程中使用的表征独立成分个数的特征量，在盲源分析中，源信号众多会给分析带来困难，特别是区域过大时，所含地物的类型很多，NVD较难确定，而NVD的大小与丰度计算误差有密切的联系。本文在对水稻像元分解前，虽然大大减少了其他地物的干扰，但水稻种植情况、水稻生长状况和稻田土壤特征在空间上都存在一定差异，导致不可能“完全分解”不同条件下的水稻像元，但能“尽可能分解”不同条件下的水稻像元。以江苏地区为例，当分别取NVD为8、9时，分解得到符合水稻生长特性的曲线都为2条(图8)，且前后两次分解得到的两条曲线之间差异很小，这就证明取NVD=8时，得到的符合水稻生长特性的独立成分已尽可能的被分解。

4.3 山区水稻面积提取误差

从丰度计算结果与统计年鉴对比中，可以看出，江苏与安徽提取的水稻面积结果准确性较高，而江西的准确性很差;对比样方验证结果却发现江西地区的样方提取的丰度值准确度较好，这说明丰度计算方法并不是误差过大的主要原因。对比2010年江西省各市统计年鉴与图1，发现江西地区土地覆盖类型数据没有很好地覆盖山区的耕地，如赣州地区多为梯田以及山间耕地，使得有山地区准确性较差，这是江西地区提取结果误差过大的主要原因。

5 结论

本文结合时序MODIS数据和样方数据，利用FastICA盲源算法对水稻的INDV时相曲线进行了丰度分解，提取水稻面积。结果显示，利用FastICA分解水稻混合像元的INDV时相曲线可以有效提取水稻面积，与统计数据、样方调查吻合较好。值得提出的是，相似性算法提取水稻像元仍有一定的局限性，而丰度的计算方法作为一种归一化的线性拉伸，也带来了一定的误差，并且山区耕地数据的缺乏也使得多山地区面积验证误差较大。我国耕地使用复杂，水稻同期作物类型繁多，在一定程度上增加了水稻遥感监测的难度，而本文的研究为大范围内动态的水稻面积监测提供了重要的参考。

图8 用FastICA提取的INDV时相曲线(黑色曲线为NVD=8时;红色曲线为NVD=9时)Fig．8 INDVtemporal profiles by FastICA(The curves are black and red when NVD=8 and NVD=9，respectively)

致谢:两位审稿专家为本文提出了宝贵的修改意见。谨致谢忱!

Boschetti M，Stroppiana D，Brivio P A，et al．2009．Multi-year monitoring of rice crop phenology through time series analysis of MODIS images［J］．Int J Remote Sens，30(18):4643-4662．

Chang C I．2003．Hyperspectral imaging:Techniques for spectral detection and classification［M］．New York:Kluwer Academic/Plenum Publishers．

陈爱军，梁学伟，卞林根，等．2012．青藏高原地区MODIS反照率的精度分析［J］．大气科学学报，35(6):664-672． Chen Aijun，Liang Xuewei，Bian Lingen，et al．2012．Assessment on the accuracy of MODIS albedos over the Tibetan Plateau［J］．Trans Atmos Scis，35(6):664-672．(in Chinese)．

邓学良，潘德炉，何冬燕，等．2009．中国海域MODIS气溶胶光学厚度检验分析［J］．大气科学学报，32(4):558-564． Deng Xueliang，Pan Delu，He Dongyan，et al．2009．Validation of MODIS aerosol optical depth retrievals over the China Sea［J］．Trans Atmos Scis，32(4):558-564．(in Chinese)．

顾晓鹤，韩立建，张锦水，等．2008．基于相似性分析的中低分辨率复合水稻种植面积测量法［J］．中国农业科学，41(4):978-985．Gu Xiaohe，Han Lijian，Zhang Jinshui，et al．2008．Monitoring of paddy rice plant area based on similar index by multi-resolution remote sensing data［J］．Scientia Agricultura Sinica，41(4):978-985．(in Chinese)．

胡茂桂，王劲峰．2010．遥感影像混合像元分解及超分辨率重建研究进展［J］．地理科学进展，29(6):747-756． Hu Maogui，Wang Jinfeng．2010．Mixed-pixel decomposition and super-resolution reconstruction of RS image［J］．Progress In Geography，29(6):747-756．(in Chinese)．

Huete A R，Liu H Q，Batchily K，et al．1997．A comparison of vegetation indices global set of TM images for EOSMODIS ［J］．Remote Sens Environ，59(3):440-451．

Lay J O，Gross M L，Zwinselman J J，et al．1983．A field ionization and collisionally activated dissociation/charge stripping study of some［C9H10］+ions［J］．Organic Mass Spectrometry，18(1):16-21．

李根，景元书，王琳，等．2014．基于MODIS时序植被指数和线性光谱混合模型的水稻面积提取［J］．大气科学学报，37(1):119-126．Li Gen，Jing Yuanshu，Wang Lin，et al．2014．Extraction of paddy planting areas based on MODIS vegetation index time series and linear spectral mixture model［J］．Trans Atmos Scis，37(1):119-126．(in Chinese)．

Ozdogan M．2010．The spatial distribution of crop types from MODIS data:Temporal unmixing using independent component analysis［J］．Remote Sens Environ，114(6):1190-1204．

Sakamoto T，Phung C V，Kotera A，et al．2009．Detection of yearly change in farming system in the vietnamese mekong delta from MODIS timeseries imagery［J］．Japan Agricultural Research Quarterly，43(3):173-185．

Wang J，Chang C I．2006．Applications of independent component analysis(ICA)in endmember extraction and abundance quantification for hyperspectral imagery［J］．Geoscience and Remote Sensing，44(9):2601-2616．

Xiao X M，Boles S，Liu J Y，et al．2005．Mapping paddy rice agriculture in southern China using multi-temporal MODIS images［J］．Remote Sens Environ，95(4):480-492．

Xiao X M，Boles S，Frolking S，et al．2006．Mapping paddy rice agriculture in South and Southeast Asia using multi-temporal MODIS images［J］．Remote Sens Environ，100(1):95-113．

杨沈斌，景元书，王琳，等．2012．基于MODIS时序数据提取河南省水稻种植分布［J］．大气科学学报，35(1):113-120． Yang Shenbin，Jing Yuanshu，Wang Lin，et al．2012．Mapping rice paddy distribution in Henan Province based on multi-temporal MODIS imagery［J］．Trans Atmos Scis，35(1):113-120．(in Chinese)．