以 MODIS 時序資料分析台灣植被物候空間分布

Volume 20, No.1, 2016, pp. 1-14 DOI：10.6574/JPRS.2016.20(1).1

1國立臺灣大學地理環境資源學系 博士後研究員 收到日期：民國 103 年 03 月 07 日

2國立臺灣大學地理環境資源學系 副教授 修改日期：民國 103 年 11 月 03 日

*通訊作者, 電話: 02-33663733, E-mail: choying@ntu.edu.tw 接受日期：民國 103 年 11 月 07 日

以 MODIS 時序資料分析台灣植被物候空間分布

張仲德

黃倬英

摘 要

地表植被物候反應生態系對於氣候變化的作用，本研究利用 2001 至 2010 年間由 MODIS (the Moderate Resolution Imaging Spectroradiometer) 波 段 反 射 率 資 料 演 算 得 到 的 月 份 光 合 作 用 植 物 (PV, photosynthetically active vegetation)資料，透過 TIMESAT 軟體分析台灣地區地表物候空間分布的特性，包 括生長季開始時間(onset)、生長季結束時間(offset)及生長季長度(length of the growing season)。結果顯示 台灣地區不同年間生長季開始時間與生長季結束時間(r = 0.67, p = 0.035)及生長季長度(r = 0.77, p = 0.009)具有顯著關係；其中生長季開始時間的延遲受到聖嬰現象(ENSO, El Niño-Southern Oscillation)影響 的春季乾旱所控制，並對生態系造成負面衝擊。TIMESAT 是相當穩定的物候分析工具，有助於釐清氣候 變化及環境擾動因素對植物生長反應的影響。

關鍵詞：中解析度成像光譜儀、物候分析、氣候變異、聖嬰現象

1. 前言

氣候變遷對於自然生態系的影響已相當顯著，

包括地表植物生長季的提前、生長季時間的延長及 物種分布的變化等物候(phenology)的改變，並進而 影響動植物生長、發展及移動等，已成為全球氣候 變遷相關議題中重要的一 環 (Parmesan & Yohe, 2003; Jeong et al., 2011)。由於陸地上的植被對氣候 改變的反應相當敏感，又與大氣及陸地之間物質、

能量交換息息相關，因此探討物候與氣候之間的關 聯性是了解及預測未來氣候變遷對全球生態系影 響的重要課題(Crabtree et al., 2009)，分析物候的季 節循環與溫度和水分獲得等氣候參數之間的連結，

以及應用不同方法決定物候的各個階段，例如生長 季開始、結束時間及生長季長度等，都受到極大的 關注。傳統物候研究在特定地點記錄特定物種各個 生長階段出現的時空特性，雖然可以提供可靠而直

接的資訊，但需要花費相當多的人力及時間，尤其 對偏遠或不易到達地區的長期監測仍存在相當大 的限制。

近二十多年來隨著衛星及遙測資訊在時間、空 間及光譜解析度的精進，有利研究者對較大尺度的 區域分析地面植物隨著不同季節生長活動所反映 出來的光譜訊號強弱或植生指標變化，評估地面植 被相對生長物候變動的情形，稱之為地表物候(land surface phenology; White & Nemani, 2006; Fisher &

Mustard, 2007)。張仲德等(2011)回顧國內外近三十 年來對植被物候現地觀測及遙測分析的發展及趨 勢，整理的結果中大致顯示寒帶及溫帶地區植被物 候主要受溫度主導，例如由中國大陸近四十年的物 候觀測來看，隨著經緯度的推移，緯度每增加一度，

開花物候延後 2.7 天；經度往東每增加一度，則延 後 0.65 天；海拔高度每增加 100 公尺，延後 1.2 天(Zheng et al., 2002)，溫度是影響物候變化最重要

的因子。而乾燥及半乾燥地區植被物候則受雨量空 間 分 布 的 變 動 所 影 響 (Nicholson et al., 1990;

Nicholson & Farrar, 1994)。相對於現地觀測記錄的 困難及空間代表性，衛星資料的分析將可提供不同 時空尺度下探討植被物候變化的可能性。

經 由 遙 測 資 訊 計 算 的 長 時 間 植 生 指 標 (vegetation index) 如 NDVI (the Normalized Difference Vegetation Index, 常態化差異植生指標) 及 EVI (the Enhanced Vegetation Index, 增揚植生指 標)，已廣泛應用於全球與區域生態環境變化的監 測、植被覆蓋變化、物候特徵萃取及作物生長量評 估 等 (Eklundh et al., 2009; Kariyeva and van Leeuwen, 2011)。植生指標所反映的植被生長時間 序列曲線應該是連續而相對平滑的，但可能受到感 測器本身問題、太陽照射角度、地表雙向反射影響、

大氣雲霧及地表冰雪等干擾因素而產生雜訊，因此 研究也經常使用最大值合成法 (MVC, maximum value composite)來降低干擾的程度(Huete et al., 2002)，但合成時間內仍會存在干擾使植生指標的 時間序列呈現不規則的變動，降低季節變化的趨勢 或者無法明確監測到植被受氣候變異所造成的影 響(Atkinson et al., 2011)。因此目前已有相當多方法 發展出來用以重建植生指標的時間序列曲線，例如 最佳斜率萃取法(BISE, best index slope extraction, Viovy et al., 1992)、門檻值法(threshold, Dall’Olmo

& Kamieli, 2002) 、 移 動 平 均 法 (moving average curve, Reed et al., 1994) 及季節中點法 (seasonal mid-point, White et al., 1997)等，其目的主要是在於 能有效消除干擾的影響，反映植被真實的物候生長 規律。從過去的研究來看，不同方法的適用性與研 究區域內所要分析的植被覆蓋類型有著緊密的關 聯，也都各具有其應用上的優點(Beck et al., 2006;

Hird & McDermid, 2009; Hufkens et al., 2012)。近年 來，TIMESAT 軟體被發展用來評估遙測時間序列 資料，能夠計算研究區域中各個像元的物候指數 (phenological metrics)，如生長季開始及結束時間、

生長季長度、變化幅度和積分等代表植被生長的規 律及變異(Jönsson & Eklundh, 2004)。也被應用到其 他的研究領域中，如環境及物候變遷空間分布的變

化(Heumann et al., 2007)、增進生態系劃分的準確 性(Tottrup et al., 2007)或提供乾燥生態系火災潛勢 的遙測及氣候指標(Verbesselt et al., 2006)等。在許 多 碳 循 環 模 式 及 生 產 力 評 估 的 研 究 中 ， 也 將 TIMESAT 視為長時期遙測資料平滑及填充的重要 分析工具(Olofsson et al., 2008; Schubert et al., 2012)，其中北美碳計劃(North American Carbon Program) 中也利用 TIMESAT 分析 MODIS (the Moderate Resolution Imaging Spectroradiometer)影 像 所 得 到 物 候 指 數 做 為 標 準 的 物 候 產 品 如 MOD09PHN 及 MOD15PHN (Jones et al., 2012)，因 此對於研究應用的進展而言提供相當多的協助。

台灣屬於熱帶及亞熱帶高山島嶼，劇烈的海拔 梯度變化造就氣候及生態多樣性的特色(Ding et al., 2008)，但不同的生態系分布範圍可能相當侷限，

面對全球氣候變遷所造成溫度及雨量的影響，即使 是微小改變對生態系的衝擊也相當大。因此地表植 被對於氣候條件變化的反應相當敏感，同時也可能 間接影響其他的生物活動和生存。台灣早期為了精 進農業栽培技術並掌握農作物生長狀態，曾對稻米 及果樹生長物候進行調查，其中果樹生長之同時線 (isochrone)隨緯度每增加一度或高度增加 100 公尺，

延遲四日；但稻米受人為灌溉影響很大，因此生長 同時線沒有差異(張月娥，1958)。對於森林地區的 物候觀察，通常是在單一試驗地從事一至三年的調 查，包括植物發芽、開花及結實等(林國銓等，1997；

楊金昌等，1999；曾喜育等，2005)，但大部分記 錄的時間長度不足以分析或探討氣候變化所造成 的影響。而衛星資料的相關研究，例如應用遙測影 像分析植被光譜季節變化，主要目的是為了區分不 同森林類型(呂明倫等，2008)，但影像來源如 SPOT (Satellite Pour l'Observation de la Terre)或 Landsat (Land satellite)的再訪週期(revisit period)較長，若受 到天氣因素的干擾並不容易得到高品質影像，對應 用衛星資訊分析植被物候產生相當大的限制。近年 來由於高時間解析度的衛星影像如 MODIS 資料取 得更加便利，包含了不同的空間解析度(250-1000 m)及光譜解析度(36 個波段)有助於不同研究項目 的進行，藉由每天取得且含蓋台灣全島的資訊能更

即時有效地分析地表植被生長活動的變化。Chang et al. (2014a, b)利用近十年 MODIS 資料探討台灣 地區不同植被類型與氣溫和雨量之間的關聯，並對 區域及大尺度氣候因素與植被生長反應有進一步 的釐清。不過相對於許多全球及大陸尺度的物候分 析，目前對熱帶及亞熱帶島嶼生態系的了解仍相當 缺乏。因此本研究將以 2001-2010 年間 MODIS 資 料為基礎，應用 TIMESAT 分析得到台灣地區地表 植被生長季開始時間、結束時間及生長季長度的變 化，以及探討影響物候參數變動的因素，此地表物 候的型態間接提供了植被覆蓋空間分布的特徵，同 時也對於植被生長的反應與年間氣候變異和關聯 提供相當有價值的資訊。

圖 1 台灣地區主要植被類型空間分布

2. 研究方法與材料

2.1 研究地區

台灣島位處歐亞大陸與太平洋交界，面積約 3

萬 6 千平方公里，在東西寛不到 150 公里的距離內，

高度的梯度變化將近 4000 公尺；北回歸線穿越台 灣中部，將台灣南北劃分成熱帶季風氣候及亞熱帶 季風氣候。年平均溫度約 22 度(攝氏溫度)，年平 均雨量約 2500 公厘，北部的月均溫變化(1 月份為 15.0 度，7 月份為 28.7 度)較南部(1 月份為 18.2 度，

7 月份為 28.4 度)更為劇烈。雨量空間分布由平地 往山地增加，東部大於西部，東北部年雨量最高可 達到 6500 公厘，而西南部最低可能低於 1500 公厘 (Chang et al., 2014a)。如此劇烈的環境梯度變化造 就台灣島嶼在氣候及生態多樣性的特性，土地覆蓋 的分布也由平原地區的稻作及旱作，到丘陵山地的 常綠闊葉林、針闊葉混淆林、針葉林及高山草原等，

森林及農業耕作的面積分別占全島的 58%及 27%，

其餘 15%為都市建成地、河流、沙灘等(圖 1, Chang et al., 2014a)，此複雜的植被覆蓋變化不只造就植 物多樣性的特性，同時也對動物的多樣性具有相當 大的貢獻(Ding et al., 2008)。

2.2 研究材料

本研究自美國國家航空暨太空總署 (NASA, National Aeronautics and Space Administration) Land Processes Distributed Active Archive Center (LP DACC, https://lpdaac.usgs.gov/get _data/

data_pool)下載取得 2001-2010 年間 8 天及 500 公尺 空間解析 度的 地表反 射率 波段資料 (MOD09A1;

surface reflectance bands 1-7)。影像經過大氣校正與 雜 訊 去 除 處 理 後 進 行 光 譜 混 合 分 析 (spectral mixture analysis)，此方法經常應用於模擬出影像中 每 一 個 像 素 (pixel) 可 行 光 合 作 用 植 物 (PV, photosynthetically active vegetation)、不可行光合作 用 植 物 (NPV, non-photosynthetically active vegetation) 與 土 壤 裸 露 地 (SRO, soils and rock outcrop)三種端元(endmember)所占的比率(Asner &

Heidebrecht, 2002) 。 此 方 法 假 設 每 一 個 像 素

(

 (  )

組合：

𝜌()_{𝑝𝑖𝑥𝑒𝑙}= ∑[𝐶𝑒𝜌()_𝑒] + 𝜀

= [𝐶𝑃𝑉𝜌()𝑃𝑉+ 𝐶𝑁𝑃𝑉𝜌()𝑁𝑃𝑉+ 𝐶𝑆𝑅𝑂𝜌()𝑆𝑅𝑂] + 𝜀

∑(𝐶

) = 1.0

式(1)中 C_e代表像素中每一個端元所佔有的比例，

而 ()e則表示在波段 時各端元的反射值， 代 表殘差；式(2)表示三個端元的總合等於 1 (Huang et al., 2013)。完成光譜混合分析後，即可抽取出每一 像素中可行光合植物(PV)覆蓋比率。由於 MODIS 的 8 天資料在台灣多雲霧覆蓋的環境仍受到相當 大的影響，因此研究中再進一步將十年間共 460 幅影像以最大值合成法(MVC, Huete et al., 2002)得 到各月份的影像(n = 120)，研究也認為以月份最大 值植生指標與氣候因素之間進行相關性分析時，相 對於其他更短時間尺度的資料能獲得更好的結果 (Suzuki et al., 2006; Prasad et al., 2007)。

另外在 TIMESAT 分析工具中提供研究者以土 地利用資料為基礎，對分析區域定義不同植被類型 空間分布，由分析人員對各主要植被類型分布區域 的時序影像資料重複測試，最後再給予不同植被類 型適當的參數設定。本研究所使用的土地利用分類 資 料 是 由 國 土 測 繪 中 心 取 得 (http://lui.nlsc.gov.tw/LUWeb/)，此資料係 2006-2008 年間利用航空照片、衛星參考資料及現地調查所建 立(Chen & Chang, 2013)。主要土地覆蓋類型中森 林及農業耕作的面積分別占全島的 58%及 27%，

其餘 15%為都市建成地、河流、沙灘等。森林地區 有 53%屬於闊葉林、21%為針葉林、19%為針闊葉 混淆林，其餘 7%為高山箭竹草原。農業耕作的面 積當中有 38%為水稻田、47%為非灌溉旱作，其他 15%為草地、散布的竹林及溼地等(圖 1, Chang et al., 2014a)。

2.3 應用 TIMESAT 進行長期資 料重建及擬合

每一時期遙測影像的植生指數代表該時間植 被生長狀況，經由影像中各個像元時間序列植被指

數的數據分析，即可了解像元內植被季節變化的訊 息(如圖 2(a)中的像元[x, y])。為了去除時間序列上 影 響 植 生 指 標 的 干 擾 ， 本 研 究 採 用 (Jönsson &

Eklundh,2002)所發展之 TIMESAT 軟體(TIMESAT v. 3.0, http:// www .nateko .lu.se /TIMESAT/

timesat.asp)先對時間序列資料進行重建及擬合，再 利用擬合後的長期資料萃取不同年度的物候參數，

其概念是利用動態門檻值方法來進行，也就是當植 生指標的增加達到當年幅度一定百分比的時間點 定義為生長季開始，當植生指標降低到當年幅度的 一定百分比的時間點定義為生長季結束時間，此幅 度通常設定為 10%，此動態門檻值的方法在於進行 多年份資料的物候分析時，不同年間的植被生長變 異不會相互影響，能得到較適合的結果(圖 2(b), Jönsson & Eklundh, 2004)。

分析中採用 Savitzky-Golay (SG)濾波函數對資 料進行平滑擬合，此方法的優點在於能保持植被生 長的動態，並降低大氣效應的影響(Chen et al., 2004)。SG 對時間序列植生指標的重建及濾波過程，

可以下列方程式表示：

𝑌̅ =

𝑁 (3)

其中𝑌̅為濾波處理後新產生之植生指標數值，𝑌_𝑗+1 代表原始植生指標數值，𝐶_𝑖為濾波係數，N 則代表 移動窗口中所包括的數據點 2m+1。應用 SG 濾波 函數來減少資料中的噪音再進行後續的分析，已是 相 當 常 見 及 可 靠 的 方 法 (Brady et al., 2007;

Kariyeva & van Leeuwen, 2011; Chang et al., 2013;

Chen et al., 2014)。

2.4 物候分析及季節參數的萃取

在 TIMESAT 中選定 SG 濾波函數對原始植生 指標 PV 時間序列進行擬合得到物候參數及空間分 布的基礎後，研究者還需要了解土地覆蓋或植被類 型的分布資訊，再進一步從分析區域中選擇代表不 同植被類型的像元，本研究中主要的植被類型包括 針葉林、針闊葉混淆林、闊葉林、農地及建成地(圖

1)，再依照要求分別設置相關的參數。(1) 植生指 標的範圍(data range)，一般的植生指標範圍介於 0 和 1 之間。(2) 噪音去除門檻方法(spike method)主 要在去除噪音所造成的離群值(spike 或 outlier)，工 具中提供兩種方式，第一種為 median filter，當時 間序列中某一個時間點的數值大於前後數值的移 動中位數時則被移除；第二種為 STL replace，主 考慮季節(Season)及趨勢(Trend)的組成，若某一時 間點的數值不適合擬合的型態，則給予較低(Low) 的權重，兩種方式的選擇由研究者視實際資料情況 來決定(Eklundh & Jönsson, 2009)。(3) 季節性參數 (seasonal parameter)的考慮通常當時間序列資料具 有一年單一生長季時設定為 1，若一年具有兩次生 長季，如台灣地區的稻作生長，設定為 0。(4) 迭 代次數(Number of envelope iterations)設定數值介 於 1 至 3，次數越大時表示降低原始曲線中較低數 值如雜訊的影響，使得擬合曲線能逼近原始曲線的 峰值。(5) 適合強度(Adaptation strength)設定數值 介於 0 至 10 之間，其目的同樣使得擬合曲線逼進 原始曲線峰值，一般的設定在 2 或 3 (Eklundh &

Jönsson, 2009)。(6) 移動平均窗口大小(SG window size)數值越高表示資料的時序變動被平滑的程度 也越高。(7) 生長季開始結束參數(season start/stop) 是用來判斷資料曲線由最低至最高生長季開始及 結束的門檻值比率，一般設定在 0.5。經由視窗化 的操作對採樣像元進行測試及擬合，判斷決定適合

該植被類型的各項參數，最後做為研究區域影像中 各個像元物候分析的基礎。利用完成物候分析產生 之結果及資料檔，研究者能從中定義所要萃取的物 候參數、萃取時間段及輸出影像格式，以利進行後 續分析。應用遙測資料分析擷取物候參數與傳統植 物觀察的物候期雖然仍會有些許落差，但能夠反應 較大區域尺度植被生長物候的關鍵訊息和相對狀 況。目前 TIMESAT 中能提供研究者擷取 11 項相 關的物候參數，圖 2(b)中的 a 點表示生長開始時間，

b 點表示生長季結束時間，而生長季長度 g 就由 a 和 b 之間的長度來決定。其中 c 和 d 是由曲線左右 兩測上升至最高值 80%的位置，生長季中點 e 則由 c 及 d 的中間位置來決定；代表植被年生長季變化 的幅度 f 是由 a-e 以 b-e 垂直變化的平均值來表示。

積分 h 為函數擬合後與左右最低值之間的面積，表 示植被季節性的活動；而積分 i 為函數擬合至植生 指標為零之間的面積，表示總植被生產力，各項物 候參數所代表之意義和研究上應用的關聯可參閱 相關說明(Eklundh & Jönsson, 2009)。在本研究中物 候分析所擷取的資訊是一般最為常用的生長季開 始時間(圖 2 中 a 點)、生長季結束時間(圖 2 中 b 點) 以及 a 與 b 之間的距離所決定的生長季長度(圖 2 中 g)，分析台灣地區地表植被生長的空間分布、不 同年間的變化和三項參數之間的關係，並探討影響 物候參數變異的原因。

圖 2 (a) 觀測資料時間序列分布及(b) 植生指數原始數值分布(細黑線)與函數擬合(function fitting, 粗 黑線)後結果及其所得到之物候參數示意圖，其中 a 及 b 分別為生長季開始和結束時間，g 表示 生長季長度。c 和 d 是由曲線左右兩測上升至最高值 80%的位置，生長季中點 e 則由 c 及 d 的中 間位置來決定；代表植被年生長季變化的幅度 f 是由 a-e 以 b-e 垂直變化的平均值來表示。積分 h 為函數擬合後與左右最低值之間的面積，表示植被季節性的活動；而積分 i 為函數擬合至植生 指標為零之間的面積，表示總植被的生產力

表 1 不同植被類型於 TIMESAT 中應用 Savitzky-Golay (SG)濾波函數進行資料擬合及物候分析時主要參數 之設定

地表種類 Spike method

Seasonal parameter

Envelope iteration

Adaptive strength

SG window

Season Start/stop

針葉林 2 1 2 3 5 0.5, 0.5

針闊葉混淆林 2 1 2 3 5 0.5, 0.5

闊葉林 2 1 2 3 5 0.5, 0.5

農地 2 0 2 2 2 0.5, 0.5

建成地 2 1 2 2 3 0.5, 0.5

3. 結果與討論

3.1 PV 時間序列平滑化結果與 分析

以不同植被類型生長季節變化來看，森林(針 葉林、針闊葉混淆林及闊葉林)相當一致，主要是 森林地區生長的季節變化有相似的趨勢，如 PV 月 份最大及最小值通常出現在 7 月(針葉林為 0.79、

針闊葉混淆林為 0.89 及闊葉林為 0.94)及 1 月(針葉 林為 0.58、針闊葉混淆林為 0.59 及闊葉林為 0.65) (圖 3(a)-(c))。台灣地區受人為耕作影響的農田具有 一年兩穫的特性，因此 PV 在季節的變動具有兩次 高峰(分別為 4 月及 9 月)，其 PV 季節變化在所有 地表植被類型當中最大(可達到 0.40)(圖 3(d))。建 成地的月份平均 PV 及季節差異是所有植被類型中 最低的(0.20-0.35) (圖 3(e))。

由擬合所需的參數來看，森林地區(針葉林、

針闊葉混淆林及闊葉林)由於植被生長的變動及季 節之間的差異較一致，因此物候分析中所需要的參 數設定也相同(表 1)。但在非森林地區，相關的參 數設定如適合強度(adaptive strength)及 SG 平均窗 口(SG window size)皆低於森林地區(表 1)，主要是 因為非森林地區 PV 的季節性差異較森林地區要低 或 者 因 為 農 地 一 年 有 兩 次 生 長 季 變 化 所 造 成 (seasonal parameter 為 0；表 1 及圖 3(d))，由於 SG 是利用濾波窗口大小及迭代次數或適合強度等過 濾干擾，較低的次數能獲得稍微平滑的結果，也能 保留異常數值；較高的次數雖可以去除異常值，但

也可能因為過度擬合而喪失植被原有的生長變化 (Chen et al., 2004)。以台灣不同植被類型分析而言，

森林及非森林地區的季節變化都相當明顯，唯有農 田地區因為人為耕作具有一年兩穫(兩次生長季) 的特性需要留意其分析時參數設定的差異，其他相 關的數值設置在不同的植被類型之間經由研究者 進行必要的調整應能得到相當良好的結果。

3.2 台灣地區地表植被覆蓋的空 間分布特性

台灣地區地表植被 PV 的空間變化由西南部的 沿海平原(0.05)往丘陵山地逐漸增加(> 0.80)，至高 海拔再略微降低(< 0.40；圖 4)；此空間上的分布特 性主要反映人為活動及氣候的影響，低海拔平原地 區氣溫及雨量皆適合植被生長，但因為都市化及農 業活動使其年平均 PV 的數值偏低，當雨量隨著海 拔的上升而增加，植被的生長及覆蓋也跟著提高。

更高海拔的地區(> 2700 公尺)由於較低的溫度及植 被類型由森林轉變為高山草原，植被覆蓋或 PV 也 跟著下降(圖 4)。雖然森林及非森林地區受人為活 動影響程度不同，但同一種植被類型的 PV 季節變 化在不同年間型態仍相當一致，其季節生長規律通 常在夏季達到最高而冬季降到最低(圖 3)，此現象 主要還是反應氣候(氣溫和雨量)季節性循環的特 性，如 Chang et al. (2014b)的分析指出不同植被類 型的 PV 與月平均溫度(r² = 0.58-0.83)及月雨量(r² = 0.50-0.78)具有相當良好的關係。而台灣全島的年 平均 PV 的變化由 2002 年的 0.63 增加至 2008 年的 0.68 (圖 4)，年間的變化也受到雨量所影響，其中 2002 年持續的乾旱使其年雨量(1520 公厘)遠低於

台灣長期平均(2500 公厘)，相對地 2008 年的則屬 於溼潤的年份(3000 公厘) (Chen et al., 2009)，這也 顯示台灣熱帶及亞熱帶的生態系對於不同年間氣 候變異或者雨量累積變化的反應相當敏感。

Chang et al. (2014b)應用主成份分析探討台灣 地區影響地表植被生長變異的主要來源，其結果指 出植被類型的空間分布特性及生長規律是地表植 被年間變異的重要組成，可以解釋 96%的生長變異；

同時植被類型各月份主成份負荷量也與月平均溫 度(r = 0.85, p< 0.01)和雨量(r = 0.56, p< 0.01)具有 相當一致的規律變化。就氣候的空間分布特徵來看，

台灣地區即使在高海拔地區不同月份的溫度分布 也甚少低於植物生長所需，但低海拔地區較高的溫 度及較低的雨量則可能造成乾旱並對植物生長形 成極大的限制(Chang et al., 2014a)。近年的研究已 顯示區域尺度的氣候有朝向乾季越加乾旱的趨勢 (Chou et al., 2013)，同時台灣地區在小雨及降雨日 數上逐漸減少，大雨頻率則明顯的增加(Chou et al., 2009)，如此可能對植物生長及發展產生負面效應，

未來對於環境生態上的衝擊仍需更多的評估才能 加以釐清。

圖 3 不同地表植被類型原始 PV 時間序列分布(黑線)及 SG 函數擬合重建結果(紅線)：(a) 針葉林 (Conifer forest)； (b) 針闊葉混淆林(Mixed forest)；(c) 闊葉林(Broadleaf forest)；(d) 農地 (Farmland)；(e) 建成地(Built-up area)

3.3 台灣地區物候分析之空間特 性

台灣地區平均之生長季開始時間(onset, Julian day, JD)分布，約由平原地區的 80 日(3 月底)逐漸 增加至中海拔地區的 160 日(6 月初) (圖 4)；生長季 結束時間的空間分布則由西南部平原的平原地區 的 240 日(8 月底)逐漸增加至台灣西南及東南部中 海拔森林地區至 340 日(12 月初) (圖 4)，平均生長 季的長度由中高海拔森林地區的 210 日減少至西 部平原農田及高海拔的 150 日(圖 4)。三個物候參 數空間分布的特徵主要也是反應植被類型的差異，

農地的生長季開始及結束時間通常較早，因此生長 季長度也較短；而森林地區生長季開始時間稍慢，

但生長季也較晚結束，因此有較長的生長季長度 (表 2)。

若以台灣地區整體平均而言，較晚的生長季開 始時間通常也會導致生長季結束時間延後(r = 0.67, p = 0.035, 圖 5(a))，但會使得生長季時間長度縮減 (r = 0.77, p = 0.009, 圖 5(b))，不過生長季結束時 間與生長季長度之間則沒有顯著的關係存在(r =

0.10, p = 0.784, 圖 5(c))，此意謂台灣地區生長季 開始時間對於地表植被的生長季長度或生長活動 具有關鍵的影響。林國銓等(1997)過去在台灣東北 部福山天然闊葉林的現地物候調查指出，通常 3-4 月為闊葉林抽芽展葉期，10 月至隔年 2 月為落葉 期；曾喜育等(2005)在台灣中部關刀溪森林對台灣 榕進行的展葉及落葉調查，指出台灣榕植株在 2-3 月開始展新葉，每年 10-11 月開始落葉，12 月至隔 年 2 月是葉片量最少的時期，同時展葉及落葉大致 與氣溫及雨量的變化具有關聯。雖然沒有確切日期 可供對比，但以本研究所呈現的空間分布來看，森 林地區生長季開始(3-4 月)及結束(10-11 月)的時間 大致能符合過去東北部及中部現地調查的結果。未 來若能應用品質良好的影像進行局部區域分析，將 有更好的機會連結現地調查及遙測分析。

由 2001-2010 年間物候參數的空間分布來看，

其中生長季開始時間延後(≥ JD 140)所佔的比例在 不同年間具有相當大的差異，由 2001 年的 33%增

加至 2002 年的 38%，再逐漸降低至 2005 年的 13%，

2006 年至 2010 年之間約佔有 25-30%的區域(圖 4)。

先前的研究已指出台灣地區生長季的延後及生產 力的減少和春季(2-3 月)雨量偏少(< 40 公厘)所造 成的乾旱具有顯著的關聯，此春季雨量的多寡與大 尺度大氣循環(ENSO, El Niño-Southern Oscillation) 具有相當緊密的連結(Chang et al., 2013)。氣候分析 也顯示就大尺度氣候循環的條件而言，前一年冬季 (11 至 2 月)若聖嬰現象(El Niño)盛行，則隔年台灣 地區的春季(2-3 月)變得相對溼潤，同時梅雨季(5-6 月)雨量也會增多；相反地，若前一年冬季反聖嬰 盛行(La Niña)，次年春季及梅雨季雨量將相對減少 (Jiang et al., 2003)。許多大陸及全球尺度的研究指 出，聖嬰及反聖嬰所造成年間氣候變異及其對植被 生長物候及生產力具有重要的擾動作用，如東北亞 地區因為聖嬰(El Niño)造成相對溫暖的春季(3-5 月)，使植被生長季開始時間提早(Park et al., 2012)。

南美洲的亞馬遜盆地在聖嬰(El Niño)盛行期間通 常會發生嚴重的乾旱，並使森林的死亡率提高及生 長受到抑制 (Nepstad et al., 2004)。對熱帶及亞熱 帶地區森林生態系而言，聖嬰/反聖嬰盛行期間所 造成的乾旱也可能導致森林結構永久性改變及生 產力的劇烈變化(Brando et al., 2008; Bastos et al., 2013)。目前氣候變遷不僅改變全球降雨強度、降 雨 季 節 分 布 並 增 加 不 同 年 間 的 變 異 (Dai, 2013;

Feng et al., 2013)，同時也將增加聖嬰及反聖嬰現象 發生的頻率及強度，因此更加頻繁持續的乾旱勢必 阻礙溼潤熱帶及亞熱帶森林生長、回復及碳收支平 衡(Marlier et al., 2013; Cai et al., 2014)，社會經濟發 展也不可避免將遭受嚴厲的衝擊。

許多研究利用植生指標衡量植被在遭受氣候 上乾旱的影響程度，主要是以乾旱事件前後的差異 或標準化方式來評估植被綠度的損失(Saleska et al., 2007; Xu et al., 2011)。但植生指標因為飽和問 題、大氣或感測器的噪音，使得綠度減損並不完全 是受到乾旱所造成；或者即使生產力已下降，但綠 度仍然沒有顯著變化而引起爭論(Asner et al., 2004;

Samanta et al., 2010; Atkinson et al., 2011)。應用其 他的技術的協助如多變量分析或物候分析等，可以

相當有效地將影響資料的雜訊去除，擷取關鍵的訊 息並與氣候變異或自然擾動等因素得到良好的連 結(Cong et al., 2013; Griffiths et al., 2013; Ma et al., 2013; Chang et al., 2014b)。本研究利用 TIMESAT 分析台灣地區地表植被物候空間分布的特性及年

間變異，其中生長季開始時間的早晚顯著影響生長 季結束及生長季長度，由於西太平洋及台灣地區春 季降雨量(2-3 月)的多寡與大尺度聖嬰現象的氣候 循環具有關聯，是影響台灣地區的地表植被物候及 生產力的關鍵因素(Chang et al., 2013)。

圖 4 由上至下分別代表 2001 年至 2010 年間年平均 PV、生長季開始時間(onset, JD, Julian date)、生長季 結束時間(offset, JD)及生長季長度(length, days)

圖 5 2001 年至 2010 年間(a) 年平均生長季開始時間(onset)與生長季結束時間(offset)關係、(b) 生長季開始 時間與生長季長度(length of the growing season)關係及(c) 生長季結束時間與生長季長度之關係

表 2 2001-2010 年間不同地表植被類型物候分析結果比較 生長季開始時間

(JD, Julian date)

生長季結束時間 (JD, Julian date)

生長季長度 (Days)

Average ( SD)

針葉林 130.6 ( 37.5) 302.3 ( 60.8) 169.8 ( 39.4) 針闊葉混淆林 122.3 ( 27.1) 305.4 ( 35.1) 183.4 ( 28.1) 闊葉林 123.9 ( 31.1) 309.9 ( 44.9) 184.7 ( 36.1) 農地 119.8 ( 40.0) 287.1 ( 78.1) 148.5 ( 60.9) 建成地 123.0 ( 39.6) 285.9 ( 77.4) 165.1 ( 54.4)

4. 結論

本研究取得 2001-2010 年間 MODIS 地表反射 率波段資料並合成月份最大值 PV 進行物候分析，

利用 TIMESAT 軟體分析得到台灣地區地表物候參 數，包括生長季開始時間、生長季結束時間及生長 季長度。結果顯示台灣地區生長季開始時間顯著影 響生長季結束時間及生長季長度，生長季開始時間 受到大尺度 ENSO 導致的春季雨量所控制，對台 灣生態系具有關鍵性的影響。在氣候變遷的情勢下，

極端氣候及大尺度氣候循環事件的頻率及強度將 更加提高，因此利用物候分析工具有助於釐清地表 植被生長反應與氣候變化之間的關聯。

致謝

本 論 文 的 完 成 主 要 感 謝 科 技 部 (NSC 98-2221-E-002-198,98-2313-B-002-062-MY2, 100-2621-B-002-001-MY3)及臺灣大學相關研究計 畫(EcoNTU:NTU-CESRP -102R7604-2)的經費補助 及支持。

參考文獻

呂明倫，葉慶龍，鍾玉龍，謝依達，2008。墾丁國 家公園四種優勢植群之地面光譜分析，航測及 遙測學刊，13(1)19-28。

林國銓，黃吳清標，劉哲政，1997。福山試驗林天 然 闊 葉 林 之 物 候 現 象 ， 台 灣 林 業 科 學 ， 12(3):347-353。

張月娥，1958。臺灣物候報告，氣象學報，4(4):6-9 頁。

張仲德，王素芬，林登秋，2011。氣候變遷與不同 尺 度 植 被 物 候 研 究 之 回 顧 ， 地 理 學 報 ， 63:1-33。

曾喜育，曾麗蓉，歐辰雄，呂福源，2005。關刀溪 森林生態系之臺灣榕物候，中華林學季刊，

38(4):377-395。

楊金昌，王亞男，姜家華，賴玉芳，1999。塔塔加 地區臺灣雲杉、臺灣鐵杉及玉山箭竹物候學之 初步研究，31(3):251-263。

Asner, G., and Heidebrecht, K.B., 2002. Spectral unmixing of vegetation, soil, and dry carbon cover in arid regions: comparing multispectral and hyperspectral observations. Int J Remote Sens, (23):3939-3958.

Asner, G., Nepstad, D., Cardinot, G., and Ray, D., 2004. Drought stress and carbon uptake in an Amazon forest measured with spaceborne imaging spectroscopy. Proc Natl Acad U S A, (101): 6039-6044.

Atkinson, P.M., Dash, J., and Jeganathan C., 2011.

Amazon vegetation greenness as measured by satellite sensors over the last decade. Geophy Res Lett, (38): L19105.

Bastos, A., Running, S.W., Gouveia, C., and Trigo, R.M., 2013. The global NPP dependence on ENSO: La Niña and the extraordinary year of 2011. J Geophys Res: Biogeosci, (118):

1247-1255.

Beck, P.S.A., Atzberger, C., Høgda, K.A., Johansen, B., and Skinmore, A.K., 2006. Improved monitoring of vegetation dynamics at very high latitudes: a new method using MODIS NDVI.

Remote Sens Environ, (100):321-334.

Brady, B.A., Joshi, R.W., Hermance, J.F., and Mustard, J.F., 2007. A curve fitting procedure to derive inter-annual phenologies from time series of noisy satellite NDVI data. Remote Sens Environ, (106): 137-145.

Brando, P.M., Nepstad, D.C., Davidson, E.A.,

Trumbore, S.E., Ray, D., and Camargo, P., 2008.

Drought effects on litterfall, wood production and belowground carbon cycling in an Amazon forest: results of a throughfall reduction experiment. Phil Trans R Soc B, 363, 1839-1848.

Cai, W., Borlace, S., Lengaigne, M., van Rensch, P., Collins, M., Vecchi, G., Timmermann, A., Santoso, A., McPhaden, M.J., Wu, L., England, M.H., Wang, G., Guilyardi, E., and Jin, F.F., 2014. Increasing frequency of extreme El Niño events due to greenhouse warming. Nature Clim Change, (4):111-116.

Chang, C.T., Wang, H.C., and Huang, C.Y., 2013.

Impacts of vegetation onset time on the net primary productivity in a mountainous island in Pacific Asia. Environ Res Lett, (8): 045030.

Chang, C.T., Wang, S.F., Vadeboncoeur, M.A., and Lin, T.C., 2014a. Relating vegetation dynamics to temperature and precipitation at monthly and annual timescales in Taiwan using MODIS vegetation indices. Int J Remote Sens, (35):

598-620.

Chang, C.T., Wang, H.C., and Huang, C.Y., 2014b.

Retrieving multi-scale climatic variations from high dimensional time-series MODIS green vegetation cover in a tropical/subtropical mountainous island. J Mt Sci, (11):407-420.

Chen, B., Zhang, X., Tao, J., Wu, J., Wang, J, Shi, P., Zhang, Y., and Yu, C., 2014. The impact of climate change and anthropogenic activities on alpine grassland over the Qinghai-Tibet Plateau.

Agr For Meteorol, 189-190, 11-18.

Chen, J., Jönsson, P., Tamura, M., Gu, Z., Matsuhita, B., and Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Goley filter.

Remote Sens Environ, (91):332-344.

Chen, S.T., Kuo, C.C., and Yu, P.S., 2009. Historical trends and variability of meteorological droughts in Taiwan. Hydrol Sci J, (54): 430-441.

Chen, X.R., and Chang, H.S., 2013. Discussion on the impact factor of agricultural land in Taiwan:

case study of Yunlin, Chiayi and Tainan.

Procedia Environ Sci, (17):28-36.

Chou, C., Huang, F., Tseng, L., Tu, J., and Tan, P.H., 2009. Annual cycle of rainfall in the Western North Pacific and East Asia Sector. J Clim, (22):2073-2094.

Chou, C., Chiang, J.C.H., Chung, C.H., Liao, Y.C., and Lee, C.J., 2013 Increase in the range between wet and dry season precipitation. Nature Geosci, (6):263-267.

Cong, N., Wang, T., Nan, H., Ma, Y.C., Wang, X., Myneni, R.B., and Piao, S., 2013. Changes in satellite-derived spring vegetation green-up date

and its linkage to climate in China from 1982-2010: a multimethod analysis. Glob Change Biol, (19):881-891.

Crabtree, R., Potter, C., Mullen, R., Sheldon, J., Huang, S., Harmsen, J., Rodman, A., and Jean, C., 2009 A modelling and spatio-temporal analysis framework for monitoring environmental change using NPP as an ecosystem indicator. Remote Sens Environ, (113):1486-1496.

Dai, A., 2013: Increasing drought under global warming in observations and models. Nature Clim Change, (3):52-58.

Dall’Olmo, G., and Kamieli, A., 2002. Monitoring phenological cycles of desert ecosystems using NDVI and LST data derived from NOAA-AVHRR imagery. Int J Remote Sens, (23)4055-4071.

Ding, T.S., Su, G., and Yuan, H.W., 2008. Breeding bird community composition in different successional vegetation in the montane coniferous forests zone of Taiwan. For Ecol Manage, (255):2038-2048.

Eklundh, L., and Jönsson, P., 2009. Timesat 3.0 software manual. Lund Univeristy, Sweden.

Eklundh, L., Johansson, T., and Solberg, S., 2009:

Mapping insect defoliation in Scots pine with MODIS time-series data. Remote Sens Environ, (113):1566-1573.

Feng, X., Porporato, A., and Rodriguez-Iturbe, I., 2013. Changes in rainfall seasonality in the tropics. Nature Clim Change, (3):811-815.

Fisher, J.I., and Mustard, J.F., 2007. Cross-scalar satellite phenology from ground, Landsat, and MODIS data. Remote Sens Environ, (109):261-273.

Griffiths, P., Kuemmerle, T., Baumann, M., Radeloff, V.C., Abrudan, I.V., Lieskovsky, J., Munteanu, C., Ostapowicz, K., and Hostert, P., 2013. Forest disturbances, forest recovery, and changes in forest types across the Carpathian ecoregion from 1985 to 2010 based on Landsat image composites. Remote Sens Environ, (151):72-88.

Heumann, B.W., Seaquist, J.W., Eklundh, L., and Jönsson, P., 2007. AVHRR derived phenological change in the Sahel and Soudan, Africa, 1982-2005. Remote Sens Environ, (108):385-392.

Hird, J.N., and McDermid, G.J., 2009. Noise reduction of NDVI time series: an empirical comparison of selected techniques. Remote Sens Environ, (113):248-258.

Huang, C., Chai, C., Chang, C., Huang, J., Hu, K., Lu, M., and Chung, Y., 2013. An integrated optical remote sensing system for environmental

perturbation research. IEEE J Sel Top Appl Earth Obs Remote Sens, (6): 2434-2444.

Huete, A.R., Didan, K., Miura, K., Rodriguez, E.P., Gao, X., and Ferreira, L.G., 2002. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens Environ, (83):195-213.

Hufkens, K., Friedl, M., Snonnenting, O., Braswell, B.H., Milliman, T., and Richardson, A.D., 2012.

Linking near-surface and satellite remote sensing measurements of deciduous broadleaf forest phenology. Remote Sens Environ, (117):

307-321.

Jeong, S.J., Ho, C.H., Gim, H.J., and Brown, M.E., 2011. Phenology shifts at start versus end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982-2008.

Glob Change Biol, (17): 2385-2399.

Jiang, Z.H., Chen, G.T.J., and Wu, M.C., 2003.

Large-scale circulation patterns associated with heavy spring rain events over Taiwan in strong ENSO and non-ENSO years. Mon Weather Rev, (131):1769-1782.

Jones, M.O., Kimball, J.S., Jones, L.A., and McDonald, K.C., 2012. Satellite passive microwave detection of North America start of season. Remote Sens Environ, (123):324-333.

Jönsson, P., and Eklundh, L., 2002. Seasonality extraction by function fitting to time-series of satellite sensor data. IEEE Trans Geosci Remote Sens, (40):1824-1832.

Jönsson, P., and Eklundh, L., 2004. TIMESAT- a program for analyzing time-series of satellite sensor data. Comput Geosci, (30): 833-845.

Kariyeva, J., and van Leeuwen, W.J., 2011.

Environment drivers of NDVI-based vegetation phenology in central Asia. Remote Sens, (3):203-246.

Ma, X., Huete, A., Yu, Q., Coupe, N.R., Davies, K., Broich, M., Ratana, P., Beringer, J., Hutley, L.B., Cleverly, J., Boulain, N., and Eamus, D., 2013.

Spatial patterns and temporal dynamics in savanna vegetation phenology across the North Australian tropical transect. Remote Sens Environ, (139): 97-115.

Marlier, M.E., DeFries, R.S., Voulgarakis, A., Kinney, P.L., Randerson, J.T., Shindell, D.T., Chen, Y., and Faluvegi, G., 2013. El Niño and health risks from landscape fire emissions in southeast Asia.

Nature Clim Change, (3):131-136.

Nepstad, D., Lefebvre, P., Lopes, U., Tomasella, J., Schlesinger, P., Solórzano, L., Moutinho, P., Ray, D., and Benito, J.G., 2004. Amazon drought and its implications for forest flammability and tree-growth: a basin-wide analysis. Glob Change

Biol, (10):704-717.

Nicholson, S.E., Davenport, M.L., and Malo, A.R., 1990. A comparison of the vegetation response to rainfall in the Sahel and East Africa, using normalized difference vegetation index from NOAA AVHRR. Clim Change, (17):209-241.

Nicholson, S.E., and Farrar, T.J., 1994. The influence of soil type on the relationships between NDVI, rainfall and soil moisture in semiarid Botswana, I.

NDVI response to rainfall. Remote Sen Environ, (50)107-120.

Olofsson, P., Lagergren, F., Lindroth, A., Lindström, J., Klemedtsson, L., Kutsch, W., and Eklundh, L., 2008. Towards operational remote sensing of forest carbon balance across Northern Europe.

Biogeosci, (5): 817-832.

Park, K.A., Bayarsaikhan, U., and Kim, K.R., 2012.

Effects of El Niño on spring phenology of the highest mountain in north-east Asia. Int J Remote Sens, (33): 5268-5288.

Parmesan, C., and Yohe, G., 2003. A globally coherent fingerprint of climate change impacts across natural systems. Nature, (421): 37-42.

Prasad, V.K., Badarinath, V.K.S., and Eaturu, A., 2007. Spatial patterns of vegetation phenology metrics and related climatic controls of eight contrasting forest types in India: analysis from remote sensing datasets. Theor Appl Climatol, (89):95-107.

Reed, G.I., Brown, J.F., VanderZee, D., Loveland, T.R., Merchant, J.W., and Ohlen, D.O., 1994.

Measuring phenological variability from satellite imagery. J Veg Sci, (5):703-714.

Saleska, S.R., Didan, K., Huete, A.R., and da Rocha, H.R., 2007. Amazon forests green-up during 2005 drought. Science, (318):612.

Samanta, A., Ganguly, S., Hashimoto, H., Devadiga, S., Vermote, Y., Knyazikhin, Y., Nemani, R.R., and Myneni, R.B., 2010. Amazon forests did not green-up during the 2005 drought. Geophys Res Lett, (37): L05401.

Schubert, P., Lagergren, F., Aurela, M., Christensen, T., Grelle, A., Heliasz, M., Klemedtsson, L., Lindroth, A., Pilegaard, K., Vesala, T., and Eklundh, L., 2012. Modeling GPP in the Nordic forest landscape with MODIS time series data- comparison with the MODIS GPP product.

Remote Sens Environ, (126):136-147.

Suzuki, R., Xu, J.Q., and Motoya, K., 2006. Global analyses of satellite-derived vegetation index related to climatological wetness and warmth. Int J Climatol, (26):425-438.

Tottrup, C., Rasmussen, M.S., Eklundh, L., and Jönsson, P., 2007. Mapping fractional forest cover across the highlands of mainland Southeast

Asia using MODIS data and regression tree modelling. Int J Remote Sens, (28):23-46.

Verbesselt, J., Jönsson, P., Lhermitte, S., van Aardt, J., and Coppin, P., 2006. Evaluating satellite and climate data-derived indices as fire risk indicators in Savanna ecosystems. IEEE Trans Geosci Remote Sens, (44):1622-1632.

Viovy, N., Arino, O., and Belward, A.S., 1992. The best index slope extraction (BISE): a method for reducing noise in NDVI time-series. Int J Remote Sens, (13):1585-1590.

White, M.A., Thornton, P.E., and Running, S.W., 1997. A continental phenology model for monitoring vegetation responses to interannual climatic variability. Glob Biogeochem Cycle, (11):217-234.

White, M.A., and Nemani, R.R., 2006. Real-time monitoring and short-term forecasting of land surface phenology. Remote Sens Environ, (104):43-49.

Xu, L., Samanta, A., Costa, M.H., Ganguly, S., Nemani, R.R., and Myneni, R.B., 2011.

Widespread decline in greenness of Amazonian vegetation due to the 2010 drought. Geophys Res Lett, (38):L07402.

Zheng, J., Ge, Q., and Hao, Z., 2002. Impacts of climate warming on plants phenophases in China for the last 40 years. Chin Sci Bull, (47):1826-1831.

1 Postdoctoral Research Fellow, Department of Geography, Received Date: Mar. 07, 2014 National Taiwan University. Revised Date: Nov. 03, 2014

2 Associate Professor, Department of Geography, Accepted Date: Nov. 07, 2014 National Taiwan University.

*Correspondent Author, E-mail: choying@ntu.edu.tw

Spatial patterns of vegetation phenology based on MODIS time-series data in Taiwan applying TIMESAT

Chung-Te Chang¹ Cho-ying Huang^2*

ABSTRACT

Land surface phenology reflects the responses of a terrestrial ecosystem to climate change. In this study, we used monthly PV (photosynthetically active vegetation) derived from MODIS (the Moderate Resolution Imaging Spectroradiometer) reflectance data between 2001 and 2010 to delineate the land surface phenology.

The phenological metrics including onset, offset, and the length of the growing season were obtained from the phenological analysis toolbox TIMESAT. Results indicated that the interannual variability of onset time was significantly correlated to offset time (r = 0.67, p = 0.035) and length of the growing season (r = 0.77, p = 0.009). In addition, the shortage of spring rainfall was influenced by ENSO events and resulted in significant delayed on vegetation onset time and cause dramatic impacts on the ecosystem. TIMESAT is a robust phenological analysis tool, which can facilitate our understanding of the vegetation responses and development associated with climate change and environmental disturbances.

Keywords: MODIS (the Moderate Resolution Imaging Spectroradiometer), phenology analysis, climatic anomaly, ENSO