2001—2015年海南島橡膠林物候時空變化特征分析

胡盈盈， 戴聲佩， 羅紅霞， 李海亮，李茂芬， 鄭倩， 禹萱， 李寧

(1.中國熱帶農(nóng)業(yè)科學(xué)院科技信息研究所/海南省熱帶作物信息技術(shù)應(yīng)用研究重點(diǎn)實(shí)驗(yàn)室，海口 571101；2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)遙感重點(diǎn)實(shí)驗(yàn)室，北京 100081； 3.中國熱帶農(nóng)業(yè)科學(xué)院環(huán)境與植物保護(hù)研究所/海南儋州熱帶農(nóng)業(yè)生態(tài)系統(tǒng)國家野外科學(xué)觀測研究站，海口 571101)

0 引言

植被物候是植物長期適應(yīng)季節(jié)性變化而形成的生長發(fā)育節(jié)律，是目前指示植物生長發(fā)育的季節(jié)特征和生理活動最常用的一個指標(biāo)[1]。森林植被作為陸地生態(tài)系統(tǒng)的重要組成部分，在全球物質(zhì)和能量循環(huán)中起著重要作用。開展森林植被物候研究不僅有助增進(jìn)森林植被對氣候變化響應(yīng)的理解，對評估植被生產(chǎn)力和碳平衡也具有重要意義[2-3]。海南島是重要農(nóng)產(chǎn)品天然橡膠保護(hù)區(qū)劃定的重點(diǎn)區(qū)域，天然橡膠是一種必不可缺的戰(zhàn)略物資，監(jiān)測海南島橡膠林物候變化特征對于重要農(nóng)產(chǎn)品保護(hù)以及生態(tài)保護(hù)具有重要意義。

長時間序列物候監(jiān)測和評估成為目前植被物候研究的主要方向，而遙感作為一門對地觀測綜合學(xué)科，由于其范圍廣、時序長、可獲取性強(qiáng)，最近幾十年來被廣泛應(yīng)用在植被物候監(jiān)測領(lǐng)域[4-6]。MODIS遙感影像由于其描述植被信息受到的干擾較小，時間序列較長和分辨率較高，許多學(xué)者利用其歸一化植被指數(shù)(normalized difference vegetation index，NDVI)對植被物候變化進(jìn)行研究[7-8]。目前，已有學(xué)者對橡膠物候進(jìn)行了研究，Razak等[9]使用Landsat NDVI數(shù)據(jù)提取了馬來西亞橡膠林物候參數(shù)，實(shí)現(xiàn)了橡膠林制圖； Kou等[10]和Zhai等[11]分別利用Landsat數(shù)據(jù)和物候觀測資料研究了西雙版納橡膠林物候特征； 陳匯林等[12]和陳小敏等[13]使用MODIS時間序列年內(nèi)均值變化及鄰域變化量來提取海南島橡膠林物候參數(shù)； 田光輝等[14]基于2001—2011年海南島橡膠林物候特征提取了橡膠林分布信息。上述研究或局限于橡膠物候的應(yīng)用，或研究物候參數(shù)不全面，對橡膠林生長季長度鮮有涉及。

海南島作為一個獨(dú)立的生態(tài)系統(tǒng)，在全球氣候變化背景下，橡膠林物候會出現(xiàn)什么樣的規(guī)律特征，需要進(jìn)一步探討。因此，本研究基于2001—2015年MODIS NDVI數(shù)據(jù)，采用Savitzky-Golay(S-G)濾波法、動態(tài)閾值法提取海南島橡膠林展葉期開始日、落葉期開始日和生長季長度，對15 a間橡膠林的物候變化進(jìn)行動態(tài)監(jiān)測，并基于橡膠林像元分析物候變化趨勢和顯著性，以期摸清海南島橡膠林物候時空變化規(guī)律，為其他研究提供基礎(chǔ)數(shù)據(jù)和參考。

1 研究區(qū)概況及數(shù)據(jù)源

1.1 研究區(qū)概況

海南島(N18°10′～20°10′，E108°37′～111°30′)地處中國最南端，是海南省的主體，是我國第二大島，全島面積約3.54萬km2。海南島是我國橡膠林的主要種植基地之一，中西部和中部地區(qū)多為橡膠林種植區(qū)，如儋州市、澄邁縣、白沙黎族自治縣、瓊中黎族苗族自治縣等[15-16]。海南島屬熱帶季風(fēng)氣候，年均氣溫為22～26 ℃，年均降水量在1 600 mm以上，地勢中間高四周低。綜合地理位置、氣候、膠林面積、膠齡、經(jīng)營主體等多種因素，在海南島西部、南部和中部分別選擇3個典型樣區(qū)。海南島典型樣區(qū)及數(shù)字高程模型(digital elevation model,DEM)如圖1所示。

1.2 數(shù)據(jù)來源與預(yù)處理

遙感數(shù)據(jù)為2001—2015年MOD13Q1 NDVI數(shù)據(jù)，空間分辨率為250 m，時間分辨率為16 d，每年有23期影像數(shù)據(jù)，數(shù)據(jù)格式為EOS-HDF。使用MODIS Reprojection Tool工具對海南島原始影像進(jìn)行鑲嵌、裁剪以及投影轉(zhuǎn)換； 然后對NDVI數(shù)據(jù)去除無效值、歸一化處理，將值控制在-1～1； 物候參數(shù)的提取在Timesat軟件中完成，需將原始數(shù)據(jù)HDF格式轉(zhuǎn)換成Timesat有效的DAT格式，歸一化處理和格式轉(zhuǎn)換均在ArcGIS軟件中完成。

海南島是天然橡膠保護(hù)區(qū)，橡膠為多年生熱帶經(jīng)濟(jì)作物，經(jīng)濟(jì)壽命長，種植至少在30 a以上[17-18]，受《森林法》保護(hù)，不可大面積砍伐。Dong等[19]研究表明基于Landsat/PALSAR橡膠林分類精度比MODIS/PLASAR高，橡膠林劃定以Chen等[18]基于Landsat /PALSAR(25 m) 繪制的2010年海南島橡膠林分布圖為準(zhǔn)，利用海南島橡膠林分布圖對NDVI時間序列數(shù)據(jù)進(jìn)行掩模處理，以便精準(zhǔn)提取橡膠林。

2 研究方法

2.1 物候提取方法

為了減少噪聲的影響，利用S-G濾波對NDVI時間序列進(jìn)行去噪平滑處理。S-G濾波[20]是利用滑動窗口對曲線進(jìn)行分段的多項(xiàng)回歸，并加權(quán)平均，擬合公式為：

(1)

式中：Y*為擬合后值；Yi+j為第i+j個像元的原始NDVI值；j為滑動窗口內(nèi)縱向原始數(shù)據(jù)下標(biāo)；Ci為第i個像元的濾波系數(shù)，由最小二乘法確定；n為窗口寬度(n=2m+1)，其中m為半個平滑窗口的寬度。采用S-G濾波對時間序列數(shù)據(jù)擬合時，移動窗口的寬度n和迭代次數(shù)d是應(yīng)用中的關(guān)鍵，n不僅決定了平滑度，也影響了生長季節(jié)開始時跟蹤曲線快速變化的能力[21]；d越高可以去除更多異常值，但是也會使結(jié)果噪聲更多。本研究依據(jù)經(jīng)驗(yàn)[17]和樣區(qū)嘗試，設(shè)定n和d都為3。

考慮到研究區(qū)橡膠林NDVI的變化特征，本研究采用動態(tài)閾值法提取橡膠林的物候參數(shù)。動態(tài)閾值法的提出者J?nsson等[21]建議植被生長季開始日期閾值為20%左右，但是對森林植被生長季開始與結(jié)束期提取時的閾值設(shè)定各有不同。根據(jù)已有經(jīng)驗(yàn)并在多次試驗(yàn)的基礎(chǔ)上，將生長季開始和生長季結(jié)束開始時間閾值分別設(shè)定為30%和60%。

2.2 物候趨勢分析

為研究海南島橡膠林每個像元點(diǎn)15 a間的物候變化趨勢，基于最小二乘法的線性回歸方法逐像元計算橡膠林物候變化率，并利用F檢驗(yàn)分析其顯著性，具體算法詳見文獻(xiàn)[22]，其公式為：

(2)

式中：Gslope為物候變化趨勢;Gi為第i年物候開始(結(jié)束)日期，N為年變化量(N=15)。當(dāng)Gslope>0，物候推遲； 當(dāng)Gslope<0，物候提前。使用F檢驗(yàn)進(jìn)行顯著性分析，當(dāng)P<0.05時，趨勢顯著。

3 結(jié)果與分析

3.1 典型樣區(qū)NDVI值擬合曲線特征

海南島3個典型橡膠林樣區(qū)Google Earth影像紋理如圖2(a)—(c)所示。本研究在典型樣區(qū)NDVI曲線重構(gòu)時，前后各添加1 a的虛擬時間，圖2(d)—(f)擬合曲線時序有69期，只選取中間時序23—46期進(jìn)行分析。儋州市兩院試驗(yàn)農(nóng)場橡膠林展葉期開始日和落葉期開始日閾值分別設(shè)置為30%和50%。圖2(d)可以看出，儋州市橡膠林在12月下旬—次年2月初明顯落葉，2月初NDVI值最低，為落葉盛期； 在2月下旬開始展葉，6月下旬進(jìn)入穩(wěn)定期； 時序38—41間NDVI值出現(xiàn)下降，經(jīng)分析該時段正處于8—9月雨季，云量較多，故NDVI值受影響下降。同理，保亭黎族苗族自治縣熱作所閾值設(shè)定為30%和60%。圖2(e)所示保亭黎族苗族自治縣橡膠林12月下旬開始落葉，2月為落葉盛期，3月下旬開始展葉，6月下旬—7月上旬NDVI值下降，可能是熱作所橡膠林地塊瑣碎(圖2(b))，像元混雜，或天氣等因素造成的NDVI值下降。瓊中黎族苗族自治縣陽江農(nóng)場閾值設(shè)定為30%和60%。圖2(f)顯示瓊中黎族苗族自治縣12月初橡膠林開始落葉，2月初—2月下旬達(dá)落葉盛期，3月中旬為橡膠林展葉期，之后NDVI值迅速上升，在5月下旬達(dá)到穩(wěn)定期。根據(jù)典型樣區(qū)試驗(yàn)，設(shè)定海南島橡膠林展葉期開始日和落葉期開始日閾值分別為30%和60%。

(a) 儋州兩院試驗(yàn)農(nóng)場四隊(duì)影像(b) 保亭熱作所十一隊(duì)影像(c) 瓊中陽江農(nóng)場老田隊(duì)影像

(d) 儋州橡膠林NDVI時間序列曲線(e) 保亭橡膠林NDVI時間序列曲線(f) 瓊中橡膠林NDVI時間序列曲線

圖2 典型樣區(qū)影像及擬合曲線

3.2 橡膠林物候結(jié)果驗(yàn)證

采用2017年儋州市(19°35′N,109°30′E)、白沙黎族自治縣(19°22′N,109°28′E)、瓊中黎族苗族自治縣(19°14′N,109°48′E)橡膠林物候觀測數(shù)據(jù)對研究結(jié)果進(jìn)行驗(yàn)證，數(shù)據(jù)來源于海南儋州熱帶農(nóng)業(yè)生態(tài)系統(tǒng)國家野外科學(xué)觀測研究站。觀測記錄數(shù)據(jù)雖然比影像時間晚，但仍可用，因?yàn)橛^測站點(diǎn)橡膠林齡為15～20 a。本研究選取展葉期開始日、落葉期開始日2種觀測數(shù)據(jù)，因生長季長度觀測站點(diǎn)未記載，根據(jù)觀測點(diǎn)記錄的第一篷葉抽發(fā)期到落葉始期之間的時間計算得出。地面觀測站點(diǎn)數(shù)據(jù)與相同站點(diǎn)15 a的平均物候數(shù)據(jù)對比，對比結(jié)果見表1。3個站點(diǎn)觀測記錄的展葉期開始日和遙感提取的結(jié)果相差不大。觀測數(shù)據(jù)的落葉期開始日均大于遙感手段監(jiān)測天數(shù)，相差6～11 d，原因是每年遙感影像最后一期(第23期)正好對應(yīng)353 d，小于全年實(shí)際天數(shù)365 d，但2種結(jié)果均表現(xiàn)出橡膠林在年末開始落葉。橡膠生長季長度觀測數(shù)據(jù)均大于遙感手段監(jiān)測數(shù)據(jù)，由于人工觀測和遙感監(jiān)測的尺度差異較大，且遙感手段監(jiān)測到的物候期和地面觀測的物候期不具有傳統(tǒng)意義上的一致性，區(qū)域內(nèi)最早展葉的植物出現(xiàn)綠色，并不一定表明植物遙感生長季的開始，因?yàn)檫b感監(jiān)測的生長季開始日期是區(qū)域內(nèi)多數(shù)植物開始展葉的日期，遙感監(jiān)測的物候期相對地面單一植株的物候期有延遲，因此遙感監(jiān)測的橡膠林生長季長度短于地面觀測記錄的生長季長度。2001—2014年海南省儋州市南豐鎮(zhèn)那王村(19°24′20″N,109°33′44″E)站點(diǎn)觀測的橡膠展葉期物候[13]平均在3月6日，本研究春季物候平均出現(xiàn)在3月4日，與14 a平均觀測數(shù)據(jù)相差不大。本文研究結(jié)果與地面觀測物候結(jié)果一致，證明了本研究結(jié)果的可靠性。

表1 不同物候測量方法的結(jié)果對比

3.3 橡膠林物候年際變化

從時間上看，海南島橡膠林展葉期開始日呈提前趨勢(圖3)，從2001年的第70天(3月上旬)，波動提前到2015年的第42天(2月中旬)，約提前28 d，平均每年提前0.94 d(R2=0.16，P>0.05)。落葉期開始日整體呈延遲趨勢，從2001年的第330天(11月下旬)延遲到2015年第356天(12月下旬)，約延遲了26 d，平均每年延遲0.84 d(R2=0.04，P>0.05)。生長季長度總體呈延長趨勢，從2001年的第246 天(8個月)波動延長到2015年的第300 天(10個月)，大致延長了54 d，平均每年延長1.79 d(R2=0.10，P>0.05)。橡膠林展葉期開始日的提前趨勢(0.94 d/a)，與1950—2015年全國木本植物春季物候變化(0.28 d/a)相比提前較多[23]，與1982—2006年東部南北樣帶亞熱帶落葉闊葉林返青期開始日變化(0.58 d/a)相比提前也較多[24]，比2001—2014年海南省儋州市南豐鎮(zhèn)那王村觀測點(diǎn)橡膠林展葉期變化(1.12 d/a)提前較少[13]，但在變化趨勢上具有一致性； 秋季落葉期開始日的延遲趨勢(0.84 d/a)，與1950—2015年全國木本植物秋季物候變化(0.19 d/a)相比延遲較多[23]，同1982—2006年東部南北樣帶亞熱帶落葉闊葉林休眠期起始日變化(0.80 d/a)相差不大[24]； 生長季長度的延長趨勢(1.79 d/a)，較1950—2015年全國木本植被生長季長度變化(0.21 d/a)延長較多[23]，且大于1982—2006年東部南北樣帶亞熱帶落葉闊葉林生長季長度延長趨勢(1.40 d/a)[24]。本研究春季展葉期開始日提前趨勢是秋季落葉期開始日延遲趨勢的1.1倍，推測生長季長度延長是由春季展葉期開始日提前引起的。

圖3 2001—2015年海南島橡膠林物候年際變化

3.4 橡膠林物候空間分布

海南島15 a間展葉期開始日自中向西和東逐漸提前，落葉期開始日則相反(圖4(a)—(b))，研究區(qū)海拔由中間向四周降低，表明物候的開始、結(jié)束日期與海拔有著密切關(guān)系。橡膠林平均展葉期開始日在第17—81天，即1月中旬—3月下旬(圖4(a))，其中第17—81天開始展葉的橡膠林占全島的81.38%，第49—81天占58.11%(圖5(a))，即全島一半左右的橡膠林在3月份開始展葉； 展葉期較早的地區(qū)集中在研究區(qū)西部和東部，均在第1—49天開始展葉。橡膠林平均落葉期開始日在第321—353天，11月中旬—12月下旬(圖4(b))，第321—次年第1天開始落葉的橡膠林有71%(圖5(b))，即絕大多數(shù)橡膠林在11月中旬—12月下旬落葉； 落葉期較早的地區(qū)集中在中部山區(qū)，落葉期較晚的地區(qū)有東方市、瓊海市和萬寧市等。橡膠林生長周期約持續(xù)209～305 d，達(dá)7～10個月左右(圖4(c))，其中生長季長度225～273 d的橡膠林占全島的54.04%(圖5(c))，即海南島大多數(shù)橡膠林生長季長度約在8～9個月； 瓊海市和萬寧市生長周期較長約在9～10個月左右。

(a) 展葉期開始日平均物候(b) 落葉期開始日平均物候(c) 生長季長度平均物候

(a) 展葉期開始日 (b) 落葉期開始日 (c) 生長季長度

3.5 橡膠林物候時空變化趨勢

從空間上看，15 a來，海南島橡膠林展葉期開始日總體上呈提前趨勢(圖6(a))，變化幅度為-2.3～-0.1 d/a，展葉期開始日顯著變化地區(qū)較少(圖6(d))，約占8.7%，變化顯著區(qū)域中約81%的橡膠林顯著提前，主要集中在白沙黎族自治縣、屯昌市、瓊中黎族苗族自治縣、萬寧市和瓊海市等地。橡膠林落葉期開始日呈延遲趨勢(圖6(b))，延遲幅度集中在0.5～2.7 d/a，顯著延遲地區(qū)占33.4%(圖6(e))。海南島橡膠林生長季長度呈增加趨勢(圖6(c))，增長趨勢為0.05～0.9 d/a，生長季長度顯著延長區(qū)域占13%(圖6(f))，主要分布在儋州市和白沙黎族自治縣等。對比橡膠林時空變化趨勢，落葉期開始日顯著延遲是橡膠林物候變化最顯著的特征。

(a) 展葉期開始日年變化趨勢(b) 落葉期開始日年變化趨勢(c) 生長季長度年變化趨勢

(d) 展葉期開始日通過(P<0.05)顯著性水平檢驗(yàn)(e) 落葉期開始日通過(P<0.05)顯著性水平檢驗(yàn)(f) 生長季長度通過(P<0.05)顯著性水平檢驗(yàn)

橡膠林春季物候的提前和秋季物候的延遲，這與近幾十年來基于地面觀測[13]和遙感監(jiān)測[25-26]的研究結(jié)果一致。此外，本研究獲取的橡膠林生長季長度是指從橡膠展葉期開始日到橡膠林落葉期開始日之間的時間累積，實(shí)際生產(chǎn)過程中，橡膠割膠時間是從篷葉發(fā)育完好后開始計算的，故本研究的橡膠林生長季長度較實(shí)際產(chǎn)膠時間偏長。就橡膠生長機(jī)理來看，該結(jié)果提取的生長季長度符合橡膠生長機(jī)理，可以準(zhǔn)確地反映橡膠林的生長季長度。

橡膠林春季物候提前，秋季物候延遲，生長季長度延長，割膠期有可能提前和延長，產(chǎn)膠量可能會相應(yīng)提高。但是，橡膠展葉期提前也會帶來一些危害的提前發(fā)生，例如第一篷葉展葉期的提前橡膠樹極大可能受到寒害，而且這個時期橡膠病蟲害極易發(fā)生，尤其是白粉病，導(dǎo)致防治時間會提前和延長。因此，隨著橡膠物候期的變化，還應(yīng)加強(qiáng)橡膠樹對極端氣象災(zāi)害和病蟲害的預(yù)防措施。

本研究中使用2010年橡膠林空間分布數(shù)據(jù)，并以此分析2001—2015年橡膠林物候時空變化特征。結(jié)合橡膠種植面積統(tǒng)計數(shù)據(jù)來看，2001—2010年和2010—2015年期間橡膠種植面積、種植區(qū)域均發(fā)生了變化。但是根據(jù)Chen等[27]研究表明，近幾十年來新增種植橡膠林主要來源于老膠園更新，其次為耕地和林地轉(zhuǎn)化，這表明橡膠林種植面積雖然有變化，但老膠園種植區(qū)的位置變化不大。從海南省各市縣橡膠熱帶作物面積統(tǒng)計數(shù)據(jù)來看，2001—2010年海南島南部、東南部內(nèi)陸地區(qū)橡膠種植面積增長幅度較小，而北部橡膠種植面積增長幅度較大，主要有儋州市、澄邁縣和臨高縣。因此，本研究結(jié)果可能由于影像時間和土地利用變化之間的差異存在誤差，其中海南島北部橡膠林物候結(jié)果差異較大，東部誤差較小。

4 結(jié)論

本研究基于MODIS NDVI數(shù)據(jù)探討了2001—2015年海南島橡膠林物候的時空變化特征，基于S-G濾波和動態(tài)閾值法提取的物候結(jié)果與地面觀測結(jié)果較一致，說明對海南島橡膠林物候監(jiān)測是可靠的，可為橡膠林生長監(jiān)測和保護(hù)提供一定的參考依據(jù)。

1)本研究設(shè)定展葉期開始日和落葉期開始日閾值分別為30%和60%，可以準(zhǔn)確地監(jiān)測海南島橡膠林物候特征。典型樣區(qū)大約在12月開始落葉，2月進(jìn)入落葉盛期，2月下旬—3月下旬陸續(xù)進(jìn)入展葉期，典型樣區(qū)物候提取方法可為局部地區(qū)物候監(jiān)測提供有效的技術(shù)參考。

2)空間上，2001—2015年海南島橡膠林從1月中旬—3月下旬開始展葉，11月中旬—12月下旬逐漸開始落葉，生長季長度長達(dá)7～10個月。同一時間段內(nèi)，海南島橡膠林物候空間分布與海拔空間分布具有較好的一致性。

3)時間上，15 a間海南島橡膠林物候變化趨勢不顯著。展葉期開始日提前約28 d，落葉期開始日約推遲26 d，生長季長度延長約54 d，展葉期開始日提前是生長季延長的主要原因。

4)物候在空間尺度上變化趨勢存在顯著空間差異。15 a間橡膠林展葉期開始日變化趨勢普遍提前，落葉期開始日變化趨勢整體推遲，生長季長度變化趨勢普遍延長，橡膠林落葉期開始日推遲成為橡膠林較顯著的物候變化特征之一。