陳 妍, 宋豫秦,*,王 偉

1 北京大學(xué)環(huán)境科學(xué)與工程學(xué)院, 北京 100871 2 中國環(huán)境科學(xué)研究院生物多樣性研究中心,北京 100012

我國是草地資源大國,但多年來的超載過牧、荒地過度開墾和濫挖亂采致使草場植被覆蓋度大幅下降,水土流失和沙漠化日趨嚴(yán)重,牧區(qū)經(jīng)濟(jì)發(fā)展及草地生態(tài)系統(tǒng)的健康因此受到了嚴(yán)重威脅[1]。研究草場變化規(guī)律是揭示草場退化驅(qū)動(dòng)力以及風(fēng)險(xiǎn)評(píng)估的基礎(chǔ),在缺乏長期監(jiān)測數(shù)據(jù)的情況下,借助衛(wèi)星遙感影像不失為有效的方法。當(dāng)前可用的遙感數(shù)據(jù)一般為氣象衛(wèi)星NOAA的1km分辨率AVHRR數(shù)據(jù)[2- 3],terra和aqua衛(wèi)星的250m分辨率MODIS數(shù)據(jù)[4- 5]和陸地衛(wèi)星Landsat系列數(shù)據(jù)。早期研究多采用分辨率較低的AVHRR數(shù)據(jù),MODIS數(shù)據(jù)雖具有較高的完整性,但仍難以滿足長時(shí)間序列研究的需求。Landsat系列衛(wèi)星自1975年開始已有7顆衛(wèi)星相繼運(yùn)行,且其30m的分辨率相較于AVHRR和MODIS來說優(yōu)勢明顯,因而被目前不少研究者所關(guān)注。

利用遙感影像研究草場變化趨勢大致可分為兩種思路,其一是通過監(jiān)督或非監(jiān)督分類,分析草場分布或面積的變化[6- 8],然而面積增減只能體現(xiàn)宏觀變化,難以反映某一具體區(qū)域在時(shí)間序列上的退化或恢復(fù)情況;其二是利用光譜信息計(jì)算植被指數(shù),通過植被指數(shù)變化直接表征生態(tài)系統(tǒng)的變化趨勢[9- 11],相對(duì)而言,這類方法的主要優(yōu)勢是將研究對(duì)象精確至像元尺度。然而干旱半干旱地區(qū)植被蓋度低,植被光譜信息弱,植被指數(shù)往往非常小,在這種情況下用其變化情況說明草場生態(tài)的改善或惡化往往意義不大。而利用植被指數(shù)反演植被蓋度不失為一種可行的替代思路。植被蓋度計(jì)算主要可通過線性光譜混合模型[12- 14]和回歸模型兩種方法實(shí)現(xiàn)。前者一般應(yīng)用于缺乏地面監(jiān)測數(shù)據(jù)的情況之下,其核心是獲得純像元下的地物光譜值。但在實(shí)際應(yīng)用中,特別是植被蓋度較低的草地生態(tài)系統(tǒng)中,典型光譜值很難獲得,且計(jì)算誤差較大[15]。因此,在具備地面監(jiān)測數(shù)據(jù)的情況下,一般通過植被指數(shù)與植被蓋度直接的相關(guān)關(guān)系建立參數(shù)模型,反演植被蓋度,其中線性回歸模型應(yīng)用最為廣泛[16- 17]。然而以往的研究多未充分考慮參數(shù)回歸苛刻的假設(shè)條件,以及多元回歸對(duì)變量間非共線性的要求,這無疑會(huì)降低反演模型的可靠性。解決此問題的途徑之一是尋求預(yù)測效果達(dá)到甚至超越參數(shù)模型的非參數(shù)方法。作為目前預(yù)測效果最好的非參數(shù)回歸模型之一,隨機(jī)森林模型與參數(shù)回歸等方法相比,無需對(duì)變量的正態(tài)性和獨(dú)立性等假設(shè)條件進(jìn)行檢驗(yàn),同時(shí)也不需要考慮多變量的共線問題[18],且運(yùn)算高效、結(jié)果準(zhǔn)確,在環(huán)境以及生態(tài)學(xué)等領(lǐng)域都有著廣泛應(yīng)用[19- 22]。在草場研究方面,隨機(jī)森林雖曾用于植被分類[23],但作為回歸模型的應(yīng)用則幾屬空白。

數(shù)據(jù)源是保證反演模型可靠性的前提,長時(shí)間序列研究會(huì)不可避免地使用到來自不同傳感器的光譜信息,而對(duì)于Landsat系列衛(wèi)星標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品而言,同一時(shí)間段內(nèi)從不同傳感器所提取的光譜信息會(huì)存在一定差異,需要對(duì)非同源數(shù)據(jù)進(jìn)行校正以及標(biāo)準(zhǔn)化處理[24- 25]。但目前草場方面已有的研究既沒有定量分析傳感器差異對(duì)反演結(jié)果造成的影響,也沒有計(jì)算校正和標(biāo)準(zhǔn)化處理對(duì)反演結(jié)果的改進(jìn)效果。

鑒此,本文利用Landsat系列衛(wèi)星TM,ETM+以及OLI的影像和植被蓋度地面監(jiān)測數(shù)據(jù)建立基于隨機(jī)森林回歸的草場植被蓋度反演模型,并將反演結(jié)果與線性回歸結(jié)果進(jìn)行比對(duì),探討其優(yōu)越性。繼之將其應(yīng)用于基于多源數(shù)據(jù)的草場植被蓋度變化趨勢分析,并探討該方法的不確定性,以期為干旱半干旱地區(qū)長時(shí)間序列的草場變化研究提供方法支持。

1 研究區(qū)域

布爾津縣位于新疆西北部阿勒泰地區(qū),屬阿爾泰山脈西南麓,介于86°25′—88°6′E,47°22′—49°11′N之間,總面積10369km2。境內(nèi)景觀異質(zhì)性強(qiáng),北部山區(qū)最高蜂海拔4374m,中部為低山丘陵、河谷地帶,南部為低平戈壁灘，海拔最低處僅為436m(圖1)。全縣生態(tài)系統(tǒng)多樣,北部山區(qū)兼具水源涵養(yǎng)與生物多樣性保護(hù)的功能,南部荒漠生態(tài)系統(tǒng)肩負(fù)防風(fēng)固沙的生態(tài)功能。草地資源豐富,各類草地類型均有分布,農(nóng)業(yè)以畜牧為主,放牧方式仍屬傳統(tǒng)游牧。當(dāng)前北部林區(qū)載畜量快速增加,林牧矛盾突出,南部以荒漠植被為主,生態(tài)環(huán)境脆弱,生態(tài)保護(hù)形勢更加嚴(yán)峻。豐富而多樣的草場分布特征,不僅可使草場植被蓋度反演模型為分析當(dāng)?shù)夭輬鰵v史變化趨勢乃至制定草場保護(hù)和畜牧業(yè)生產(chǎn)政策奠定基礎(chǔ),同時(shí)也可為其他地區(qū)的各類草場研究提供一定借鑒。

2 數(shù)據(jù)獲取與預(yù)處理

為了建立基于隨機(jī)森林回歸的草場植被蓋度反演模型,本文所用之?dāng)?shù)據(jù)主要包括基于地面調(diào)查的草場蓋度數(shù)據(jù)和衛(wèi)星遙感影像。

植被蓋度地面調(diào)查數(shù)據(jù)來源于布爾津縣畜牧獸醫(yī)局草原站提供的2010、2011、2013、2014、2015年6—9月的草本、半灌木及矮小灌木草原樣方調(diào)查表。調(diào)查樣地涵蓋全縣8個(gè)典型草場類型,即低平地草甸、高寒草甸、山地草甸、溫性草甸草原、溫性草原化荒漠、溫性荒漠以及溫性荒漠草原。依據(jù)典型性原則,樣方在能夠代表整個(gè)樣地草原植被、地形及土壤等特征的區(qū)域隨機(jī)布設(shè),監(jiān)測樣地共計(jì)205個(gè),分布及類型如圖1所示,每個(gè)樣地設(shè)置3個(gè)樣方,樣方大小為1m×1m,監(jiān)測結(jié)果用3個(gè)樣方的平均值表示。部分監(jiān)測點(diǎn)位雖然超越了縣界,但并未超越本文所用的遙感影像范圍,因此本文將此類監(jiān)測點(diǎn)也納入分析范圍。植被蓋度是由目測法測定的樣方內(nèi)植物地上部分垂直投影面積占樣方面積的比率,監(jiān)測頻次為每年一次。

遙感數(shù)據(jù)為下載自美國地質(zhì)調(diào)查局(U.S. Geological Survey)網(wǎng)站(http://glovis.usgs.gov/)的Landsat系列衛(wèi)星遙感影像,整個(gè)縣域需要兩景影像覆蓋,行帶號(hào)分別為144/26、144/27,影像時(shí)間及傳感器信息如表1所示。

本文所使用的每一景影像都包含有目前此類研究通用的一級(jí)標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品(簡稱L1數(shù)據(jù))以及經(jīng)過校正和標(biāo)準(zhǔn)化處理的地表反射率數(shù)據(jù)產(chǎn)品(簡稱SR數(shù)據(jù)),用以定量研究傳感器差異對(duì)反演結(jié)果的影響。

本文的研究對(duì)象為草場,因此需要對(duì)影像進(jìn)行監(jiān)督分類,提取出草場范圍。首先,利用ENVI 5.1軟件對(duì)2015年Landsat8 OLI的兩景 L1影像(144/27、144/27)分別進(jìn)行剪裁、輻射定標(biāo)、融合以及FLAASH大氣校正等預(yù)處理。其次初步將影像定義為林地、草地、荒地、水體、冰雪區(qū)五大類型,分別劃定訓(xùn)練樣本,通過最大似然法進(jìn)行分類。經(jīng)過分類后處理,以Google Earth Pro.的高分辨率影像為基礎(chǔ),劃定驗(yàn)證樣本,通過混肴矩陣計(jì)算分類精度。兩景影像(144/27、144/27)的總體分類精度分別為97.7705%和98.8749%,Kappa系數(shù)分別為0.9692和0.9642。合并草地、荒地兩種類型的區(qū)域,通過目視解譯去除其中的耕地和建設(shè)用地區(qū)域,得到布爾津縣草場分布范圍(圖2)。

表1 遙感影像信息列表

圖1 研究區(qū)地形及監(jiān)測點(diǎn)位分布圖 Fig.1 Elevation and monitoring points′ distribution of study area

圖2 草場范圍及抽樣點(diǎn)分布圖 Fig.2 Distribution of grassland and sampling points in study area

3 研究方法

3.1 隨機(jī)森林模型的建立與驗(yàn)證

3.1.1 植被指數(shù)選擇

歸一化植被指數(shù)(NDVI)是草場退化研究中使用最為廣泛的植被指數(shù)[26- 28],本文首先選擇其作為模型的解釋變量之一。由于研究區(qū)南部大片區(qū)域?yàn)橹脖簧w度較低的溫性草原化荒漠、溫性荒漠以及溫性荒漠草原,僅使用NDVI作為解釋變量,反演效果可能并不能達(dá)到預(yù)期,必須考慮土壤背景的影響[29],因此進(jìn)一步引入可用于修正土壤背景敏感性以及氣溶膠散射影響三個(gè)植被指數(shù)——土壤調(diào)節(jié)植被指數(shù)(SAVI)、修正土壤調(diào)節(jié)植被指數(shù)(MSAVI)、增強(qiáng)植被指數(shù)(EVI)。上述植被指數(shù)的計(jì)算方法如式1—式4所示,其中參數(shù)的選取依據(jù)了USGS發(fā)布的植被指數(shù)產(chǎn)品相關(guān)指南[30]。

(1)

(2)

(3)

(4)

3.1.2 模型建立

圖1中監(jiān)測點(diǎn)位的植被蓋度數(shù)據(jù)分別來自2010、2011、2013、2014、2015年6—9月,經(jīng)統(tǒng)計(jì)分析,在6—9月期間,植被蓋度并無顯著區(qū)別,因此在對(duì)應(yīng)年份的6—9月中選擇云量最少的影像加以利用,即利用2010—2011年8月的TM影像和2013—2015年7月的OLI影像計(jì)算相應(yīng)點(diǎn)位的4個(gè)植被指數(shù)。

利用R軟件在205組數(shù)數(shù)據(jù)(植被蓋度及其對(duì)應(yīng)植被指數(shù))中隨機(jī)抽取150組作為訓(xùn)練數(shù)據(jù)用于模型建立。

隨機(jī)森林是由多棵相互沒有關(guān)聯(lián)的決策樹組成的集成決策樹,用{h(X,Θk),k=1, …}表示,其中X為輸入向量,{Θk}為獨(dú)立同分布隨機(jī)向量[31]。與經(jīng)典回顧分析類似,隨機(jī)森林可以解釋多個(gè)自變量(X1、X2…Xk)對(duì)因變量Y的作用。此處將植被蓋度的實(shí)地監(jiān)測值作為因變量Y,NDVI、SAVI、MSAVI和EVI作為自變量。隨機(jī)森林模型的建立通過調(diào)用R語言中“randomForests”程序包[32]來實(shí)現(xiàn)。該方法首先完成兩個(gè)隨機(jī)采樣過程,即通過自助法重采樣技術(shù)有放回的在150組訓(xùn)練數(shù)據(jù)中重復(fù)隨機(jī)抽取150個(gè)訓(xùn)練樣本,未被抽取到的數(shù)據(jù)被稱為“袋外”(out of bag)數(shù)據(jù)。再從NDVI、SAVI、MSAVI和EVI的4個(gè)變量中隨機(jī)選擇若干個(gè)變量(以mtry表示)建立決策樹,模型中mtry的省缺值一般為總變量的1/3,本文有4個(gè)變量,則mtry值為1。最后重復(fù)n次上述過程，生成由n棵決策樹組成的隨機(jī)森林,n值越大預(yù)測越好。隨著n值的增大,袋外數(shù)據(jù)誤差在顯著降低后會(huì)基本保持穩(wěn)定。為了節(jié)省計(jì)算時(shí)間,n值保證袋外數(shù)據(jù)誤差穩(wěn)定即可,本文n值取150。

模型輸出所有決策樹計(jì)算結(jié)果的均值,并通過計(jì)算“解釋方差百分比”來評(píng)定模型預(yù)測能力。對(duì)于各個(gè)自變量對(duì)因變量的影響程度,用方差增量以及節(jié)點(diǎn)純度增量兩個(gè)指標(biāo)來定性表征。前者指將某一變量替換成隨機(jī)變量后對(duì)預(yù)測結(jié)果造成的影響,若用于替換的隨機(jī)變量顯著改變了方差,則認(rèn)為原變量重要性很高;后者從同質(zhì)性增加的角度去表征變量的重要性[33]。隨機(jī)森林模型的具體算法請參照Biau等[34]的文章。

為了將隨機(jī)森林回歸方法與傳統(tǒng)的線性回歸進(jìn)行比較,本文采用同樣的訓(xùn)練數(shù)據(jù)同時(shí)建立了植被蓋度與NDVI、SAVI、MSAVI和EVI的一元線性回歸及多元線性回歸模型。

3.1.3 模型驗(yàn)證

利用205組數(shù)據(jù)經(jīng)過隨機(jī)抽樣剩余的55組數(shù)據(jù)用于模型驗(yàn)證。應(yīng)用前文建立的隨機(jī)森林模型和各個(gè)線性回歸模型分別計(jì)算55個(gè)點(diǎn)位的植被蓋度預(yù)測值,并將實(shí)測值與預(yù)測值相減,根據(jù)差值分布情況比較模型的預(yù)測能力。此外,通過繪制殘差散點(diǎn)圖和正態(tài)Q-Q圖來檢驗(yàn)線性回歸模型方差齊性以及殘差的正態(tài)性,并通過Kappa系數(shù)來判斷變量的共線性。

3.2 模型的應(yīng)用及其不確定性分析

3.2.1 模型應(yīng)用

利用ArcMap 10軟件在研究區(qū)域內(nèi)隨機(jī)生成1000個(gè)相互間隔大于1km的點(diǎn)位,去除草場范圍之外的點(diǎn)位,得到草場抽樣點(diǎn)。為了研究模型在基于多源數(shù)據(jù)的植被蓋度反演領(lǐng)域的應(yīng)用,本研究以2005—2008年ETM+影像、2009—2011年TM影像、2013—2015年OLI影像為基礎(chǔ),分析2005—2015年間(除2012年),植被蓋度的變化趨勢。

首先,剔除研究時(shí)段內(nèi)落在云層和云影中的草場抽樣點(diǎn),最終得到的356個(gè)點(diǎn)位用于本小節(jié)的抽樣分析(圖2)。其次,分別利用上述影像的L1和SR數(shù)據(jù)計(jì)算各個(gè)點(diǎn)位的NDVI、SAVI、MSAVI和EVI值,將其輸入3.1中建立的隨機(jī)森林模型,分別得到監(jiān)測點(diǎn)位基于L1數(shù)據(jù)和SR數(shù)據(jù)的植被蓋度值。最后,將每年356個(gè)監(jiān)測點(diǎn)位的植被蓋度預(yù)測值求平均,繪制出基于L1數(shù)據(jù)和SR數(shù)據(jù)的植被蓋度變化趨勢圖。

3.2.2 不確定性分析

由于2005—2015年間遙感影像源自3個(gè)不同的傳感器,而傳感器的差異會(huì)給反演結(jié)果帶來一定影響,因此本研究進(jìn)一步應(yīng)用下載自美國地質(zhì)調(diào)查局地球資源觀察和科學(xué)中心(Earth Resources Observation and Science Center)的MODIS—MOD13Q1 V005數(shù)據(jù)集建立隨機(jī)森林模型,反演2005—2015年的植被蓋度。由于MODIS數(shù)據(jù)不存在傳感器差異問題,因此可將基于此數(shù)據(jù)反演出的植被蓋度變化趨勢作為參考,將3.2.1中得到的基于L1和SR數(shù)據(jù)的趨勢曲線與之相比較,變化趨勢越接近基于MODIS數(shù)據(jù)的反演結(jié)果,則說明該數(shù)據(jù)中不同傳感器給反演結(jié)果造成的誤差越小。

前文提到SR數(shù)據(jù)是L1數(shù)據(jù)經(jīng)過校正和標(biāo)準(zhǔn)化處理的地表反射率數(shù)據(jù)產(chǎn)品,可在一定程度上克服傳感器的差異對(duì)結(jié)果帶來的影響。為了定量分析SR數(shù)據(jù)對(duì)反演結(jié)果的改進(jìn)效果及不確定性,本研究進(jìn)一步比對(duì)同一時(shí)期內(nèi),基于不同傳感器的反演結(jié)果之間存在的差異。

研究區(qū)范圍草場生長期內(nèi),同時(shí)具有TM和ETM+影像的時(shí)段是2009年8月,同時(shí)具有ETM+和OLI影像的時(shí)段是2015年7月和2016年8月。在草場抽樣點(diǎn)中去除09、15、16年3年內(nèi)云層、云影中的點(diǎn)位以及ETM+數(shù)據(jù)中壞條帶中的點(diǎn)位后,剩余663個(gè)點(diǎn)位用于不確定分析。利用隨機(jī)森林模型分別計(jì)算同一時(shí)段基于不同傳感器的反演結(jié)果并相減,分別繪制基于L1數(shù)據(jù)和SR數(shù)據(jù)的TM—ETM+ 、OLI—ETM+的蓋度反演結(jié)果差異概率密度圖。

4 結(jié)果與討論

4.1 擬合結(jié)果

隨機(jī)森林模型的計(jì)算結(jié)果顯示,基于L1數(shù)據(jù)和SR數(shù)據(jù)而建立的隨機(jī)森林回歸模型輸出的解釋方差百分比分別為69.82%和72.42%,因此后者的預(yù)測效果優(yōu)于前者;重要性方面,各個(gè)植被指數(shù)并未表現(xiàn)出明顯的差別(表2)。

表2 隨機(jī)森林回歸模型變量重要性

由于L1數(shù)據(jù)未經(jīng)過標(biāo)準(zhǔn)化處理,TM和OLI數(shù)據(jù)計(jì)算的EVI值不在同一數(shù)量級(jí),無法用于模型的建立；NDVI：歸一化植被指數(shù),Normalized Difference Vegetation Index;SAVI：土壤調(diào)節(jié)植被指數(shù),Soil Adjusted Vegetation Index;MSAVI：修正土壤調(diào)節(jié)植被指數(shù),Modified Soil Adjusted Vegetation Index;EVI：增強(qiáng)植被指數(shù),Enhanced Vegetation Index

本研究一共建立線性回歸模型23個(gè),其中P值顯著的模型及相應(yīng)參數(shù)如表3所示。由表3可知,對(duì)于一元線性回歸模型來說L1數(shù)據(jù)和SR數(shù)據(jù)的R2差別不大;對(duì)于多元回歸來說所有的L1數(shù)據(jù)都無法得到顯著結(jié)果;從R2的數(shù)值上來看SR數(shù)據(jù)的多元回歸的效果要明顯優(yōu)于一元回歸。

4.2 模型比較

為了進(jìn)一步分析上述模型的預(yù)測能力,本文將55個(gè)驗(yàn)證點(diǎn)位的監(jiān)測值和預(yù)測值作差,結(jié)果的分布情況如圖3所示。首先,雖然L1和SR數(shù)據(jù)的一元回歸效果在表3中無法通過R2判斷,但通過比較紅色和綠色箱線圖可明顯看出,SR數(shù)據(jù)實(shí)測值和預(yù)測值的差異明顯小于L1數(shù)據(jù);第二,隨機(jī)森林回歸模型的預(yù)測能力略優(yōu)于多元回歸模型,且這兩者明顯優(yōu)于一元線性回歸;第三,所有線性回歸的結(jié)果較真實(shí)值來說普遍偏小,但隨機(jī)森林的預(yù)測偏差相對(duì)均衡。

在本文所建立的所有模型中,預(yù)測效果最好的是基于L1和SR數(shù)據(jù)的隨機(jī)森林模型,以及由變量SAVI-MSAVI建立的SR數(shù)據(jù)二元回歸模型。

表3 線性回歸模型結(jié)果

由于L1標(biāo)準(zhǔn)化處理,TM和OLI數(shù)據(jù)計(jì)算的EVI值不在同一數(shù)量級(jí),無法用于模型的建立

圖3 植被蓋度的實(shí)測值與預(yù)測值差異Fig.3 The difference between monitoring values and predicted valuesNDVIL1：基于L1數(shù)據(jù)NDVI指數(shù)的線性回歸模型,Linear model based on NDVI from level 1 standard product;SAVIL1：基于L1數(shù)據(jù)SAVI指數(shù)的線性回歸模型,Linear model based on SAVI from level 1 standard product; MSAVIL1：基于L1數(shù)據(jù)MSAVI指數(shù)的線性回歸模型,Linear model based on MSAVI from level 1 standard product;NDVISR：基于SR數(shù)據(jù)NDVI指數(shù)的線性回歸模型,Linear model based on NDVI from surface reflectance product;SAVISR：基于SR數(shù)據(jù)SAVI指數(shù)的線性回歸模型,Linear model based on SAVI from surface reflectance product;MSAVISR：基于SR數(shù)據(jù)MSAVI指數(shù)的線性回歸模型,Linear model based on MSAVI from surface reflectance product;EVISR：基于SR數(shù)據(jù)EVI指數(shù)的線性回歸模型,Linear model based on EVI from surface reflectance product;SRSM：基于SR數(shù)據(jù)SAVI及MSAVI指數(shù)的線性回歸模型,Linear model based on SAVI and MSAVI from surface reflectance product;SRSE：基于SR數(shù)據(jù)SAVI及EVI指數(shù)的線性回歸模型,Linear model based on SAVI and EVI from surface reflectance product;SRNSE：基于SR數(shù)據(jù)NDVI、SAVI及EVI指數(shù)的線性回歸模型,Linear model based on NDVI,SAVI and EVI from surface reflectance product;rfL1：基于L1數(shù)據(jù)的隨機(jī)森林回歸模型,Random forest model based on level 1 standard product;rfSR：基于SR數(shù)據(jù)的隨機(jī)森林回歸模型,Random forest model based on surface reflectance product

圖4為SRSM方差齊性和殘差正態(tài)性檢驗(yàn)結(jié)果。理想狀態(tài)下,殘差-擬合值關(guān)系圖上的點(diǎn)應(yīng)隨機(jī)分布,而正態(tài)QQ圖中的點(diǎn)應(yīng)分布于虛線之上,但SRSM模型并不能滿足其基本假設(shè)條件,特別是殘差的正態(tài)分布性。不僅如此,本研究所有其他的線性回歸模型也都不能很好的滿足其基本假設(shè)。此外,SRSM的各變量應(yīng)滿足非共線性,即Kappa值小于100,但經(jīng)過計(jì)算該模型的Kappa值高達(dá)5570。從對(duì)模型的診斷結(jié)果看,SRSM模型雖然預(yù)測效果較好,但應(yīng)用于植被蓋度的反演存在較大的缺陷。

圖4 SRSM模型的假設(shè)檢驗(yàn) Fig.4 Hypotheses testing of SRSM modelSRSM：基于SR數(shù)據(jù)SAVI及MSAVI指數(shù)的線性回歸模型,Linear model based on SAVI and MSAVI from surface reflectance product

綜上,本文認(rèn)為隨機(jī)森林模型相較于線性回歸模型來說不但不受諸多假設(shè)條件和非共線性要求的制約,并能達(dá)到更好的預(yù)測效果,因此更適合作為植被蓋度的反演模型。

4.3 2005—2015年草場植被蓋度變化趨勢及不確定性分析

圖5 基于MODIS、SR及L1數(shù)據(jù)的抽樣點(diǎn)蓋度平均值 Fig.5 Mean vegetation cover of sampling points based on MODIS, SR and L1 data L1： 一級(jí)標(biāo)準(zhǔn)數(shù)據(jù)產(chǎn)品Level1 Standard Data Product;SR: 經(jīng)過校正和標(biāo)準(zhǔn)化處理的地表反射率數(shù)據(jù)產(chǎn)品Surface Reflectance Data Products,

基于L1、SR以及MODIS數(shù)據(jù)分別得到的草場蓋度平均值變化趨勢如圖5所示。除2006年之外,SR和L1數(shù)據(jù)反演得到的植被蓋度變化情況基本保持了一致的變化趨勢,但2005—2009年,也即數(shù)據(jù)源為ETM+的時(shí)間段內(nèi),基于L1的反演結(jié)果要明顯小于SR數(shù)據(jù)的反演結(jié)果。SR數(shù)據(jù)與不存在傳感器差異問題的MODIS數(shù)據(jù)反演所得的植被蓋度值相比,差異小于5%,但L1數(shù)據(jù)的偏差較大,ETM+數(shù)據(jù)的預(yù)測結(jié)果普遍偏小,2006年的植被蓋度均值甚至偏小超過20%。

為了分析傳感器差異給反演結(jié)果造成的影響,本文進(jìn)一步對(duì)比了同一時(shí)期內(nèi)源自不同傳感器光譜信息反演得到的植被蓋度,對(duì)比結(jié)果如圖6所示。該圖顯示了2009年ETM+數(shù)據(jù)與TM數(shù)據(jù),以及2015—2016年ETM+數(shù)據(jù)與OLI數(shù)據(jù)分別反演得到的植被蓋度的差異。

圖6 同時(shí)段不同數(shù)據(jù)源的預(yù)測結(jié)果差異 Fig.6 Difference in predicted values based on different dada sources in the same periodSRrfETM_TM：SR數(shù)據(jù)基于ETM+傳感器與TM傳感器隨機(jī)森林反演結(jié)果的差異;SRrfETM_OLI：SR數(shù)據(jù)基于ETM+傳感器與OLI傳感器隨機(jī)森林反演結(jié)果的差異;L1rfETM_TM：L1數(shù)據(jù)基于ETM+傳感器與TM傳感器隨機(jī)森林反演結(jié)果的差異;L1rfETM_OLI：L1數(shù)據(jù)基于ETM+傳感器與OLI傳感器隨機(jī)森林反演結(jié)果的差異

由圖6的概率密度曲線可知,同一時(shí)段L1數(shù)據(jù)不同傳感器對(duì)計(jì)算結(jié)果差異巨大,特別是ETM+數(shù)據(jù)的預(yù)測結(jié)果普遍小于TM數(shù)據(jù),上述結(jié)果與圖5的結(jié)論一致。從具體計(jì)算的數(shù)值來看：對(duì)于SR數(shù)據(jù),幾乎90%的測試點(diǎn)位上不同傳感器對(duì)計(jì)算結(jié)果產(chǎn)生的影響都能落在(-10%,11%)的區(qū)間之中;但對(duì)于L1數(shù)據(jù),90%測試點(diǎn)位的ETM+與OLI的差異落在(-30%,11%)范圍內(nèi),而ETM+與TM的差異更是落在了(-40%,9%)之中。由此可見,不同傳感器的差異在L1數(shù)據(jù)上表現(xiàn)的尤為突出,該數(shù)據(jù)不適于基于多元數(shù)據(jù)的草場植被蓋度變化趨勢研究,即便是經(jīng)過校正以及標(biāo)準(zhǔn)化處理的SR數(shù)據(jù),反演結(jié)果也存在約±10%的不確定性。

4.4 討論

就反演方法而言,以往類似的研究多基于相關(guān)分析線性回歸以及多項(xiàng)式模型[16- 17,35- 37],且模型的篩選通常僅僅通過P值以及決定系數(shù)R2來評(píng)判，而不對(duì)模型預(yù)測結(jié)果做進(jìn)一步驗(yàn)證。從本文的線性回歸結(jié)果看,L1數(shù)據(jù)和SR數(shù)據(jù)的一元線性模型都具有顯著性且R2水平相當(dāng),在以MSAVI作為解釋變量時(shí),L1數(shù)據(jù)的擬合效果甚至要優(yōu)于SR數(shù)據(jù)(表3)。但從圖3中可以清楚地看到前者的誤差要遠(yuǎn)遠(yuǎn)大于后者。可見,P值以及R2并未幫助研究者對(duì)比出SR數(shù)據(jù)的優(yōu)越性。更為重要的是,上述模型都屬于參數(shù)估計(jì)類,若不對(duì)其基本假設(shè)進(jìn)行驗(yàn)證,就無法保證所得結(jié)論的可靠性。本研究顯示,所有的線性模型都未能通過方差齊性和殘差正態(tài)性檢驗(yàn)。而遺憾的是,迄今大多數(shù)研究都忽略了這種統(tǒng)計(jì)學(xué)的基本要求。

需要強(qiáng)調(diào)的是,除了模型本身之外,數(shù)據(jù)源也是影響預(yù)測效果的重要因素。由圖5中基于MODIS數(shù)據(jù)反演得到的變化趨勢可知,近10年來草場蓋度的平均值總的變化幅度不足20%左右,而2006年L1數(shù)據(jù)反演所得之蓋度平均值與MODIS數(shù)據(jù)計(jì)算結(jié)果的偏差就已經(jīng)超過了20%。該問題在圖6中體現(xiàn)得更加突出。這便意味著,傳感器差異帶來的影響甚至可能超過草場變化本身。雖然本研究所用的SR數(shù)據(jù)與傳統(tǒng)研究常用的L1數(shù)據(jù)相比,已經(jīng)在一定程度是減少了這類影響[38],但仍未將其完全消除。因此,在涉及長時(shí)間序列草場植被變化以及驅(qū)動(dòng)力的研究中,對(duì)這種影響進(jìn)行量化就顯得尤為必要。

本文的局限性主要體現(xiàn)在以下幾個(gè)方面：首先,使用的實(shí)地監(jiān)測數(shù)據(jù)為1m×1m樣方的測量值,而遙感影像為30m×30m精度,如果監(jiān)測點(diǎn)周圍數(shù)10m范圍內(nèi)不能保持均一性,擬合結(jié)果會(huì)因此受到影響;其次,本研究所涉及的草場種類眾多,但植被指數(shù)計(jì)算涉及的參數(shù)均采用了USGS網(wǎng)站提供產(chǎn)品指南中的推薦值,后續(xù)研究中應(yīng)對(duì)不同類型草場分別設(shè)定參數(shù),以提高結(jié)果的可靠程度;再次,本研究初步估算了2005—2015年草場蓋度的平均值,今后可進(jìn)一步應(yīng)用隨機(jī)森林模型分析草場植被蓋度的空間變化特征;最后,本文僅初步研究了傳感器差異給反演結(jié)果造成的影響,如何通過歸一化處理等方法減小這種差異將是未來研究中需要探索的問題。

5 結(jié)論

1) 本文通過隨機(jī)森林回歸方法建立了基于Landsat系列衛(wèi)星數(shù)據(jù)的草場植被蓋度反演模型,并應(yīng)用該模型分析了近10年來布爾津縣草場蓋度的變化趨勢。

2) 通過與線性回歸反演模型的對(duì)比可知,隨機(jī)森林回歸方法不但能規(guī)避線性回歸苛刻的假設(shè)條件,同時(shí)可超越絕大多數(shù)線性回歸模型的預(yù)測能力。

3) 模型不確定性方面,Landsat ETM+的標(biāo)準(zhǔn)數(shù)據(jù)反演得到的植被蓋度較之TM和OLI數(shù)據(jù)普遍偏小;地表反射率數(shù)據(jù)雖然可以大幅降低Landsat系列不同傳感器對(duì)反演結(jié)果的影響,但反演得到的植被蓋度仍存在(-10%,11%)的不確定性。

致謝：北京大學(xué)生命科學(xué)學(xué)院王昊老師提供MODIS數(shù)據(jù),中國環(huán)境科學(xué)研究院張玉波博士和付剛在遙感影像處理方面給予幫助,大理大學(xué)東喜馬拉雅資源與環(huán)境研究所任國鵬老師為本研究提出建議,特此致謝。

參考文獻(xiàn)(References):

[1] 樊江文, 鐘華平, 陳立波, 張文彥. 我國北方干旱和半干旱區(qū)草地退化的若干科學(xué)問題. 中國草地學(xué)報(bào), 2007, 29(5): 95- 101.

[2] 黃敬峰, 桑長青. 天山北坡中段天然草場牧草產(chǎn)量遙感動(dòng)態(tài)監(jiān)測模式. 干旱區(qū)資源與環(huán)境, 1993, 7(2): 53- 59.

[3] 李京, 陳晉, 袁清. 應(yīng)用NOAA/AVHRR遙感資料對(duì)大面積草場進(jìn)行產(chǎn)草量定量估算的方法研究. 自然資源學(xué)報(bào), 1994, 9(4): 365- 368.

[4] Hmimina G, Dufrêne E, Pontailler J Y, Delpierre N, Aubinet M, Caquet B, De Grandcourt A, Burban B, Flechard C, Granier A, Gross P, Heinesch B, Longdoz B, Moureaux C, Ourcival J M, Rambal S, Saint André L, Soudani K. Evaluation of the potential of MODIS satellite data to predict vegetation phenology in different biomes: an investigation using ground-based NDVI measurements. Remote Sensing of Environment, 2013, 132: 145- 158.

[5] Geerken R, Batikha N, Celis D, DePauw E. Differentiation of rangeland vegetation and assessment of its status: field investigations and MODIS and SPOT VEGETATION data analyses. International Journal of Remote Sensing, 2005, 26(20): 4499- 4526.

[6] Chen S B, Rao P. Land degradation monitoring using multi-temporal Landsat TM/ETM data in a transition zone between grassland and cropland of northeast China. International Journal of Remote Sensing, 2008, 29(7): 2055- 2073.

[7] Feng Y, Lu Q, Tokola T, Liu H, Wang X. Assessment of grassland degradation in Guinan county, Qinghai Province, China, in the past 30 years. Land Degradation & Development, 2009, 20(1): 55- 68.

[8] Geerken R, Zaitchik B, Evans J P. Classifying rangeland vegetation type and coverage from NDVI time series using Fourier Filtered Cycle Similarity. International Journal of Remote Sensing, 2005, 26(24): 5535- 5554.

[9] Cai H Y, Yang X H, Xu X L. Human-induced grassland degradation/restoration in the central Tibetan Plateau: the effects of ecological protection and restoration projects. Ecological Engineering, 2015, 83: 112- 119.

[10] Wessels K J, Prince S D, Carroll M, Malherbe J. Relevance of rangeland degradation in semiarid northeastern South Africa to the nonequilibrium theory. Ecological Applications, 2007, 17(3): 815- 827.

[11] 魏小琴. 阿勒泰地區(qū)NDVI變化及其主要驅(qū)動(dòng)因子分析[D]. 烏魯木齊: 新疆農(nóng)業(yè)大學(xué), 2015.

[12] R?der A, Udelhoven T, Hill J, Del Barrio G, Tsiourlis G. Trend analysis of Landsat-TM and -ETM+ imagery to monitor grazing impact in a rangeland ecosystem in Northern Greece. Remote Sensing of Environment, 2008, 112(6): 2863- 2875.

[13] Li J Y, Yang X C, Jin Y X, Yang Z, Huang W G, Zhao L N, Gao T, Yu H D, Ma H L, Qin Z H, Xu B. Monitoring and analysis of grassland desertification dynamics using Landsat images in Ningxia, China. Remote Sensing of Environment, 2013, 138: 19- 26.

[14] 楊強(qiáng), 王婷婷, 陳昊, 王運(yùn)動(dòng). 基于MODIS EVI數(shù)據(jù)的錫林郭勒盟植被覆蓋度變化特征. 農(nóng)業(yè)工程學(xué)報(bào), 2015, 31(22): 191- 198.

[15] 呂長春, 王忠武, 錢少猛. 混合像元分解模型綜述. 遙感信息, 2003, (3): 55- 58, 60- 60.

[16] 郭輝, 黃粵, 李向義, 包安明, 宋洋, 孟凡浩. 基于多尺度遙感數(shù)據(jù)的塔里木河干流地區(qū)植被覆蓋動(dòng)態(tài). 中國沙漠, 2016, 36(5): 1472- 1480.

[17] 馬中剛, 孫華, 王廣興, 林輝, 佘宇晨, 鄒琪. 基于Landsat 8-OLI的荒漠化地區(qū)植被覆蓋度反演模型研究. 中南林業(yè)科技大學(xué)學(xué)報(bào), 2016, 36(9): 12- 18.

[18] 李欣海. 隨機(jī)森林模型在分類與回歸分析中的應(yīng)用. 應(yīng)用昆蟲學(xué)報(bào), 2013, 50(4): 1190- 1197.

[19] Cutler D R, Edwards Jr T C, Beard K H, Cutler A, Hess K T, Gibson J, Lawler J J. Random forests for classification in ecology. Ecology, 2007, 88(11): 2783- 2792.

[20] 金宇, 周可新, 方穎, 劉欣. 基于隨機(jī)森林模型預(yù)估氣候變化對(duì)動(dòng)物物種潛在生境的影響. 生態(tài)與農(nóng)村環(huán)境學(xué)報(bào), 2014, 30(4): 416- 422.

[21] 王盼, 陸寶宏, 張瀚文, 張巍, 孫銀鳳, 季妤. 基于隨機(jī)森林模型的需水預(yù)測模型及其應(yīng)用. 水資源保護(hù), 2014, 30(1): 34- 37, 89- 89.

[22] Rodríguez-Veiga P, Wheeler J, Louis V, Tansey K, Balzter H. Quantifying forest biomass carbon stocks from space. Current Forestry Reports, 2017, 3(1): 1- 18.

[23] Mansour K, Mutanga O. Classifying increaser species as an indicator of different levels of rangeland degradation using WorldView- 2 imagery. Journal of Applied Remote Sensing, 2012, 6(1): 063558.

[24] Chander G, Markham B L, Helder D L. Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO- 1 ALI sensors. Remote Sensing of Environment, 2009, 113(5): 893- 903.

[25] Vogelmann J E, Gallant A L, Shi H, Zhu Z. Perspectives on monitoring gradual change across the continuity of Landsat sensors using time-series data. Remote Sensing of Environment, 2016, 185: 258- 270.

[27] Pattison R R, Jorgenson J C, Raynolds M K, Welker J M. Trends in NDVI and tundra community composition in the arctic of NE Alaska between 1984 and 2009. Ecosystems, 2015, 18(4): 707- 719.

[28] 毛志春, 宋宇, 李蒙蒙. 基于MODIS反演數(shù)據(jù)的河套地區(qū)荒漠化研究. 北京大學(xué)學(xué)報(bào): 自然科學(xué)版, 2015, 51(6): 1102- 1110.

[29] Liu B, You G Y, Li R, Shen W S, Yue Y M, Lin N F. Spectral characteristics of alpine grassland and their changes responding to grassland degradation on the Tibetan Plateau. Environmental Earth Sciences, 2015, 74(3): 2115- 2123.

[30] USGS. Landsat surface reflectance-derived spectral indices, Version 3.3. 2016. [2016- 12- 13] https://landsat.usgs.gov/sites/default/files/documents/si_product_guide.pdf

[31] Breiman L. Random forests. Machine Learning, 2001, 45(1): 5- 32.

[32] Liaw A, Wiener M. Classification and regression by randomForest. R News, 2002, 2(3): 18- 22.

[33] Echeverry-Galvis M A, Peterson J K, Sulo-Caceres R. The social nestwork: tree structure determines nest placement in Kenyan weaverbird colonies. PLoS One, 2014, 9(2): e88761.

[34] Biau G, Scornet E. A random forest guided tour. TEST, 2016, 25(2): 197- 227.

[35] Brinkmann K, Dickhoefer U, Schlecht E, Buerkert A. Quantification of aboveground rangeland productivity and anthropogenic degradation on the Arabian Peninsula using Landsat imagery and field inventory data. Remote Sensing of Environment, 2011, 115(2): 465- 474.

[36] 萬紅梅, 李霞, 董道瑞. 基于多源遙感數(shù)據(jù)的荒漠植被覆蓋度估測. 應(yīng)用生態(tài)學(xué)報(bào), 2012, 23(12): 3331- 3337.

[37] 董建軍, 牛建明, 張慶, 康薩如拉, 韓芳. 基于多源衛(wèi)星數(shù)據(jù)的典型草原遙感估產(chǎn)研究. 中國草地學(xué)報(bào), 2013, 35(6): 64- 69.

[38] Claverie M, Vermote E F, Franch B, Masek J G. Evaluation of the Landsat- 5 TM and Landsat- 7 ETM+ surface reflectance products. Remote Sensing of Environment, 2015, 169: 390- 403.