作物空间分布信息的及时获取对农业生产管理、农业可持续发展及国家粮食安全具有重要意义，是区域作物产量估测、农作物种植结构调整和优化的主要依据(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傅里叶变换后的谐波特征作为不同作物识别的数据源，选取训练样本通过在线字典学习算法构建判别字典，根据字典计算待测样本的稀疏系数，通过最小重构误差确定作物种植类型，从而识别冬小麦、春玉米、夏玉米、棉花和果园种植区域，并对识别结果的位置精度进行验证。

华北平原位于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)。

通过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) Figure 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)。

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

式中， $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$ 。

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

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

假设有 $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)

式中， ${{{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}}$ 构成的集合中

式中， ${{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$ 。

稀疏表示要求字典具有过完备性，即字典中的原子个数要远大于样本的特征维度( $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$ 次迭代的稀疏系数

式中， $\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}$ 即为所求。

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

式中， ${\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}$ 。

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

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

利用上包络线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曲线 Figure 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级谐波曲线及谐波叠加曲线 Figure 3 Time series curves of S-G filtered LAI, harmonic terms 0—5 and their sum

表 1(Table 1)

表 1 研究区域主要作物的谐波参数 Table 1 Parameters of harmonic function for main crops in study area 种植模式 谐波 振幅 初始相位 种植模式 谐波 振幅 初始相位 冬小麦—夏玉米 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

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

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

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

每一类随机选取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)

表 2 混淆矩阵 Table 2 Confusion matrix 分类数据 验证数据 用户精度/% 冬小麦— 夏玉米 冬小麦—夏玉米& 其他地物 春玉米 早播夏玉米 棉花 果园 其他 总计 冬小麦—夏玉米 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

表 2 混淆矩阵 Table 2 Confusion matrix

本研究采用的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年的样本数据通过在线字典学习算法构建判别字典，根据判别字典批量化处理多年时间序列数据，进而提取多年作物种植区域，提高了数据处理效率。

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

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