肖京格，喬彥友，王成波，夏昊，付東

(1.中國科學院遙感與數(shù)字地球研究所，北京 100094；2. 中國科學院大學，北京100049)

0 引言

土地覆蓋分類為全球碳預算、森林資源管理、生態(tài)環(huán)境調(diào)查和可持續(xù)發(fā)展研究等提供了重要基礎[1-3]，對研究人類與自然相互作用具有重要意義。衛(wèi)星遙感數(shù)據(jù)因其空間覆蓋范圍廣、持續(xù)觀測時間長和成本低等優(yōu)勢成為土地覆蓋分類的重要數(shù)據(jù)源。Landsat數(shù)據(jù)自2008年免費開放以來[4]，因其易獲取性、顯著表達地物特征的能力和逐步完善的數(shù)據(jù)生產(chǎn)體系[5]，被越來越多地應用到大面積土地覆蓋分類研究中。然而，常用的Landsat多光譜數(shù)據(jù)只有30 m空間分辨率，且可用光譜波段不多，即使是最新投入使用的陸地成像儀(operational land imager，OLI)，也只能獲取到9個不低于30 m空間分辨率的波段[6]，導致其生產(chǎn)的土地覆蓋分類產(chǎn)品精度有限。通常遙感數(shù)據(jù)的評價與分析多從空間、光譜和時間三個維度展開[7]。由于同一地點一年內(nèi)可以獲得數(shù)十景Landsat衛(wèi)星影像，具備時序特點，為從時間維度提高分類精度提供了可能。

基于多時相影像數(shù)據(jù)的類別特征構造多以差值法和均值法為主[8]，前者通過對多時相數(shù)據(jù)的特定波段或光譜指數(shù)做差，以突出目標地物；后者通過對數(shù)據(jù)堆疊并計算像元均值來達到集成多時相信息的目的。這兩種方法的共性問題是對遙感數(shù)據(jù)時相的選擇依賴性太強，時空移植性差；同時由于粗放地對數(shù)據(jù)進行差值和均值計算，并沒有充分挖掘和利用多時相數(shù)據(jù)中所蘊含的信息[9]，后來有學者嘗試從時間序列數(shù)據(jù)中挖掘信息以構造分類特征。如Julien等通過將時序LST和NDVI與一種土地覆蓋動力學方法結合，來對不同農(nóng)作物覆蓋類型進行分類[10]；Joseph等首先將時序數(shù)據(jù)構造成季節(jié)性三元組，然后用纓帽(tasseled-cap)亮度、濕度和歸一化植被指數(shù)(normalized difference vegetation index，NDVI)執(zhí)行數(shù)據(jù)分類[11]；Li M等則將地形特征如海拔、坡度和坡向等信息與時序光譜信息結合開展土地覆蓋分類[12]。然而，Landsat時間序列數(shù)據(jù)所固有的復雜性也對相關研究結果造成了不利影響。Landsat重訪周期較長(16 d)，此外，增強型專題繪圖儀(enhanced thematic mapper plus，ETM+)的掃描線校正器(scan line corrector，SLC)在2003年失效，再加上云覆蓋的影響，使大面積土地覆蓋監(jiān)測中，常難以獲得足夠理想的時序影像。如何在現(xiàn)有數(shù)據(jù)質(zhì)量的前提下構造更有效的特征，提高分類精度，儼然成為一個難題。

對于分類器的選擇，傳統(tǒng)遙感影像分類中有最大似然、決策樹、支持向量機、隨機森林等[13]；2006年后，隨著深度學習的興起，深度置信網(wǎng)絡(deep belief network，DBN)[14]、卷積神經(jīng)網(wǎng)絡(convolutional neural network，CNN)[15]、自編碼器(autoencoder)[16]等被嘗試應用到遙感影像分類中。這些方法大多源自于傳統(tǒng)圖像分析與自然語言處理等領域。由于遙感數(shù)據(jù)的多光譜和地物的多空間尺度、異質(zhì)性，以及訓練樣本有限、物候變化等問題，導致部分方法并不能取得理想效果，有待探索更具實用性和準確性的分類算法。

為了解決應用Landsat數(shù)據(jù)進行大面積土地覆蓋遙感監(jiān)測中特征構造與分類器選擇方面所存在的問題，本文提出一種基于年度時間序列數(shù)據(jù)的分類方法(LandUTime)。該方法更多關注于不同土地覆蓋類型物理性質(zhì)的年度變化特征，從不同的變化模式中提取信息，實現(xiàn)類別劃分，克服數(shù)據(jù)質(zhì)量、復雜地表、季候變化等造成的干擾。根據(jù)所生成的時序特征，設計了一種更具適用性的集成分類算法(UniBagging)，以獲得更高精度的土地覆蓋分類結果。

1 研究理論

該方法通過特征構造和分類2個階段實現(xiàn)(圖1)。特征構造階段，以對時間序列數(shù)據(jù)進行回歸分析為主，首先預處理時序數(shù)據(jù)，然后將獲得的原始特征層(如各個波段地表反射率值、NDVI等光譜指數(shù))在時間維上進行數(shù)據(jù)疊加，在像元級別上構建時序軌跡，利用回歸分析，對軌跡點進行擬合，對不同土地覆蓋類型的物理性質(zhì)變化進行建模，提取不同像元位置的模型參數(shù)，構造特征集；分類階段中，分類器集成與訓練是關鍵步驟，首先結合新構造的時序特征生成訓練樣本集和測試樣本集，之后根據(jù)原始特征層的選擇，將新生成的特征空間劃分為不同子空間，按子空間集成分類器，接著在各個子空間內(nèi)進行分類器的訓練，最后加權整合不同子空間中基分類器的預測值，輸出最終分類結果。

圖1 LandUTime原理流程

1.1 時間序列模型

不同地物具有不同的光譜值和光譜變化模式。LandUTime旨在基于Landsat時間序列數(shù)據(jù)對地表變化進行建模，并從構建的模型中提取特征，根據(jù)不同的變化模式區(qū)分土地覆蓋類型。本文中時間序列模型的意義則在于從復雜、無組織的數(shù)據(jù)中提取有用信息，用于模式分類。以NIR波段為例，選擇典型土地覆蓋區(qū)域的單個像元，取其一年中的時序像元值(云遮蓋部分被掩膜剔除)繪制散點圖(圖2(a))，若僅根據(jù)數(shù)據(jù)點的當前狀態(tài)來區(qū)分不同土地覆蓋類型顯得非常困難?；貧w分析作為一種幫助人們了解自變量變化時因變量如何變化的方法，在本文中被用于表征表面反射率和時間之間的關系，獲得回歸曲線(圖2(b))，這時根據(jù)不同曲線來區(qū)分地物類型則變得相對容易。

圖2 時序數(shù)據(jù)信息提取示意圖

對于像元位置 (m，n)，構造2個向量(1)和(2)，其中，xm，n是位置(m，n)所有可用時序像元的獲取時間所構成的向量，為方便計算，本文采用儒略日表示；ym，n代表對應的地表反射率值。通常情況下，一年內(nèi)土地覆被類型的變化可以被分解為兩部分[9]：季候因素引起的周期性變化，植被生長/退化、氣候或人為因素引起的漸緩式變化。因此，構建一種同時包含兩種變化類型的時間序列模型(3)。利用最小二乘法，基于剔除云遮蓋部分后剩余的時序清晰像元，進行回歸分析。

xm，n=(xm，n，1，xm，n，2，xm，n，3，…xm，n，U)

(1)

ym，n=(ym，n，1，ym，n，2，ym，n，3，…ym，n，U)

(2)

(3)

為了對本文所提出時間序列模型的有效性進行比較評價，特引入常用的多項式模型作為參照。出于保證多項式模型精度和防止過擬合的考慮，本文使用四階多項式模型(4)作為參照。

(4)

1.2 基分類器集成算法(UniBagging)

為了充分利用時間序列模型中蘊含的特征信息，本文基于集成學習的理念設計了一種專用于時序特征的分類方法，如圖3(a)所示。對于一個像元位置，該方法首先提取所有時間序列模型的系數(shù)(b0，b1，b2，b3)和擬合優(yōu)度評價指數(shù)(R2)作為初始時序特征。前文模型的回歸分析在特定原始特征層完成，此處特征塊需根據(jù)不同模型進行構造。特征塊將初始時序特征空間劃分為不同子空間，將分類器組織和集成到各個子空間中。在圖3(b)中，橢圓表示基分類器，矩形表示不同的特征子空間。像元位置(m，n)的原始時序特征經(jīng)分塊重組后用λ表示，v是特征塊，i 代表不同原始特征層的索引，則得到：

λ=(v1，v2，v3，…，vi，…)

(5)

(6)

之后，根據(jù)劃分后的子空間，將樣本映射為不同子集，輸入分類器進行訓練和驗證。具體算法實現(xiàn)流程見表1。

圖3 UniBagging的特征劃分與分類器子空間集成原理

表1 UniBagging算法實現(xiàn)流程

2 實驗與分析

2.1 實驗設計

為進行驗證與分析，本文共設置5組實驗，見圖4。實驗1采用Landsat 6個常用波段(BLUE，GREEN，RED，NIR，SWIR 1，SWIR 2)的地表反射率值作為特征輸入，實驗2采用由模型(4)提取的多項式特征，實驗3、實驗4、實驗5則采用本文所提出的時間序列特征，但使用了不同的分類算法。通過比較實驗1、實驗2、實驗3來評價時序特征的有效性，通過比較實驗3、實驗4、實驗5來評價分類方法UniBagging的有效性。為了方便比較分析，減少不確定因素的干擾，本文采用KNN為實驗5的基分類器，采用與實驗1相同的6個常用波段作為LandUTime的原始輸入特征層。結合模型(3)可知，初始時序特征數(shù)為5×6=30個，每個像元位置獲得的特征塊為6個。

圖4 方法比較的實驗設置

2.2 實驗區(qū)與數(shù)據(jù)

研究區(qū)位于中國安徽省西南部，所屬WRS行列號為039/121，占地面積4 107.87 km2，覆蓋Landsat圖像2 401×1 901個像元。該研究區(qū)處在亞熱帶與北溫帶過渡帶，四季分明，氣候溫和，雨量充沛，光照充足，無霜期長，適宜于暖溫帶和亞熱帶各種林木生長。有各種典型的植被(包括比例相當?shù)某＞G林與落葉林)和土地覆蓋類型。森林覆蓋率高，且人為因素影響下的土地覆蓋(建設用地、耕地等)區(qū)域也比較明顯，具有比較強的區(qū)域代表性。

以2014年該區(qū)域的土地覆蓋情況為研究對象。實驗1中的影像獲取于2014年5月1日，這是該年可以獲取到的距離生長季最近且目標區(qū)域無云覆蓋的影像。實驗2～實驗5使用了相同的原始數(shù)據(jù)，下載了該區(qū)域云覆蓋率小于80%的所有Landsat ETM+/OLI影像，共24 景，其中Landsat-7 ETM+影像7 景，Landsat-8 OLI影像17 景。本研究中使用的所有影像數(shù)據(jù)均為“ESPA On Demand Product”產(chǎn)品(https：∥espa.cr.usgs.gov/)。影像數(shù)據(jù)經(jīng)過了輻射定標、大氣校正、幾何校正和正射校正等處理。其中，Landsat-7 ETM +影像使用LEDAPS系統(tǒng)進行預處理，Landsat-8 OLI數(shù)據(jù)用LaSRC系統(tǒng)進行預處理。為了消除云覆蓋對分類結果的干擾，采用Fmask算法逐一去除每期影像上的云和云陰影像元。統(tǒng)計分析顯示，大部分剩余像元均較清晰，幾乎不存在云覆蓋與云陰影像元，表明去除效果較徹底。時序數(shù)據(jù)集中像元總數(shù)為109 543 224個(2 401×1 901×24)，剔除像元47 754 867個，剩余數(shù)據(jù)中各空間位置的平均時序像元數(shù)為10，足以保證時序建模。

根據(jù)實驗區(qū)地面實際情況，參照美國國家土地覆蓋數(shù)據(jù)庫(national land cover database，NLCD)2011版的標準[18]，確定了7種土地覆蓋類型：裸地、耕地、落葉林地、常綠林地、混交林地、建設用地和水體，每種類型的訓練與驗證樣本數(shù)見表2。在驗證階段，通過分層抽樣，每個類別選取700個像元構成驗證樣本集。樣本是基于實地調(diào)查和高分2號融合影像(分辨率0.8 m)人工解譯獲得。

表2 訓練樣本與驗證樣本數(shù)量表

2.3 實驗結果與分析

為了驗證和評價LandUTime土地覆蓋分類方法的有效性，本文進行了5組實驗。圖5顯示了在目標區(qū)內(nèi)的分類結果。表3顯示了5組實驗得到的分類精度。

1)時序特征有效性評價。從表3可以看出，在實驗1、實驗2、實驗3中，實驗3獲得的總體精度(OA)和Kappa系數(shù)(Kappa)最高，表明了本文所提出的時序特征對提高分類精度的有效性，而實驗1中6個地表反射率波段在區(qū)分不同土地覆蓋類型方面則顯得不夠理想。這主要由地物光譜復雜性所致，同物異譜和同譜異物現(xiàn)象普遍存在。如在生長季，常綠林地、落葉林地和耕地則很難精確區(qū)分。

圖5 土地覆蓋分類結果

表3 實驗精度與算法運行時間

注：OA為總體精度；Kappa為Kappa系數(shù)；Tt為分類器訓練所需時間，單位為s；Tp為分類器完成實驗區(qū)土地覆蓋分類所需時間，單位為s。

然而，可以發(fā)現(xiàn)實驗2的結果也比較差，為進一步分析原因，將多項式模型提取到的特征與本文所提出模型進行比較，特將二者進行回歸分析時獲得的R2(用于評價擬合優(yōu)度)進行比較，對研究區(qū)內(nèi)所有像元位置的R2值做直方圖統(tǒng)計，繪制直方圖變化曲線如圖6。其中L1是由本文所提出模型回歸的R2進行統(tǒng)計而繪制，L2是根據(jù)多項式模型回歸的結果繪制。其中L1所代表數(shù)據(jù)集的均值為0.56，方差為0.19；L2所代表數(shù)據(jù)集的均值為0.53，方差為0.18。很明顯，4階多項式對時序數(shù)據(jù)的回歸擬合精度要比本文所提出的模型高，但其分類精度卻較低。這主要是因Landsat數(shù)據(jù)問題與多項式模型自身局限性所致。由于云覆蓋和Landsat-7數(shù)據(jù)條帶的影響，大部分時序軌跡都存在數(shù)據(jù)點缺失的問題。對于這種情況，多項式模型的回歸，只能做到對現(xiàn)有數(shù)據(jù)擬合誤差的最小化，而忽略了時序軌跡所代表的地物物理性質(zhì)變化，導致所獲取的特征不能對實際地物變化進行較為準確的建模，從而加重誤分類。

圖6 不同模型R2統(tǒng)計直方圖變化曲線

本文所提出的模型具有周期性，符合土地覆蓋類型年度物候變化特性。通過對時序數(shù)據(jù)進行回歸分析，可以對不同土地覆蓋類型進行科學建模。而且該模型主要由正弦、余弦和一元線性成分組成，數(shù)據(jù)依賴性低，較少的幾個關鍵數(shù)據(jù)點即可擬合出可靠的模型參數(shù)。這些因素極大降低了像元值缺失問題對特征提取環(huán)節(jié)的干擾，使得該方法對云覆蓋和ETM+數(shù)據(jù)條帶問題的影響表現(xiàn)得更魯棒。

2)分類算法差異性分析。在使用本文所提出的時序特征的3組實驗(實驗3、實驗4、實驗5)中，隨機森林得到了最低的總體精度和Kappa值。雖然實驗3中KNN獲得了比較好的分類結果，但在處理高維數(shù)據(jù)時，KNN的性能受到了嚴重影響，實驗2和實驗3中，KNN完成目標區(qū)域土地覆蓋分類分別用了5 226.73和2 426.67 s，比其他組實驗耗時都要長。實驗5中，UniBagging得到了最高的準確性，同時也能很好地處理高維數(shù)據(jù)。因為采用KNN作為基本分類器時，每個KNN只接受一個子空間中的變量作為輸入，這使得當特征維度升高時，UniBagging的性能依然不受影響。

隨機森林之所以會得出比較差的分類結果，主要是由于其采用了隨機子空間的樣本采樣方法，而這種方法并不適用于本文中的時序特征，因為來自不同模型的系數(shù)和R2不能被混合。如果將來自藍波段的b0、b3與來自紅波段的b1、b2和R2被隨機采樣出來，它們并不能構成任何有意義的曲線，而隨機子空間恰恰會出現(xiàn)這一狀況。只有相應的模型系數(shù)才可以正確重建表征地物光譜年度變化的曲線。

3)不同覆蓋類別的精度差異。為了分析LandUTime對土地覆蓋類型的分類精度，計算了實驗1和實驗5的混淆矩陣和各類別用戶精度與生產(chǎn)精度。如表4所示，LandUTime在區(qū)分不同覆蓋類型方面的優(yōu)勢非常明顯，尤其對于耕地、常綠林地、落葉林地等，取得了非常好的效果。若使用單一時相影像，無法達到這樣的精度。同時可以看出，實驗1和實驗5對混交林地的分類準確度相對較低，主要是因為該種地類異質(zhì)性較強，實際自然界中混交林、落葉林和常綠林三者的邊界本來就很模糊，與空間尺度極其相關。還有是裸地、建設用地和水體這3種類型間存在少量錯分，主要是由于LandUTim方法依據(jù)地物年度變化模式的不同來區(qū)分地類，故具有顯著年度季節(jié)變化的對象(如耕地，落葉林地)與年度變化較小的地物(如建設用地、水體等)很容易區(qū)分，但變化模式相近的地物類型則可能出現(xiàn)些許混淆。

表4 實驗1和實驗5的混淆矩陣與類別精度

注：BL為裸地；CC為耕地；DF為落葉林地；EF為常綠林地；MF為混交林地；DL為建設用地；OW為水體；PA為生產(chǎn)者精度；UA為用戶精度。

4 結束語

本文面向大面積土地覆蓋分類問題，提出了一種基于Landsat年度時間序列數(shù)據(jù)的方法，該方法主要包括兩部分，首先是利用回歸分析的時序特征構造策略，然后是基于時序特征子空間的集成分類算法。實驗結果證明，該方法取得了非常好的分類效果，它既充分利用了Landsat數(shù)據(jù)在時間維度中所蘊含的信息，又克服了云覆蓋以及Landsat-7中掃描線校正器(SLC)故障等造成的時序數(shù)據(jù)缺失的問題；同樣基于集成學習的理念，該方法取得了比隨機森林更高的精度，克服了隨機子空間法所造成的特征結構被破壞的缺陷；隨著特征維度的提升，該方法的性能依然穩(wěn)健。

本文研究表明從遙感數(shù)據(jù)的時間維挖掘特征信息對于提高土地覆蓋分類精度具有重大潛力?！巴V異物，同物異譜”的問題導致常規(guī)從空間維和光譜維提取到的特征具有一定局限性，從時間維著手可以克服現(xiàn)有問題，為實現(xiàn)分類精度提升提供了的新可能。

LandUTime具有很強的可擴展性。為了驗證該方法的有效性，文中只采用了6個常用Landsat波段，還可把光譜指數(shù)(如NDVI、NDMI、EVI等)或光譜指數(shù)與Landsat波段組合作為原始特征層進行時間序列回歸分析。同時，該方法在其他區(qū)域環(huán)境的普適性探索也將成為未來研究的重點。