陳 黔,李曉松,修曉敏,楊廣斌

1 貴州師范大學(xué) 地理與環(huán)境科學(xué)學(xué)院,貴陽 550001 2 中國科學(xué)院 遙感與數(shù)字地球研究所 數(shù)字地球重點實驗室,北京 100094

與喬木林相比,灌木林在整個生態(tài)系統(tǒng)中所占的比重較小,但由于其對水分、養(yǎng)分條件要求不高,一直以來在干旱、半干旱區(qū)有著大量的分布[1]。因此,灌木對于干旱生態(tài)系統(tǒng)的穩(wěn)定與平衡起著十分重要的作用[2]。準確掌握沙地灌木分布及動態(tài)具有重要的生態(tài)學(xué)意義。傳統(tǒng)植被分布調(diào)查主要以地面觀測為主,該方法具有一定的局限性,例如在人力、物力上的花費都比較大,樣地的外推也面臨著較大的不確定性[3]。隨著遙感對地觀測技術(shù)的不斷發(fā)展,使得大區(qū)域、大尺度的植被監(jiān)測成為了可能,而植被類型分類與覆蓋度估算一直是植被遙感監(jiān)測領(lǐng)域的熱點,已在全球不同尺度、不同地區(qū)得到了廣泛的應(yīng)用[4-8]。目前,絕大部分植被類型遙感分類以像元或若干像元的組合(對象)為基本單元,將其劃歸為某一指定植被類型[9-10],并未考慮像元中不同植被類型混合存在的情況,而這一問題對于中低分辨率遙感數(shù)據(jù)來說相當普遍,特別是在地表空間異質(zhì)性較強的區(qū)域[11-12]。意識到這一問題,估算每一像元特定植被類型覆蓋度的研究得到了開展,其中最為成熟、應(yīng)用最為廣泛的是全球MODIS FVC樹木覆蓋度產(chǎn)品[13],并且Sexton等人[14]已通過尺度轉(zhuǎn)換方法生產(chǎn)了全球30m空間分辨率樹木覆蓋度產(chǎn)品。值得指出的是,FVC產(chǎn)品中樹木覆蓋度主要反映的是高于5m的樹木的覆蓋度,并將灌木與草本等定義為其他植被覆蓋度,顯然無法滿足灌木信息提取的需求,另外其精度在以稀疏森林為主的干旱地區(qū)也存在較大疑問[15]。針對干旱地區(qū)的特殊情況,木本植被覆蓋度(包括喬木/灌木)的估算研究得到了開展[16],基于時間序列MODIS數(shù)據(jù)提取的季節(jié)參數(shù)估算了整個薩赫勒地區(qū)的木本覆蓋度,Higginbottom等[17]基于Landsat數(shù)據(jù)利用機器學(xué)習(xí)的方法估算了南非較小區(qū)域的木本覆蓋度。與FVC產(chǎn)品不同,這些研究將喬木與灌木信息混合在了一起,也無法滿足灌木覆蓋度提取的需求。因此,發(fā)展一個大尺度、中高分辨率、適用于干旱地區(qū)沙地、可與高大喬木與草本區(qū)分的灌木覆蓋度遙感估算方法具有重要的意義。

植被覆蓋度估算最為常用的就是植被指數(shù)法,一般而言植被指數(shù)與植被覆蓋度之間具有較高的相關(guān)性,但單純利用植被指數(shù)往往無法區(qū)分不同植被類型,并且受到荒漠地區(qū)植被稀疏及地表異質(zhì)性較強的影響,植被指數(shù)對于以根莖生長為主的灌木探測能力不足[11]?；旌舷裨纸夥ㄒ驳玫搅藦V泛應(yīng)用,該方法的優(yōu)勢是不需要采集大量地面樣本,然而合理的端元選取對估算精度影響較大,導(dǎo)致該方法的估算結(jié)果可能存在較大誤差[18-19]。與上述方法相比,機器學(xué)習(xí)方法運算速度更快、穩(wěn)定性更高且不容易產(chǎn)生過擬合,通過構(gòu)建樣本數(shù)據(jù)集,結(jié)合地形、氣候、植被指數(shù)、光譜信息等眾多變量,往往能構(gòu)建出更為穩(wěn)定的估算模型[14,16,20]?？紤]到大區(qū)域、多時相中高分辨率光學(xué)遙感數(shù)據(jù)的可用性及機器學(xué)習(xí)算法的計算需求,Google Earth Engine(GEE)遙感分析云平臺已成為開展相應(yīng)研究的最佳工具之一。GEE提供了校正后的全球Landsat系列數(shù)據(jù),并提供了機器學(xué)習(xí)算法接口[21]。研究人員已經(jīng)利用GEE開展了國家尺度植被動態(tài)監(jiān)測[22]、洲際尺度耕地制圖[23-24]、全球地表溫度估算[25]等等。然而,尚沒有針對沙地灌木覆蓋度估算的相關(guān)研究。

為此,本研究提出了充分利用GEE遙感云平臺,通過Collect Earth工具構(gòu)建地面灌木覆蓋度樣本庫,以Landsat-8影像及其他輔助數(shù)據(jù)為變量,通過機器學(xué)習(xí)方法開展大尺度時間序列灌木覆蓋度估算的技術(shù)體系,并以中國北方毛烏素沙地為示范區(qū)開展了應(yīng)用,并對其應(yīng)用效果進行了評價,以期為大尺度沙地灌木覆蓋度估算提供參考。

1 研究區(qū)概況

毛烏素沙地位于鄂爾多斯高原東南部和陜北黃土高原以北的洼地,覆蓋內(nèi)蒙古自治區(qū)的鄂爾多斯南部、陜西省榆林市的北部風沙區(qū)和寧夏回族自治區(qū)鹽池縣東北部小片區(qū)域[26],如(圖1)所示。地理位置為 37°27′30″—39°22′30″N,107°20′—111°30′E,地處中國東部季風的尾閭地區(qū),屬典型大陸性半干旱氣候,年降水量在250—440mm之間,面積約為4.22萬km2,在中國北方第二大沙地。毛烏素沙地流動沙地、半固定沙地分布較為廣泛,荒漠化進程迅速,因此,一直以來都是我國荒漠化治理的重點地區(qū)[27]。土壤在毛烏素沙地分布具有明顯的地帶性變化特征。在梁地、沙地地區(qū),土壤主要以栗鈣土和棕鈣土為主；在灘地中土壤主要呈草甸土—鹽化草甸土—鹽土過渡的形式；在河谷階地中,土壤由淡栗鈣土過渡為鹽化草甸土。整體來看,毛烏素沙地土壤偏沙性,其中栗鈣土、棕鈣土、灰鈣土、風沙土所占面積較大,養(yǎng)分含量低[28]。灌木植被是毛烏素沙地最具代表性,同時也是分布最為廣泛的植物群落。西鄂爾多斯荒漠地區(qū)大量分布超旱生灌木,如低矮灌木狹葉錦雞兒(Caraganastenophylla)、駝絨黎(Ceratoideslatens)等,中部及中東部干草原地區(qū)植被物種豐富多樣,并保有大量的原生灌木、半灌木植被,這其中就包括油蒿(Artemisiaordosica)、白沙蒿(Artemisiasphaerocephalakrasch)等蒿類植被以及楊柴(Hedysarummongolicum)、檸條(Caraganakorshinskii)等豆科植被[29]。隨著三北防護林、京津冀風沙源等一系列生態(tài)恢復(fù)工程的推進,據(jù)第五次全國荒漠化和沙化監(jiān)測最新數(shù)據(jù)顯示,毛烏素沙地作為重點治理區(qū)域,其植被覆蓋情況,尤其是作為優(yōu)勢植被的灌木覆蓋情況發(fā)生了顯著變化[30]。

2 研究方法

整體研究流程如(圖2)所示。首先,數(shù)據(jù)準備工作包括兩個部分,其一是對多源遙感數(shù)據(jù)進行選取并進行預(yù)處理,其二是通過兩種采樣法獲取樣本并構(gòu)建樣本數(shù)據(jù)集。其次,提取樣方對應(yīng)多源數(shù)據(jù),包括光譜參數(shù)、植被指數(shù)、物候因子以及地形因子,并對因子進行相關(guān)性分析。最后,構(gòu)建CART和SVM兩種估算模型并進行精度驗證,最終產(chǎn)出灌木覆蓋度產(chǎn)品。

圖2 沙地灌木覆蓋度估算技術(shù)流程圖Fig.2 Flowchart of methodology for shrub coverage estimation in sandy land

2.1 輸入因子的獲取與預(yù)處理

(1)衛(wèi)星遙感數(shù)據(jù)

本研究基于GEE平臺,以陸地數(shù)據(jù)Landsat-5 TM以及陸地數(shù)據(jù)Landsat-8 OLI地表反射率影像為數(shù)據(jù)源,利用GEE平臺強大的運算能力對影像進行質(zhì)量篩選、云掩膜及平均值計算,獲取毛烏素沙地2000—2017年(2012、2016除外)每年5—9月共80期地表反射率影像。在此基礎(chǔ)上,通過NDVI最大合成的方法生成了各年度生長季地表反射率影像。

除反射率外,本研究計算了歸一化植被指數(shù)(NDVI)、土壤調(diào)節(jié)植被指數(shù)(SAVI)、修正的土壤調(diào)節(jié)植被指數(shù)(MSAVI)、差值植被指數(shù)(DVI)[11,32]4個植被指數(shù)作為輸入?yún)?shù),以更好估算灌木覆蓋度。

(1)

(2)

(3)

DVI=NIR-R

(4)

式中,NIR為近紅外波段反射率,R為紅光波段反射率,L為降低土壤背景影響而設(shè)定的一個土壤調(diào)節(jié)參數(shù),該參數(shù)的值隨植被的密度變化而變化,Huete[33]認為在大面積稀疏植被覆蓋區(qū)域L取值為0.5時對消除土壤背景影響效果最佳。

(2)地形數(shù)據(jù)獲取與處理

地形數(shù)據(jù)調(diào)用GEE平臺上的SRTM數(shù)字高程模型,地形因子提取包括高程、坡度、坡向3個因子。

(3)氣候數(shù)據(jù)獲取與處理

氣候數(shù)據(jù)調(diào)用TerraClimate全球陸地表面月平均氣候數(shù)據(jù)集[31]。氣候因子提取包括降水、風速、帕爾默干旱指數(shù)、蒸散量4個因子作為輸入?yún)?shù)。

輸入變量匯總?cè)绫?所示。四類因子共包含17個變量。

2.2 樣本庫構(gòu)建

樣本數(shù)據(jù)集構(gòu)建是本研究的關(guān)鍵步驟。在此,我們運用了Collect Earth樣本收集器,這是一款由聯(lián)合國糧農(nóng)組織(FAO)開發(fā)的專業(yè)影像解譯軟件,主要用于陸地數(shù)據(jù)的收集、管理及分析。該軟件可以輕松地訪問多個免費的高時空分辨率衛(wèi)星影像存檔,利用這些數(shù)據(jù)的協(xié)同作用,實現(xiàn)了一種全新的土地監(jiān)測方法,稱為增強視覺解譯法[34]。該方法可提升空間采樣精度,特別是在地表異質(zhì)性較強的干旱、半干旱地區(qū),可解決植被稀疏所導(dǎo)致的采樣困難問題。本研究對毛烏素沙地進行樣本庫構(gòu)建,樣方以10km為間隔均勻分布在研究區(qū)內(nèi),共415個大小為100m×100m的樣方,每個樣方中等距分布有5×5,共25個采樣點,每個采樣點對應(yīng)4%的地塊,采用等距取樣法對每個采樣方分別進行覆蓋度提取,如(圖3)所示,通過這種方式實現(xiàn)覆蓋度的量化。根據(jù)實際情況對樣本進一步篩選,將影像時效性較差、影像被云遮蓋、影像分辨率較低的樣方進行去除,最終通過Collect Earth樣本收集器獲取376個有效樣方。

表1 輸入變量及特性統(tǒng)計表Table 1 Input factors and characteristics

2017年8月,課題組依據(jù)標準一致的野外采樣方案,對毛烏素沙地開展灌木覆蓋度調(diào)查,共采集大小為30m×30m的樣方24個。上述共400個樣方數(shù)據(jù)為毛烏素沙地灌木覆蓋度遙感估算打下了重要基礎(chǔ)。

圖3 基于Collect Earth樣本收集器樣點采集Fig.3 Illustration of sampling design and project properties

2.3 灌木覆蓋度估算方法

為提升數(shù)據(jù)的時間匹配度,本研究運用GEE平臺,分別對Landsat-8及Landsat-5影像進行時間序列選取、云量篩選、生長季NDVI最大值合成以及鑲嵌等處理,對2010至2015年5年間(2012年除外)獲取樣方對應(yīng)影像波段光譜值、植被指數(shù)、地形因子以及氣候因子進行提取,在此基礎(chǔ)上分別構(gòu)建分類回歸樹、支持向量機兩種估算模型,并對模型估算精度進行評價。

2.3.1CART(Classification And Regression Tree)模型

CART算法是Breiman[35]于1984年提出的一種基于決策樹構(gòu)建的算法,實質(zhì)上也是歸納學(xué)習(xí)算法中的一種。其主要原理是從已有的樣本庫中歸納出整體的某些特征規(guī)律,生成樹形分類結(jié)構(gòu)從而用于除樣本之外的事例分類?；驹硎峭ㄟ^訓(xùn)練數(shù)據(jù)集構(gòu)建預(yù)測模型,采用二分遞歸分割方法,使得生成的非葉子節(jié)點都有兩個分支。將多個分支預(yù)測模型進行擬合,通過預(yù)測連續(xù)變量獲得樹形結(jié)構(gòu)模型。該算法既可用于分類,也可用于連續(xù)變量的預(yù)測,可有效處理大區(qū)域、大樣本量和髙維度的數(shù)據(jù)。

CART算法在樹形結(jié)構(gòu)生長過程中,采用經(jīng)濟學(xué)中基尼(Gini)系數(shù)作為選擇最優(yōu)測試變量和分割閾值的準則,其數(shù)學(xué)定義如下所示：

(5)

(6)

式中,P(j|h)是從訓(xùn)練樣本集中隨機抽取的一個樣本,當某一測試變量值為h時屬于第j類的概率,njh為訓(xùn)練樣本中該測試變量值為h時屬于第j類的樣本個數(shù),n(h)為訓(xùn)練樣本中該測試變量值為h的樣本個數(shù),j為類別個數(shù)。

在本次研究中,目標變量是灌木覆蓋度,測試變量是用于構(gòu)建模型的各種輔助因子,但如果僅按照上述方法生成的樹形結(jié)構(gòu)模型往往會出現(xiàn)“過擬合”的現(xiàn)象,這是由于完整的樹形結(jié)構(gòu)模型對于樣本數(shù)據(jù)集特征的描述“過于精確”其原因是受到了噪聲影響,失去了一般代表性從而無法對新數(shù)據(jù)進行下一步精確的分類,針對該問題本研究采用K折交叉驗證的方法進行葉子結(jié)點修剪,在驗證過程中引入一個“可調(diào)錯誤率”的概念,即對某一樹枝的全部葉子節(jié)點增加一個懲罰因子,若該樹枝仍然保持較低錯誤率,則予以保留,否則給予剪除。最終得出一顆適用于本研究且兼顧高復(fù)雜度與低錯誤率的最優(yōu)二叉樹[36]。

2.3.2SVM(Support Vector Machine)模型

支持向量機(SVM)通過引入核函數(shù)低維空間中線性不可分的問題映射到高維空間中,其基本思想是通過使用位于類域邊緣的支持向量,在特征空間中的類之間構(gòu)造分離超平面,在高維空間中最大化不同類別之間的間隔尋求最佳超平面。

1997年,Vnpiak等人[37]基于回歸函數(shù)構(gòu)建了用于預(yù)測連續(xù)變量的支持向量回歸(SVR)模型,即給定訓(xùn)練樣本集T={x(x1,y1),(x2,y2),...(xN,yN)},y1∈R,目的是學(xué)習(xí)到一個f(x)使得其與y盡可能的接近,其數(shù)學(xué)定義如下所示：

f(x)=wTx+b(w,x∈Rd)

(7)

式中，w和b為回歸函數(shù)f(x)的斜率和偏移,d為維度輸入空間。

選擇合適的核函數(shù)對于SVR模型的關(guān)系構(gòu)建非常重要,前人的研究中較為常見的核函數(shù)有4種：1)線性核函數(shù)：K(xi,xj)=xi×xj；2)多項式核函數(shù)：K(xi,xj)=(xi×xj+1)d；3)Sigmoid核函數(shù)：K(xi,xj)=tanh[a(xi×xj)+b]；4)徑向基核函數(shù)(Radial Basis Function,RBF)：K(xi,xj)=exp(-g‖xi-xj‖2)。在不同核函數(shù)的作用下,數(shù)據(jù)在高維特征空間中的分布結(jié)構(gòu)將不同,從而直接決定了SVM與核方法的最終性能[38-39]。

本研究選擇應(yīng)用最廣、普適性最強的RBF作為構(gòu)建模型的核函數(shù),以基于樣本點提取的各種輔助因子作為f(x)函數(shù)中的輸入因子,灌木覆蓋度作為輸出結(jié)果,借助臺灣大學(xué)林智仁(Lin Chih-Jen)教授開發(fā)的Libsvm工具箱[40]進行參數(shù)選擇,該工具箱具有參數(shù)調(diào)節(jié)相對較少以及提供很多默認參數(shù)的特點,通過網(wǎng)格搜索法對RBF核函數(shù)中g(shù)參數(shù)和C(懲罰)參數(shù)進行選擇,主要步驟包括：(1)按照Libsvm軟件包要求的格式準備數(shù)據(jù)集；(2)對數(shù)據(jù)集進行歸一化處理；(3)選用RBF函數(shù)作為核函數(shù)并運用網(wǎng)格搜索算法；(4)采用交叉驗證方法選擇最佳參數(shù)C與g；(5)采用最優(yōu)參數(shù)懲罰因子C和參數(shù)g構(gòu)建支持向量機模型。

2.4 精度評價

本研究樣本數(shù)據(jù)庫共包含400個樣方,每塊樣方均包含了喬木、灌木、草地及沙地四種覆蓋類型。所有樣方采集時間與影像信息提取時間基本保持一致,最大限度的保證數(shù)據(jù)的時間匹配度。隨機選取25%的樣方對模型精度進行驗證,通過均方根誤差(RMSE)、確定系數(shù)R2和估算精度(Estimation Accuracy,EA)[41]3個指標來評價不同模型估算灌木覆蓋度的精度,計算公式為：

(8)

(9)

(10)

3 結(jié)果與討論

3.1 毛烏素沙地灌木覆蓋情況

本研究采用Collect Earth樣本收集器結(jié)合野外樣地調(diào)查的方式構(gòu)建樣本數(shù)據(jù)集,共獲取400個樣方。其中,0—20%覆蓋度的樣方共205個,占總樣方數(shù)的51.3%；20%—40%覆蓋度的樣方共86個,占總樣方數(shù)的21.5%；40%—60%覆蓋度的樣方共63個,占總樣方數(shù)的15.8%；60%—80%覆蓋度的樣方共37個,占總樣方數(shù)的9.2%；80%—100%覆蓋度的樣方共9個,占總樣方數(shù)的2.3%。

圖4 樣本數(shù)據(jù)集灌木覆蓋度地理位置分布圖Fig.4 The distribution of different shrub coverage degree sampling points in Mu Us Sandy Land

從不同覆蓋度樣方空間分布情況來看,如(圖4)所示,覆蓋度較高的區(qū)域主要集中在毛烏素沙地的東北緣,有較多覆蓋度大于60%的樣方分布,中東部地區(qū)主要分布覆蓋度40%—60%的樣方,中部及中西部地區(qū)主要分布覆蓋度低于20%的樣方,零星分布有覆蓋度在20%—40%的樣方。從整體上來看,毛烏素沙地東部地區(qū)樣方平均覆蓋度要高于西部地區(qū)。

3.2 估算變量選擇

本研究涉及到了Landsat-5與Landsat-8兩種數(shù)據(jù),盡管Landsat系列數(shù)據(jù)具有較好地延續(xù)性,但比較GEE平臺上大氣校正后Landsat-5與Landsat-8數(shù)據(jù)一致性仍具有重要意義。因此,我們選取一塊未曾變化的典型區(qū)域(包括喬木、灌木、草地、沙地),分別利用2011年和2013年最大NDVI合成影像進行相關(guān)性分析(共計11746個像元)。通過相關(guān)性分析得到Landsat-5和Landsat-8得到的各對應(yīng)波段及植被指數(shù)決定系數(shù)分別為0.94,0.92,0.88和0.92,如(圖5)所示,均達到極顯著水平(P<0.001),表明兩種傳感器的反射率存在顯著的相關(guān)[42]。其中近紅外和短波紅外波段的斜率趨近于1,紅光波段在4種波段中斜率相對較小但也達到了0.84?？紤]到數(shù)據(jù)年份不一致,上述結(jié)果說明兩種傳感器對應(yīng)波段的反射率具有較高的一致性,用于長時間序列樣本光譜信息提取具有可行性。

圖5 Landsat-5與Landsat-8影像對應(yīng)波段相關(guān)性分析Fig.5 Correlation analysis of corresponding channels between Landsat-5 and Landsat-8

不同類型變量與灌木覆蓋度相關(guān)性分析結(jié)果如(圖6)所示。光譜參數(shù)各波段均在P=0.01水平上顯著相關(guān),其中近紅外波段相關(guān)性最強；各種植被指數(shù)均在P=0.01水平上顯著相關(guān),其中NDVI、MSAVI、SAVI相關(guān)性一致,DVI相關(guān)性相對較低；物候因子中蒸散量以及降水在P=0.01水平上顯著相關(guān),而風速與干旱指數(shù)相關(guān)性低；地形因子中坡度及坡向P=0.05水平上顯著相關(guān),而高程相關(guān)性低?？紤]到覆蓋度與遙感影像之間關(guān)系復(fù)雜,單一變量不能滿足估算要求,因此對四類參數(shù)因子進行綜合選擇,分別選取各類因子中與覆蓋度相關(guān)性較高的因子。針對植被指數(shù),NDVI、SAVI、MSAVI 3個指數(shù)與覆蓋度相關(guān)性一致,并且三者之間自相關(guān)性強,為避免波段信息重疊對模型產(chǎn)生影響,因此只選其中一個(即NDVI)作為模型參數(shù)。同時,相關(guān)研究表明DVI對低覆蓋植被具有較強探測能力,因此也將其選作模型參數(shù)[3]。最終選取近紅外波段、NDVI、DVI、Precipitation、AET以及Slope作為構(gòu)建模型的自變量。

圖6 覆蓋度與17個變量的相關(guān)系數(shù)矩陣Fig.6 Correlation coefficient matrix between the coverage and 17 variablesElevation：高程；Slope：坡度；Aspect：坡向；PDSI：帕爾默干旱指數(shù)，Palmer Drought Severity Index；Precipitation：降水量；Wind speed：風速；AET：實際蒸散量，Actual Evapotranspiration；NDVI：歸一化植被指數(shù)，Normalized Difference Vegetation Index；MSAVI：修正的土壤調(diào)節(jié)植被指數(shù)，Modified Soil-Adjusted Vegetation Index；SAVI：土壤調(diào)節(jié)植被指數(shù)，Soil-Adjusted Vegetation Index；DVI：差值植被指數(shù)，Difference Vegetation Index；Blue：藍光波段；Green：綠光波段；Red：紅光波段；NIR：近紅外波段，Near Infrared；SWIR1：短波紅外波段1，Shortwave Infrared 1；SWIR2：短波紅外波段2，Shortwave Infrared 2；Coverage：覆蓋度

3.3 灌木覆蓋度估算模型構(gòu)建及精度驗證

3.3.1CART灌木覆蓋度估算模型

利用CART對由上述的自變量和因變量組成的樣本進行分析,(圖7)展示了CART樹形結(jié)構(gòu)的一部分,為了防止過擬合,對模型進行K-fold交叉驗證(K=5),樹形結(jié)構(gòu)為：首先以凹凸點NIR≤2309.00為切割點構(gòu)造二叉樹,當NIR≤2309.00時,通過Precipitation構(gòu)造下一顆二叉樹；當NIR>2309.00時,通過NDVI構(gòu)造下一顆二叉樹,樹形結(jié)構(gòu)不斷生長,最終通過剪枝得到一顆最優(yōu)二叉樹。

(表2)統(tǒng)計了CART模型每個輸入自變量對于構(gòu)建模型的貢獻度,可以看出貢獻度最大的是NIR波段,其次是NDVI和DVI,貢獻度分別達到64.60和59.36,AET和Precipitation貢獻度均為13.51而Aspect的貢獻度最低僅為10.97。這個排名在一定程度上體現(xiàn)了每個自變量在構(gòu)建模型中的權(quán)重。其中NIR波段貢獻度最高,原因是灌木屬于木本植物,主要生長方式是通過根和莖因增粗生長從而形成大量的木質(zhì)結(jié)構(gòu),因此,與葉綠素含量較高的草本植物相比,灌木對于紅光波段的吸收能力有限,反射的近紅外光也越少。由此看來,NIR波段是區(qū)分沙地草本植物和灌木植物的最佳因子。

圖7 回歸樹模型結(jié)構(gòu)示意圖Fig.7 The structure of Regression Tree model

表2 自變量對構(gòu)建模型貢獻度Table 2 The contribution degree of different factors

NIR：近紅外波段，Near Infrared；NDVI：歸一化植被指數(shù)，Normalized Difference Vegetation Index；DVI：差值植被指數(shù)，Difference Vegetation Index；AET：實際蒸散量，Actual Evapotranspiration；Precipitation：降水量；Precipitation：降水量

3.3.2SVM灌木覆蓋度估算模型

圖8 SVM模型參數(shù)尋優(yōu)Fig.8 The SVM model parameter optimization

將選出的自變量代入SVM模型中,調(diào)用Libsvm工具箱使用交叉驗證的方法確定支持向量機模型參數(shù)gamma和c。如(圖8)所示,研究中將參數(shù)gamma和c的變化范圍設(shè)置在-10—10之間,本研究選取樹脂復(fù)雜度較小的RBF核函數(shù),它能將樣本非線性地映射到一個更高維的空間,與線性核不同,RBF核函數(shù)能夠處理屬性間的非線性關(guān)系。調(diào)用網(wǎng)格搜索參數(shù)算法并設(shè)置K-fold交叉驗證(K=5),最終得出最優(yōu)懲罰系數(shù)c為512,gamma為1。

3.3.3模型精度驗證

圖9 兩種機器學(xué)習(xí)模型估算和實測覆蓋度對比散點圖Fig.9 Measured and estimated fraction based on two different machine-learning models

模型驗證結(jié)果如(圖9)所示。從圖中可知,兩種機器學(xué)習(xí)模型均可實現(xiàn)一定精度的灌木覆蓋度估算,CART與SVM模型的R2分別達到了0.73和0.72,RMSE分別為13.66%和13.73%,估算精度分別為61.8%和61.6%。盡管兩種模型估算精度總體相差不大,對特定區(qū)間灌木覆蓋度估算效果卻有所差異。當灌木覆蓋度低于50%時,SVM模型估算精度較高,而當覆蓋度大于50%時CART模型估算精度較高。

3.4 時間序列灌木覆蓋度變化分析

圖10 2000—2015年毛烏素沙地灌木覆蓋度分布圖Fig.10 Map of shrub coverage distribution of Mu Us Sandy from 2000—2015

本研究采用CART及SVM算法對2000—2015年間毛烏素沙地灌木覆蓋度進行計算(2012年除外),部分年份估算結(jié)果如(圖10)所示。從總體分布來看,可以發(fā)現(xiàn)兩種模型的估算結(jié)果一致性較強。對2015年估算結(jié)果進行分析表明,灌木覆蓋度高于60%的區(qū)域主要分布在毛烏素沙地的南北兩緣以及東部地區(qū),西部地區(qū)主要以中低覆蓋度為主,中部及西南地區(qū)以流動沙地為主。統(tǒng)計結(jié)果顯示,CART模型及SVM模型估算毛烏素沙地平均灌木覆蓋度分別為28.6%和32.9%,其中,覆蓋度低于25%的區(qū)域分別占總面積的53.6%和36.6%；覆蓋度在25%—50%之間的區(qū)域分別占總面積的31.4%和34.3%；覆蓋度在50%—75%之間的區(qū)域分別占總面積的11.5%和25.7%；覆蓋度大于75%的區(qū)域分別占總面積的3.5%和3.4%。從2000—2015年間時間動態(tài)來看,毛烏素沙地灌木覆蓋度呈整體增長趨勢,其中南北兩緣及東部地區(qū)增長趨勢明顯,中部及中西部地區(qū)長期處于較低覆蓋水平。

為說明本研究所發(fā)展產(chǎn)品特點,我們以毛烏素沙地為例比較了灌木覆蓋度與樹木覆蓋度產(chǎn)品。美國馬里蘭大學(xué)的Hansen團隊發(fā)布了首套全球30m分辨率樹木覆蓋度產(chǎn)品[43],但該產(chǎn)品計算了高于5m的樹木的覆蓋度,該產(chǎn)品目前已成為區(qū)域木本植被豐富性的主要數(shù)據(jù)。然而,在以灌木覆蓋為主的中國北方沙地,該數(shù)據(jù)對大范圍存在的灌木分布無法反映(圖11),考慮到該類地區(qū)在全球及中國的廣泛分布,開展灌木覆蓋度的科學(xué)估算及動態(tài)分析具有重要的應(yīng)用意義。本研究提出的灌木覆蓋度估算技術(shù)體系,具備全球推廣生成相應(yīng)產(chǎn)品的潛力,可以對樹木覆蓋度產(chǎn)品形成有效補充,可為干旱地區(qū)生態(tài)系統(tǒng)管理、科學(xué)研究提供重要支撐。

圖11 灌木覆蓋產(chǎn)品與樹木覆蓋產(chǎn)品對比Fig.11 Comparison of shrubcover and treecover remote sensing imagery products

4 結(jié)論

本文提出了一種利用Collect Earth進行地面樣本采集、通過Google Earth Engine云平臺與機器學(xué)習(xí)方法估算中高分辨率灌木覆蓋度的方法,并以毛烏素沙地為例開展了應(yīng)用研究。主要結(jié)論如下：

(1)本研究采用Collect Earth樣本收集器進行等距采樣的方式構(gòu)建樣本數(shù)據(jù)集,解決了樣本數(shù)據(jù)集偏少或樣本提取精度較低的問題,為今后更大尺度的估算提供了新的手段。

(2)機器學(xué)習(xí)算法對毛烏素沙地灌木覆蓋度估算精度均達到61%以上,其中,分類回歸樹模型確定系數(shù)R2為0.73,均方根誤差為13.66%,支持向量機模型確定系數(shù)R2為0.72,均方根誤差為13.73%。兩種模型估算精度基本持平,但對中低和中高覆蓋度估算能力各有所長,分類回歸樹模型對中高覆蓋度區(qū)域估算精度高于支持向量機模型,而在中低覆蓋度區(qū)域支持向量機模型估算精度高于分類回歸樹模型。

(3)GEE云平臺存儲了公開可用的全球尺度近40年的遙感影像以及其他遙感產(chǎn)品,通過Google數(shù)據(jù)中心眾多超級計算機的快速運算將過去需要花費幾天甚至幾周時間的數(shù)據(jù)處理流程在幾小時甚至分秒間完成。BP級的數(shù)據(jù)存儲提供了更多的數(shù)據(jù)選擇性,大大提高了數(shù)據(jù)質(zhì)量,不同的數(shù)據(jù)源也拉長了數(shù)據(jù)的時間序列,使得今后更大尺度上的干旱、半干旱區(qū)稀疏植被估算成為了可能。