基于諧波分析和線性混合模型的河北平原區土地覆被分類研究

張 連 翀，王 衛

(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