基于多源遙感數(shù)據(jù)橡膠林識(shí)別及其種植適宜區(qū)研究

蔣永泉,史正濤,高書鵬

(云南師范大學(xué) 地理學(xué)部,昆明 650500)

橡膠因其天然獨(dú)特的物化特性而成為我國(guó)重要的戰(zhàn)略基礎(chǔ)物資之一[1]。橡膠樹于20世紀(jì)初引入云南省種植[2],目前,已成為我國(guó)最大的天然橡膠生產(chǎn)基地。20世紀(jì)60年代云南省瑞麗市開始大面積種植橡膠林,截至2017年,全市共種植橡膠7 053.55hm2,總產(chǎn)值高達(dá)6 000萬(wàn)元余[3]。長(zhǎng)久以來,橡膠林種植多以犧牲天然林為主,導(dǎo)致了當(dāng)?shù)匚锓N多樣性減少、景觀連通性降低、水土流失等[4]諸多問題。因此,準(zhǔn)確監(jiān)測(cè)橡膠林的種植狀況對(duì)于平衡經(jīng)濟(jì)發(fā)展和生態(tài)環(huán)境保護(hù)具有重要意義。目前,遙感技術(shù)在橡膠林識(shí)別中取得了長(zhǎng)足的進(jìn)展,MODIS數(shù)據(jù)是識(shí)別橡膠林使用最廣的遙感數(shù)據(jù)源之一[5-6],但由于MODIS空間分辨率低,識(shí)別結(jié)果較差,很難將斑塊較小的橡膠林識(shí)別出來。隨著遙感技術(shù)的發(fā)展,多源遙感數(shù)據(jù)識(shí)別地物已成為遙感應(yīng)用的趨勢(shì)和熱點(diǎn)。與單一數(shù)據(jù)源相比,多源數(shù)據(jù)可充分利用各種數(shù)據(jù)的高光譜、高時(shí)空分辨率構(gòu)建各種特征,提高識(shí)別精度以滿足應(yīng)用需求。高書鵬等[7]基于時(shí)空數(shù)據(jù)融合算法獲取的ETM+，OLI，Sentinel-2A時(shí)序數(shù)據(jù)提取西雙版納的橡膠林,Chen等[8]和Kou等[9]都基于雷達(dá)數(shù)據(jù)和Landsat數(shù)據(jù)分別提取海南省和西雙版納的橡膠林,結(jié)果顯示,同時(shí)運(yùn)用2種數(shù)據(jù)源均比單一數(shù)據(jù)源的提取精度高。因此,多源遙感數(shù)據(jù)可以彌補(bǔ)單一數(shù)據(jù)存在的缺陷,實(shí)現(xiàn)不同數(shù)據(jù)源之間的優(yōu)勢(shì)互補(bǔ),將多源遙感數(shù)據(jù)進(jìn)行融合,在一定程度上優(yōu)化了分類結(jié)果,提高了識(shí)別精度[10],在橡膠林監(jiān)測(cè)管理中具有較好的應(yīng)用前景。

原種植于亞馬遜河流域的橡膠樹,要求高溫、多雨、靜風(fēng)、光照充足的氣候條件,而我國(guó)屬于非傳統(tǒng)植膠區(qū),適宜種植區(qū)非常有限[11]。瑞麗市作為我國(guó)北緣橡膠分布區(qū),最冷月平均氣溫僅有12.5℃,是世界上種植橡膠氣溫最低的地區(qū),經(jīng)常受到寒潮和寒害的影響[12]。1970—1971年、1975—1976年以及1982—1983年瑞麗市橡膠林都分別遭受了不同程度的寒害,使當(dāng)?shù)亟?jīng)濟(jì)受到重創(chuàng)[13-14]。隨著市場(chǎng)需求量擴(kuò)大和膠價(jià)上漲,2004年后,膠農(nóng)開始大面積砍伐森林種植橡膠林并向不適宜種植區(qū)擴(kuò)張,導(dǎo)致了橡膠生長(zhǎng)速度慢、產(chǎn)膠量低[15-16]。因此,掌握橡膠林的空間分布特征及其種植適宜區(qū)對(duì)規(guī)劃管理具有重要的意義。本文綜合利用Sentinel-2和MODIS-NDVI多源遙感數(shù)據(jù)識(shí)別橡膠林,該方法數(shù)據(jù)易于獲取,可實(shí)現(xiàn)大范圍的監(jiān)測(cè),基于識(shí)別結(jié)果,對(duì)瑞麗市橡膠林空間分布特征及其種植適宜區(qū)分析,可為相關(guān)部門決策和生態(tài)恢復(fù)提供參考。

1 研究區(qū)概況

研究區(qū)位于云南省西部瑞麗市,地處橫斷山脈高黎貢山余脈的向南延伸部分,地勢(shì)西北高東南低,位于北緯23.38′～24.14′,東經(jīng)97.31′～98.02′之間,其西北、西南、東南三面與緬甸接壤[17]。該區(qū)屬于南亞熱帶季風(fēng)性氣候區(qū),冬無嚴(yán)寒,夏無酷暑,全年分旱雨兩季,年平均氣溫21℃,年降水量1 394.8mm,年平均日照2 330h[18]。瑞麗市總面積為1 020km2,其中森林面積占總面積的46.67%,主要森林類型有熱帶山地雨林,亞熱帶季風(fēng)常綠闊葉林,亞熱帶山地落葉闊葉林,針葉林和竹林。高等植物類型豐富,共有高等植物284科,943屬,3 733種,其中有國(guó)家級(jí)保護(hù)植物45種,省級(jí)保護(hù)植物24種[19]。全市植被覆蓋度較高,其中橡膠林是當(dāng)?shù)胤N植面積最大的人工林[20]。

2 研究方法

2.1 數(shù)據(jù)準(zhǔn)備及預(yù)處理

1) 選用的Sentinel-2遙感影像數(shù)據(jù)可通過歐空局(European Space Agency,ESA)網(wǎng)站(https://scihub.copernicus.eu/dhus/#/home)免費(fèi)下載獲取。大氣校正通過Sen2Cor[21]小插件完成,其核心算法是大氣輻射傳輸模型libRadtran[22],目的是將大氣上層表觀反射率轉(zhuǎn)化為大氣下層地表反射率。根據(jù)研究需要提取出10m的紅光波段、綠光波段、藍(lán)光波段和近紅外波段,以及4個(gè)20m的植被紅邊波段并將其重采樣為10m,然后裁剪出研究區(qū)大小的區(qū)域以備后續(xù)使用。

2) MOD13Q1是由美國(guó)國(guó)家航空航天局(NASA)提供的16d合成的250m空間分辨率NDVI產(chǎn)品,可在NASA網(wǎng)站(https://search.earthdata.nasa.gov/search)免費(fèi)獲取。獲取的原始MODIS數(shù)據(jù)格式為HDF,投影類型與Sentinel-2數(shù)據(jù)不一致,因此本文采用MRT工具做投影轉(zhuǎn)換,將其地理坐標(biāo)和投影類型分別轉(zhuǎn)為WGS84和UTM,然后裁剪出研究區(qū)大小的區(qū)域。為了得到10m×10m像元大小的NDVI數(shù)據(jù),將23個(gè)波段的NDVI數(shù)據(jù)用立方卷積重采樣處理,使其與Sentinel-2多光譜數(shù)據(jù)具有相同的像元大小。

3) 樣本數(shù)據(jù)包括訓(xùn)練樣本和驗(yàn)證樣本,本文基于Google Earth高清影像,對(duì)研究區(qū)進(jìn)行目視解譯,將土地利用類型分為5類,分別為橡膠林、自然林、耕地、建設(shè)用地和水體。根據(jù)各地類的空間分布情況、復(fù)雜程度以及比例,一共隨機(jī)選取了2 563個(gè)點(diǎn),其中橡膠林527個(gè)、自然林764個(gè)、耕地681個(gè)、建設(shè)用地337、水體254個(gè)。并將各地類按1∶1的比例分開分別作為訓(xùn)練樣本和驗(yàn)證樣本。

2.2 地類的NDVI時(shí)序變化特征和光譜特征分析

為分析瑞麗市橡膠林與其他地類在不同時(shí)段的生長(zhǎng)變化差異情況,選用2019年的MOD13Q1數(shù)據(jù),分析5個(gè)地類在23個(gè)時(shí)相的NDVI響應(yīng)關(guān)系(圖1)。從總體上看,不同地類對(duì)NDVI的響應(yīng)不同,橡膠林在整個(gè)時(shí)序變化過程中與其他地類存在明顯的差異。第49d(2月)后,橡膠林開始萌芽進(jìn)入返青期,NDVI值上升；第129 d(5月)—177d(6月)橡膠林的NDVI值急劇降低,這是由于夏季研究區(qū)內(nèi)多云雨天氣,云量較大NDVI值偏低引起的;在第241—289 d(9—10月)NDVI達(dá)到最大值;第305 d(11月)后,橡膠林的NDVI值開始緩慢下降,此時(shí)橡膠林開始落葉。

不同地類對(duì)不同波段的吸收和反射光譜存在一定的差異,圖2為各地類光譜反射率均值曲線圖。在藍(lán)光波段、綠光波段和紅光波段橡膠林易與水體混淆；在植被紅光波段1,橡膠林的光譜值快速升高,且可與其他地類較好的分離；在植被紅邊2和3波段,橡膠林易與建設(shè)用地混淆；在植被紅邊波段4,橡膠林的光譜值達(dá)到最高,與其他有較好的可分性。因此,本文選用植被紅邊波段1和4用于區(qū)分橡膠林和其他地類。

圖1 地類NDVI時(shí)間序列曲線圖

注:“B2”藍(lán)光波段；“B3”綠光波段；“B4”紅光波段；“B5”植被紅邊波段1；“B6”植被紅邊波段2；“B7”植被紅邊波段3；“B8A”植被紅邊波段4

2.3 地物分類特征提取

基于Sentinel-2多光譜遙感影像提取分類特征,包括植被指數(shù)、紋理特征和光譜特征共23個(gè)特征,以及MOD13Q1時(shí)間序列數(shù)據(jù)的23個(gè)NDVI值(每16d一個(gè)),如表1所示。 歸一化植被指數(shù)(NDVI)是反映植被生長(zhǎng)狀態(tài)的重要因子[23],與植被的覆蓋度有密切的相關(guān)性。NDVI的計(jì)算公式如下[24]:

(1)

式中:ρNIR和ρred分別表示近紅外波段和紅波段的表光反射率。

灰度共生矩陣(GLCM)是20世紀(jì)70年代初由Haralick等[25]提出的,是通過灰度的空間特性來描述紋理的方法,能用來反映地物類別在空間特征的差異[26]?；叶裙采仃嚬灿?4種特征指標(biāo)來定量描述紋理,本文選用了最常用的5種,分別為均值(Mean)、同質(zhì)性(Homogeneity)、相關(guān)性(Correlation)、熵(Entropy)和對(duì)比度(Contrast)。

(2)

(3)

(4)

(5)

(6)

式中:ux,uy,σx,σy,分別為行和列的均值和標(biāo)準(zhǔn)差；n為共生矩陣的階數(shù)；i,j為共生矩陣的坐標(biāo)；P為(i,j)處的共生矩陣數(shù)值。

表1 地物分類特征表

3 結(jié)果與分析

3.1 橡膠林識(shí)別結(jié)果

本文采用支持向量機(jī)分類器進(jìn)行分類,分別對(duì)Sentinel-2多光譜遙感影像、MOD13Q1-NDVI時(shí)序數(shù)據(jù)以及多源遙感數(shù)據(jù)進(jìn)行分類,結(jié)果如圖3所示。從總體上看,3種不同數(shù)據(jù)源得到的瑞麗市橡膠林分布在空間上趨于一致。僅用單一的MOD13Q1-NDVI分類(圖3(a)),不能較準(zhǔn)確地識(shí)別出一些面積較小且分布零散的橡膠林。僅用單一Sentinel-2數(shù)據(jù)分類(圖3(b)),雖然能將面積較小的橡膠林識(shí)別出來,但得到的結(jié)果比較破碎。圖3(c)為綜合運(yùn)用多源數(shù)據(jù)的分類結(jié)果,由于綜合了兩種數(shù)據(jù)源在空間和時(shí)間上的優(yōu)勢(shì),總的分類結(jié)果得到了很大的改善,能在準(zhǔn)確識(shí)別出一些破碎地塊的同時(shí),又改善了Sentinel-2數(shù)據(jù)提取橡膠林結(jié)果的破碎程度,有效地避免了單一數(shù)據(jù)源的分類錯(cuò)誤。

圖3 橡膠林識(shí)別結(jié)果

本文基于混淆矩陣分別計(jì)算分類的制圖精度(Producer Accuracy,PA)、用戶精度(User Accuracy,UA)、總體精度(Overall Accuracy,OA)以及Kappa系數(shù)來評(píng)價(jià)結(jié)果的可行性。從表2可以看出,綜合應(yīng)用多源遙感數(shù)據(jù)能夠?qū)崿F(xiàn)瑞麗市橡膠林的高精度識(shí)別,優(yōu)于僅用單一遙感數(shù)據(jù)源分類方法,聯(lián)合多源遙感數(shù)據(jù)能夠有效提高分類精度。綜合運(yùn)用多源數(shù)據(jù),橡膠林的制圖精度PA為89.73%、用戶精度UA為73.98%、總體精度OA為91.41%、Kappa系數(shù)為0.76；與僅用MODIS數(shù)據(jù)相比,PA,UA,OA和Kappa系數(shù)分別提高了18.25%,11.52%,6.10%和0.19；與僅用Sentinel-2數(shù)據(jù)相比,PA,UA,OA和Kappa系數(shù)分別提高了16.73%,3.39%,3.21%和0.12。

表2 橡膠林識(shí)別精度評(píng)價(jià)對(duì)比表

3.2 橡膠林空間分布特征

從分類結(jié)果可以看出,瑞麗市橡膠林主要分布在的西南和東北方向上,對(duì)結(jié)果進(jìn)行聚類分析處理去除破碎的錯(cuò)分對(duì)象,基于30m的DEM數(shù)據(jù)分析橡膠林空間分布特征。瑞麗市DEM數(shù)據(jù)的海拔高度范圍為682～1 997m,將其劃分為6個(gè)等級(jí):682～900m,900～1 100m,1 100～1 200m,1 200～1 300m,1 300～1 400m,≥1 400m；根據(jù)我國(guó)水利部頒發(fā)《土壤侵蝕分類分級(jí)標(biāo)準(zhǔn)》[27]坡度范圍,將瑞麗市坡度劃分為6個(gè)等級(jí):0～5°,5～8°,8～15°,15～25°,25～35°,≥35°。得到表3、表4的統(tǒng)計(jì)結(jié)果。

表3 瑞麗市橡膠林在海拔高度上的分布特征

表4 瑞麗市橡膠林在坡度上的分布特征

截至2019年底,瑞麗市橡膠總種植面積為7 693.83hm2,占土地總面積的7.54%,其中有6 433.11hm2種植在海拔為682～900m高度范圍內(nèi),其次是種植在900～1 000m范圍內(nèi),分別占橡膠林總種植面積的83.61%和15.26%；主要種植的坡度為5～15°,占橡膠林總種植面積的83.61%。

3.3 橡膠種植適宜區(qū)

中華人民共和國(guó)農(nóng)業(yè)部發(fā)布的《橡膠樹種植基地建設(shè)標(biāo)準(zhǔn)》(NY/T221-2016)[28]根據(jù)橡膠樹對(duì)溫度、水分、風(fēng)速和光照等膠樹生長(zhǎng)要求以及產(chǎn)膠潛力,將云南省西部海拔≥900m高度以上的地區(qū)劃分為不宜作為常規(guī)膠園宜林地。依據(jù)此標(biāo)準(zhǔn)的海拔高度變化將瑞麗市劃分為橡膠適宜種植和不適宜種植2個(gè)區(qū)域,然后對(duì)2019年瑞麗市植橡膠林適宜區(qū)進(jìn)行統(tǒng)計(jì)分析(圖4),瑞麗市有6 433.11hm2橡膠種植在適宜種植區(qū),占橡膠林總種植面積的83.61%,已有1 260.72hm2橡膠林向不適宜種植區(qū)擴(kuò)張,占橡膠林總種植面積的16.39%。

圖4 瑞麗市橡膠林適宜種植區(qū)

4 結(jié)論

1) 綜合Sentinel-2多光譜影像和MOD13Q1 NDVI時(shí)間序列數(shù)據(jù)的分類方法,可充分發(fā)揮2種遙感數(shù)據(jù)源的不同特征優(yōu)勢(shì),與僅用單一數(shù)據(jù)的結(jié)果相比,有效提高了橡膠林的識(shí)別精度,實(shí)現(xiàn)了瑞麗市橡膠林高精度識(shí)別,綜合運(yùn)用多源遙數(shù)據(jù)橡膠林識(shí)別結(jié)果的制圖精度PA為89.73%、用戶精度UA為73.98%、總體精度OA為91.41%,Kappa系數(shù)為0.76。

2) 截至2019年底,瑞麗市橡膠種植面積為7 693.83hm2,有83.61%的橡膠林種植在海拔為684～900m高度范圍內(nèi),占比最大；其次是種植在900～1 100m高度范圍內(nèi),占橡膠林總種植面積大的15.26%。主要種植的坡度為5～15°,占橡膠林總種植面積的83.61%。

3) 根據(jù)農(nóng)業(yè)部發(fā)布的《橡膠樹種植基地建設(shè)標(biāo)準(zhǔn)》將瑞麗市的土地按海拔高度劃分為適宜種植區(qū)和不適宜種植區(qū),截至2019年底,瑞麗市有6 433.11hm2橡膠種植在適宜區(qū),占橡膠林總面積的83.61%,已有1 260.72hm2橡膠林向不適宜種植區(qū)擴(kuò)張,占橡膠林總面積的16.39%。