台東太麻里溪集水區地景變遷之研究

全文

(1)國立臺灣師範大學地理學系第二十二屆博士論文 Department of Geography College of Liberal Arts National Taiwan Normal University Doctoral Dissertation. 台東太麻里溪集水區地景變遷之研究 Analysis of Landscape Changes in the Taimali Watershed of Eastern Taiwan. 指導教授：廖學誠 Advisor: Shyue-Cherng Liaw, Ph.D.. 研究生：葉春國 Chun-Kuo Yeh. 中華民國 103 年 6 月.

(2) 中文摘要土地覆蓋變遷分析對於提供集水區土地經營、監測及規劃等資訊是極為重要的。本研究呈現 2005～2011 年間東台灣太麻里集水區大規模地景變遷與其頻繁自然干擾間的關係。為了解土地覆蓋變遷與自然干擾間之關係，我們採取複合式的地景分析途徑，應用了地景指標、馬可夫鏈模式、改變軌跡及邏輯斯迴歸模式。組成性地景指標分析結果顯示崩塌地面積擴張 7.5 倍（15.52 km2），而森林面積縮減 10.6%（20.90 km2）。空間型態地景指標分析結果顯示森林區塊面積變得較小、形狀不規則及空間上破碎；然而，崩塌地區塊面積變得較大及空間上較為聚集。因此，研究區土地覆蓋變遷主要發生於森林覆蓋的減少和破碎，以及崩塌地數目和面積的增加。在研究期間，地景中最顯著的面積改變是林地轉變為崩塌地和河道，分別約為 15.82 km2 和 6.73 km2。就轉移機率而言，人為區塊轉變為林地和河道的機率為最大。透過馬可夫鏈的分析，本研究產生了三種各類土地覆蓋比例未來（2017～2035）的預測狀況。三種預測狀況分別被視為對於下游居民生命及產財具有高、中、低的風險性。 FL、FC 及 VR 三個改變軌跡共占整個集水區全部改變面積的 75.65%，他們被視為主要的改變趨勢，而且代表太麻里集水區的整體地景（OL）變遷趨勢。改變分析結果顯示地景變遷易發生於畢祿山層地質、鄰近斷層和河道、坡度大於 35°及東向坡的條件下。就模式驗證而言，FC 模式具有最高的 AUC 值。然而， RCI 指標值顯示變遷機率的預測皆能高度符合實際的改變軌跡。因此，本研究得出的四個邏輯斯迴歸模式有助於地景變遷的預測。 . 關鍵詞：自然干擾、地景指標、馬可夫鏈、改變軌跡、邏輯斯迴歸。. i .

(3) . ABSTRACT Analysis of land cover changes is fundamental for providing information in relation to watershed land management, monitoring, and planning. This study reveals large-scale land cover transformation in relation to frequent natural disturbances within the Taimali watershed in eastern Taiwan during 2005–2011. To understand land cover changes in relation to natural disturbances, a combined landscape analysis approach using landscape metrics, the Markov chain model, change trajectories, and the logistic regression model is used. Results of composition metrics analysis show that areas of landslides within the region had expanded 7.5 times (by 15.52 km2), but areas of forest had shrunk by 10.6% (20.90 km2). Spatial configuration metrics analysis indicates that patches of forest are becoming smaller, more irregular, and more spatially fragmented, but that landslide patches are expanding and becoming more spatially aggregated. It is considered that land cover changes in the area have occurred mainly through loss and fragmentation of forest cover, and an increase in the number and area of landslides. During 2005–2011, the most noticeable area changes related to transitions from forest cover to landslides (15.82 km2) and channels (6.73 km2). Results show the transition probabilities of human-made patches changing into forest cover and channel corridors are the greatest. Through Markov chain analysis, three future (2017–2035) projections of the proportions of each land cover type are produced. These three predictive statuses are regarded as posing a high, moderate, and low risk, respectively, to life and property downstream. Three change trajectories, FL, FC, and VR, covering 75.65% of the entire changed area, are considered main changes in trend and represent overall landscape (OL) changes in the Taimali watershed. Results of change analysis indicate that the. ii .

(4) . occurrence of landscape change is subject to the geologic condition of Pilushan Formation, the vicinities of faults and rivers, a gradient greater than 35°, and the eastward slope. As for model validation, the FC model has the highest AUC value. RCI values suggest that the prediction of change probabilities could correspond with the actual change trajectories very well. Therefore, the four logistic regression models are useful tools in the landscape change prediction.. Keywords: natural disturbance; landscape metrics; Markov chain; change trajectory; logistic regression.. iii .

(5) . CONTENTS. 1. INTRODUCTION ------------------------------------------------------------------------- 1 1.1 Background ----------------------------------------------------------------------------- 1 1.2 Objectives ------------------------------------------------------------------------------- 3 1.3 Organization ---------------------------------------------------------------------------- 3. 2. LITERATURE REVIEW ----------------------------------------------------------------- 6 2.1 Landscape ecology and the applications of landscape metric-------------------- 6 2.2 Applications of Markov chains ------------------------------------------------------ 7 2.3 Land cover change trajectory -------------------------------------------------------- 9 2.4 Logistic regression ------------------------------------------------------------------- 11. 3. MATERIALS AND METHODS ------------------------------------------------------- 13 3.1 Study area and disturbance events ------------------------------------------------- 13 3.2 Landscape image and classification ----------------------------------------------- 16 3.3 Landscape metrics ------------------------------------------------------------------- 17 3.4 Markov chains ------------------------------------------------------------------------ 20 3.5 Method of trajectory calculation --------------------------------------------------- 21 3.6 Environmental variables in relation to change trajectories --------------------- 22 3.6.1 Environmental variables of the study area --------------------------------- 22 3.6.2 Change analysis in association with environmental variables ----------- 29 3.7 Logistic regression model ----------------------------------------------------------- 30 3.7.1 Mathematics of logistic regression ------------------------------------------ 30 3.7.2 Data sampling and multicollinearity diagnostics -------------------------- 31 3.7.3 Analysis of model fit ---------------------------------------------------------- 32 3.7.4 Model Validation -------------------------------------------------------------- 34 iv . .

(6) . 4. RESULTS --------------------------------------------------------------------------------- 36 4.1 Landscape metrics analysis --------------------------------------------------------- 36 4.1.1 Landscape composition changes -------------------------------------------- 36 4.1.2 Spatial configuration changes ----------------------------------------------- 39 4.1.3 Summary ----------------------------------------------------------------------- 42 4.2 Markov chains analysis -------------------------------------------------------------- 43 4.2.1 Area change matrix ----------------------------------------------------------- 43 4.2.2 Transition probability matrix ------------------------------------------------ 44 4.2.3 Predictions of Markov chain ------------------------------------------------- 46 4.2.4 Summary ----------------------------------------------------------------------- 48 4.3 Analysis of land cover change trajectories --------------------------------------- 50 4.3.1 Calculation and classification of land cover change trajectories -------- 50 4.3.2 Change analysis in connection with environmental variables ----------- 54 4.3.3 Summary ----------------------------------------------------------------------- 64 4.4 Logistic regression analysis --------------------------------------------------------- 66 4.4.1 Multicollinearity diagnostics ------------------------------------------------ 66 4.4.2 Model results ------------------------------------------------------------------- 71 4.4.3 Model validation --------------------------------------------------------------- 81 4.4.4 Summary ----------------------------------------------------------------------- 90. 5. DISCUSSION ----------------------------------------------------------------------------- 92 5.1 Assumptions and predictions of Markov chain ------------------------------ 92 5.2 Selection of independent variables for logistic models --------------------- 93. 6. CONCLUSIONS ------------------------------------------------------------------------- 95 REFERENCES --------------------------------------------------------------------------- 99. v .

(7) . List of Figures Figure. Page. 1.1 Flowchart of the study process---------------------------------------------------------- 5 3.1 Location of the study area with typhoon tracks ------------------------------------- 14 3.2 Lithology map --------------------------------------------------------------------------- 26 3.3 Distance to faults map ------------------------------------------------------------------ 26 3.4 Average annual rainfall map ----------------------------------------------------------- 26 3.5 Distance to rivers map ------------------------------------------------------------------ 27 3.6 Elevation map --------------------------------------------------------------------------- 27 3.7 Slope map -------------------------------------------------------------------------------- 27 3.8 Aspect map ------------------------------------------------------------------------------- 28 3.9 Curvature map --------------------------------------------------------------------------- 28 4.1 Classified land cover maps of the study area from 2005 to 2011 ----------------- 38 4.2 Changes in landscape metrics from 2005 to 2011 ---------------------------------- 41 4.3 Prediction of land cover changes during 2017–2035 (%) ------------------------- 48 4.4 Total trajectories of land cover in the Taimali watershed ------------------------- 53 4.5 Main land cover change trajectories in the Taimali watershed ------------------- 53 4.6 Binary map of land cover trajectories ------------------------------------------------ 54 4.7 Map of RCI and percentage of the entire area of change for lithology----------- 55 4.8 Map of RCI and percentage of the entire area of change for distance to faults - 56 4.9 Map of RCI and percentage of the entire area of change for rainfall ------------ 57 4.10 Map of RCI and percentage of the entire area of change for distance to rivers -------------------------------------------------------------------------------------------- 59 4.11 Scatter plot between distance to rivers and RCI value ---------------------------- 59 4.12 Map of RCI and percentage of the entire area of change for elevation --------- 60. vi .

(8) . 4.13 Map of RCI and percentage of the entire area of change for slope ------------- 61 4.14 Map of RCI and percentage of the entire area of change for aspect ------------ 63 4.15 Map of RCI and percentage of the entire area of change for curvature -------- 64 4.16 Map of probability of change for OL change trajectories ------------------------ 73 4.17 Map of probability of change for FL change trajectories ------------------------ 76 4.18 Map of probability of change for FC change trajectories ------------------------ 78 4.19 Map of probability of change for VR change trajectories ------------------------ 81 4.20 The ROC curve for change trajectories of OL ------------------------------------- 82 4.21 Map of RCI and percentage of the entire area of change for OL change trajectories ----------------------------------------------------------------------------- 83 4.22 The ROC curve for change trajectories of FL ------------------------------------- 84 4.23 Map of RCI and percentage of the entire area of change for FL change trajectories ----------------------------------------------------------------------------- 85 4.24 The ROC curve for change trajectories of FC ------------------------------------- 86 4.25 Map of RCI and percentage of the entire area of change for FC change trajectories ----------------------------------------------------------------------------- 87 4.26 The ROC curve for change trajectories of VR ------------------------------------- 88 4.27 Map of RCI and percentage of the entire area of change for VR change trajectories ----------------------------------------------------------------------------- 89. vii .

(9) . List of Tables Table. Page. 3.1 Earthquake occurrence in the vicinity of the study area --------------------------- 14 3.2 The results of accuracy assessment (%) ---------------------------------------------- 17 3.3 Class-level landscape metrics used in this study ------------------------------------ 19 3.4 Independent variables used in the logistic regression model --------------------- 25 4.1 The area distribution of landscape elements from 2005 to 2011 ------------------ 37 4.2 Matrices of area change (km2) and transition probability for 2005–2008 ------- 45 4.3 Matrices of area change (km2) and transition probability for 2008–2011 ------- 45 4.4 Matrices of area change (km2) and transition probability for 2005–2011 ------- 45 4.5 The top 7 trajectories in the Taimali watershed ------------------------------------- 52 4.6 The percentages of land cover change trajectories --------------------------------- 52 4.7 The area distribution for lithological classes ---------------------------------------- 55 4.8 The area distribution for classes of distance to faults ------------------------------ 56 4.9 The area distribution for rainfall classes --------------------------------------------- 57 4.10 The area distribution for classes of distance to rivers ----------------------------- 58 4.11 Linear correlation analysis between RCI and distance to rivers ----------------- 59 4.12 The area distribution for elevation classes ----------------------------------------- 60 4.13 The area distribution for slope classes ---------------------------------------------- 61 4.14 The area distribution for slope classes ---------------------------------------------- 62 4.15 The area distribution for curvature classes ----------------------------------------- 64 4.16 Statistics of multicollinearity diagnostics of all independent variables for change trajectories of OL ----------------------------------------------------------- 67 4.17 Statistics of multicollinearity diagnostics of seven independent variables for change trajectories of OL ----------------------------------------------------------- 67. viii .

(10) . 4.18 Statistics of multicollinearity diagnostics of all independent variables for change trajectories of FL ----------------------------------------------------------- 68 4.19 Statistics of multicollinearity diagnostics of seven independent variables for change trajectories of FL ----------------------------------------------------------- 68 4.20 Statistics of multicollinearity diagnostics of all independent variables for change trajectories of FC ----------------------------------------------------------- 69 4.21 Statistics of multicollinearity diagnostics of seven independent variables for change trajectories of FC ----------------------------------------------------------- 69 4.22 Statistics of multicollinearity diagnostics of all independent variables for change trajectories of VR ----------------------------------------------------------- 70 4.23 Statistics of multicollinearity diagnostics of seven independent variables for change trajectories of VR ----------------------------------------------------------- 70 4.24 Summary of model fit with seven independent variables for overall landscape change --------------------------------------------------------------------------------- 72 4.25 Classification table based on the logistic regression model for overall landscape change --------------------------------------------------------------------------------- 72 4.26 Regression coefficients and Wald statistics for overall landscape change ----- 73 4.27 Summary of model fit with six independent variables for forest-landslide change --------------------------------------------------------------------------------- 74 4.28 Classification table based on the logistic regression model for forest-landslide change --------------------------------------------------------------------------------- 75 4.29 Regression coefficients and Wald statistics for forest-landslide change ------- 75 4.30 Summary of model fit with five independent variables for forest-channel change --------------------------------------------------------------------------------- 77 4.31 Classification table based on the logistic regression model for forest-channel change --------------------------------------------------------------------------------- 77 ix .

(11) . 4.32 Regression coefficients and Wald statistics for forest-channel change --------- 78 4.33 Summary of model fit with five independent variables for vegetation recovery ------------------------------------------------------------------------------------------- 79 4.34 Classification table based on the logistic regression model for vegetation recovery -------------------------------------------------------------------------------- 80. 4.35 Regression coefficients and Wald statistics for vegetation recovery ------------- 80 4.36 Characteristics of the five classes of probability of change for OL change trajectories ----------------------------------------------------------------------------- 83 4.37 Characteristics of the five classes of probability of change for FL change trajectories ----------------------------------------------------------------------------- 85 4.38 Characteristics of the five classes of probability of change for FC change trajectories ----------------------------------------------------------------------------- 87 4.39 Characteristics of the five classes of probability of change for VR change trajectories ----------------------------------------------------------------------------- 89. x .

(12) 1. INTRODUCTION. 1.1 Background Landscape changes are acknowledged as being a part of global environmental change. However, whether such changes are natural or manmade, they are a subject of concern. Landscapes in Taiwan are susceptible to natural disturbances in relation to the country’s unique global geographical position. Taiwan is situated on the plate boundary between the Philippine Sea plate and the Eurasian plate, where two arc–trench systems have been interacting since the Miocene (Wu et al. 2006b). Large earthquakes are therefore a major hazard, and their strong motion results in landslides, rock falls, ground displacement, and damages to buildings in the vicinity of the epicenter (Wu et al. 2006b; Wu et al. 2006c). Large earthquakes occur more often in eastern Taiwan in relation to the active arc–continent collision that has a convergence rate of 7–8 cm/year (Wu et al. 2006c). Taiwan sits on the border between the Tropics and the Subtropics, and as such lies in the way of typhoons as they track their way across the Northwestern Pacific Ocean. Typhoons occur in the summer and autumn every year, bringing extreme rainfall that triggers debris flow, landslides, and floods, and which poses a substantial threat to homes, infrastructure, and life (Chang et al. 2011). Typhoons often strike Taiwan along the eastern coast, and consequently the area is more prone to natural disturbances than other regions. The Taimali watershed is representative of an area in eastern Taiwan that is susceptible to natural disturbances; and hence is considered to be suitable for an investigation and analysis which will enhance our understanding of land cover change dynamics under the influence of frequent natural disturbances. Studies of land cover. 1 .

(13) changes can contribute to monitoring changes on temporal and spatial scales, and can be used to assess the status of natural resources and predict future changes (Munsi et al. 2010). In addition, land managers can benefit from land cover change mapping, as this provides useful information in relation to watershed land management, monitoring, and planning. In Taiwan, previous studies of land cover changes under the effects of natural disturbances have focused on central Taiwan, and in particular on the Chenyulan watershed, in relation to the 1999 Mw 7.6 Chi-Chi earthquake and several subsequent large typhoons (Chu et al. 2009; Hong et al. 2010; Lin et al. 2006; Lin et al. 2009; Lin et al. 2010). These studies detected land cover changes using landscape metrics, spatial sampling, simulation, and spatial analysis. However, only limited research has been performed in eastern Taiwan, despite the occurrence of frequent natural disturbances. In addition, few studies have attempted to use the combined landscape analysis approach by applying landscape metrics, Markov chains, change trajectories, and logistic regression to delineate land cover changes, and have only computed the transition matrix, but not predicted future changes (Malaviya et al. 2010; Munsi et al. 2010). To assess land cover changes in the Taimali watershed, this research therefore applies combined landscape analysis and uses three FORMOSAT-2 satellite images, from 2005, 2008, and 2011. The present study provides a new dimension in the understanding of land cover changes within a region experiencing frequent natural disturbances.. 2 .

(14) 1.2 Objectives This study provides a comprehensive analysis and projection of the dynamics of landscape change induced by the natural disturbances in the Taimali watershed of eastern Taiwan from 2005 to 2011. The landscape change predominantly occurred in the loss and fragmentation of forest cover and the increment in the number and area of landslides. We focus on three land cover types (landslide, forest, and channel), since they provide important information useful to disaster risk assessment at the Taimali watershed. The primary contributions of this paper are to: (1) detect land cover changes that occurred on a spatio-temporal scale during 2005–2011. (2) predict the future (2017–2035) possible distribution of land cover according to the use of three different disturbance regimes. (3) discern the main land cover change trajectories and connect them with environmental variables. (4) model the probabilities of occurrence of change using logistic regression for four land cover change trajectories.. 1.3 Organization This dissertation is composed of six chapters. The first chapter provides an introduction to the background, objectives, and organization. A review of the literature is presented in detail in Chapter 2 and the review topics contain the theory and application of landscape ecology, Markov chains, land cover change trajectory, and logistic regression. The uses of the materials and methods are demonstrated in Chapter 3. The materials include three landscape images from 2005, 2008, and 2011, three large typhoon events between 2005 and 2011, and historic records of earthquake from 2001 to 2011. The methods employed in the present study contain the landscape 3 .

(15) metrics analysis, the Markov chains model, the change trajectories analysis, and the logistic regression model. The results are shown in Chapter 4. The topics in this section include: (1) changes in the composition and spatial configuration of land cover; (2) matrices of the area change and transition probability, and predictions of the Markov chains; (3) detection of major change trajectories and their relationship with environmental variables; and (4) modelling the probabilities of occurrence of change trajectories and validation of logistic models. Finally, discussion of the results are presented and some conclusions are drawn in Chapter 5 and 6, respectively. The flowchart of the study process is illustrated in Fig. 1.1. First, three FORMOSAT-2 satellite images are used to perform image classification. Through accuracy assessment and combination, four landscape elements are detected. Land cover data plus other eight environmental variables are applied to carry out four major research topics using landscape metrics, Markov chains, pixel history, and logistic regression, respectively. All these processes are to achieve the assessment of landscape change in the Taimali watershed of eastern Taiwan.. 4 .

(16) Fig. 1.1 Flowchart of the study process. 5 .

(17) 2. LITERATURE REVIEW To enhance our understanding about earthquakes and typhoons as a disturbance process on landscape pattern, we take a combined landscape analysis approach that applies landscape metrics, the Markov chain model, change trajectories, and the logistic regression model. These techniques are utilized to assess the effects of typhoons on the characteristics of land cover change, to predict possible trends for the future, and to map probabilities of change for the main change trajectories.. 2.1 Landscape ecology and the applications of landscape metrics Landscape ecology is simply the ecology of landscapes and a landscape is described as a kilometer-wide mosaic where the mix of local ecosystems or land uses is repeated in a similar form (Forman 1995). Nevertheless, when the emphasis is placed on spatial heterogeneity, the concept of landscape as a spatial mosaic can be at various spatial scales from the beetle’s lifeworld within a few square meters to nitrogen dynamics in a watershed (Turner et al. 2001). A landscape mosaic is composed of three types of landscape elements that are relatively homogeneous units recognized at the landscape scale and refer to patches, corridors, and a background matrix (Forman 1995). The subjects of study in landscape ecology are the three landscape characteristics, inclusive of structure (form), function, and change (Forman 1995; Zonneveld 1990). Structure refers to the numbers, sizes, shapes, kinds, and spatial pattern or arrangement of landscape elements; function refers to the flows of energy, materials, and organisms through the structure; and change refers to the dynamics or alteration in the structure and function over time (Dramstad et al. 1996; Turner 1989). Therefore, developing quantitative methods are necessary to describe. 6 .

(18) the landscape structure, compare different landscapes, identify significant changes over time, and relate landscape patterns to ecological function (Turner 1989). Landscape metrics, derived from information theory, describe and quantify landscape patterns that are linked to ecological processes and can predict them on a broad scale (Turner, 1989; Turner et al. 2001). Many studies have demonstrated the importance of the application of landscape metrics in quantifying landscape patterns, their changes over time, and their effects on ecological processes (Antwi et al. 2008; Frondoni et al. 2011; Kintz et al. 2006; Paudel and Yuan, 2012). These metrics of landscape status can be applied to correlate with specific aspects of ecosystem function; therefore changes in metrics are indicators of changes in the ecological condition of the landscape (O’Neill et al. 1997). Besides, on account of the interplay between spatial patterns and ecological processes, landscape metrics can inform environmental planners about landscape functions that are usually difficult or impossible to measure directly (Botequilha Leitão et al. 2006). Olsen et al. (2007) further indicated that the use of landscape metrics is important in building an effective environmental monitoring system.. 2.2 Applications of Markov chains A Markov chain is a stochastic process, and represents a system of elements undergoing transitions from one state to another over time. It can be parameterized by empirically estimating transition probabilities between discrete states in the observed system (Balzter, 2000). Tattoni et al. (2011) indicated that Markov chains were useful mathematical tools for modeling a process with a finite number of states and known probabilities of moving from each state to the others. Reddy et al. (2008) also suggested that the Markov chain analysis was a convenient tool for modeling land use change because it was difficult to describe changes and processes in the landscape. 7 .

(19) Apart from quantifying land use changes, the Markov model procedure can reveal obscure data trends (Muller and Middleton 1994). Therefore, the Markov chain analysis can describe land cover changes from one time point to another, and it further employs this information as the basis to predict future changes. Performing the Markov chain analysis first requires the computation of the transition probabilities that are reported in a matrix showing the rates of transition between states (Tattoni et al. 2011). Markov chains have been widely applied to model changes in land use and land cover from nature-dominated to human-dominated landscapes on a variety of spatial scales. With respect to natural landscapes, research has mostly been related to deforestation, vegetation recovery, forest succession, etc. Through calculating transition probability, Reddy et al. (2009) found that extensive deforestation of Nawarangpur in India resulted from forest converting into agriculture land. In the study of landslide areas in central Taiwan, the Markov chain model was applied for vegetation recovery assessment and simulating future land cover changes was based on the transitive probability matrix (Chuang et al. 2011). In addition, from the examples of Erzgebirge in Germany, the results showed that Markov chain models were an efficient tool to organize the existing knowledge of forest succession into a system of quantitative predictions (Benabdellah et al. 2003). In relation to human-modified landscapes, research has mostly focused on urban expansion. In a case study of Beijing which is a mega city of China, Wu et al. (2006a) found that urban that growth engulfed the surrounding agriculture lands in suburbs and they projected land use changes for the next 20 years with the application of Markov chains and regression analyses. The same examples that croplands were replaced with urban development could be found in the studies of Muller and Middleton (1994) and Weng (2002). 8 .

(20) The disturbance regime of a landscape refers to the spatio-temporal dynamics of disturbances over a prolonged period of time, using characteristics such as spatial distribution, frequency, size, and intensity. It remains a current topic of active research that aims to improve our quantitative understanding and ability to predict the effects of changing disturbance regimes on landscape patterns (Turner et al. 2001). In the Taimali watershed, each typhoon and earthquake event causes a varying degree of change to the landscape, due to the differing size, intensity, and frequency of each event. This research uses transition probabilities calculated from different time periods as a proxy for variable disturbance regimes, with the aim of predicting future land cover changes under varying statuses of disturbance regimes.. 2.3 Land cover change trajectory From a methodological point of view, there are four different ways to quantify land cover changes, including: (1) simple measures of change, like the proportion of change, (2) annual rates of change, (3) changes in the spatial arrangement or structure of land cover, and (4) spatially-explicit representations of change at the pixel level using the concept of pixel histories (Mena 2008). The fourth approach described above is namely the change trajectory that tracks pixel histories across the time series (Rindfuss et al. 2004). Land cover change trajectories refer to a succession of land cover types for a given sampling unit in the raster images of time series over more than two observation years (Mertens and Lambin 2000; Zhou et al. 2008a; Zhou et al. 2008b). Land cover change can be combined to create a high diversity of change trajectories that are determined by the local conditions, regional contexts, government policies, and external influences (Mertens and Lambin 2000; Verburg et al. 2010). By coupling the change history with a specific region, trajectories not only allow the illustration of change patterns in temporal dimension, but also consider their 9 .

(21) continuity and direction (Ruiz and Domon 2009). The notion of trajectory, targeting at identifying the paths by which landscapes change over time, has recently drawn a lot of attention in studies of landscape change analysis (Ruiz and Domon 2009). By means of classification of land cover change trajectories, researchers could find the main driving forces of the landscape changes, which were dominated by natural forces or human activities. A small case study in the Northern Intensive Study Area of Northern Ecuadorian Amazon showed that agricultural trajectories which included pasture and large and small scale agriculture were the dominant change trajectories causing the loss of forest (Mena 2008). Zhou et al. (2008a and 2008b) demonstrated a research on land cover change trajectories in the fragile ecosystem of China’s arid environment and the results indicated that natural forces such as floods were dominant in environmental change; however, human impact, such as cultivation and dam/reservoir construction, has increased dramatically in recent years. In a study of highly dynamic landscape in southern Chile, Carmona and Nahuelhual (2012) found that the trajectories of deforestation and forest degradation, which converted from old growth forest to scrubland and secondary forest, respectively, were the two most important groups of change trajectories. Obtaining comprehensive knowledge of the dynamics of land cover changes could be beneficial to reconstruct past land cover changes, discern their main driving forces, predict future changes, and thus help develop sustainable management practices pointing at preserving fundamental landscape functions (Domon and Bouchard 2007; Hietel et al. 2004). This study tries to show that the analysis of land cover change trajectories can be spatially explicit in response to the dynamics of natural disturbances. In addition, for a case study reconstituting the landscape of Godmanchester (Quebec, Canada), Domon and Bouchard (2007) emphasized that the research on landscape dynamics must not only allow the recording of the 10 .

(22) transformations and their consequences, but also attempt to explain them by identifying the relationships between the biophysical, anthropic and technological factors. Since the study area is a mountainous watershed and the area of human activity is less than 1.84% of the entire watershed, this paper also attempts to reveal the relationships between change trajectories and physical environmental conditions.. 2.4 Logistic regression One of our research objectives is to develop a spatially-explicit model for predicting the probability of land cover change in the Taimali watershed based on geological, meteorological, topographical, and distance-related factors selected according to their importance to the landscape change and data availability. Land cover trajectories can be divided into two categories: stable versus dynamic trajectories (Mena 2008). Stable trajectories are characterized by pixel histories that are unchanging across the time-series while dynamic ones are described by pixel histories with varied classes over time. For each pixel’s trajectory type, it is either stable or dynamic; therefore the dependent variable is a binary variable. Due to the characteristics of the dependent variable, logistic regression was found to be the most suitable method for analyzing the probability of land cover change and hence employed in the present study. Logistic regression is a method for estimating the probability of occurrence of an event from dichotomous data, using a set of variables which may be continuous, discrete, or a combination of both (Bui et al. 2011). The logistic regression has two main advantages. One is that independent variables do not necessarily have normal distributions; the other is that the predicted values are expressed as probabilities lying in the range between 0 and 1 giving the risk of event occurrence in an area (Kleinbaum and Klein 2010). Many researches simultaneously applied the concept of pixel history and the 11 .

(23) method of logistic regression analysis to estimate the probability of land cover change for every pixel while they often focused on landslide hazard assessment. Taking a case study in Central Taiwan for instance, the logistic regression results revealed that the frequency of landslide occurrence played an important role in assessing landslide hazards (Lin et al. 2010). Besides, analyzing historical and recent landslides in the northeastern Ohio, U.S.A., Nandi and Shakoor (2009) suggested that the multivariate logistic regression model, compared with bivariate statistical analysis, was the better model in predicting landslide susceptibility. In fact, relatively few studies were found to assess the probability of change on every pixel for land cover change trajectories (Braimoh and Vlek 2005; Mena 2008). An example from Northern Ecuadorian Amazon (NEA), Mena (2008) mapped the probability of landscape transitions using logistic regression analysis. The results demonstrated that the landscape of the NEA was more dynamic at its core rather than in the periphery. Therefore, binary logistic regression was utilized to model the probability of change, i.e., risk of occurrence of change trajectories, for a given set of independent variables.. 12 .

(24) 3. MATERIALS AND METHODS. 3.1 Study area and disturbance events The study area used in this research is the Taimali watershed, a region susceptible to earthquake and typhoon occurrence in eastern Taiwan (Fig. 3.1). Located in Taitung County, it is a typical mountainous watershed. The Taimali stream moves in a direction from west to east before finally flowing into the Pacific Ocean. The watershed rises to a height of 3000 m, with an average elevation of 935.4 m, (almost 35% of the entire watershed area is over 1000 m). It measures about 27 km in length, is west–east oriented, and stretches from the Taimali stream estuary on the eastern boundary to Mount Dawu on the western border. This watershed also has a total area of 211.5 km2 and an average gradient of 57.6%. Based on historic records from 2001–2011 provided by the Central Weather Bureau (CWB) of Taiwan, three earthquakes with Moment magnitudes (Mw) greater than 6 (and seismic intensities estimated as 5) have occurred in the vicinity of the Taimali watershed (Table 3.1). The Chengkung earthquake occurred in eastern Taiwan in 2003, causing landslides, rock falls, and construction damage in the epicentral area (Wu et al. 2006c). The 2006 Taitung earthquake produced significant ground displacement, and a large numbers of aftershocks occurred throughout the following month (Wu et al. 2006b). The 2010 Jiashian earthquake in SW Taiwan caused moderate damage, but no surface rupture was observed, implying a deep source, which is a relatively rare occurrence in western Taiwan (Huang et al. 2013). However, although no deaths were reported, 96 people were injured, and two high speed trains were disrupted at the time of quake due to ground shaking (Rau et al. 2012).. 13 .

(25) Fig. 3.1 Location of the study area with typhoon tracks. Table 3.1 Earthquake occurrence in the vicinity of the study area Name. Chengkung earthquake Taitung earthquake Jiashian earthquake. Date. GMT+08:00. Latitude. Longitude. Largest intensity. Depth Mw. at Taimali. (Taiwan). (°N). (°E). (km). 10/12/2003. 12:38:13.5. 23.07. 121.40. 17.7. 6.8. 5. 01/04/2006. 18:02:19.5. 22.88. 121.08. 7.2. 6.1. 5. 04/03/2010. 08:18:52.1. 22.97. 120.71. 22.6. 6.3. 5. Moment magnitude (Mw) based on Wu et al. (2006b), Wu et al. (2006c), and Huang et al. (2013); all other information obtained from the CWB. 14 .

(26) Typhoons usually strike Taiwan along the eastern coast, and as the mouth of the Taimali watershed faces east and is adjacent to the Pacific Ocean the area is susceptible to typhoons. Information related to typhoon occurrence in eastern Taiwan was obtained from the typhoon database provided by the CWB of Taiwan. During 2005–2011, a total of 23 typhoons hit Taiwan, and of these three fierce typhoons struck the Taimali watershed (Fig. 3.1). Firstly, the center of Typhoon Haitang swept across northern Taiwan from the east towards the west on July 18, 2005, with a maximum wind speed of 55.0 m/s and a radius of 280 km (CWB 2005). This typhoon caused heavy rainfall and subsequent significant flooding, which destroyed infrastructure and houses in the Taimali watershed. Typhoon Morakot then occurred on August 7–8, 2009; its center moved from east to northwest with a maximum wind speed of 40.0 m/s and a radius of 250 km (CWB 2009). Typhoon Morakot brought record-breaking torrential rainfall that triggered serious flooding, debris flow, and landslides in eastern, central, and southern parts of Taiwan; causing a tremendous amount of damage to property and infrastructure in the Taimali watershed. On September 19–20, 2010, the center of Typhoon Fanapi then passed through central Taiwan from east to west, with a maximum wind speed of 45.0 m/s and a radius of 200 km (CWB 2010). This typhoon produced heavy rainfall in the eastern and central regions of Taiwan, and the flooding that occurred in the Taimali stream flushed away several buildings and approximately 120 m of rail track. The typhoons described here caused both the destruction and creation of land cover, thereby resulting in land cover change.. 15 .

(27) 3.2 Landscape image and classification To obtain land cover data, three FORMOSAT-2 satellite images from 2005, 2008, and 2011, with spatial resolutions of 8 m and minimal cloud cover, were acquired from the Spatial Information Research Center (SIRC) of the National Taiwan University. At our request, three images underwent preprocessing of geo-registration, rectification, and projection into the Transverse Mercator coordinating system by the SIRC. ERDAS IMAGINE version 9.2 software (Leica Geosystems, Inc.) was employed to perform the image classification. Field investigation and land cover maps of the Taimali watershed produced by previous studies were used as reference information to successfully accomplish the land cover classification and the subsequent accuracy assessment. Six types of land cover were identified (including landslide, forest, farmland, buildup, channel sand, and water) using the supervised classification method and the parametric rule of maximum likelihood. In the present study, “forest” is defined as a mixture of grassland and woodland. Each land cover type was recognized using several training samples, and the number of training samples was based on the size and degree of the spatial dispersion of a specific class. Table 3.2 demonstrates the accuracy assessment of land cover classification evaluated from 256 randomly generated reference points obtained from using stratified random sampling (Xiao and Ji 2007). All the overall accuracies and Kappa coefficients are shown to be greater than 83.90% and 0.80, respectively. Since “farmland” and “buildup” were seen to cover a smaller area, both of these land cover types were combined and then categorized as “human-made patches” to simplify the information. The land cover types, “water” and “channel sand” were also combined as “channel corridors” for the same purpose. Therefore, four landscape elements: landslide patches, forest matrix, human-made patches, and channel corridors were used for the analysis of landscape changes. In addition, to explicitly 16 .

(28) detect changes in channel width, channel corridors were digitized using ArcGIS version 9.3.1 in three satellite images. To retain the connectivity of channel corridors, all classified images originally had a pixel size of 8 m, but were converted to a pixel size of 4 m to facilitate further analysis. The resulting three sets of land cover data were then used to calculate landscape metrics in 2005, 2008, and 2011, respectively, along with their transition probabilities.. Table 3.2 The results of accuracy assessment (%) 2005. 2008. 2011. 77.08/92.50. 82.22/90.24. 68.97/80.00. 98.21/100.00. 92.98/98.15. 96.45/95.32. Farmland. 92.86/95.12. 83.72/90.00. 76.47/92.86. Buildup. 93.94/77.50. 89.47/42.50. 100.00/54.55. Channel sand. 81.58/77.50. 68.63/85.37. 74.07/80.00. Water. 97.44/95.00. 90.24/92.50. 87.50/70.00. Overall accuracy. 90.23. 83.98. 89.45. Kappa coefficient. 0.882. 0.807. 0.801. Landslide Forest. 77.08/92.50 indicates Producer’s accuracy/User’s accuracy.. 3.3 Landscape metrics A large number of landscape metrics or indices have been developed and used to characterize a wide variety of spatial patterns (Gustafson 1998; O’ Neill et al. 1997). Spatial patterns have an influence on a wide variety of ecological processes, which 17 .

(29) then in turn influence the spatial patterns (Botequilha Leitão et al. 2006). Landscape metrics are able to measure the composition and spatial configuration of landscape elements; in this context, composition refers to features associated with the variety and abundance of patch types within the landscape, without reference to spatial attributes, and spatial configuration refers to the spatial character and arrangement, position, or orientation of patches within the class or landscape (McGarigal et al. 2012). In our research, we utilized nine class-level metrics calculated by FRAGSTATS version 4.0, including Class Area (CA), Percentage of Landscape (PLAND), Mean Patch Area (AREA_MN), Edge Density (ED), Mean Shape Index (SHAPE_MN), Mean Patch Fractal Dimension (FRAC_MN), Number of Patches (NP), Mean Euclidean Nearest-Neighbor Distance (ENN_MN), and Aggregation Index (AI) (Table 3.3). The first two metrics belong to composition metrics and the others belong to spatial configuration metrics. We calculated class-level metrics because each land cover type provides different information in relation to watershed management at Taimali. These metrics were selected because they were able to assess the variation in aspects of the landscape’s characteristics. They measure the aspects of the patches, which includes proportional abundance, area, edge, shape complexity, subdivision, isolation, and aggregation. Subdivision refers to the degree to which a patch type is broken up into separate patches. Isolation deals with the degree to which patches are spatially isolated from each other, and aggregation refers to the tendency of patch types to be spatially aggregated (and it is an umbrella term used to describe several related concepts, including subdivision and isolation) (McGarigal et al. 2012).. 18 .

(30) Table 3.3 Class-level landscape metrics used in this study Metrics. Abbreviation. Unit. Description. Class Area. CA. ha. Measuring how much of the landscape is comprised of a particular patch type. Percentage of. PLAND. %. Landscape Mean Patch Size. Quantifying the proportional abundance of each patch type in the landscape. AREA_MN. ha. Measuring landscape configuration by the average condition of patch size. Edge Density. ED. m/ha. Standardizing edge to a per unit area basis that facilitates comparisons among landscapes of varying size. Mean Shape Index. SHAPE_MN. none. Measuring the complexity of patch shape compared to a standard shape (square) of the same size. Mean Patch. FRAC_MN. none. Fractal Dimension Number of. dimension of the object NP. none. Patches Mean Euclidean. Characterizing patch shapes via the fractal. Measuring the extent of subdivision or fragmentation of the patch type. ENN_MN. m. Nearest-Neighbor. Measuring patch isolation by using simple Euclidean geometry. Distance Aggregation Index AI. %. Measuring aggregation based on the ratio of the observed number of like adjacencies to the maximum possible number of like adjacencies. The computation and description of metrics are based on FRAGSTATS 4.0 (McGarigal et al. 2012).. 19 .

(31) 3.4 Markov chains A Markov chain, defined as a stochastic process, shows that the value of a process at time t,. , depends only on its value at time. values. ,. ,…,. 1,. , and not on the sequence of. that the process passed through before arriving at. (Weng 2002). It can be expressed as follows (Wu et al. 2006a):. |. ,. ,...,. (3.1). |. |. The conditional probabilities probabilities of order. from state. , to state. are called transition. for all indices 0 ≤ s < t, with. 1 ≤ i, j ≤ k, and are denoted as the transition probability matrix P (Balzter 2000; Logofet and Lesnaya 2000). For K states, the first-order Markov transition matrix P takes the form:. (3.2). P. The future distribution of states distribution of states. can be calculated by multiplying the initial. with the state transition probability matrix P having. elements pij that represent transition probabilities from state. to state , and can be. expressed as (Baasch et al. 2010):. (3.3). In practical application of Markov chain analysis, researchers are initially required 20 .

(32) to establish a transition probability matrix of land cover changes from one time point to another. The matrix shows the character of change at each land cover type, and can then serve as a basis for projecting the proportion of land cover at any considered future time period. In the present study, future projections (2017–2035) of land cover changes under different disturbance regimes were performed using the Markov chain.. 3.5 Method of trajectory calculation Our objective is to acquire a categorical map that demonstrated the pixel history, namely land cover trajectories, at the pixel level. The categorical map can be called the trajectory layer. Each land cover class is described as a code in the land cover raster layer. In the present study, codes in the land cover raster layer from 1 to 4 represent the landslide patch, forest matrix, human-made patch, and channel corridor, respectively. To make the codes in the land cover layers of year 2005, 2008, and 2011 discriminative to each other, the three land cover layers are multiplied by 100, 10, and 1, respectively, and then they are summed to obtain a trajectory layer. The trajectory layer shows a sequence of codes for every pixel, which can be called trajectory codes. When there are less than ten classes of land cover, the trajectory codes can be calculated by the following formula (Mena 2008; Wang et al. 2012; Wang et al. 2013):. TC. C. 10. C. 10. C. 10. (3.4). where TCij is the trajectory code of pixel at row i and column j in the trajectory layer; n is the quantity of the time points; C. , C. , and C. are the codes in land. cover layers of different time points at the given pixel, and they have to be ranked by time points in ascending order. 21 .

(33) 3.6 Environmental variables in relation to change trajectories To realize the relationship between the changed area and the environmental variables, we first convert the map of the main land cover change trajectories into a binary map composed of the changed area and unchanged area. Second, we need to prepare the distribution maps of related environmental variables, and then overlay the binary map with them to calculate the percentage and intensity of land cover change in each subregion of the distribution map of environmental variables.. 3.6.1 Environmental variables of the study area The environmental variables are listed in Table 3.4 and are described respectively as follows: (1) Lithology Lithology is considered an underlying driving factor in the research exploring the relationship between landslides and environmental factors (Bai et al. 2010; Lin et al. 2010; Bui et al. 2011). Since landslide occurrences are the most significant event of land cover change in the study area, the lithology variable is included in the analysis. Geological formations of the study area are composed of three strata, inclusive of Pilushan Formation, Lushan Formation, and Tananao Schist (Fig. 3.2). Lithological characteristics are identified according to Ho (1986). Pilushan Formation is an Eocene system and mainly consists of slate and phyllite. Lushan Formation, a Miocene system, primarily contains charcoal graywacke interbeded with black-charcoal argillite, slate, and phyllite. Tananao Schist, a Pre-Tertiary metamorphic complex, is made up of various schists and marble, both of which are metamorphosed from sedimentary and volcanic rocks.. 22 .

(34) (2) Distance to faults There are two faults in the study area and their spatial distribution is similar to the letter Y shape (Fig. 3.3). One passes across the central part of Taimali watershed in a north-south direction and the other is elongated in a north-east to south-west direction at the northern part of Taimali watershed. The distance to faults is calculated by measuring the straight-line distance to faults with ArcGIS Spatial Analyst tools and is classified into five classes, including 0-2000 m, 2000-4000 m, 4000-6000 m, 6000-8000 m, and ＞8000 m (Fig. 3.3).. (3) Rainfall Since there is no rainfall station in the study area, kriging interpolation is applied to estimate the average annual rainfall of the study area. Average annual rainfall map is produced by utilizing the ordinary Kriging method with the spherical semivariogram model to interpolate a surface from the neighboring twelve rainfall stations. Average annual rainfall is divided into four classes, inclusive of 2000-3000 mm, 3000-4000 mm, 4000-5000 mm, and 5000-6000 mm (Fig. 3.4).. (4) Distance to rivers Distance to rivers is calculated by measuring the straight-line distance to Taimali stream in 2005 with ArcGIS Spatial Analyst tools and is classified into five classes, including 0-300 m, 300-600 m, 600-900 m, 900-1200 m, and ＞1200 m (Fig. 3.5).. (5) Elevation Elevation data are taken from a digital terrain model (DTM) provided by the Aerial Survey Office, Forestry Bureau, Council of Agriculture, R.O.C. The resolution of DTM is 40 × 40 m. Elevation is categorized into six classes, inclusive of 0-500 m, 23 .

(35) 500-1000 m, 1000-1500 m, 1500-2000 m, 2000-2500 m, and ＞2500 m (Fig. 3.6).. (6) Slope The slope gradient is extracted from the DTM and then divided into six slope gradient categories, including 0°-5°, 5°-15°, 15°-25°, 25°-35°, 35°-45°, and ＞45° (Fig. 3.7).. (7) Aspect The slope aspect is also obtained from the DTM and the resulting aspect image is classified and labeled with 9 identifiers based on Miliaresis (2008) (Fig. 3.8). In the distribution map of aspect, azimuth N stands for aspects in the range 337.5° to 22.5°, NE for 22.5° to 67.5°, E for 67.5° to 112.5°, SE for 112.5° to 157.5°, S for 157.5° to 202.5°, SW for 202.5 ° to 247.5°, W for 247.5° to 292.5°, NW for 292.5° to 337.5°, and F for flat terrain that is an undefined aspect.. (8) Curvature Curvature is mathematically defined as the change in slope angle along a very small arc of the curve, and can be expressed as the inverse of the radius of a circle which is tangent over a small arc of the curve (Ohlmacher 2007). A positive curvature indicates that the surface is a convex landform at that cell (e.g., hill and ridge), a negative curvature indicates that the surface is concave landform at that cell (e.g., depression and valley), and a curvature of zero indicates that the surface is flat (Florinsky 2012). The distribution map of the curvature is presented in Fig 3.9.. 24 .

(36) Table 3.4 Independent variables used in the logistic regression model Independent. Type. Unit. Range. Source. variables Central geological Lithology. Categorical –. –. survey, MOEA, R.O.C. Central geological. Distance to Continuous Meters. 0–14295.73. faults. survey, MOEA, R.O.C. Water Resources. Rainfall. Continuous Millimeters. 2457.08–5641.97. Agency, MOEA, R.O.C.. Distance to. Digital Terrain Continuous Meters. 0–2863.56. rivers. Model Digital Terrain. Elevation. Continuous Meters. 3–3090. Model Digital Terrain. Slope. Continuous Degrees. 0–82.60. Model Digital Terrain. Aspect. Categorical Degrees. 0-360 Model Digital Terrain. Curvature. Continuous 1/hectometer -105.94–56.75. 25 . Model.

(37) Fig. 3.2 Lithology map. Fig. 3.3 Distance to faults map. Fig. 3.4 Average annual rainfall map 26 .

(38) Fig. 3.5 Distance to rivers map. Fig. 3.6 Elevation map. Fig. 3.7 Slope map. 27 .

(39) Fig. 3.8 Aspect map . Fig. 3.9 Curvature map . 28 .

(40) 3.6.2 Change analysis in association with environmental variables (1) Percentage of the entire area of change The changed area in each subregion of the distribution map of the environmental variables could be computed using the“Tabulate Area＂of spatial analyst tools in ArcGIS. Then the changed area in each subregion divided by that of the entire watershed gives rise to the percentage of the entire area of change, which reveals what the percentage of change in each subregion is (Wang et al. 2012).. (2) Relative change intensity (RCI) index To explicitly quantify the concentration of land cover change in a specific subregion of the entire watershed, the relative change intensity (RCI) index is applied in our study. The concept of RCI, derived from the research of Wang et al. (2012), is analogous to the Location Quotient that is a tool for comparing regional industry composition. The RCI index can be used to understand change intensity in a subregion as compared with that in other subregions or the entire watershed. The RCI value can be calculated by the following equations (Wang et al. 2012):. RCI. Rsc/Rec. (3.5). Rsc. Asc/As. (3.6). Rec. Aec/Ae. (3.7). where RCI can be expressed by the proportion of the change ratio in a subregion Rsc to that in the entire watershed Rec; Rsc equals the proportion of changed area in the subregion Asc to the entire area of the subregion As; Rec equals the proportion of changed area in the entire watershed Aec to the entire area of the watershed Ae. The RCI index has a value of 1 when the change ratio in the subregion equals that in the 29 .

(41) entire watershed, and it increases when the change ratio in the subregion is larger than that in the entire watershed.. 3.7 Logistic regression model 3.7.1 Mathematics of logistic regression In the present study, we observe some independent variables, such as elevation, slope, aspect, etc., on every pixel for which we have also determine the status of change, as either 1 if “with change” or 0 if “without change”. We would like to use these independent variables to estimate the probability that the land cover change will occur for every pixel during a defined study period. Since the dependent variable is a binary variable being 1 versus 0, it would be proper to apply the logistic regression model to describe the probability of change. The logistic regression model, based on the logistic function, can be expressed in the formula (Bai et al. 2010; Kleinbaum and Klein 2010):. 1. (3.8). where P is probability of occurrence of land cover change for every pixel, which ranges between 0 and 1 on a S-shape curve; z can be substituted for the linear sum expression of β plus β1 times X1 plus β2 times X2, and so on to βk times Xk:. z. β. β X. β X. β X. (3.9). where β , a constant term, is the intercept of the model; βi (i = 1, 2,…,k) are the slope coefficients of the model; and Xi (i = 1, 2,…,k) are independent variables of interest.. 30 .

(42) The terms β. and βi in the model representing unknown parameters could be. estimated based on the data of variables Xi and Y (status of change) of the pixels. To find the best fitting set of parameters, maximum likelihood (ML) estimation, a preferred estimation method for logistic regression, is used to estimate the parameters in the model.. 3.7.2 Data sampling and multicollinearity diagnostics To maintain the equal ratio of change presence (1)/change absence (0) and decrease the effect of spatial autocorrelation, the sample data of pixels are selected by applying a stratified random sampling technique. The sample sizes are determined by the dependent variables with different sizes of changed area. Three different sample sizes are used for the logistic regression analysis, including ten percent of the total pixels for overall landscape change trajectory (OL), two percent of the total pixels for forest-landslide trajectory (from forest to landslide, FL) and forest-channel trajectory (from forest to channel, FC), and one percent of the total pixels for vegetation recovery (from landslide to forest,VR), respectively. Model fitting via logistic regression is sensitive to collinearities among the independent variables in the model, which is the same case in linear regression (Hosmer and Lemeshow 2000). Multicollinearity results from strong correlation between independent variables and may result in incorrect signs and magnitudes of regression coefficients for dependent variables, thereby, leading to inaccurate inferences about relations between independent and dependent variables (Frans 2008). The Pearson correlation coefficient is ordinarily used to assess correlation between two variables, but it would not discern the situation where multiple independent variables may be interdependent. However, the Tolerance (TOL) and Variance Inflation Factor (VIF) could be applied to discern multicollinearity in this situation, 31 .

(43) and they are two important indexes widely used for multicollinearity diagnosis. The TOL is defined as 1－R2, where R2 is the coefficient of determination calculated by regressing each variable on all the other independent variables (Allison 1999). The lower the TOL value, the higher the multicollinearity. The TOL value below 0.2 demonstrates the presence of multicollinearity between independent variables, and the value below 0.1 reveals serious muliticollinearity among them (Menard 2002). The VIF is obtained from the reciprocal of the tolerance, and shows how inflated the variance of the coefficient is, compared to what it would be in the absence of collinearity (Allison 1999). There are no formal cutoff values for TOL and VIF to determine whether multicollinearity is present between variables. Allison (1999) suggests that a TOL value below 0.4, i.e. VIF value above 2.5, may reveal the presence of multicollinearity. In our study, variables with TOL<0.4 and VIF>2.5 are rejected from the logistic regression analysis.. 3.7.3 Analysis of model fit (1) Pseudo-R2 statistics The R2 statistic, also called the coefficient of determination, which represents the proportion of variation in the data accounted for by a linear regression model, is not appropriate for use with the logistic regression model (Hilbe 2009). The pseudo-R2 statistics are devised by statisticians to be more appropriate for logistic models and have similar properties to the R2 statistic of linear regression model. Cox & Snell R2 and Nagelkerke R2, both of which are a kind of estimate of the pseudo-R2 value, could indicate how much the variability of the dependent variable may be explained by all included predictor variables (George and Mallery 2010). Ranging from 0 to 1, the higher the Cox & Snell R2 and Nagelkerke R2 values, the better the fit of the model. Both pseudo-R2 statistics are widely used by many researchers (Bai et al. 2010; Bui et 32 .

(44) al. 2011; Lin et al. 2010) and thus are employed by the present study to assess the model fit.. (2) Classification table The classification table, summarizing the results of a fitted model, is one of the original fit tests used by statisticians with logistic regression (Hilbe 2009). The classification table is the result of cross-classifying the outcome variable, y, with a dichotomous variable whose values are derived from estimated logistic probabilities (Hosmer and Lemeshow 2000). Each estimated probability is compared to a defined cutoff point whose most commonly used value is 0.5. If the estimated logistic probabilities are greater than 0.5, then we categorize the cases as the value 1 for the independent variable; and if the probabilities are less than 0.5, then we categorize the cases as the value 0. Based on the 2 × 2 classification table, three kinds of probabilities could be estimated. One is what proportion of "changed" sites are correctly classified (sensitivity), another is what proportion of "unchanged" sites are correctly classified (specificity), and the other is what overall percentage is correctly classified. For a study evaluating the discrimination performance of the model derived from logistic regression, Pearce and Ferrier (2000) uses a 2 × 2 classification table to examine the agreement between predictions and actual observations. Therefore, crosstabulating the observed values with the predicted values for the dependent variable helps us determine how well the model identifies changes in the landscape of Taimali watershed.. 33 .

(45) 3.7.4 Model Validation Implementing a validation for the prediction results is regarded as one of the important and essential study procedures in prediction modeling, and the prediction model and image will have hardly any scientific significance without some kind of validation (Chung and Fabbri 2003). In addition, model validation may be especially important when the fitted model is applied to predict the outcome for future subjects (Hosmer and Lemeshow 2000). To evaluate data mining models, it is important to separate the data into two sets: one with 70 percent of the sample data, for training the model, and one with 30 percent of the sample data, for testing the model. The accuracy of logistic regression analysis in modeling the probability of land cover change in the Taimali watershed is validated by evaluating the performance of the classification results obtained from the testing data sets, and by computing the AUC values, RCI values, and percentage of the area of observed change trajectories in various classes of probability of change.. (1) ROC analysis A more detailed description of classification accuracy could be acquired from the area under the ROC (Receiver Operating Characteristic) curve that provides a measure of the logistic regression model’s ability to classify cases into one of two groups (Hosmer and Lemeshow 2000). The ROC curve is constructed by calculating multiple sensitivity/specificity pairs and drawing a plot of sensitivity (true positive rate) on the vertical axis versus 1-specificity (false positive rate) on the horizontal axis for all possible cutoff points (Hilbe 2009). The area under the ROC curve (AUC) ranges from 0 to 1 with higher values indicative of better fit. As a general rule, an AUC value above 0.7 is considered acceptable while a value exceeding 0.9 is regarded as an outstanding discrimination (Hosmer and Lemeshow 2000). AUC 34 .

(46) values of 0.95 and larger are highly unlikely (Hilbe 2009). Since the ROC curve provides a useful way to evaluate the performance of classification results and compare the different models, it is applied to implement the goodness-of-fit assessment of the model. In the present study, the ROC curves are plotted according to the true positive rate of identified change presence and false positive rate of identified change presence, as the values of the cutoff point varied.. (2) Relative change intensity (RCI) index Relative change intensity (RCI) index could measure the concentration of land cover change in a subregion of the whole watershed; hence the index is also utilized in this section to validate the model performance. The calculation method refers to equation 3.5-3.7.. 35 .

(47) 4. RESULTS. 4.1 Landscape metrics analysis Landscape metrics are used to measure the landscape composition and spatial configuration, and thus analysis of the metrics at any different time point enables the recognition of dynamic changes in the landscape pattern. This study used nine landscape metrics to quantify the landscape pattern, and demonstrates that natural disturbances have exerted a great influence on landscape changes in the Taimali watershed.. 4.1.1 Landscape composition changes Three land cover maps of the study area from the period 2005–2011 were obtained from satellite image classification (as shown in Fig. 4.1). Based on the primary objective of this study, four major land cover types were detected and delineated. By calculating the two class-level composition metrics for each landscape element, both the class area and percentage distribution during 2005–2011 were produced (Table 4.1). As listed in Table 4.1, the CA of the landslide patches rapidly enlarged from 2.08 km2 in 2005 to 10.68 km2 in 2008, and further extended to 17.60 km2 by 2011. In addition, the CA of the channel corridors showed a similar tendency to that of the landslide patches, and extended from 7.51 km2 in 2005 to 15.05 km2 in 2011. In contrast, however, the changes in CA for the forest matrix and human-made patches showed a substantial decline. For example, forest covered an area of 198.04 km2 in 2005, but by 2008 had been reduced to 187.78 km2, with a further reduction to 177.14 km2 in 2011. The CA of human-made patches decreased from 3.90 km2 in 2005 to 1.75 km2 in 2011.. 36 .

(48) For the percentage land area (PLAND), forest matrix was found to be the most dominant land cover type in the Taimali watershed over the study period, and its PLAND was maintained at over 83% at all times. Although the forest type was the landscape matrix, its PLAND have decreased to a great extent from 93.62% in 2005 to 88.77% in 2008, and further to 83.74% in 2011. The PLAND of human-made patches also decreased from 1.84% in 2005 to 0.83% in 2011. In contrast, however, the PLAND for both landslide patches and channel corridors showed a considerable increase, particularly for landslide patches, which increased significantly increased from 0.99% in 2005 to 8.32% in 2011. The PLAND of channel corridors increased from 3.55% to 7.11% during the period 2005–2011. In total, the CA of landslide patches and channel corridors expanded 7.5 times (15.52 km2) and 1.0 times (7.54 km2), respectively, between 2005 and 2011. However, the CA of forest matrix and human-made patches shrank by 10.6% (20.90 km2) and 55.1% (2.15 km2), respectively.. Table 4.1 The area distribution of landscape elements from 2005 to 2011 Elements. Landslide patches. 2005. 2008. CA. PLAND. CA. PLAND. CA. PLAND. (km2). (%). (km2). (%). (km2). (%). 2.08. 0.99. 10.68. 5.05. 17.60. 8.32. 198.04. 93.62. 187.78. 88.77. 177.14. 83.74. Human-made patches. 3.90. 1.84. 2.48. 1.17. 1.75. 0.83. Channel corridors. 7.51. 3.55. 10.60. 5.01. 15.05. 7.11. Forest matrix. 37 . 2011.

(49) 2005. 2008. 2011. Fig. 4.1 Classified land cover maps of the study area from 2005 to 2011. 38 .

(50) 4.1.2 Spatial configuration changes By computing seven class-level spatial configuration metrics for each landscape element, the trend in landscape change between 2005 and 2011 was obtained (Fig. 4.2). The patch area is perhaps the single most important and useful piece of information contained in the landscape, and has a large amount of ecological utility in its own right (McGarigal et al. 2012). In addition, the calculation of patch edge provides a valuable measure of the extent of spatial pattern dissection (Cardille and Turner 2002). Edge metrics are even more meaningful when interpreted in combination with other metrics, such as mean patch area (Olsen et al. 2007). NP is also a fundamental metric, and is able to measure the degree of subdivision or fragmentation of a particular patch type (McGarigal et al. 2012). The changes in metrics of AREA_MN, ED, and NP are demonstrated in Fig. 4.2 a–c. In terms of forest matrix, AREA_MN displayed a significant decrease, while ED and NP showed a considerable increment. This suggests that forest fragmentation has created a number of smaller patches with a greater number of dissected edges. In contrast with the forest matrix, channel corridors showed the opposite tendency. Due to their connectivity, the NP of channels was always maintained at 1 throughout the study period. In addition, the increasing AREA_MN and decreasing ED for channels indicates that their patch area expanded and their edges slightly decreased, owing to the expansion of channels in relation to flooding. However, the NP, AREA_MN, and ED for landslide patches all considerably increased over the study period, showing that landslide patches expanded in amount, size, and edge. Landslides, considered a major sediment source, can supply a large amount of unstable sediment to a basin (Koi et al. 2008). The expansion of landslide patches is related to a series of natural disturbances, and it is considered that this occurrence should receive attention in relation to hazard management. In comparison with values 39 .

(51) of landslides patches, the values NP, AREA_MN, and ED for human-made patches all decreased gradually, thereby indicating that human-made patches shrank in amount, size, and edge. Both SHAPE_MN and FRAC_MN measure the complexity of a patch shape. SHAPE_MN has a value of 1 when the patch shape is square and highly compact, and its value increases without an upper limit as the patch shape becomes more irregular (McGarigal et al. 2012). The value of FRAC_MN ranges from 1 to 2, approaching 1 for patch shapes with simple perimeters such as squares or circles, and approaching 2 for patch shapes with highly convoluted perimeters (Mandelbrot 1977). As shown in Fig. 4.2 d and e, both these metrics show a similar trend in change, and the values of both show a decreasing trend for all land cover types except for forest. For forest matrix, SHAPE_MN increased from 1.56 in 2005 to 1.58 in 2011, and FRAC_MN increased from 1.064 in 2005 to 1.072 in 2011. These results suggest that the shape of the forest land cover type became more irregular, but the shape of the other land cover types transformed into more compact and simple shapes. ENN_MN and AI are applied to measure the extent of isolation and aggregation, respectively, for the patch types. ENN_MN approaches 0 when the distance to the nearest patch of the same type is reduced, and AI ranges from 0 to 100, reaching 100 when the focal patch type is maximally aggregated into a single and compact patch (McGarigal et al. 2012). As shown in Fig. 4.2 f and g, the decreasing ENN_MN but increasing AI for landslide patches indicates that these areas have increasingly become patches of the same type and spatial aggregation. Conversely, increasing ENN_MN and decreasing AI for forest matrix and human-made patches reveals that these two patch types have become more spatially isolated and dispersed. As channel corridors represent only one patch, the ENN_MN of channel corridors could not be defined, and has not been shown in Fig. 4.2 f. However, the AI of channel corridors 40 .