基于MODIS-Landsat時空融合的陜北黃土高原植被覆蓋變化研究

劉詠梅，馬 黎,黃 昌,凱 楠

(1.西北大學 城市與環境學院, 陜西 西安 710127；2.水利部 黃土高原水土流失過程與控制重點實驗室，河南 鄭州 450003)

植被作為氣候和土壤之間的聯系在全球氣候變化和物質能源流動的交換中發揮著重要作用[1],植被覆蓋變化監測已成為生態環境研究的重心[2-5]。陜北黃土高原地表破碎、溝壑縱橫,生態環境脆弱[6],植被恢復對該區域的水土保持起著重要作用[7]。自1999年退耕還林生態建設工程實施以來,該區域的植被狀況發生了很大變化,因此,開展陜北黃土高原植被覆蓋的動態監測對黃土高原生態環境建設和退耕還林工程效果評價具有重要意義。

多位學者對陜北黃土高原植被覆蓋變化及影響因子開展了系列研究,李登科等[6]利用GIMMS，SPOTVEGETATION兩種數據對陜北黃土丘陵溝壑區不同時期的植被覆蓋變化及其與氣候的關系進行了分析;宋富強等[8]運用MODIS數據,從不同土地利用類型和不同坡度NDVI植被指數變化兩方面分析了退耕還林對陜北黃土高原植被變化的影響;湛青青等[9]運用Landsat TM數據分別計算1990，2000，2010年的植被覆蓋度,分析該區植被覆蓋的空間變化情況。上述研究開展時必須在時空尺度選擇上進行取舍:MODIS，SPOTVEGETATION等數據可以較好地表達植被覆蓋的時序變化,但250～1 000m的低空間分辨率對于陜北黃土高原的復雜地形來說，在空間細節的表達上不夠精確[10];Landsat等傳感器獲取的30m多光譜數據空間分辨率較高但有較長的重訪周期,使得植被監測缺少時間上的連續性。時空數據融合綜合了Landsat數據高空間分辨率和MODIS數據高時間分辨率的優點[11],Zhu等[12]學者提出了一種考慮像元反射率時間變化特征的增強型自適應時空融合模型ESTARFM(enhanced spatial and temporal adaptive reflectance fusion model),模型提高了對非均勻異質景觀的預測精度,在復雜地形條件ESTARFM模型的預測結果與真實影像間的相關性最高[13]。

本文基于ESTARFM模型,采用Landsat和MODIS數據融合獲得2008—2016年6～8月高時空分辨率的NDVI時序數據,并利用像元二分模型估算植被覆蓋度,結合土地利用分類數據和氣溫、降水等氣象數據,分析陜北黃土高原植被覆蓋時空變化規律及影響因素,探討時空分辨率融合影像在植被動態監測中的適用性,以期為該區域水土保持、環境保護等相關應用提供基礎數據和科學參考。

1 研究區與數據

1.1 研究區概況

陜北黃土高原地處于陜西省北部,地理位置在35°21′～39°34′N,107°15′～111°14′E之間,海拔高度為800m～1 800m。東隔黃河與山西省毗鄰,南邊與陜西省渭南市、銅川市、咸陽市相接,西以子午嶺為界與寧夏回族自治區、甘肅省接壤,北接內蒙古自治區[14]，如圖1所示。氣候上處在暖溫帶大陸性季風半濕潤氣候向溫帶半干旱氣候的過渡地帶[8],四季變化交替明顯,年平均溫度介于7.5～12.3℃之間,年降水量為350～660mm[15]。陜北黃土高原北部屬于風沙高原區,南部屬于黃土高原區,地表包括農田、森林、草地、荒漠等類型,地形復雜破碎,溝壑縱橫[8]。區內植被類型多樣,且分布差異較明顯,主要分布有溫帶針葉林、落葉闊葉林、溫帶落葉灌叢、溫帶中生草甸、溫帶草原和暖溫帶草甸草原等植被類型,并呈現由東南向西北的落葉闊葉林向溫帶草原過渡的特征[16]。

圖1 研究區地理位置示意圖Fig.1 The geographical location of the study area

1.2 數據源與預處理

本研究采用的陜北黃土高原Landsat數據分別為2008—2012年的Landsat5 TM和2013—2016年的Landsat8 OLI影像(http://www.gscloud.cn)。MODIS數據采用500m空間分辨率8天合成的MOD09A1數據產品(http://www.nasa.gov),數據獲取時間與Landsat數據相對照且間隔不超過3天。土地利用數據采用基于Landsat TM/OLI影像解譯的2015年全國土地利用數據。氣象數據來自寒區旱區科學數據中心(http://westdc.westgis.ac.cn)的SWAT模型中國大氣同化驅動數據集(CMADS V1.0),空間分辨率為1/3°×1/3°的逐天格網數據,時間尺度為2008—2016年。

利用ENVI5.1軟件對Landsat影像進行輻射定標、大氣校正和幾何校正。運用MRT(MODIS Reprojection Tool)對MOD09A1數據進行投影轉換,重采樣到Landsat的30m空間分辨率。對兩種影像分別拼接后，用研究區矢量邊界裁剪。由于影像數據來自不同的傳感器,需調整影像的波段順序相一致(表1)。為確保數據大小、影像柵格行列號等完全一致，裁剪時需在統一的圖像處理軟件下裁剪。

表1 Landsat和MODIS影像波段對照表Tab.1 Corresponding bands between Landsat and MODIS

2 研究方法

2.1 ESTARFM時空融合模型

ESTARFM時空融合模型旨在將Landsat影像與MODIS影像融合,從而得到具有高空間高時間分辨率特征的數據。對于指定區域,假設在同一日期由不同傳感器獲得的遙感數據可以相互比較,并在數據預處理后相互間有關聯。ESTARFM的主要思想是利用關聯來融合多源數據,同時最小化系統偏差。模型以預測像元為中心設置一定大小的滑動窗口,對窗口內像元利用權重函數進行卷積運算,得到中心像元預測值[11]?；瑒哟翱谠谡跋裆现鹨灰苿?從而得到需要預測的影像(式(1))[12]。

RM(xi,yi,tt,B))

(1)

式中:RL和RM分別為Landsat和MODIS影像數據;w為滑動窗口大小;(xw/2,yw/2)為待測像元的中心位置;(xiyi)為第i個相似像元的位置;tt為影像獲取時間;B為影像波段;N為相似像元的數量;Wi是由空間、時間、光譜的距離共同決定的權重大小;Vi為轉換系數[12]。選擇ta和tb時的兩對MODIS和Landsat數據分別通過公式(1)結合tp時的MODIS數據預測tp時的Landsat數據,結果記為RLa(xw/2,yw/2,tp,B)和RLb(xw/2,yw/2,tp,B)。結合這兩種預測結果,預測的中心像元反射率會更準確。以越靠近預測時間則權重更大為準則,該權重計算為式(2)[12]。最后中心像元反射率的預測值為公式(3)[12]。

βt=

(t=a,b)

(2)

(3)

ESTARFM模型可以通過BI(Blend-then-Index)和IB(Index-then-Blend)兩種方式實現[17],BI即將原始遙感影像輸入到模型中,融合得到影像后再進行后續指數的計算[17];IB通過原始遙感影像計算研究所需的指數,再將計算好的指數輸入到模型中,以模擬預測日期的指數[17]。由于運用ESTARFM模型時IB方式具有更高的準確性且計算時間較短[17],所以本文通過ESTARFM-IB方式進行歸一化植被指數的模擬。即采用經過ta和tb時刻的兩對MODIS NDVI和Landsat NDVI數據,結合tp時刻的MODIS NDVI數據預測同一時刻的高空間分辨率的Landsat NDVI數據。

2.2 植被覆蓋度計算

經前期處理共獲得陜北黃土高原2008—2016年6～8月共27期高時空分辨率NDVI數據,采用像元二分模型估算植被覆蓋度,其計算公式如式(4)，(5)[18]。

NDVI=(NIR-R)/(NIR+R)，

(4)

Fc=(NDVI-NDVIsoil)/(NDVIveg-NDVIsoil)

(5)

式中:NIR為近紅外波段反射率,R為紅波段反射率;NDVIveg為純植被覆蓋時的NDVI值,NDVIsoil為純土壤地表的NDVI值[19]。根據整幅影像的NDVI灰度分布頻率累計表,選擇5%和95%的置信度區間內上下限累計頻率的值分別為NDVIsoil和NDVIveg。

3 結果與分析

3.1 ESTARFM模型融合結果與分析

由于Landsat數據存在云量大和數據缺失等問題,依據相關研究成果[10-13]和實際檢驗，本研究中采取預測日期前或后的兩對MODIS和Landsat數據進行融合。例如2010年6月27日和2011年7月17日的兩組MODIS NDVI和Landsat NDVI數據作為基礎數據,通過ESTARFM-IB融合方式得到的2011年6月17日的ESTARFM NDVI數據,用2011年6月17日真實的Landsat NDVI數據作為參考為例,驗證融合數據結果的準確和可靠性。

圖2(a)～(b)顯示整體上融合影像與真實影像的相似程度較高,NDVI像元值在空間分布上基本一致,清晰再現了30m空間分辨率的空間細節信息,數據紋理結構明顯。相關性分析表明(圖2(c)),ESTARFM NDVI與真實Landsat NDVI的散點大部分分布在y=x附近,斜率K為0.874 3,兩者的相關性系數達到0.972 8(P<0.01)。因此融合影像能夠較好的反映同期Landsat NDVI的影像光譜信息,可以作為后續計算研究區植被覆蓋度的基礎數據。

圖2 預測NDVI與真實NDVI對比分析Fig.2 Comparison between predicted NDVI and real NDVI

3.2 陜北黃土高原植被覆蓋時空變化特征

3.2.1 植被覆蓋年際變化特征 計算陜北黃土高原2008—2016年每年6～8月的平均覆蓋度并繪制其變化曲線圖。圖3顯示,2008—2016陜北黃土高原整體植被覆蓋變化呈顯著增加的趨勢,平均植被覆蓋度由2008年的0.402增加到2016年的0.639,增速為0.031 2/a,增幅達58.9%,線性增長趨勢呈顯著水平。研究區植被覆蓋曲線變化呈現波動變化上升的狀態,但其年際間波動較大。表明研究區植被覆蓋狀況得到了明顯改善,也顯示出當地生態環境的脆弱性。

圖3 2008—2016陜北黃土高原植被覆蓋變化曲線Fig.3 The variation of annual FVC of Loess Plateau in Northern Shaanxi Province during 2008—2016

3.2.2 植被覆蓋空間變化特征 根據陜北黃土高原的植被生長狀況將植被覆蓋度分為極低覆蓋度(0～20%)、低覆蓋度(20～40%)、中覆蓋度(40～60%)、中高覆蓋度(60～80%)、高覆蓋度(80%～100%)5個等級[9]。2008，2012和2016年該區域的植被覆蓋度空間分布如圖4所示,呈現較為明顯的區域差異,由東南向西北遞減,高值區主要分布在陜北黃土高原東南部地區,低值區分布在西北部地區。2008年北部定邊、靖邊、橫山、榆陽、神木、子洲、綏德等縣區主要為極低、低、中覆蓋度植被區域,南部甘泉、富縣、洛川、黃陵、黃龍、宜川等縣區為中高、高植被覆蓋度。陜北黃土高原北部,植被稀少,稀樹灌木叢草原較多;南部特別是橋山、黃龍山林區保有大量天然次生林,植被繁茂[2]。2012年全區極低植被覆蓋度區域面積占比縮小,中覆蓋度面積擴大,約占整個研究區的60%。到2016年中高覆蓋度面積占比明顯增多,主要增多的區域位于研究區中部、中東部和東北部各縣區。陜北黃土高原植被覆蓋度等級結構有所好轉。

為研究近10年來陜北黃土高原植被覆蓋變化趨勢的空間差異,對每個像元的植被覆蓋度與時間進行一元線性回歸分析。圖5顯示陜北黃土高原大部分區域的植被覆蓋度均表現增加趨勢,東北部府谷、神木、佳縣和中部米脂、綏德到中南部的清澗、延長等縣的植被覆蓋度多呈中度和明顯增加的趨勢,該區域經過近10年退耕還林等生態工程的持續實施,林地、草地等面積的快速增加使得植被恢復效果較好。南部的黃龍、黃陵、富縣和中北部的定邊、榆陽等縣區的部分區域多呈現基本不變或略微減少的趨勢。南部地區由于天然次生林較多,植被生長狀況較為穩定,植被變化不明顯;中北部部分地區屬于長城沿線風沙區,在退耕還林工程實施中多采用封禁措施[8]，加之降水稀少使得植被的恢復效果相對較低。對比圖4,植被覆蓋增加趨勢較為明顯的區域多為中、中高植被覆蓋度,基本不變或降低的區域多為高覆蓋度。運用自然間斷點分級法將變化斜率分為7類(表2),由統計結果看來2008—2016年研究區78%的土地面積植被覆蓋度在逐年增加。

圖4 2008，2012，2016年陜北黃土高原植被蓋度空間分布圖Fig.4 Spatial distribution of vegetation coverage of Loess Plateau in Northern Shaanxi Province in 2008, 2012 and 2016

圖5 2008—2016陜北黃土高原植被覆蓋度變化斜率分級圖Fig.5 The classified changing slope of FVC of Loess Plateau in Northern Shaanxi Province(2008—2016)

Tab.2 Statistical results of vegetation cover variation trend of Loess Plateau in Northern Shaanxi Province(2008—2016)

植被覆蓋變化斜率范圍變化程度面積/km2百分比/% slope<-0.045明顯減少1760.22 -0.045≤slope<-0.027中度減少3600.45 -0.027≤slope<-0.009略微減少1 6802.10 -0.009≤slope<0.009基本不變15 28019.1 0.009≤slope<0.027略微增加33 49541.8 0.027≤slope<0.045中度增加20 26425.3 0.045≤slope明顯增加8 72010.9

圖6 不同土地利用類型植被蓋變化曲線Fig.6 The change of vegetation coverage trend of different types of land use types

3.2.3 不同土地利用類型植被覆蓋變化 目前針對陜北黃土高原地區不同植被類型覆蓋度變化特征的研究相對較少。本文對該區域的主要植被類型:有林地、灌木林、疏林地、其他林地、高覆蓋度草地、中覆蓋度草地和低覆蓋度草地的植被覆蓋變化進行分析可知(圖6),2008—2016年陜北黃土高原7種主要植被類型6～8月的蓋度均值介于0.372～0.744之間,依次為有林地>灌木林>高覆蓋度草地>疏林地>其他林地>中覆蓋度草地>低覆蓋度草地,其年均數值分別為0.744 9，0.695 5，0.695 2，0.643 6，0.502 0，0.473 3，0.372 8。各類型平均蓋度的曲線變化特點基本一致,均呈波動增加的趨勢。計算各類型年均植被蓋度與年份的相關系數可知灌木林和疏林地與年份的相關性最高(P<0.05),表明這兩種類型的植被狀況恢復較好,其余各類型與年份的相關性分析均沒有通過顯著性檢驗。有林地和高覆蓋度草地的標準差較小,分別為0.028和0.021,表明在9a間這兩種植被類型的覆蓋變化較為穩定和平衡,沒有明顯的增加或減低的趨勢;疏林地和其他林地的標準差較高,約為0.05,這兩種類型屬于植被分布稀疏的區域,主要由各種園地、苗圃構成,受人類活動影響而植被蓋度年際波動較大。

3.3 氣溫和降水對植被覆蓋變化的影響

植被的變化受自然和人為因素的共同影響,其中自然因素中氣溫和降水是對植被變化影響較大的兩個因素[8]。對研究區71個站氣象資料的分析表明(圖7),2008—2016年該區6～8月平均氣溫呈微弱的下降趨勢(slope=-0.021℃/a),6～8月降水量呈增加趨勢(slope=12.85mm/a),同時氣溫和降水的年際變化較大。

圖7 2008—2016年6～8月陜北黃土高原降水、氣溫的變化趨勢Fig.7 Trends of precipitation and temperature of loess plateau in Northern Shaanxi Province from june to august(2008—2016)

本文從空間分布的角度進一步分析該區域6～8月平均植被蓋度與同期氣溫、降水的關系,以像元為單位計算植被覆蓋度與氣溫和降水的相關系數(圖8,表3)。植被蓋度與降水間呈顯著正相關的區域主要分布在陜北黃土高原中部的定邊至吳堡和東北部的佳縣、府谷等縣區(P<0.05),占全區總面積的10.1%。中部個別縣區的部分地區相關系數呈極顯著相關(P<0.01),占全區總面積的6.35%。中北部區域的主要植被類型為草地[8],在降水豐沛的時段長勢較好,地表植被覆蓋度就會增加,反之亦然。與降水相關系數為負的區域主要分布在研究區南部的黃陵、黃龍、宜川各縣,其中0.33%的區域呈極顯著負相關、1.79%呈顯著負相關。這些區域主要是針闊混交林[8],南部林區有大面積的水源涵養林,這些森林植被短期內受降水的影響較小。植被蓋度與氣溫呈顯著和極顯著正相關的區域主要分布在南部各縣區,僅占全區總面積1.36%。植被蓋度與氣溫呈顯著和極顯著負相關的區域極少,僅占0.78%,零星分布于北部地區。相比研究區植被蓋度與同時期溫度變化的關系,研究區植被蓋度對降水因子的響應更敏感。

圖8 植被覆蓋度與降水、氣溫的相關性分析圖Fig.8 Correlation coefficient maps of vegetation coverage with precipitation and temperature

Tab.3 The vegetation coverage correlation significance percentage with temperature and precipitation of Loess Plateau in Northern Shaanxi Province(2008—2016)

相關性植被覆蓋度與降水/%植被覆蓋度與氣溫/% 極顯著負相關0.330.15 顯著負相關1.790.63 不顯著負相關25.3047.83 不顯著正相關56.1450.04 顯著正相關10.101.18 極顯著正相關6.350.18

4 結 論

1)ESTARFM模型得到的黃土高原地區NDVI時序數據整體和局部細節上與真實數據相似性程度很高,獲得了較高精度的光譜與空間信息,表明通過ESTARFM模型得到的高時空分辨率NDVI數據能夠更好地應用于陜北黃土高原地形復雜地區的植被覆蓋變化監測研究。

2)2008—2016年陜北黃土高原地區植被覆蓋度整體呈顯著增加的趨勢,空間上呈現東南高西北低的特點。植被覆蓋結構整體上好轉,趨勢分析表明近10年來研究區植被蓋度整體呈現穩定的增加狀態。該區域7種主要植被類型的蓋度均呈波動增加趨勢,灌木林和疏林地的增長趨勢最為顯著;有林地和高覆蓋度草地的波動變化較為穩定,疏林地和其他林地的植被覆蓋變化年際波動大。

3)2008—2016年陜北黃土高原6～8月的降水量呈上升趨勢,而氣溫則呈現微弱的下降趨勢?？臻g相關性分析表明,整體上6～8月的降水和氣溫對植被覆蓋度的影響區域差異較為明顯。植被覆蓋度與降水呈顯著相關的面積為18.57%,與氣溫呈顯著相關的面積僅為2.04%,整體上植被蓋度變化與氣候因素的相關性較弱。黃土高原實施退耕還林還草工程以來,植被恢復狀況較好,也使得植被覆蓋對氣溫、降水等自然因子的影響敏感度減弱。

5 討 論

在黃土高原地區植被覆蓋變化的相關研究中, 多采用MODIS NDVI 250m/16d， SPOT VGT NDVI 1km/10d， GIMMS NDVI 8km/15d等低空間分辨率的時序數據。 其中MODIS， SPOT數據多用在2000年以后的黃土高原相關區域植被覆蓋變化的研究中[2,8,23,25], 而GIMMS數據多用在較長時間尺度的研究[1,6,24,26]中。 上述研究數據在空間精度和時間連續性上并不能達到統一, 因而研究結論多偏向區域整體植被覆蓋的變化或在研究時段的選擇上有所限制。 而對于陜北黃土高原特殊復雜的地形特點來說,基于ESTARFM模型獲取的高時空分辨率的遙感影像數據相對較好的解決了以上問題, 這種數據同時具備高度的空間細節表達力和快速時序變化能力, 為在陜北黃土高原開展高精度的植被動態監測研究提供了支撐。 同時研究發現連續的時空分辨率合成Landsat數據提高了MODIS等時序數據的空間分辨率, 對于陜北黃土高原復雜的地形特點來說, 可以更加精確的得到每一個柵格像元對應的NDVI數值, 確保了植被蓋度趨勢演變過程的可靠性。 ESTARFM模型為綜合多源遙感數據在黃土丘陵溝壑地形復雜地區開展高精度的植被覆蓋遙感動態監測提供了較優的數據源。本文針對研究區植被覆蓋的整體時空變化、各等級植被覆蓋變化和不同土地利用類型植被覆蓋變化進行了探討,今后在對高時空分辨率數據的挖掘和細化方面需要做進一步研究。