出版日期: 2019-09-25
点击次数：
下载次数：
DOI: 10.11834/jrs.20197391
2019 | Volumn23 | Number 5
上一篇  |  下一篇
技术方法 
基于时间序列叶面积指数稀疏表示的作物种植区域提取
expand article info 王鹏新1,2 , 荀兰1,2 , 李俐1,2 , 王蕾1,2 , 孔庆玲1,2
1. 中国农业大学 信息与电气工程学院，北京 100083
2. 农业农村部 农业灾害遥感重点实验室，北京 100083
收稿日期: 2017-09-18; 预印本: 2018-05-27
基金项目: 国家重点研发计划(编号：2016YFD0300603-3)
第一作者简介: 王鹏新，1965年生，男，教授，研究方向为定量遥感及其农业应用。E-mail：wangpx@cau.edu.cn
通信作者简介: 李俐，1976年生，女，副教授，研究方向为微波遥感及其农业应用。E-mail：lilixch@cau.edu.cn

摘要

以华北平原黄河以北地区为研究区域，以时间序列叶面积指数LAI(Leaf Area Index)傅里叶变换的谐波特征作为不同作物识别的数据源，利用稀疏表示的分类方法识别2007年—2016年冬小麦、春玉米、夏玉米等主要农作物种植区域。首先利用上包络线Savitzky-Golay滤波分别对2007年—2016年的时间序列MODIS LAI曲线进行重构，进而对重构的年时间序列LAI进行傅里叶变换，以0—5级谐波振幅、1—5级谐波相位作为作物识别的依据，基于各类地物的训练样本，通过在线字典学习算法构建稀疏表示方法的判别字典，对每个待测样本利用正交匹配追踪算法求解稀疏系数，从而计算对应于各类地物的重构误差，根据最小重构误差判定待测样本的作物类型，并对作物识别结果的位置精度进行验证。结果表明，2007年—2016年作物识别的总体精度为77.97%，Kappa系数为0.74，表明本文提出的方法可以用于研究区域主要作物种植区域的提取。

关键词

遥感, 叶面积指数, 稀疏表示, Savitzky-Golay滤波, 华北平原, 农作物, 识别

Extraction of planting areas of main crops based on sparse representation of time-series leaf area index
expand article info WANG Pengxin1,2 , XUN Lan1,2 , LI Li1,2 , WANG Lei1,2 , KONG Qingling1,2
1.College of Information and Electrical Engineering, China Agricultural University, Beijing 100083, China
2.Key Laboratory of Remote Sensing for Agri-Hazards, Ministry of Agriculture and Rural Affairs, Beijing 100083, China

Abstract

Crop mapping is an important component of agriculture monitoring. Accurate information on crop area coverage is vital for food security and the agricultural industry, and the demand for timely crop mapping is high. Previous research indicated that remote sensing technology is a practical and feasible method for agricultural crop area extraction. In this study, the north area of the Yellow River in the North China Plain is chosen as the study area, where the main crops are winter wheat, maize, cotton, and soybean. To obtain the distribution information of crops, the yearly four-day composite MODIS time-series Leaf Area Index (LAI) with 500 m spatial resolution is collected. A total of 92 MODIS LAI images obtained yearly from 2007 to 2016 are used to build time-series LAI curves. To avoid the edge effect of the time-series LAI caused by the Savitzky–Golay filter, the last two phases of LAI images in the last year and the first two phases of LAI images in the next year are added to build the time-series LAI in a year. The Savitzky–Golay filter is then applied on the yearly time-series LAI pixel by pixel to minimize effects of anomalous values caused by atmospheric haze, cloud contamination, and so on. Fourier transform method based on reconstructed LAI is further employed to extract the key parameters. The 11 parameters, including the amplitudes of 0–5 terms and the phases of 1–5 terms, are taken as the features for crop identification. The training samples and verification samples of various crops are obtained through ground investigation and Google Earth images. On the basis of the training samples of various crops, online dictionary learning algorithm is applied to construct the dictionary used to identify the crops. With the dictionary, the orthogonal matching pursuit algorithm is further applied on samples under testing to obtain the sparse representation coefficient. Then the crops are identified according to the minimum reconstruction error, which can be calculated by the dictionary and the coefficient. Therefore, the areas planting winter wheat, spring maize, summer maize, cotton, and orchard from 2007 to 2016 are extracted in the study area. Lastly, the accuracy of the identification results is evaluated yearly by a confusion matrix. Results show that the reconstructed time-series LAI curves are smooth and consistent with crop growth and development characteristics. Overall identification accuracy reaches 77.97% with a Kappa coefficient of 0.74 from 2007 to 2016. User accuracies for individual crops are as follows: winter wheat and summer maize, 90.60%; spring maize, 73.40%; early summer maize, 81.80%; cotton, 69.40%; and orchard, 81.60%. Annual overall accuracies from 2007 to 2016 range between 70.57% and 83.71% and Kappa coefficients range from 0.66 to 0.81. In conclusion, combining the harmonic characteristics of the time-series LAI with the sparse representation can effectively identify the areas for planting different crops. The approach developed in this study is feasible for extracting information on main crop distribution in the study area.

Key words

remote sensing, Leaf Area Index, sparse representation, Savitzky-Golay filter, North China Plain, crops, identification

1 引　言

作物空间分布信息的及时获取对农业生产管理、农业可持续发展及国家粮食安全具有重要意义，是区域作物产量估测、农作物种植结构调整和优化的主要依据(Belgiu和Csillik，2018；Whelen和Siqueira，2017)。作物类型的遥感识别主要依据作物的光谱特征差异和物候期差异(Conrad 等，2011；Pan 等，2012；Jia 等，2014)。相比于单一时相的遥感数据，多时相遥感数据可以反映不同作物整个生育期的生长状况，具有更强的区分不同作物的能力。近年来，利用时间序列归一化植被指数NDVI(Normalized Difference Vegetation Index)、增强型植被指数EVI(Enhanced Vegetation Index)数据进行作物种植区域遥感提取的研究较多，然而在高密度和多层植被情况下 NDVI 会出现饱和现象。时间序列叶面积指数LAI(Leaf Area Index)能够描述作物年内或年际生长变化特征，体现其生长变化过程，然而，基于遥感数据获取大区域的LAI需要大量的地面实测LAI数据。MODIS卫星数据已经进行多年的观测，由于其重复周期短、覆盖面积广等特点，适用于进行长期、大面积的作物监测，被广泛应用于作物提取工作中(Wardlow和Egbert，2008；Arvor 等，2011；Skakun 等，2017)。通过对时间序列MODIS LAI进行分析，可以实现不同作物类型的识别(方墨人和田庆久，2004)。现有研究中，大多采用监督分类和非监督分类的方法提取作物种植区域，且需要利用地面实测数据作为训练样本和确定输出类别的依据，然而在进行多年作物种植区域提取的研究时，地面实测数据的获取耗时耗力，因此，在保证分类精度的基础上，如何减少对地面实测数据的依赖、提高数据处理效率是至关重要的(Massey 等，2017)。

近年来，稀疏表示理论成为国内外的研究热点。Mallat和Zhang(1993)基于小波分析提出了利用冗余字典对信号进行稀疏分解的思想，并将其引入匹配追踪算法，使分解得到的表达结果十分简洁，也为后续理论的发展奠定了基础。稀疏表示的本质是用尽可能少的非零系数高效、简洁地表示待求解向量，即从过完备字典当中找到一种求解系数存在尽可能多的零的向量表示方式，从而利用稀疏的表达形式揭示信号的主要特征和结构，简化信号处理过程(Candes和Wakin，2008)。对于绝大部分元素为零的矩阵，稀疏表示方法只需要处理矩阵的非零元素便可以揭示信号的内在属性，大大降低了矩阵运算的难度和计算量。因此，将稀疏表示应用于遥感图像处理，可以提高数据处理效率。

稀疏表示在遥感图像处理中的应用主要集中在图像融合(尹雯 等，2013)、图像去噪(夏琴 等，2016)、超分辨率遥感图像重构(刘帅 等，2015)、高光谱遥感图像分类(宋相法和焦李成，2012；宋琳 等，2012)等方面。Wright等(2009)提出基于稀疏表示的分类算法SRC(Sparse Representation based Classification)并用于人脸识别，SRC的主要思想是利用训练样本的一种线性组合拟合测试样本，然后根据拟合结果和测试样本计算重构误差，最后将测试样本划分为最小重构误差对应的类别中。宋相法和焦李成(2012)利用高光谱遥感图像数据集构造学习字典，根据学习字典计算每个像素的稀疏系数，从而获得像素的稀疏表示特征，最后根据稀疏表示特征和光谱信息分别构造随机森林，通过投票机制获得分类结果，结果表明，该方法的分类总体精度和Kappa系数均高于光谱信息方法和稀疏表示特征方法。唐晓晴等(2016)在稀疏表示分类模型中引入局部二值模式LBP(Local Binary Pattern)特征，使用LBP对遥感图像进行特征提取，获得遥感图像的局部纹理特征，根据LBP直方图训练结构化字典，建立了基于稀疏表示的遥感图像分类模型，结果表明，与传统的基于LBP的支持向量机和基于LBP的K最近邻方法相比，该方法能够有效提高分类精度。

本文以华北平原黄河以北地区为研究区域，以时间序列MODIS LAI傅里叶变换后的谐波特征作为不同作物识别的数据源，选取训练样本通过在线字典学习算法构建判别字典，根据字典计算待测样本的稀疏系数，通过最小重构误差确定作物种植类型，从而识别冬小麦、春玉米、夏玉米、棉花和果园种植区域，并对识别结果的位置精度进行验证。

2 材料与方法

2.1 研究区域概况

华北平原位于30°00′N—40°24′N、112°48′E—122°45′E，北起燕山，西沿太行山、伏牛山，南抵淮河干流及苏北灌溉总渠，东邻渤海与黄海，涵盖北京、天津、河北、山东、河南、江苏和安徽五省二市的大部或部分地区，主要由黄河、淮河、海河及其支流冲积而成(莫兴国 等，2011)。该地区海拔多低于50 m，多年平均气温为10.0—14.2℃，年均降水量为500—900 mm，属于温带大陆性季风气候，是我国重要的粮食作物和经济作物产区，适宜种植小麦、玉米、棉花等多种农作物(王鹏新 等，2017；王学 等，2015；唐鹏钦 等，2011)。本文以华北平原黄河以北地区(34°48′N—40°27′N、113°11′E—119°49′E)为研究区域进行作物识别(图1)。

2.2 数据来源及预处理

通过NASA网站(http://reverb.echo.nasa.gov/reverb/[2017-09-18])获取研究区域2007年—2016年的MCD15A3H产品(轨道号分别为h26v04、h26v05、h27v04和h27v05)，包括LAI和光合有效辐射，其空间分辨率为500 m，时间分辨率为4 d，每年共92个时相的LAI影像。对MODIS LAI影像进行影像镶嵌和投影转换等预处理，根据研究区域矢量图进行裁剪，构建研究区域每年92个时相的时间序列LAI数据。为避免利用Savitzky-Golay滤波对时间序列LAI进行处理时造成的边缘效应，在每年时间序列LAI的首、尾分别增加前一年第91和第92时相、后一年第1和第2时相的LAI，从而获得研究区域每年96个时相的LAI影像，即每年每个像素对应包含96个LAI的变化曲线，对其进行Savitzky-Golay滤波后采用当年92个时相的时间序列LAI进一步进行傅里叶变换分析(王鹏新 等，2017)。


图 1 研究区域MODIS影像(NIR/R/G)
Fig. 1 MODIS image of study area (NIR/R/G)

根据中国农业信息网提供的农作物种植信息，确定冬小麦、春玉米、夏玉米、棉花和果园为本文主要识别的作物。根据播种时间的早晚，将夏玉米分为早播夏玉米与轮作夏玉米(即冬小麦—夏玉米轮作)。早播夏玉米在白地(未种植作物的耕地)播种或在蔬菜/豌豆收获后播种。根据地面调查发现，在MODIS 500 m像素覆盖范围内，当冬小麦—夏玉米分布较少，混有其他地物(如秋收作物、村庄、厂房等)较多时，时间序列LAI曲线仍呈现为双峰，但第一个波峰不显著，因此，本文将一年两季作物中第一个峰值LAI小于1.5的像素对应区域定义为冬小麦—夏玉米&其他地物种植区域。最终确定冬小麦—夏玉米、冬小麦—夏玉米&其他地物、春玉米、早播夏玉米、棉花、果园、其他为本文的主要识别地物类型，城市、林地、草地、水体和荒地区域通过结合MCD12Q1产品及MCD15A3H产品的填充值区域进行识别。结合野外实地调查数据和研究区域Google Earth影像，根据LAI曲线的峰值个数、峰值出现时间、曲线斜率最大值出现时间及不同作物的物候特征确定700个样点，每类选取50个样点作为训练样本用于字典学习，其余50个样点用于分类精度验证，采用混淆矩阵的方式对本文的作物识别结果进行精度验证，评价指标包括总体精度、用户精度、制图精度和Kappa系数4项(王鹏新 等，2017)。

2.3 时间序列LAI的Savitzky-Golay滤波

由于云和大气等噪声的影响，时间序列LAI曲线呈锯齿状的不规则波动，导致作物LAI曲线反映在季节上的变化趋势不明显，影响了信息提取的精度(吴文斌 等，2009)。Savitzky和Golay(1964)提出的Savitzky-Golay滤波是一种移动窗口的加权平均算法，其加权系数通过在滑动窗口内对给定高阶多项式的最小二乘拟合得出。利用Savitzky-Golay滤波对遥感时间序列数据进行重构可以去除云和大气等噪声的影响，本文将其用于时间序列LAI的重构，表示为

$ {\rm{LAI}}_t^a = \sum\limits_{h = - g}^{h = g} {\frac{{{C_h}{\rm{LAI}}_{t + h}^0}}{{2g + 1}}} $ (1)

式中， $t$ 表示LAI影像的时相，取值范围为1—92； $g$ 为滑动窗口的半径，设为2； ${\rm{LAI}}_{t + h}^0$ 为滑动窗口内第h个原始LAI； ${C_h}$ 为窗口内第h个原始LAI的滤波系数，通过查表(Savitzky和Golay，1964)获得； ${\rm{LAI}}_t^a$ 为S-G滤波后的LAI(王鹏新 等，2017)。

本文采用上包络线S-G滤波(Chen 等，2004)通过多次循环迭代对时间序列MODIS LAI进行平滑处理，获得MODIS LAI的上包络线 ${\rm{LAI}}_t^b$ 。

2.4 基于傅里叶变换的特征提取

傅里叶变换可以把时间序列数据从时域空间转换到频域空间，将时间序列LAI分解为一系列由不同频率的余弦曲线构成的谐波曲线，且时间序列曲线的大部分信息集中在前几级谐波中，谐波由振幅和初始相位确定(梁守真 等，2011)。对S-G滤波后的时间序列 ${\rm{LAI}}_t^b$ 进行傅里叶变换

$ {F_{\textit{z}}} = \frac{1}{{92}}\sum\limits_{t = 0}^{91} {{\rm{LAI}}_t^b\exp \left(- \frac{{2{\text{π}} i{\textit{z}}t}}{{92}}\right)} $ (2)

式中， ${F_{\textit{z}}}$ 为 ${\rm{LAI}}_t^b$ 的第z个傅里叶系数，z取值范围为0—91；i为虚数单位，i2=−1。

2.5 基于稀疏表示的作物识别

假设有 $M$ 类待提取的作物种植类型，第 $m$ 类作物有 ${N_m}$ 个训练样本 ${\left\{ {{{d}}_j^m} \right\}_{j = 1,\cdots,{N_m}}}$ ，其中， ${{d}}_j^m$ 表示由第 $m$ 类作物中第 $j$ 个训练样本的 $w$ 个谐波特征参数构成的列向量。设 ${{y}}$ 为由 $w$ 个谐波特征参数表示的待测样本，若 ${{y}}$ 属于第 $m$ 类作物，则 ${{y}}$ 可以由这些训练样本的线性组合近似表示(Chen 等，2011)

$\begin{split} {{y}} \approx & {{x}}_1^m{{d}}_1^m + {{x}}_2^m{{d}}_2^m + \cdots + {{x}}_{{N_m}}^m{{d}}_{{N_m}}^m = \\ & \left[ {{{d}}_1^m,\;{{d}}_2^m,\; \cdots\!,\; {{d}}_{{N_m}}^m} \right]{\kern 1pt} {\kern 1pt} {\left[ {{{x}}_1^m,\;{{x}}_2^m,\; \cdots\!,\; {{x}}_{{N_m}}^m} \right]^{\rm{T}}}= {{{D}}^m}{{{x}}^m} \end{split}$ (3)

式中， ${{{D}}^m}$ 为 $w \times {N_m}$ 维的子字典，由第 $m$ 类作物的训练样本构成，每一列 ${{d}}_j^m$ 为 ${{{D}}^m}$ 的一个原子； ${{{x}}^m}$ 为 ${N_m} \times 1$ 维的稀疏向量，它的每一个元素为 ${{{D}}^m}$ 中对应原子的权重。

因此，待测样本 ${{y}}$ 位于由 $M$ 个子字典 ${\left\{ {{{{D}}^m}} \right\}_{m = 1,\cdots,M}}$ 构成的集合中

$\begin{split} {{y}} \approx & {{{D}}^1}{{{x}}^1} + {{{D}}^2}{{{x}}^2} + \cdots + {{{D}}^M}{{{x}}^M} = \\ & \left[ {{{{D}}^1},\; \cdots,\; {{{D}}^M}} \right]{\left[ {{{{x}}^1},\; \cdots,\; {{{x}}^M}} \right]^{\rm{T}}} = {{Dx}} \end{split}$ (4)

式中， ${{D}}$ 为 $w \times N$ 维的结构化字典，包含N个训练样本， $N = \sum\limits_{m = 1}^M {{N_m}} $ ； ${{x}}$ 为 $N \times 1$ 维的稀疏向量，由 ${\left\{ {{{{x}}^m}} \right\}_{m = 1,\cdots,M}}$ 构成，理想情况下，若 ${{y}}$ 为第 $m$ 类作物，则 ${{{x}}^j} = 0$ ， $\forall j = 1,\cdots,m,{\kern 1pt} j \ne m$ 。

2.5.1 过完备字典的构造

稀疏表示要求字典具有过完备性，即字典中的原子个数要远大于样本的特征维度( $N \gg w$ )。本文通过调用Matlab的Spams工具箱，应用在线字典学习ODL(Online Dictionary Learning)算法(Mairal 等，2010)进行字典构建，该算法的具体过程为

(1)初始化 ${{{A}}_0} = 0({{{A}}_0} \in {{\bf{R}}^{k \times k}})$ ， ${{{B}}_0} = 0({{{B}}_0} \in {{\bf{R}}^{w \times k}})$ ，基于2015年的样本数据，每一类别随机选取50个训练样本，利用训练样本所在像素的时序LAI傅里叶变换的 $w$ 个谐波特征参数构成的列向量初始化字典 ${{{D}}_0} \in {{\bf{R}}^{w \times N}}$ ，令 ${{D}} = {{{D}}_0}$ ，将迭代次数 $q$ 初始化为1；

(2)利用最小角回归LARS(Least Angle Regression)算法(Efron 等，2004)计算经过 $q$ 次迭代的稀疏系数

$ {{\hat{ x}}_q}\buildrel \Delta \over = \mathop {\arg \min }\limits_{{{x}} \in {{\bf{R}}^{1 \times N}}} \frac{1}{2}\left\| {{{{y}}_q} - {{{D}}_{q - 1}}{{x}}} \right\|_2^2 + \gamma {\left\| {{x}} \right\|_1} $ (5)

式中， $\gamma \in {\bf{R}}$ ，为正则化参数，设置为0.15； ${\left\| {{x}} \right\|_1}$ 为 ${{x}}$ 的 ${l_1}$ −范数；

(3)计算 ${{{A}}_q} = {{{A}}_{q - 1}} + {{\hat{ x}}_q}{{\hat{ x}}_q}^{\rm{T}}$ ， ${{{B}}_q} = {{{B}}_{q - 1}} + {{{y}}_q}{{\hat{ x}}_q}^{\rm{T}}$ ；

(4)更新经过 $q$ 次迭代的字典 ${{{D}}_q}$

1) 初始化 $j = 1$ ；

2) 更新第 $j$ 列字典原子 ${{{d}}_j}$ ： ${{{u}}_j} = \dfrac{1}{{{{A}}\left[ {j,j} \right]}}({{{b}}_j} - {{D}}{{{a}}_j}) + $ $ {{{d}}_j}$ ， ${{{d}}_j} = \dfrac{1}{{\max ({{\left\| {{{{u}}_j}} \right\|}_2},1)}}{{{u}}_j}$ ， ${{{a}}_j}$ 、 ${{{b}}_j}$ 分别表示 ${{{A}}_q}$ 、 ${{{B}}_q}$ 中的第 $j$ 列；

3) 若 $j < N$ ，令 $j = j + 1$ ，返回步骤2)，否则，计算 ${{{D}}_q}\buildrel \Delta \over =\mathop {\arg \min }\limits_{{{D}} \in {{C}}} \dfrac{1}{q}\sum\limits_{s = 1}^q {(\left\| {{{{y}}_s} - {{D}}{{{\hat{ x}}}_s}} \right\|_2^2 + {\textit{λ}} {{\left\| {{{{\hat{ x}}}_s}} \right\|}_1})} =\mathop {\arg \min }\limits_{{{D}} \in {{C}}} \dfrac{1}{q}$ $(\dfrac{1}{2}{\rm{Tr}}({{{D}}^{\rm{T}}}{{D}}{{{A}}_q}) - {\rm{Tr}}({{{D}}^{\rm{T}}}{{{B}}_q})) $ ，其中， ${{C}} \buildrel \Delta \over = \{ {{{D}}_q} \in {{\bf{R}}^{\omega \times N}}\;$ ${\rm{ s.t.}}\;\forall \;j = 1,\cdots,N,\;\;{{d}}_j^{\rm{T}}{{{d}}_j} \leqslant 1\} $ ， ${\rm{Tr}}$ 为矩阵的际，即矩阵的对角线元素之和；

(5)设置最大迭代次数 ${q_{{\rm{max}}}}$ 为200，若 $q < 200$ ，令 $q = q + 1$ ，重复步骤(2)—(4)进行循环迭代；当 $q = 200$ 时，退出循环，字典 ${{{D}}_q}$ 即为所求。

2.5.2 稀疏系数的求解

对于待测样本 ${{y}}$ ，计算稀疏系数

$ {\hat{ x}} = \mathop {\arg }\limits_{} \min {\left\| {{{Dx}} - {{y}}} \right\|_2}, \;{\rm{s.t.}}\;{\left\| {{x}} \right\|_0} \leqslant {K_0} $ (6)

式中， ${\left\| {{x}} \right\|_0}$ 为 ${{x}}$ 的 ${l_0}$ −范数，即向量 ${{x}}$ 中非零元素的个数； ${K_0}$ 为稀疏度上界。

然而上述问题为NP-hard(Nondeterministic Polynomial-time hard)问题，利用传统的方法很难得到精确的解，本文通过调用Matlab的Spams工具箱，利用贪婪算法中的正交匹配追踪OMP(Orthogonal Matching Pursuit)算法(Tropp和Gilbert，2007)求解。

OMP算法已知 $w \times N$ 维字典 ${{D}} = \left[ {{{{d}}_1},\;{{{d}}_2},\;\cdots,\;} \right.$ $\left.{{{d}}_N}\right] $ ，求解 $w \times 1$ 维待测样本 ${{y}}$ 对应于 ${{D}}$ 的稀疏表示系数 ${\hat{ x}}$ ，具体过程为

(1)初始化残差 ${{{r}}_0} = {{y}}$ ，索引集 ${{{\varLambda }}_0} = {\rm{\phi }}$ ，将迭代次数 $q$ 初始化为1，将 ${K_0}$ 设置为10；

(2)选取经过 $q$ 次迭代时与残差最为匹配的原子索引 ${{{\textit{λ}}}_q} = \arg \mathop {\max }\limits_{v = 1,\cdots,N} \left| {\left\langle {{{r}}_{q - 1}^{},{{{d}}_v}} \right\rangle } \right|$ ；

(3)更新索引集和原子集： ${{{\varLambda }}_q} = {{{\varLambda }}_{q - 1}} \cup \left\{ {{{\textit{λ}}_q}} \right\}$ ， ${{{D}}_q} = \left[ {{{{D}}_{q - 1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{{d}}_{{{\textit{λ}}_q}}}} \right]$ ， ${{{D}}_0}$ 为空矩阵；

(4)利用LARS算法(Efron 等，2004)计算稀疏系数 ${{\hat{ x}}_q}$ ，同式(5)；

(5)更新残差 ${{{r}}_q} = {{y}} - {{{D}}_q}{{\hat{ x}}_q}$ ；

(6)若 $q < {K_0}$ ，令 $q = q + 1$ ，重复步骤(2)—(5)进行循环迭代；当 $q = {K_0}$ 时，退出循环，输出 $N$ 维稀疏系数 ${{\hat{ x}}_q}$ 。

2.5.3 作物类型识别

待测样本 ${{y}}$ 所代表的作物类型可以直接根据OMP算法所求解的稀疏系数 ${\hat{ x}}$ 确定，定义 ${{y}}$ 被判别字典中第 $m$ 类作物对应的字典原子构成的集合 ${{{D}}^m}$ 所重构的误差为

$ {r^m}({{y}}) = {\left\| {{{y}} - {{{D}}^m}{{{\hat{ x}}}^m}} \right\|_2},\;\;m = 1,2,\cdots,M $ (7)

式中， ${{\hat{ x}}^m}$ 为与 ${{{D}}^m}$ 相对应的稀疏系数。 ${{y}}$ 对应于重构误差 ${r^m}({{y}})$ 最小的作物类型。

3 结果与分析

3.1 主要作物的时间序列LAI曲线分析

利用上包络线S-G滤波分别对研究区域2007年—2016年每一像素的时序MODIS LAI曲线进行平滑处理，获得重构的年时间序列LAI(以某一冬小麦—夏玉米地面调查样点所在像素的2016年时序LAI曲线为例，如图2(a))。分析2016年各类待识别作物的重构LAI曲线(图2(b))，可以发现，冬小麦—夏玉米种植区域的时间序列LAI曲线有2个显著峰值，分别对应于冬小麦抽穗期和夏玉米抽雄期，第290天时LAI开始呈现小幅度增长，至第330天形成一个微小的波峰，对应于冬小麦的分蘖期。冬小麦—夏玉米&其他地物种植区域的时间序列LAI曲线有2个峰值，且LAI峰值较小，第一个峰值不显著，根据地面调查，原因可能在于该点所在的 500 m×500 m 像素覆盖范围内有厂房、河道或其他秋收作物(如春玉米)。冬小麦—夏玉米与春玉米交叉种植时，LAI曲线的第一个峰值较小，第二个峰值较大，第二个峰值反映了春玉米、夏玉米的抽雄期。春玉米、夏玉米、棉花种植区域的时间序列LAI曲线均有1个显著峰值，且春玉米的峰值出现时间早于夏玉米，棉花的峰值持续时间长于春玉米和夏玉米。果园种植区域的LAI在4月底至9月初一直保持一年中的较高值，时间跨度较长。因此，根据不同作物类型的时间序列LAI可以进行作物种植区域的提取。


图 2 研究区域主要作物的时间序列LAI曲线
Fig. 2 Time series LAI curves of main crops in study area

利用快速傅里叶变换分别对研究区域2007年—2016年S-G滤波后的LAI曲线进行分解，获得各级谐波曲线。以图2(b)中各类作物的LAI曲线为例，进行傅里叶变换分析(图3，表1)。可以发现，在不同作物种植模式下，0级谐波振幅均较大，反映了年时间序列LAI的均值，其中，冬小麦—夏玉米种植区域的0级谐波振幅最大。在一年一季作物种植区域，随着谐波级数增加，谐波振幅整体上呈现为逐渐减小的趋势。春玉米的抽雄期早于早播夏玉米，该物候差异反映为二者的1级谐波初始相位不同，分别为2.68和2.57。在冬小麦—夏玉米和冬小麦—夏玉米&其他地物种植区域，1级和3级谐波振幅较大，2级谐波振幅较小。0—5级谐波叠加曲线与S-G滤波后的曲线较为接近，可以反映时间序列LAI的曲线特征。因此，选取傅里叶变换后的0—5级谐波振幅、1—5级谐波初始相位作为作物识别的特征参数。


图 3 S-G滤波LAI、傅里叶变换的0—5级谐波曲线及谐波叠加曲线
Fig. 3 Time series curves of S-G filtered LAI, harmonic terms 0—5 and their sum

表 1 研究区域主要作物的谐波参数
Table 1 Parameters of harmonic function for main crops in study area

下载CSV 
种植模式 谐波 振幅 初始相位 种植模式 谐波 振幅 初始相位
冬小麦—夏玉米 0 1.81 — 冬小麦—夏玉米&其他地物 0 1.13 —
1 0.81 −3.08 1 0.54 2.71
2 0.19 −2.58 2 0.20 −1.41
3 0.78 0.70 3 0.27 0.86
4 0.51 −2.28 4 0.13 −2.74
5 0.15 1.55 5 0.02 1.28
春玉米 0 1.14 — 早播夏玉米 0 1.08 —
1 0.81 2.68 1 0.71 2.57
2 0.50 −0.92 2 0.48 −1.37
3 0.30 1.92 3 0.41 1.14
4 0.18 −1.26 4 0.28 −2.45
5 0.14 1.90 5 0.19 0.28
棉花 0 0.87 — 果园 0 1.19 —
1 0.58 2.53 1 0.61 2.83
2 0.31 −1.36 2 0.04 −1.05
3 0.16 1.14 3 0.09 −0.56
4 0.05 −2.35 4 0.01 −2.18
5 0.01 2.99 5 0.04 2.38

3.2 主要作物的空间分布

分析研究区域2007年—2016年主要作物的空间分布特征(图4)可以发现，研究区域主要为一年一季和一年两季的作物种植模式，且一年两季作物种植区域大于一年一季作物种植区域，并存在明显的空间分异规律。研究区域中部及南部以一年两季的种植模式为主，研究区域中北部以一年一季的种植模式为主。小麦和玉米为研究区域的主要粮食作物，2007年—2016年研究区域中冬小麦和玉米的种植面积呈现波动变化。棉花种植区域总体呈现为减少的趋势，这与近年来国家政策引导棉花种植向新疆集中的状况一致。河北省赵县、藁城市、晋州市一带以及深州市的果园种植面积较为稳定，这与地面调查中果园的分布情况一致，该地区有悠久的果树种植历史，赵县一带主要种植梨树，深州市主要种植桃树。


图 4 2007年—2016年主要农作物提取和地表分类结果
Fig. 4 Classification results of land use, including main crops from 2007 to 2016

3.3 精度验证

每一类随机选取50个样点采用混淆矩阵的方式分别对2007年—2016年作物提取结果的位置精度进行验证，结果表明，每年的分类总体精度介于70.57%—83.71%，Kappa系数介于0.66—0.81。将10年的混淆矩阵进行整合(表2)，2007年—2016年的分类总体精度为77.97%，Kappa系数为0.74，说明利用遥感时序数据结合稀疏表示的方法在较长时间尺度内提取作物的精度较为稳定，提取结果可较为客观地反映区域作物的空间分布信息。冬小麦—夏玉米种植区域的用户精度最高，为90.60%。早播夏玉米和果园种植区域的用户精度较高，分别为81.80%和81.60%，果园种植区域的制图精度达到最高，为97.61%。春玉米、冬小麦—夏玉米&其他地物和棉花种植区域的用户精度较低，分别为73.40%、71.60%和69.40%。

表 2 混淆矩阵
Table 2 Confusion matrix

下载CSV 
分类数据 验证数据 用户精度/%
冬小麦—
夏玉米 		冬小麦—夏玉米&
其他地物 		春玉米 早播夏玉米 棉花 果园 其他 总计
冬小麦—夏玉米 453 32 0 6 0 2 7 500 90.60
冬小麦—夏玉米&其他地物 37 358 2 27 9 0 67 500 71.60
春玉米 9 31 367 13 25 0 55 500 73.40
早播夏玉米 13 24 5 409 9 0 40 500 81.80
棉花 4 26 14 12 347 1 96 500 69.40
果园 17 9 7 0 1 408 58 500 81.60
其他 8 42 5 33 18 7 387 500 77.40
总计 541 522 400 500 409 418 710 3500
制图精度/% 83.73 68.58 91.75 81.80 84.84 97.61 54.51
总体精度/% 77.97 Kappa系数 0.74

4 讨　论

本研究采用的MODIS LAI影像空间分辨率为500 m，这可能导致很多破碎耕地地块提取不出来，且耕地周边的村庄、厂房、道路、河道、树木、杂草等会对同一像素的时间序列LAI 曲线造成干扰，以致影响作物提取精度，将其用于农作物种植面积数据的精确获取尤为不足。对于生长期接近的不同作物，时间序列LAI曲线特征较相似，容易造成混淆。另外，小宗作物(如大豆、花生)的种植相对较为破碎，会混淆在大宗作物类别中，导致冬小麦、玉米的遥感提取面积可能被高估。表2中冬小麦—夏玉米种植区域的用户精度最高，主要原因是作物类型的识别精度主要取决于待测样本的谐波特征与判别字典中各类作物训练样本谐波特征的差异，冬小麦—夏玉米种植区域的时间序列LAI曲线有两个显著波峰，与其他作物类型相比，谐波特征差别较大。由于MODIS LAI 500 m空间分辨率的影响，冬小麦—夏玉米&其他地物所在像素内不仅有冬小麦—夏玉米，可能还种植有春玉米、早播夏玉米、棉花等作物，LAI曲线特征混淆相对较为严重，因此冬小麦—夏玉米&其他地物的制图精度较低。

本文基于LAI上包络线傅里叶变换的谐波特征进行作物识别，经S-G滤波重构后的LAI上包络线较好地反映了研究区域农田种植多熟制度下的耕地LAI年内动态变化趋势以及局部的突变信息，在此基础上，提取11个谐波特征参数，一定程度上减少了数据冗余。若直接基于LAI上包络线进行作物识别，则增加了对数据存储空间的需求，同时降低了数据处理效率。若直接基于原始MODIS LAI傅里叶变换后的谐波特征进行作物识别，云、大气等噪声的影响会导致LAI曲线特征不能更好地反映，不利于进行时间序列LAI曲线的变化趋势分析和信息提取。

现有研究中，稀疏表示在遥感分类中的应用大多基于高光谱影像(宋相法和焦李成，2012；宋琳 等，2012)，基于遥感时序数据结合稀疏表示方法进行分类的研究较少。与传统的监督分类法相比，本文提出的作物种植区域提取方法，不需要在每次分类时重新选择训练样本，可以利用某一年的数据选择待提取地物的训练样本用于字典构建，本文基于2015年的样本数据通过在线字典学习算法构建判别字典，根据判别字典批量化处理多年时间序列数据，进而提取多年作物种植区域，提高了数据处理效率。

5 结　论

本文利用稀疏表示的方法提取华北平原黄河以北地区2007年—2016年冬小麦—夏玉米、冬小麦—夏玉米&其他地物、春玉米、早播夏玉米、棉花和果园种植区域，分类总体精度为77.97%，Kappa系数为0.74，表明基于时间序列LAI傅里叶变换后的谐波特征利用稀疏表示的方法提取研究区主要作物的种植区域具有可行性。对于多年作物种植区域的提取，本文的方法相比于传统的监督分类法更具有优越性。

在今后的研究工作中需要进一步探索MODIS数据混合像素问题的有效解决方法，以降低混合像素的影响，从而在稀疏表示方法的基础上，进一步提高区域作物的识别精度。本文的研究对象是华北平原黄河以北地区的主要农作物，提出的作物种植区域提取方法有待于在整个华北平原地区内进行试验。

参考文献(References)

  • Arvor D, Jonathan M, Meirelles M S P, Dubreuil V and Durieux L. 2011. Classification of MODIS EVI time series for crop mapping in the state of Mato Grosso, Brazil. International Journal of Remote Sensing, 32 (22): 7847–7871. [DOI: 10.1080/01431161.2010.531783]
  • Belgiu M and Csillik O. 2018. Sentinel-2 cropland mapping using pixel-based and object-based time-weighted dynamic time warping analysis. Remote Sensing of Environment, 204 : 509–523. [DOI: 10.1016/j.rse.2017.10.005]
  • Candes E J and Wakin M B. 2008. An introduction to compressive sampling. IEEE Signal Processing Magazine, 25 (2): 21–30. [DOI: 10.1109/msp.2007.914731]
  • Chen J, Jönsson P, Tamura M, Gu Z H, Matsushita B and Eklundh L. 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the savitzky–golay filter. Remote Sensing of Environment, 91 (3/4): 332–344. [DOI: 10.1016/j.rse.2004.03.014]
  • Chen Y, Nasrabadi N M and Tran T D. 2011. Hyperspectral image classification using dictionary-based sparse representation. IEEE Transactions on Geoscience and Remote Sensing, 49 (10): 3973–3985. [DOI: 10.1109/tgrs.2011.2129595]
  • Conrad C, Colditz R, Dech S, Klein D and Vlek P L G. 2011. Temporal segmentation of MODIS time series for improving crop classification in Central Asian irrigation systems. International Journal of Remote Sensing, 32 (23): 8763–8778. [DOI: 10.1080/01431161.2010.550647]
  • Efron B, Hastie T, Johnstone I and Tibshirani R. 2004. Least angle regression. The Annals of Statistics, 32 (2): 407–499. [DOI: 10.1214/009053604000000067]
  • Fang M R and Tian Q J. 2004. The study of land type classification technology based on time-series LAI of MODIS. Remote Sensing for Land and Resources, 16 (3): 5–7, 22. [DOI: 10.6046/gtzyyg.2004.03.02] ( 方墨人, 田庆久. 2004. 基于MODIS的LAI时间序列谱的地物分类方法研究. 国土资源遥感, 16 (3): 5–7, 22. [DOI: 10.6046/gtzyyg.2004.03.02] )
  • Jia K, Liang S L, Wei X Q, Yao Y J, Su Y R, Jiang B and Wang X X. 2014. Land cover classification of Landsat data with phenological features extracted from time series MODIS NDVI data. Remote Sensing, 6 (11): 11518–11532. [DOI: 10.3390/rs61111518]
  • Liang S Z, Xing Q G, Shi P and Zhou D. 2011. Harmonic analysis on NDVI time series of typical land covers in Shandong Province. Chinese Journal of Ecology, 30 (1): 59–65. [DOI: 10.13292/j.1000-4890.2011.0006] ( 梁守真, 邢前国, 施平, 周迪. 2011. 山东省典型地表覆被NDVI时间序列谐波分析. 生态学杂志, 30 (1): 59–65. [DOI: 10.13292/j.1000-4890.2011.0006] )
  • Liu S, Zhu Y J and Xue L. 2015. Remote sensing image super-resolution method using sparse representation and classified texture patches. Geomatics and Information Science of Wuhan University, 40 (5): 578–582. [DOI: 10.13203/j.whugis20130385] ( 刘帅, 朱亚杰, 薛磊. 2015. 一种结合稀疏表示和纹理分块的遥感影像超分辨率方法. 武汉大学学报(信息科学版), 40 (5): 578–582. [DOI: 10.13203/j.whugis20130385] )
  • Mairal J, Bach F, Ponce J and Sapiro G. 2010. Online learning for matrix factorization and sparse coding. The Journal of Machine Learning Research, 11 : 19–60.
  • Mallat S G and Zhang Z F. 1993. Matching pursuits with time-frequency dictionaries. IEEE Transactions on Signal Processing, 41 (12): 3397–3415. [DOI: 10.1109/78.258082]
  • Massey R, Sankey T T, Congalton R G, Yadav K, Thenkabail P S, Ozdogan M and Meador A J S. 2017. MODIS phenology-derived, multi-year distribution of conterminous U.S. crop types. Remote Sensing of Environment, 198 : 490–503. [DOI: 10.1016/j.rse.2017.06.033]
  • Mo X G, Liu S X, Lin Z H and Qiu J X. 2011. Patterns of evapotranspiration and GPP and their responses to climate variations over the North China Plain. Acta Geographica Sinica, 66 (5): 589–598. [DOI: 10.11821/xb201105002] ( 莫兴国, 刘苏峡, 林忠辉, 邱建秀. 2011. 华北平原蒸散和GPP格局及其对气候波动的响应. 地理学报, 66 (5): 589–598. [DOI: 10.11821/xb201105002] )
  • Pan Y Z, Li L, Zhang J S, Liang S L, Zhu X F and Sulla-Menashe D. 2012. Winter wheat area estimation from MODIS-EVI time series data using the Crop Proportion Phenology Index. Remote Sensing of Environment, 119 : 232–242. [DOI: 10.1016/j.rse.2011.10.011]
  • Savitzky A and Golay M J E. 1964. Smoothing and differentiation of data by simplified least squares procedures. Analytical Chemistry, 36 (8): 1627–1639. [DOI: 10.1021/ac60214a047]
  • Skakun S, Franch B, Vermote E, Roger J C, Becker-Reshef I, Justice C and Kussul N. 2017. Early season large-area winter crop mapping using MODIS NDVI data, growing degree days information and a Gaussian mixture model. Remote Sensing of Environment, 195 : 244–258. [DOI: 10.1016/j.rse.2017.04.026]
  • Song L, Cheng Y M and Zhao Y Q. 2012. Hyper-spectrum classification based on sparse representation model and auto-regressive model. Acta Optica Sinica, 32 (3): 0330003 [DOI: 10.3788/aos201232.0330003] ( 宋琳, 程咏梅, 赵永强. 2012. 基于稀疏表示模型和自回归模型的高光谱分类. 光学学报, 32 (3): 0330003 [DOI: 10.3788/aos201232.0330003] )
  • Song X F and Jiao L C. 2012. Classification of hyperspectral remote sensing image based on sparse representation and spectral information. Journal of Electronics and Information Technology, 34 (2): 268–272. [DOI: 10.3724/SP.J.1146.2011.00540] ( 宋相法, 焦李成. 2012. 基于稀疏表示及光谱信息的高光谱遥感图像分类. 电子与信息学报, 34 (2): 268–272. [DOI: 10.3724/SP.J.1146.2011.00540] )
  • Tang P Q, Wu W B, Yao Y M and Yang P. 2011. New method for extracting multiple cropping index of North China Plain based on wavelet transform. Transactions of the CSAE, 27 (7): 220–225. [DOI: 10.3969/j.issn.1002-6819.2011.07.039] ( 唐鹏钦, 吴文斌, 姚艳敏, 杨鹏. 2011. 基于小波变换的华北平原耕地复种指数提取. 农业工程学报, 27 (7): 220–225. [DOI: 10.3969/j.issn.1002-6819.2011.07.039] )
  • Tang X Q, Liu Y Z and Chen J L. 2016. Improvement of remote sensing image classification method based on sparse representation. Computer Engineering, 42 (3): 254–258, 265. [DOI: 10.3969/j.issn.1000-3428.2016.03.046] ( 唐晓晴, 刘亚洲, 陈骏龙. 2016. 基于稀疏表示的遥感图像分类方法改进. 计算机工程, 42 (3): 254–258, 265. [DOI: 10.3969/j.issn.1000-3428.2016.03.046] )
  • Tropp J A and Gilbert A C. 2007. Signal recovery from random measurements via orthogonal matching pursuit. IEEE Transactions on Information Theory, 53 (12): 4655–4666. [DOI: 10.1109/tit.2007.909108]
  • Wang P X, Xun L, Li L, Xie Y and Wang L. 2017. Extraction of planting areas of main crops based on Fourier transformed characteristics of time series leaf area index products. Transactions of the Chinese Society of Agricultural Engineering, 33 (21): 207–215. [DOI: 10.11975/j.issn.1002-6819.2017.21.025] ( 王鹏新, 荀兰, 李俐, 解毅, 王蕾. 2017. 基于时间序列叶面积指数傅里叶变换的作物种植区域提取. 农业工程学报, 33 (21): 207–215. [DOI: 10.11975/j.issn.1002-6819.2017.21.025] )
  • Wang X, Li X B, Tan M H and Xin L J. 2015. Remote sensing monitoring of changes in winter wheat area in North China Plain from 2001 to 2011. Transactions of the Chinese Society of Agricultural Engineering, 31 (8): 190–199. [DOI: 10.3969/j.issn.1002-6819.2015.08.028] ( 王学, 李秀彬, 谈明洪, 辛良杰. 2015. 华北平原2001-2011年冬小麦播种面积变化遥感监测. 农业工程学报, 31 (8): 190–199. [DOI: 10.3969/j.issn.1002-6819.2015.08.028] )
  • Wardlow B D and Egbert S L. 2008. Large-area crop mapping using time-series MODIS 250 m NDVI data: an assessment for the U. S. Central Great Plains. Remote Sensing of Environment, 112 (3): 1096–1116. [DOI: 10.1016/j.rse.2007.07.019]
  • Whelen T and Siqueira P. 2017. Use of time-series L-band UAVSAR data for the classification of agricultural fields in the San Joaquin Valley. Remote Sensing of Environment, 193 : 216–224. [DOI: 10.1016/j.rse.2017.03.014]
  • Wright J, Yang A Y, Ganesh A, Sastry S S and Ma Y. 2009. Robust face recognition via sparse representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 31 (2): 210–227. [DOI: 10.1109/TPAMI.2008.79]
  • Wu W B, Yang P, Tang H J, Zhou Q B, Ryosuke S and Zhang L. 2009. Comparison of two fitting methods of NDVI time series datasets. Transactions of the CSAE, 25 (11): 183–188. [DOI: 10.3969/j.issn.1002-6819.2009.11.033] ( 吴文斌, 杨鹏, 唐华俊, 周清波, Ryosuke S, 张莉. 2009. 两种NDVI时间序列数据拟合方法比较. 农业工程学报, 25 (11): 183–188. [DOI: 10.3969/j.issn.1002-6819.2009.11.033] )
  • Xia Q, Xing S, Ma D Y, Mo D L, Li P C and Ge Z X. 2016. An improved K-SVD-based denoising method for remote sensing satellite images. Journal of Remote Sensing, 20 (3): 441–449. [DOI: 10.11834/jrs.20165149] ( 夏琴, 邢帅, 马东洋, 莫德林, 李鹏程, 葛忠孝. 2016. 遥感卫星影像K-SVD稀疏表示去噪. 遥感学报, 20 (3): 441–449. [DOI: 10.11834/jrs.20165149] )
  • Yin W, Li Y X, Zhou Z M and Liu S Q. 2013. Remote sensing image fusion based on sparse representation. Acta Optica Sinica, 33 (4): 0428003 [DOI: 10.3788/AOS201333.0428003] ( 尹雯, 李元祥, 周则明, 刘世前. 2013. 基于稀疏表示的遥感图像融合方法. 光学学报, 33 (4): 0428003 [DOI: 10.3788/AOS201333.0428003] )
数据专区
  • 目录 目录
  • 文章信息 文章信息
  • 引用格式 引用格式
  • 地图 地图
  • 评论 评论
  • 图 图
  • 表 表
  • 二维码 二维码
  • 参考文献 参考文献
  • 多媒体 多媒体