牛海鹏 王占奇 肖东洋

(1.河南理工大学测绘与国土信息工程学院， 焦作 454000；2.河南理工大学黄河生态文明与高质量发展研究院， 焦作 454000)

0 引言

水稻是世界第三大粮食作物[1]。中国是世界上最大的水稻生产国，水稻产量在中国农业生产中占有重要地位[2]。黄河两岸地区是河南省水稻种植区。随着自然、社会和经济的发展，原阳县水稻种植面积从约30 000 hm2下滑至约18 000 hm2。因此，适时、准确监测水稻种植区，对粮食政策制定、粮食安全以及农业发展具有重要意义[3]。

随着卫星技术的发展，利用遥感技术监测水稻种植面积已经得到广泛应用。目前最常用的提取水稻等农作物空间分布的遥感卫星为Landsat数据和MODIS数据[4-7]。文献[8]使用多时段MODIS影像提取了南亚和东南亚地区的水稻分布信息，文献[9]基于多时相MODIS数据提取了河南省水稻种植区的分布。文献[10]利用Landsat数据并结合农户调查，提取了翻阳湖平原水稻主产区的水稻种植区。Landsat数据具有较高的空间分辨率，但是时间分辨率较低，并且容易受云雨天气影响，无法获取作物关键生长期的时序数据[11]，因此只根据Landsat数据进行农作物区提取会受到数据缺失的限制。MODIS数据具有36个光谱波段，时间分辨率较高，但其空间分辨率低，存在大量混合像元，对地面的细节表征能力不强[12]，使基于MODIS数据的农作物提取只能在大尺度区域上取得较好的效果。

近年来，国内外学者针对遥感影像时空融合展开了研究，并取得一些成果[13-15]。文献[16]根据像元的空间相关性和光谱相似性，提出一种自适应遥感影像融合模型(Spatial and temporal adaptive reflectance fusion model, STARFM)，该模型能较准确预测出影像的反射率。文献[17]提出了反射率变化的时空自适应融合算法(Spatial temporal adaptive algorithm for mapping reflectance change, STAARCH)。文献[18]在STARFM模型的基础上提出ESTARFM (Enhanced spatial and temporal adaptive reflectance fusion model)模型，该模型能通过时间和光谱分离理论分析反射率的变化趋势，提高融合模型预测Landsat数据的精度。文献[19]基于类别反射率和类内像元反射率的时间变化特征，提出了STDFM (Spatial and temporal data fusion model)模型。文献[20]对STDFM模型进行改进，提出了ESTDFM(Enhanced spatial and temporal data fusion model)模型。文献[21]对5种主流遥感数据时空融合模型进行对比实验，结果表明，ESTARFM模型的融合效果最好，更能反映地物的细节特征。

利用时空融合模型融合Landsat与MODIS数据对水稻种植面积进行提取的方法值得研究。本文以河南省沿黄稻区原阳县为例，基于ESTARFM模型，对Landsat和MODIS影像进行融合，构建完整的时序Landsat数据，计算得到完整的时序NDVI曲线，采用决策树分类的方法进行水稻种植区域的提取，以期在县域尺度上应用ESTARFM模型提取水稻种植面积取得理想的结果，为综合利用高空间分辨率Landsat数据和高时间分辨率MODIS数据进行农作物遥感识别和监测提供一定的实验基础。

1 研究区与数据

1.1 研究区概况

原阳县隶属河南省新乡市，地处豫北平原，位于113°36′～114°15′E，34°55′～35°11′N之间，区域总面积1 022 km2。南邻黄河，北面是余河通道，地势从西南至东北呈逐渐降低趋势，地貌属于黄河冲积平原。农业以小麦、水稻、玉米为主，其中水稻种植以一季稻为主，每年5月上旬育苗播种，10月中上旬收割。

1.2 数据与预处理

空间分辨率为30 m的Landsat数据(表1)来源于美国地质调查局(United States Geological Survey，USGS)，下载2015年品质较好、云量覆盖小的影像。Landsat8 OLI影像已经基于地形数据进行几何校正处理，因此本次实验中不再对Landsat数据进行几何校正。投影坐标系为UTM-WGS84坐标系，利用ENVI 5.3软件对影像进行辐射定标后，使用FLAASH模型进行大气校正。

表1 遥感数据类型及获取日期Tab.1 Remote sensing data types and acquisition date

全年共23期的空间分辨率为500 m的MODIS13Q1数据(表1)来源于NASA网站(https://ladsweb.nascom.nasa.gov/search)。对MODIS影像利用MRT工具(MODIS reprojection tool)重投影为UTM-WGS84投影，把HDF格式转换为GEOTIFF格式，采用双线性内插法(Bilinear)将影像重采样为30 m分辨率。处理后，影像与Landsat具有相同的投影方式、分辨率和像元尺寸。对Landsat与MODIS的对应波段按顺序一一对应，如表2所示。

表2 Landsat8 OLI与MODIS数据波段设置及对应关系Tab.2 Band settings of MODIS and Landsat8 images

2 研究方法

本文采用ESTARFM模型将Landsat8 OLI数据与MODIS13Q1数据进行融合，从而重构出完整时序的融合Landsat数据；计算得到融合影像的NDVI数据，通过TIMESAT软件中savitzky-Golay(S-G)滤波法进行平滑降噪；得到6种主要地物的NDVI曲线，根据6种地物的曲线特征建立专家决策树进行分类，进而提取出水稻种植区域。具体工作流程如图1所示。

图1 技术流程图Fig.1 Flow chart of technical structure

2.1 ESTARFM时空数据融合模型

本文采用的ESTARFM模型是在STARFM模型基础上改进的。假设在理想情况下：Landsat和MODIS影像之间的反射率差异只存在系统偏差，并且2个时间的影像没有产生较大的差异，此时ESTARFM模型预测融合影像的公式为

F(x,y,tp,B)=F(x,y,t0,B)+a(C(x,y,tp,B)-
C(x,y,t0,B))

(1)

式中F——Landsat影像的反射率

(x,y)——像元位置坐标

B——影像参与计算的波段

C——MODIS影像的反射率

tp、t0——影像获取的时间

a——进行线性回归时的转换系数，由传感器之间的系统偏差确定

由于MODIS影像中含有大量的混合像元，此时式(1)不成立。文献[22]提出假设将混合像元的反射率建模表示为该像元中存在的不同土地覆盖成分反射率的线性组合，其权重由其像元内混合土地覆盖的区域覆盖率决定。此时模型预测公式为

F(x,y,tp,B)=F(x,y,t0,B)+v(x,y)(C(x,y,tp,B)-
C(x,y,t0,B))

(2)

式中v(x，y)——对Landsat和MODIS影像同一像元进行线性回归得到的转换系数

考虑到相邻的同类像元具有相似的反射率变化，文献[16]提出一种移动窗口的方法，移动窗口用于搜索窗口内的相似像元，然后将相似像元的信息集成到所需预测的Landsat反射率计算中，此时模型预测公式为

F(xw/2,yw/2,tp,B)=F(xw/2,yw/2,t0,B)+

(3)

式中 (xw/2,yw/2)——搜索窗口中心坐标

N——包括中央预测像元的相似像元的数目

Wi——第i个相似像元的权重

Vi——第i个相似像元的转换系数

(xi,yi)——第i个相似像元的坐标

选择2个不同时刻m和n的Landsat和tp时刻的MODIS影像进行预测融合，用m时刻的Landsat影像和在P时刻的MODIS影像预测P时刻的Landsat影像反射率，用n时刻的观测数据预测P时刻的Landsat影像反射率。结果分别记为Fm(xw/2,yw/2,tp,B)和Fn(xw/2,yw/2,tp,B)。通过2种预测结果的加权组合，可以获得更精确的P时刻反射率。因此在这种情况下，为Landsat影像输入设置更大的时间权重是合理的。因此模型公式为

(4)

最终模型预测P时刻的Landsat影像反射率计算公式为

F(xw/2,yw/2,tp,B)=TmFm(xw/2,yw/2,tp,B)+
TnFn(xw/2,yw/2,tp,B)

(5)

2.2 融合影像数据的时序NDVI曲线

归一化植被指数(NDVI)是用于表达植被生长状况、覆盖情况的参数。本文利用2015年已有的质量较好的Landsat影像和ESTARFM融合预测出其他对应时间的Landsat融合影像，通过ENVI 5.3的波段运算，可以得到全年23期的NDVI数据。由于原始NDVI时序数据受噪声和空值的影响较大[23]，造成NDVI时序曲线出现不规则的波动，所以需要对NDVI时间序列数据进行平滑降噪处理。本实验利用TIMESAT软件对构建的NDVI时序数据进行滤波，该软件有非对称高斯函数拟合、双Logistic函数拟合和S-G滤波法，在本次研究中，通过多次实验对比发现，S-G滤波方法基于滑动窗口的平均滤波，以局部拟合为主，有较强的细节拟合能力，能有效地去除时间序列数据中的噪声，保留原始时序曲线的有效信息[24-25]。因此本文采用S-G滤波方法对23期的NDVI数据进行滤波得到时序NDVI曲线。

图2 真实Landsat影像与ESTARFM模型融合影像以及STARFM模型融合影像对比Fig.2 Comparisons of real Landsat image with ESTARFM model fusion image and STARFM model fusion image

3 结果分析

3.1 ESTARFM时空融合模型结果与精度评价

根据实验数据，对ESTARFM模型运行结果进行精度评价，并与STARFM模型运行结果进行对比。把第209天和第257天的Landsat数据和MODIS13Q1数据、第225天的MODIS13Q1数据作为ESTARFM模型的输入数据；把第257天的Landsat数据和MODIS13Q1数据、第225天的MODIS13Q1数据作为STARFM模型的输入数据；融合预测出第225天的ESTARFM融合数据和STARFM融合数据；以此融合数据与实际第225天的Landsat数据做对比，以显示ESTARFM模型的预测精度。影像数据以R(Blue)G(Red)B(NIR)假彩色合成。

选取其中3块对比明显且具有代表性的区域(建筑物、水体、植被、农作物等)进行对比，如图2所示。在目视解译的基础上，ESTARFM模型时空融合的效果相比STARFM模型效果更好，各种地物大体上纹理清晰，轮廓特征比较准确，与原始影像在空间上基本保持一致。在图2中标记区域可以看出，相对于STARFM模型，ESTARFM模型的融合结果更加准确，与真实Landsat影像更具有一致性。只有小部分纹理比较粗糙，这是由于在构建融合数据时，使用不同时刻Landsat数据和MODIS数据存在时相的差异，这会造成融合影像与真实影像间的偏差。

为进一步验证ESTARFM模型的融合预测效果，选取融合出的第225天的影像，分3个波段与原始影像进行对比，获得散点图(图3)，从图3看出散点基本都在1∶1趋势线的周围，计算真实影像与融合影像在蓝波段、红波段和近红外波段的决定系数R2分别为0.824 6、0.920 3、0.940 2。决定系数R2都在0.82以上，且在红波段和近红外波段达到0.92以上，表明真实反射率与融合预测的反射率相关程度较高，因此证实了ESTARFM模型融合预测精度较高，能取得较好的融合效果，可以作为研究区水稻面积提取的数据依据。

3.2 时序NDVI曲线

结合实地调研与Google Earth影像目视解释，将原阳县的主要地物分为6类：建筑用地、裸地、水体、林地、小麦-水稻以及小麦-玉米。对计算得到的2015年内23个时相的NDVI数据，利用TIMESAT软件中的S-G滤波函数对时序NDVI数据进行平滑去噪处理。利用实地调研样点，得到每种地物的NDVI时序曲线，如图4所示。

图4 主要地物时序NDVI曲线Fig.4 Time series NDVI curves of main features

从图4可以看出：主要作物为小麦-水稻轮作与小麦-玉米轮作，一年有2个波峰，在5月中下旬小麦成熟，此时收割小麦造成在儒略日第145天左右NDVI出现下降至波谷，并且由于水稻处于生长初期，稻田中含有大量的水，因此稻田的NDVI相比玉米的NDVI呈现出更低的现象；随后水稻和玉米逐渐生长，在7月下旬至8月上旬NDVI呈现波峰。因此说明ESTARFM模型融合影像得到的时序NDVI曲线与实际作物NDVI变化情况相一致，能够较好地反映地物的季节性变化情况，为后续提取水稻及其他地物信息提供数据基础。

3.3 主要地物类别的NDVI阈值设置

根据各个地物NDVI曲线显示不同的特征，本文利用设置NDVI阈值的方法，建立分层决策树的分类方法对主要地类进行提取。具体算法如图5所示。

图5 分类决策树Fig.5 Classification decision tree

(1)由图4的NDVI时序曲线可知，植被和非植被的曲线特征差别比较明显。非植被地物全年趋于平缓，且全年的NDVI相对于植被的NDVI明显呈现较低的水平。在第97天时植被与非植被差异最明显，具有较好的区分性。因此利用第97天的NDVI小于0.45作为植被与非植被用地的阈值。

(2)水体的NDVI一年中大部分都为负值，虽然有少数水体像元在第145天呈现出正值，但都小于0.10，因此用第145天的NDVI小于0.10作为水体和建筑用地、裸地的阈值。

(3)建筑用地和裸地的时序NDVI曲线大体趋势相似，都呈现为较低值的水平。但是裸地在春夏交替时期，呈现出增长的趋势，原因是裸地在这个时期会生长出一些野生绿色植物。因此可以用第161天的NDVI是否小于0.35作为判别建筑用地和裸地的阈值。

(4)林地与农作物的NDVI曲线相比，林地的峰值持续时间较长，呈现峰值后，NDVI没有出现较大的波动，与水稻、玉米等农作物的时序NDVI曲线有明显的区别，主要农作物水稻和玉米的NDVI在5月下旬至6月上旬呈现最低值，在7月下旬至8月上旬呈现峰值，因此利用波峰与波谷的差值区分开林地与农作物。分别选取100个林地与农作物的样点，求取NDVI平均值，计算第209天与第161天的NDVI差值，经反复实验，对比发现，利用第209天与第161天的NDVI差值小于0.22时分类效果最好。

(5)由图4可以看出，水稻、玉米在生长初期的NDVI较低，并且水稻在生长初期田间有大量的水，NDVI比玉米更低，在7月下旬至8月上旬同时呈现峰值，且水稻的NDVI峰值比玉米略高，并且在9月下旬至10月上旬，水稻和玉米成熟，NDVI迅速下降，但下降后，水稻的NDVI比玉米的高，玉米的NDVI基本都在0.20以下，因此利用第225天与第289天NDVI的差值大于0.56作为判别玉米和水稻的阈值，第225天与第289天NDVI的差值小于0.20的判别为其他作物。

3.4 分类结果及精度验证

根据研究区内水稻与其他主要地物的时序NDVI曲线特征，利用专家决策树分类的方法，对研究区进行土地利用分类，提取出水稻种植区域。具体的分类结果如图6所示。

图6 决策树分类结果Fig.6 Decision tree classification results

利用实地采集和在Google Earth上选取的各个地类的样点共2 560个，通过分类后的混淆矩阵对分类结果进行定量分析(表3)，总体分类精度为92.23%，Kappa系数为0.904 3。水稻、玉米、建筑用地的制图精度和用户精度均在90%以上，且水稻的制图精度达到96.73%，用户精度也达到了93.51%。由于林地分布比较分散且与其他地类相连紧密，造成林地的混合像元问题严重，制图精度和用户精度都小于90%。

表3 分类结果精度评价Tab.3 Classification result accuracy evaluation

根据原阳县2015年统计年鉴对分类结果进行定性分析，并且向当地农户访问，原阳县大面积种植水稻的区域集中在祝楼乡、太平镇和葛埠口乡，与分类结果一致。统计年鉴记录原阳县2015年的水稻种植面积为20 468 hm2，分类结果中水稻面积为18 556.38 hm2，提取水稻面积与统计数据的一致性为90.66%。

4 讨论

不同地物的时序NDVI曲线具有不同的特征，因此可以利用水稻的时序NDVI将其与其他地物区分。中高空间分辨率的遥感影像在农作物识别提取中作为重要的数据源，但是由于受云雨天气的影响，经常无法获取连续的中高空间分辨率影像，因此利用ESTARFM模型等数据融合技术，可以解决遥感数据缺失的问题，作为一种对农作物识别提取的有效手段。但是，时空融合模型构建数据的品质受基准Landsat和MODIS影像数据的影响，获取相近时刻并且高品质的Landsat和MODIS影像是构建高品质时序融合数据的关键，但是这在部分地区是比较困难的。对于零碎分散地块，融合得到的30 m分辨率的Landsat数据依然存在明显的混合像元问题，导致提取精度不高。另外本次研究只选取了单一的时序NDVI指数进行分类提取，可能对有些地物反映不明显。

因此，在以后的研究中，利用ESTARFM模型融合更高空间分辨率影像和更高时间分辨率影像，进一步解决混合像元严重的问题，并且能监测出短期内地物发生快速变化的情况。另外，可以开展多种指数共同对地物提取，如EVI/NDWI指数等。探讨多种指数对水稻和其他农作物识别能力的适用性及准确性。

5 结论

(1)ESTARFM时空融合模型能较好地融合出缺失时相的Landsat影像，结果显示，空间细节信息比较清晰，在Red波段和NIR波段，与真实Landsat的Red波段和NIR波段一致性程度较高，决定系数R2分别达到0.920 3和0.940 2，因此NDVI的一致性同样较高。说明时空融合模型能够较好地重构高空间分辨率数据，有效解决多云多雨地区受天气影响的Landsat数据缺失问题。

(2)利用TIMESAT软件对研究区时序影像进行S-G滤波平滑处理，得到主要地类的时序NDVI曲线，判别出水稻以及其他主要地类的不同特征，按照一定的阈值范围设计决策树进行分类，有较好的分类效果。总体分类精度为92.23%，Kappa系数为0.904 3；水稻的制图精度为96.73%，用户精度为93.51%；与统计数据的一致性为90.66%。分类精度较高，在技术上具有可行性。