基于谐波分析和线性混合模型的河北平原区土地覆被分类研究

张 连 翀，王 卫

(1．中国科学院遥感与数字地球研究所，北京 100094；2．中国科学院大学，北京 100049；3．河北师范大学资源与环境科学学院河北省环境演变与生态建设重点实验室，河北 石家庄 050024)



基于谐波分析和线性混合模型的河北平原区土地覆被分类研究

张 连 翀1,2,3，王 卫3*

(1．中国科学院遥感与数字地球研究所，北京 100094；2．中国科学院大学，北京 100049；3．河北师范大学资源与环境科学学院河北省环境演变与生态建设重点实验室，河北 石家庄 050024)

土地覆被分类是研究土地利用/覆被变化的基础数据和关键环节。以16 d合成MODIS-EVI时间序列数据为主要数据源，采用谐波分析方法分析不同土地覆被类型的季节性变化规律和物候特征差异，引入谐波特征值构建线性混合模型，提取不同端元的丰度值。从土地覆被类型较齐全、谐波特征具有代表性的石家庄地区高空间分辨率影像上选择训练样本，确定MODIS纯净像元和混合像元的划分阈值，对河北平原区进行土地覆被分类制图。结果表明，与河北省县级土地调查统计数据对比，一年两熟耕地、一年一熟耕地、园地及有林地、自然陆地表面的总量精度分别为90.19%、86.17%、85.96%和77.82%，平均总量精度为85.03%；与石家庄地区9个县(市)一年两熟耕地和一年一熟耕地基于TM的分类结果对比，平均面积相对误差分别为10.25%、13.98%。受粗空间分辨率和合成周期、水热条件以及种植模式破碎化限制，混合像元主要集中在河北平原中东部地区，一年两熟耕地、一年一熟耕地、园地及有林地混合面积比重较大。

谐波分析；线性混合模型；MODIS时间序列；混合像元；土地覆被分类

0 引言

土地覆被分类是区域土地利用/覆被变化驱动力分析、净初级生产力、土壤呼吸、蒸散量估算的基础数据和关键环节。植被指数(Vegetation Index，VI)时间序列所具有的季节规律是植被季相变化特征的综合反映，因此植被指数时间序列分析在研究区域乃至全球范围的物候现象、改进土地覆盖分类效果、研究气候变化效应等方面发挥着重要作用[1]。

利用长时段卫星数据获取地表植被物候信息，具有周期短、覆盖范围广、不受地理条件限制、成本低和信息量大等特点，尤其是数据序列的时空一致性好[2]。针对不同植被的物候期差异，确定其土地覆被类型，国内外学者已经取得了大量研究成果。Jakubauskas等[3]基于NOAA/AVHRR-NDVI数据对欧洲作物进行分类；候光雷等[4]基于SPOT/VGT-NDVI数据通过神经网络分类法提取东北地区耕地资源；张霞等[5]基于MODIS-EVI数据通过决策树法实现华北地区土地覆盖分类；陆广勇等[6]基于MODIS-NDVI数据通过决策树法获取黄土丘陵区土地覆被分类结果。研究表明，同一区域内不同植被类型在物候信息上表现出明显的生长差异特征，且与其空间分布格局存在相关性。因此，利用植被指数时间谱蕴涵的物候周期特征对植被覆盖分布具有较为显著差异性的地区进行分类研究，具有较高的可行性。

本文基于MODIS-EVI时间序列数据对河北平原区进行土地覆被分类研究：利用谐波分析法对原始数据进行降噪重建，在降低云覆盖、大气环境变化和双向反射特性等噪声污染的同时，又能突出不同植被的物候信息差异；引入谐波特征值构建线性混合模型对河北平原区典型土地覆被类型进行丰度解混，使不同端元对像元信息的贡献度得到有效表达；结合石家庄地区TM数据选择训练样本并确定划分阈值，得到亚像元级的土地覆被分类结果。

1 数据源和数据处理

1.1 研究区域

河北平原区位于河北省域东南部，总面积约8.16万km2，占全省面积的43.4%。自太行山山麓至滨海地带，地势低平，是我国华北平原的主要组成部分。一年两熟耕地土地熟化度高，田块平整且规模连片种植，以冬小麦、夏玉米为主，粮食高产县多集中于京广铁路沿线的山麓平原地带；一年一熟耕地主要连片分布在黑龙港低平原等商品棉生产基地，其次零散分布在冀东平原等地区。裸地、休闲耕地主要位于滨海地区，土壤盐碱化严重。果园及落叶阔叶林是该区最主要的木本植被类型。

1.2 研究数据

MODIS数据由NASA提供(http://reverb.echo.nasa.gov)，包括16 d合成增强型植被指数(EVI)数据(MOD13Q1，250 m)和年合成地表覆盖类型数据(MOD12Q1，500 m)，其中MOD13Q1数据选择时段为2011年1月1日至2011年12月31日，共计23个时相的观测数据。站点数据选自中国农作物产量资料旬值数据集(http://cdc.nmic.cn)，包括研究区内38个农业气象观测站的冬小麦、夏玉米和37个农业气象观测站的棉花、春玉米等作物名称、发育期名称、发育期日期、发育程度、发育期距平等数据。这些资料均通过调查大田作物获得，因此对大面积的作物物候具有很好的面上代表性。本文主要利用站点地理位置信息，为选择不同熟制土地覆被类型的训练样本提供依据。辅助数据中的遥感数据选自2011年不同季节云量较少的16景Landsat TM石家庄地区9个县(市)30 m空间分辨率影像(http://glovis.usgs.gov)，用于选择不同土地覆被类型的训练样本、确定纯净像元和混合像元的划分阈值和面积相对误差评价；统计数据选自2010年河北省土地调查统计数据，用于和分类结果进行精度比较。

1.3 数据预处理

首先，利用MODIS数据处理软件MRT(MODIS Reprojection Tool)对研究区同一时相MODIS13Q1的所有Tile进行EVI数据提取，并对其进行批量拼接、WG84坐标系和UTM投影转换等预处理操作；其次，将MOD12Q1-IGBP分类系统中不符合季节性变化特征的水体、城市和建成区作为掩膜对MODIS-EVI数据进行裁切，最终获得研究区时间序列数据。

2 分类体系和研究方法

2.1 分类体系

本研究在比较分析IGBP、UMD、CAS1990、GLC2000四种土地覆被分类体系的基础上，结合河北平原区自然、社会经济和耕地制度特征等实际情况，将研究区分为一年两熟耕地、一年一熟耕地、园地及有林地(果园、乔木、灌丛等自然林地)、自然陆地表面(裸地、休闲耕地、沙地)、人工陆地表面(居民点、交通、工矿用地)、水体(陆地水域、沿海滩涂(含盐田)、沟渠等用地)、一年一熟耕地/一年两熟耕地、一年一熟耕地/园地及有林地、一年两熟耕地/园地及有林地等9种土地覆被类型。

2.2 时间序列谐波分析

谐波分析法(Harmonic Analysis of Time Series，HANTS)由Jakubauskas提出[3]，即将时间序列数据按傅立叶级数展开为不同阶频率谐波叠加的形式。它将植被物候特征信息集中在低阶谐波，云覆盖、大气环境变化和双向反射特性等噪声信息则被分配到高阶谐波。将高阶谐波剔除后，选择若干能够反映地表植被动态的低阶谐波进行叠加拟合，从而实现原始时间序列数据的高质量重构。该方法揭示了像元时间特征维的周期性规律，定量化地表达植被动态变化信息，为植被空间分布格局的分析和提取提供了有效手段。

谐波分析法的核心是二维傅里叶变换和最小二乘法拟合，特征值包括谐波余项、谐波幅值、谐波相位等。其中，谐波余项表示数据的年内均值，反映了植被综合生产力状况；谐波幅值表示数据在不同周期模式内的波动范围，幅值越大说明植被对应周期变化的模式越明显；谐波相位表示数据峰值出现时间，反映了植被的最大生长期[7-9]。

设f(t)是周期为N的离散时间序列信号，即：

(1)

它的傅里叶级数展开式为：

(2)

(3)

(4)

(5)

(6)

(7)

(8)

式中：A0为谐波余项，k为谐波阶数，Ak、ωk、φk为第k阶谐波的幅值、频率、相位。ak、bk为第k阶谐波的傅里叶系数,采用最小二乘法拟合：

(FT×F)×C=FT×Y

(9)

式中：C为系数矩阵，F为傅立叶矩阵，FT为傅立叶矩阵的转置矩阵。

分类前，首先要验证不同地表类型谐波特征值存在显著差异。受空间分辨率和地物复杂多样性限制，研究区内像元混合严重，在图像中很难找到只包含一种土地覆被的纯净像元。因此，本文将具有相同或相似生长周期的土地覆被类型(如具有相同复种指数的棉花、春玉米)归为一类，结合研究区农业气象观测站点数据和MOD12Q1年合成地表覆盖类型数据，选取一年两熟耕地、一年一熟耕地、园地及有林地、自然陆地表面等进行谐波特征值分析。

2.2.1 一年两熟耕地 冬小麦从2月底开始返青，在5月初进入抽穗阶段，形成第1个峰值，6月中旬成熟收割后形成第1个波谷；夏玉米此时开始生长出苗，8月中下旬进入抽穗阶段形成第2个峰值，至10月初成熟收割后形成第2个波谷；10月中下旬开始播种冬小麦。此阶段出现峰值原因是研究区内光温条件能够保证冬小麦越冬前仍具有一定的植被覆盖度，但并不代表其1个生长季。此后冬小麦进入越冬期直至次年2月底，形成第3个波谷(图1a)。在谐波分量特征值中(图1b)，1、3阶谐波幅值较高，2阶谐波幅值非常小。1阶谐波幅值反映了作物的全年生长整体情况，3阶谐波幅值反映了作物的全年生长过程。3阶谐波相位与时间序列信号峰值时间一致，3个波谷分别对应冬小麦的初始生长期、冬小麦的收割期、夏玉米的收割期，3个波峰分别对应冬小麦的生长高峰、夏玉米的生长高峰、冬小麦的出苗高峰。1、3阶谐波特征值包含了该类型绝大部分植被信息，拟合曲线能够代表其季变规律和物候特征。

图1 一年两熟耕地时间序列谐波分析

Fig.1 Harmonic analysis of time series of two crops a year

2.2.2 一年一熟耕地 一年一熟耕地作物类型多样，生育关键期迭和，16 d的合成周期反映了该熟制作物的综合物候信息。鉴于棉花在该区域内集中分布，选择其作为代表进行分析。棉花从4月中旬开始播种，7月底达到峰值，至10月中下旬成熟收割，冬季达到最低值(图2a)。在谐波分量特征值中(图2b)，1阶谐波幅值最大，说明该类型是以全年为生长周期的单峰变化模式，峰值在7月底至8月中旬附近；2阶谐波的较大幅值主要受春季杂草生长的影响，该分量峰值并不代表实际生长过程。1、2阶谐波特征值包含了一年一熟耕地的绝大部分植被信息，拟合曲线能够代表其季节性变化规律和物候特征。

图2 一年一熟耕地时间序列谐波分析

Fig.2 Harmonic analysis of time series of one crop a year

2.2.3 园地及有林地 园地及有林地植被生长过程基本一致。植被在3月中旬开始发芽，8月底达最大值，果实成熟后开始落叶，EVI值逐渐下降，冬季达最低值，生长周期长(图3a)。在谐波分量特征值中(图3b)，1阶谐波幅值远大于其他阶谐波幅值，说明该类型是以全年为生长周期的单峰变化模式。1阶谐波特征值包含园地及有林地的绝大部分植被信息，拟合曲线代表了其季节性变化规律和物候特征。

图3 园地及有林地时间序列谐波分析

Fig.3 Harmonic analysis of time series of garden woodland

2.2.4 自然陆地表面 该类型在时间序列数据上呈现出与植被生长相似的变化趋势。这是由于植被与非植被的混合像元普遍存在，MODIS-EVI值是两者在像元尺度上的线性加权和，随像元内植被观测值的变化而变化(图4a)。在谐波分量特征值中(图4b)，1阶谐波幅值远远大于其他阶谐波幅值，说明该类型是以全年为生长周期的单峰变化模式。1阶谐波特征值包含了自然陆地表面的绝大部分植被信息，拟合曲线能够代表像元中植被的季节性变化规律和物候特征。

图4 自然陆地表面时间序列谐波分析

Fig.4 Harmonic analysis of time series of natural land

2.3 线性混合模型

不同土地覆被类型时间序列谐波特征值差异明显，为构建线性混合模型提供了可分性依据。MODIS时间序列数据被认为是不同土地覆被观测特征和植被生长动态的综合表达。设r是影像中任意N维时间谱向量，M=[m1，m2，…，mp]是大小为N×P的端元光谱矩阵，s=[s1，s2，…，sp]T是P维系数矢量，其各分量元素对应于像元r中端元mi的组分丰度，n为N维噪声项，则线性混合像元模型表示为：

(10)

式中：组分向量s满足“和为1约束”，即s1+s2+…+sp=1，以及“非负约束”，即si≥0。考虑影像中所有像元，式(10)可以改写为如下矩阵形式：

R=MS+N

(11)

式中：组分矩阵S满足“和为1约束”和“非负约束”。

首先，将前3阶谐波幅值、相位做最小噪声分离变换，用前5个特征值计算纯净像素指数，设置阈值在100～200之间的像元为感兴趣区，获得具有最纯像元的多维散点图；其次，根据最纯净像元应在数据云团最外侧拐角，其内部点是这些最纯净像元线性组合的原理，进一步选取纯净端元；再次，将端元平均值作为该端元值，在丰度值总和为1的限制性条件下求解线性混合方程，得到初始丰度值(图5)；最后用均方根误差(RMSE)评价分解精度。

(12)

式中：p表示端元数目，εi是第i种端元的残余误差。

RMSE值越小说明分类精度越高。初始RMSE最大值为0.017855，最小值为0.000091，平均值为0.000637。重新提取误差较大区域的端元，反复运行线性混合模型，至RMSE中没有较高误差且所有丰度值非负、总和小于等于1时完成迭代。最终误差呈随机分布，95%以上的RMSE小于0.000178。

3 结果与分析

3.1 土地覆被分类结果

图5 线性混合模型分解结果

Fig.5 Results from linear mixture model

选择土地覆被类型较齐全、谐波特征具有代表性的石家庄地区作为训练样本。1)对该地区Landsat TM影像按分类体系进行最大似然度分类；2)将结果重采样到250 m分辨率后与MODIS影像配准；3)从分类结果上对各土地覆被类型随机选取100个纯净像元，统计对应MODIS像元土地覆被端元的丰度频率累积值，考虑到配准及试验误差，取频率为5%的丰度值作为区分MODIS纯净像元和混合像元的阈值(表1)，如果端元丰度值大于或等于阈值，为纯净像元，反之为混合像元；4)将混合像元内的各端元丰度值两两相加，将最大加和值对应的端元组合确定为该混合像元的土地覆被类型，最终得到研究区土地覆被分类结果(图6，见封3)。

表1 各端元纯净像元阈值

Table 1 Threshold of pure pixel of endmembers

一年两熟耕地一年一熟耕地园地及有林地自然陆地表面阈值0．88430．90950．89020．9240

分类结果表明，一年两熟耕地主要分布于冀中南地区的淀西清北平原、滹滏平原、滏西平原等地区；一年一熟耕地主要分布于黑龙港平原、运东平原、滦河及冀东平原等地区；园地及有林地集中分布于赵县、深州、沧县、昌黎等地区；人工陆地表面一部分位于人口数量多、经济水平发达的城市地区，另一部分位于铁路、公路沿线呈带状分布，平原区村庄众多也使得该类型呈点状零散分布；自然陆地表面主要分布于黄骅、海兴、盐山、孟村、文安、大城、任丘等市县所辖区域。混合像元主要集中在河北平原区中部地区，包括廊坊、保定、沧州、邢台、邯郸等部分地区，特别是渤海沿岸，该地区在省内地势最低，地下水质量差，既有旱、涝、碱和水资源缺乏等障碍因素，又遭受大风、海啸等自然灾害，导致种植业生产受到限制，休闲耕地或裸地较多，其中，一年两熟耕地与一年一熟耕地、一年两熟耕地与园地及有林地、一年一熟耕地与园地及有林地组合所占面积比重较大。

3.2 分类结果验证

研究区内MODIS像元端元面积等于该端元丰度值乘以像元面积，对所有像元的端元面积求和得到研究区端元的总面积。公式为：

(13)

式中：Si为端元组分的总面积，Ai为第i个像元中端元组分丰度值。

将2010年河北省县级土地调查统计数据中土地覆被类型按端元重分类后面积S0与研究区各端元面积Si对比，得研究区总量精度Kr(表2)。公式为：

(14)

表2 县级土地覆被分类面积与统计面积总量精度比较

Table 2 A county-level comparison of land cover area from MODIS data and statistical data

类别总量精度(%)平均精度(%)一年两熟耕地90．19一年一熟耕地86．17园地及有林地85．96自然及人工陆地表面77．8285．03

选择石家庄地区9个县(市)不同熟制耕地面积与基于TM的分类结果进行对比，一年两熟和一年一熟耕地的平均面积相对误差分别为10.25%、13.98%，平均相对误差为12.11%，与全省精度分析检验结果基本一致(表3)。

表3 石家庄市所辖9县(市)不同熟制耕地面积提取相对误差

Table 3 Relative error of different cropping system extraction of 9 counties in Shijiazhuang

耕地类型相对误差(%)藁城晋州辛集新乐深泽无极赵县高邑栾城平均误差(%)一年两熟9．397．174．6212．2112．5410．7513．5311．2910．7510．25一年一熟15．858．8722．5816．0714．8517．218．715．1916．5213．98

4 结论与讨论

本文利用16 d合成MODIS-EVI时间序列数据，对河北平原区不同土地覆被类型进行谐波分析，实现了植被季节性变化规律物候特征的参数化表达；将谐波特征值用于构建线性混合模型，实现了一年两熟耕地、一年一熟耕地、园地及有林地、自然陆地表面的丰度提取，获得了研究区土地覆被分类结果；与河北省县级土地调查统计数据对比，一年两熟耕地、一年一熟耕地、园地及有林地、自然陆地表面的总精度分别为90.19%、86.17%、85.96%、77.82%，平均总精度为85.03%；与石家庄地区9个县(市)一年两熟耕地和一年一熟耕地统计基于TM的分类结果对比，平均面积相对误差分别为10.25%、13.98%。

受粗空间分辨率和16 d合成周期的限制，谐波特征值表达的是地物综合物候信息，细节特征被忽略；一年一熟耕地作物类型多样，生育关键期迭和，导致“异物同谱”现象严重；以均方根误差、统计数据和同时期遥感数据对端元丰度和端元总量进行精度评价，缺乏其他多源数据的辅助参考。考虑到单一数据进行土地覆盖分类的局限性，今后考虑将不同高空间分辨率和高时间分辨率的遥感数据进行时空融合，并辅以其他遥感特征参量(如地表温度LST)的支持，以提高分类精度。

[1] 张霞，李儒，岳跃民，等．谐波改进的植被指数时间序列重建算法[J]．遥感学报，2010，14(3)：437-447．

[2] 林忠辉，莫兴国．NDVI时间序列谐波分析与地表物候信息获取[J]．农业工程学报，2006，22(12)：138-144．

[3] JAKUBAUSKAS M E，LEGATES D R，KASTENS J H．Crop identification using harmonic analysis of time-series AVHRR NDVI data[J]．Computers and Electronics in Agriculture，2002，37：127-139．

[4] 侯光雷，张洪岩，王野乔，等．基于时间序列谐波分析的东北地区耕地资源提取[J]．自然资源学报，2010，25(9)：1607-1617．

[5] 张霞，焦全军，张兵，等．利用MODIS_EVI图像时间序列提取作物种植模式初探[J]．农业工程学报，2008，24(5)：161-165．

[6] 陆广勇，杨勤科，王海江．基于MODIS-NDVI时序数据的黄土丘陵区土地覆被分类研究[J]．水土保持研究．2011，18(2)：112-120．

[7] ZHANG S W,LEI Y P,WANG L P,et al.Crop classification using MODIS NDVI data denoised：A case study in Hebei Plain，China by Wavelet[J].Chinese Geographical Science，2011，21(3)：322-333．

[8] ZHANG M W,ZHOU Q B,CHEN Z X,et al.Crop discrimination in Northern China with double cropping systems using Fourier analysis of time-series MODIS data[J].International Journal of Applied Earth Observation and Geoinformation，2008，10(4):476-485．

[9] 梁守真，邢前国，施平，等.山东省典型地表覆被NDVI时间序列谐波分析[J].生态学杂志，2011(1):59-65.

ChineseAcademyofSciences，Beijing100049；3．HebeiKeyLaboratoryofEnvironmentEvolutionandEcologyConstruction，SchoolofResourcesandEnvironmentalScience，HebeiNormalUniversity，Shijiazhuang050024，China)

Land Cover Classification of Hebei Plain Based on Harmonics Analysis and Linear Mixture Model

ZHANG Lian-chong1,2,3，WANG Wei3

(1．InstituteofRemoteSensingandDigitalEarth，ChineseAcademyofSciences，Beijing100094；2．Universityof

In this paper，the time series data of 16-days composite MODIS Enhanced Vegetation Index (EVI)were used to examine the land use/cover classification in Hebei Plain.The harmonic analysis of time series (HANTS)algorithm based on Fourier transformation was applied to predigest and compress the MODIS-EVI time series images,at the same time,the phenological features of two crops a year,one crop a year,garden woodland and natural land were extracted.The first to third eigenvalues were used to establish the linear mixture model and calculate the abundance of different endmembers.Shijiazhuang district,which had all types and characteristics of land cover and eigenvalue,was used to set the thresholds of the pure and mixed pixel based on the correlation between MODIS and Landsat images.The land cover mapping was completed by these thresholds．The result indicated that the overall accuracy of two crops a year,one crop a year,garden woodland and natural land was 90.19%，86.17%，85.96%,77.82% with the average was 85.03% when compared with official statistics at county-level.The average area error of two crops a year and one crop a year in Shijiazhuang were 10.25% and 13.98%.Mixed pixels are mainly concentrated in middle east of Hebei Plain and larger proportion is two crops a year，one crop a year，garden woodland.The error was caused mostly by the coarse resolution and composite period of MODIS，the water and heat conditions and the great spatial variability．

harmonic analysis；linear mixture model；MODIS time series；mixed pixel；land cover classification

2014-08-18；

2014-12-18

国家自然科学基金项目(41471091)；中国科学院“一三五”规划项目(Y3SG0300CX)；国家863计划项目(2013AA12A301)；河北省高校重点学科建设项目

张连翀(1985-)，男，博士研究生，主要从事高性能地学计算研究。*通讯作者E-mail：wangwei@mail.hebtu.edu.cn

10.3969/j.issn.1672-0504.2015.03.019

F301

A

1672-0504(2015)03-0098-06