Spatial Autocorrelation Analysis of Multi-Scale Damaged Vegetation in the Wenchuan Earthquake-Affected Area , Southwest China

Major earthquakes can cause serious vegetation destruction in affected areas. However, little is known about the spatial patterns of damaged vegetation and its influencing factors. Elucidating the main influencing factors and finding out the key vegetation type to reflect spatial patterns of damaged vegetation are of great interest in order to improve the assessment of vegetation loss and the prediction of the spatial distribution of damaged vegetation caused by earthquakes. In this study, we used Moran’s I correlograms to study the spatial autocorrelation of damaged vegetation and its potential driving factors in the nine worst-hit Wenchuan earthquake-affected cities and counties. Both dependent and independent variables showed a positive spatial autocorrelation but with great differences at four aggregation levels (625 × 625 m, 1250 × 1250 m, 2500 × 2500 m, and 5000 × 5000 m). Shrubs can represent the characteristics of all damaged vegetation due to the significant linear relationship between their Moran’s I at the four aggregation levels. Clustering of similar high coverage of damaged vegetation occurred in the study area. The residuals of the standard linear regression model also show a significantly positive autocorrelation, indicating that the standard linear regression model cannot explain all the spatial patterns in damaged vegetation. Spatial autoregressive models without spatially autocorrelated residuals had the better goodness-of-fit to deal with damaged vegetation. The aggregation level 8 × 8 is a scale threshold for spatial autocorrelation. There are other environmental factors affecting vegetation destruction. Our study provides useful information for the countermeasures of vegetation protection and conservation, as well as the prediction of the spatial distribution of damaged vegetation, to improve vegetation restoration in earthquake-affected areas.


Introduction
Major earthquakes are geologically disastrous and can trigger serious ecological degradation.In addition to causing massive human casualties, such earthquakes induce vegetation destruction [1], heavy economic losses [2], biodiversity reduction [3], aggravated sedimentation [4], and landscape fragmentation [5].Ecological restoration in earthquake-affected areas is a slow process due to the unstable, crude, and nutrient-poor soil conditions left in the wake of secondary geo-hazards [1,3,6].To improve the restoration process and to hasten vegetation recovery in earthquake-degraded ecosystems, many researchers have opted to report the control factors of post-earthquake vegetation cover and plant species compositions using remote sensing data and field surveys [6][7][8][9].However, the effects of topographic and seismic factors on vegetation destruction in a major earthquake have not been reported, leaving a great need for insight into these factors when assessing vegetation loss and predicting the spatial distribution of damaged vegetation after earthquakes worldwide.Hence, during the design and implementation of restoration programs, it is necessary to elucidate the main factors influencing the degree of vegetation damage.
The catastrophic 8.0 Ms Wenchuan earthquake struck Sichuan Province, China on May 12, 2008.It killed at least 68,000 people and caused an estimated 845.1 billion Ren Min Bi (RMB) in direct economic losses [10].Wild and horticultural vegetation were also seriously damaged and subsequently buried by debris flow, landslides, and formerly dammed lakes [1].It is estimated that across the Sichuan province, 32.867 × 10 4 ha of vegetation destruction and 2098.63 × 10 4 m 3 of stocking volume loss were caused by the earthquake and subsequent geo-hazards [11].
Vegetation can naturally recover from large disturbances via succession, but in the wake of severe disturbances, such a recovery takes longer without human facilitation [12,13].Therefore, vegetation recovery programs are of great importance in the years after catastrophic earthquakes and are a regular restoration method applied to degraded lands [14].Several studies have assessed post-earthquake vegetation recovery and its effects on soil erosion control [15][16][17].These studies used remote sensing images to analyze the recovery potential of damaged vegetation in earthquake-affected areas.While such studies certainly improved vegetation recovery through informing restoration protocols [17], there is still a need to evaluate which topographic and seismic factors significantly influenced vegetation destruction, particularly with respect to specific countermeasures for different vegetation types.However, such studies are lacking because unaided natural recovery occurred in over half of the total restoration area [17].Our study aims to fill in the gaps and contribute to the growing body of knowledge on regionalized countermeasures for the vegetation and spatial distribution prediction of vegetation destruction of earthquake-affected areas by presenting a spatial autocorrelation analysis of the damaged vegetation and its influencing factors.
Conventional statistical methods based on linear and logistic regressions assume the data and random distribution of their residuals to be statistically independent and identically distributed [18].However, the spatial land use data have a spatial dependency, a phenomenon known as spatial autocorrelation, which contains useful information but requires appropriate statistical methods to analyze [19].Spatial autocorrelation can be used to measure the degree of spatial association of random variables with nearby variables across a geo-referenced space [18].Though the spatial autocorrelation analysis could be seen as a methodological disadvantage [19], applying it has successfully described the spatial variability of land use [19][20][21].Previous studies show that the released energy of the Wenchuan earthquake caused a sudden dislocation in the approximate 300 km-long Yingxiu-Beichuan fracture along the faults of the Longmengshan fault system [1,22], indicating that spatial dependency may exist in the destruction of earthquake-affected areas.However, to our knowledge, the spatial autocorrelation analysis has rarely been used to characterize the spatial structure of damaged vegetation, especially at different aggregation levels.In addition, our study area was the main component of the restoration areas where the Chinese government has carried out a nearly 160 billion USD recovery plan [23].
Our research provides meaningful targeted information for future vegetation recovery, vegetation protection and conservation, and the prediction of vegetation destruction caused by earthquakes.
Nine severely damaged cities and counties-including Jiangyou, Mianzhu, Shifang, Beichuan County, Anxian County, Maoxian County, Pingwu County, Qingchuan County, and Wenchuan County-were chosen as our study area, where we collected data on topographic factors (including slope angles, orientations, distance to river, and elevation), seismic factors (including seismic intensity and distance to fault), human accessibility (including distances to road and residential area), and damaged vegetation.To thoroughly understand the relationship between the damaged vegetation and influencing factors in the Wenchuan earthquake-affected area, we used a spatial autocorrelation analysis to reveal the regional distribution, scale effect, and main influencing factors of damaged vegetation.
The main objectives of this study are (1) to characterize the spatial variability of damaged vegetation at different aggregation levels in the Wenchuan earthquake-affected area, (2) to identify the representative components of damaged vegetation, and (3) to clarify the main influencing factors on the spatial variability of damaged vegetation.

Study Area and Geological Setting
As seen in Figure 1, the nine cities and counties chosen for our study (N 30  2 and are located at the eastern margin of the Tibetan Plateau in the most severely earthquake-affected areas of west Sichuan Province, China.They are characterized by semiarid hot and subtropical humid monsoon climates, with average annual precipitations of about 506 mm and 1503 mm, respectively.The mean annual temperatures are 13.8 and 15.0 °C respectively, and the mean annual evapotranspirations are 1392 mm and 1156 mm.From north to east, the elevation gradually decreases from 5230 to 462 m above sea level (a.s.l.), with a slope from 0 • to 76 • .The dominant soil types include brown soil, valley brown soil, and yellow soil.
The Wenchuan earthquake's energy source was the famous active thrust Longmenshan fault, located at the conjunction of the Sichuan Basin western margin and the eastern edge of the Tibetan Plateau [24,25].The Longmenshan fault initially formed in the Mesozoic and further developed throughout the Cenozoic from the collision of the Indian Plate onto the Eurasian Plate, resulting in the rapid uplift and tectonic deformation of the Longmenshan Mountain range [26].The Longmenshan fault includes three regions: the Back fault, the Central fault, and the Front fault.As reverse faults, the Back and the Front faults consist of a series of reverse imbricate faults with a NE 250-450 orientation.They extend from Wenchuan County to Maoxian County and from Pengxian County to Guanxian County, respectively.Unlike the Back and the Front fault, the Central fault is a thrust fault with a NE 350-450 orientation, extending from the town of Yingxiu to Beichuan County.The Central fault and the Front fault both extend for 200 km and show dextral strike-slip movements along the eastern margin of the Tibetan Plateau, with higher slip rates (1-10 mm•yr −1 ) than throw rates (<1 mm•yr −1 ) since the Epipleistocene [25].

Data Collection
To investigate the spatial variation of damaged vegetation in the Wenchuan earthquake-affected area, information on the vegetation damage (Table 1) was extracted from an overlay analysis between the geo-hazard interpretations of remote sensing images from the China-Brazil Earth Resources Satellite (CBERS-02B) (Figure 1) and pre-earthquake vegetation distribution supplied by the Sichuan Forest Inventory and Planning Institution (SFIPI) [27].The CBERS-02B images have a swath of 113 km and 5 wave bands, as seen in Table 2. To avoid the negative effects of cloud and haze on remote sensing images, 49 CBERS-02B CCD images containing 5 spectral bands (including red, green, blue, NIR1, and NIR2) with a pixel resolution of 19.5 m were acquired from May 14 to 28, 2008, excluding May 21 and 27 (Table 3).This was the period during which the Chinese government reset the orbit parameter of CBERS-02B to monitor the geo-hazards daily following the Wenchuan earthquake.Due to spatial incompatibility between band-5 and the other bands in the CBERS-02B CCD images, we used the methods described by Li et al. [28] to re-project band-5 and then generated composite images using the re-projected band-5 and the other four bands.We applied the geometric correction using a second-order polynomial model with 24 ground control points (GCP) to ensure accuracy, based on primary scale Chinese national topographic maps (1:50,000).The spatial inconsistency errors were bounded within 1 pixel in this study.Subsequently, we interpreted the distribution of geo-hazards (including rock avalanches, landslides, debris flows, and landslide-dammed lakes) and used the GIS overlay analysis of pre-earthquake vegetation distribution, interpreted by SFIPI and the post-earthquake geo-hazards distribution, to determine the distribution of damaged vegetation, which was defined as a certain proportion of the pre-earthquake vegetation in a pixel, removed or destroyed by geo-hazards.landslides, debris flows, and landslide-dammed lakes) and used the GIS overlay analysis of pre-earthquake vegetation distribution, interpreted by SFIPI and the post-earthquake geo-hazards distribution, to determine the distribution of damaged vegetation, which was defined as a certain proportion of the pre-earthquake vegetation in a pixel, removed or destroyed by geo-hazards.In this study, the digital elevation model (DEM) of the 25 m × 25 m cell size Sichuan Province study area (Figure 2), provided by the Institute of Mountain Hazards and Environment, was used to calculate and extract the topographic factors, including river systems, elevation, slope, and aspect (Figure 2), with the spatial analyst tools in ARC/INFO 9.2.The fundamental geographic information features of 1:250,000 for the Sichuan Province, including roads and residential areas, were chosen as proxies for human accessibility.The fault data were extracted from the map of horizontal peak ground acceleration (PGA) in the Wenchuan earthquake-affected area, created by Tang et al. [29].In addition, the seismic intensity zone data were extracted from the Ms 8.0 Wenchuan earthquake seismic intensity map created by Xu et al. [30].All data were combined with the damaged vegetation data using an overlay analysis.In this study, a set of 19 potential drivers of earthquake damage degree (Table 4) were selected for the damaged vegetation at four different aggregation levels.The highest resolution in the data set is 1 × 1 (625 × 625 m) grid cells.With the accumulative data of 2 × 2, 4 × 4, and 8 × 8 grid cells, a total of four aggregation levels were created to simulate different scales using the methods described in Overmars et al. [19].
Forests 2018, 9, x FOR PEER REVIEW 6 of 21 In this study, the digital elevation model (DEM) of the 25 m × 25 m cell size Sichuan Province study area (Figure 2), provided by the Institute of Mountain Hazards and Environment, was used to calculate and extract the topographic factors, including river systems, elevation, slope, and aspect (Figure 2), with the spatial analyst tools in ARC/INFO 9.2.The fundamental geographic information features of 1:250000 for the Sichuan Province, including roads and residential areas, were chosen as proxies for human accessibility.The fault data were extracted from the map of horizontal peak ground acceleration (PGA) in the Wenchuan earthquake-affected area, created by Tang et al. [29].In addition, the seismic intensity zone data were extracted from the Ms 8.0 Wenchuan earthquake seismic intensity map created by Xu et al. [30].All data were combined with the damaged vegetation data using an overlay analysis.In this study, a set of 19 potential drivers of earthquake damage degree (Table 4) were selected for the damaged vegetation at four different aggregation levels.The highest resolution in the data set is 1 × 1 (625 × 625m) grid cells.With the accumulative data of 2 × 2, 4 × 4, and 8 × 8 grid cells, a total of four aggregation levels were created to simulate different scales using the methods described in Overmars et al. [19].

Spatial Autocorrelation Analysis
Spatial autocorrelation is defined as the phenomenon of geographic proximity's potential to determine the values of random variables over distances that are similar or dissimilar to randomly associated pairs of observations [31].Researchers have used spatial autocorrelation to describe the dependency and covariance of variables within a spatial neighborhood [32][33][34].

Global Spatial Autocorrelation Statistics
Correlograms calculated with univariate (Moran's I or Geary's c) and multivariate data (Mantel correlogram) can be used to quantify the spatial dependency per distance class [35][36][37].In this study, we used correlograms of the Moran's I as described in Overmars et al. [19].The formula is expressed as follows [31]: Moran's I for h = i: where y h and y i are the values of the observed variable at sites h and i, respectively; the values of w hi are weight matrices; and W in a (n × n) weight matrix is the sum of the weights w hi for a given distance class, i.e., the number of pairs used to calculate the coefficient.Moran's I can be calculated using a free software, Geoda [38], in which the value of Moran's I generally ranges from −1 to 1. High positive Moran's I values indicate a positive autocorrelation, which implies the clustering of similar values; low negative values, by contrast, represent a negative autocorrelation, which implies the clustering of dissimilar values.A value close to zero indicates no autocorrelation [39][40][41].

Local Spatial Autocorrelation Statistics
Global spatial autocorrelation may not pick up the aberrant local spatial pattern that causes some local patterns to be opposite to the global spatial trend [41].To assess the influence of aberrant local spatial patterns on the magnitude of the global statistic, Anselin [42] outlined a class of local indicators of spatial association (LISA) (1) to assess the extent of significant local spatial clustering around an individual location, (2) to indicate the local pockets of spatial non-stationarity, and (3) to identify outliers or spatial regimes.The formula expressed in the form of a LISA is as follows [42]: Local Moran's I can also be calculated by "GeoDa" with some visualization results-the significance map of Moran's I and spatial clustering [38].

Spatial Autoregressive Models
To avoid the effects of the biased estimation of error variance, t-test significance levels, and the overestimation of R 2 caused by conventional statistical methods when dealing with spatial autocorrelation, a commonly used spatial autoregressive model is Equation (3) created by Anselin [43].From Equation (3), we can derive several specific models by imposing restrictions [19,44].Setting W 2 × 0 produces a spatial lag model (SLM) or mixed regressive-spatial autoregressive model (Equation ( 4)) and expresses how the magnitude of a decision variable for a geographical grid depends on the magnitudes of the decision variables set by other grids [44,45].This model can explain the additional variation y over the spatial sample via additional explanatory variables in the matrix X [43].Setting W 1 = 0 produces a spatial error model (SEM) (Equation ( 5)) and describes the spatial autocorrelation in the disturbances [44].
In Equations ( 3)-( 5), y is an n × 1 vector of cross-sectional dependent variables and X contains an n × k matrix of explanatory variables with an associated k × 1 regression coefficient vector β.W 1 and W 2 are n × n spatial weight matrices describing the interconnections between different locations, ρ is the coefficient of the spatially lagged dependent variable, λ is the coefficient of the spatially correlated errors, and ε is a vector or random error term [19,46].
In this study, two spatial autoregressive models (Equations ( 4) and ( 5)) are used.The maximum likelihood estimation is used to estimate the parameters in GeoDa, since an ordinary least squares estimation for spatial autoregressive models is biased [38].We use the pseudo R 2 (defined as the ratio of variance of the predicted values over the variance of the observed values for the dependent variable) and the value of the maximized log likelihood (LIK) to assess the fit in the spatial autoregressive models and all other models [19,44].In addition, we select the ρ value and its significance as the percentage of the prediction that is obtained using the spatial part and the significance of spatial autocorrelation [19].

Data Analysis
GeoDa was used to calculate the Moran's I statistic.The weight matrices were established based on distance [47] and were equal within a lag.The lag sizes chosen were 625 m, 1250 m, 2500 m, and 5000 m, respectively at the 1 × 1, 2 × 2, 4 × 4, and 8 × 8 aggregation levels.Correlograms were calculated and compared for all combinations of the four aggregation levels, driving factors, and ten vegetation types.Within a same lag, we used correlograms pairs between the Moran's I of the percentage of all damaged vegetation and ten vegetation types to create regression models to determine their relationships.
To calculate SLM and SEM, GeoDa was used as well.We used the independent variables selected method of de Koning et al. [48] to remove insignificant variables one by one until a model with only significant variables remained.

Relationships between Correlograms of All Destructed Vegetation and Ten Vegetation Types at the Four Cell Levels
As seen in Figure 3, the correlograms of the Moran's I (distances as 625 m, 1250 m, 2500 m, and 5000 m of a lag respectively) of the percentage of all the damaged vegetation and ten vegetation types are presented in 1 × 1, 2 × 2, 4 × 4, and 8 × 8 cells.The correlograms showed large differences in Moran's I between vegetation types.All the damaged vegetation and ten vegetation types, represented in Figure 3, showed a significantly positive spatial autocorrelation (P < 0.001), which decreased gradually and still had a significant autocorrelation with distance at all four aggregation levels.All damaged vegetation, economic forests, evergreen coniferous forests, meadows, and shrubs had a comparatively high Moran's I across the lag distances at all four aggregation levels.Moran's I of ten vegetation types were significantly related to that of all the destructed vegetation across the lag distances at all four aggregation levels (Table 5), indicating that the components of all the destructed vegetation could represent the characteristics of its spatial autocorrelation.However, unlike the other vegetation types having exponential function or power function relationships with all destructed vegetation, shrubs had linear function relationships with all the destructed vegetation at all four aggregation levels (Table 5).To calculate SLM and SEM, GeoDa was used as well.We used the independent variables selected method of de Koning et al. [48] to remove insignificant variables one by one until a model with only significant variables remained.As seen in Figure 3, the correlograms of the Moran's I (distances as 625 m, 1250 m, 2500 m, and 5000 m of a lag respectively) of the percentage of all the damaged vegetation and ten vegetation types are presented in 1 × 1, 2 × 2, 4 × 4, and 8 × 8 cells.The correlograms showed large differences in Moran's I between vegetation types.All the damaged vegetation and ten vegetation types, represented in Figure 3, showed a significantly positive spatial autocorrelation (P < 0.001), which decreased gradually and still had a significant autocorrelation with distance at all four aggregation levels.All damaged vegetation, economic forests, evergreen coniferous forests, meadows, and shrubs had a comparatively high Moran's I across the lag distances at all four aggregation levels.4.

Relationships between correlograms of all destructed vegetation and ten vegetation types at the four cell levels
Table 5.The regression models between the Moran's I of ten vegetation types and those of all the destructed vegetation across the lag distance of 625 m, 1250 m, 2500 m, and 5000 m comparatively at four corresponding aggregation levels.The vegetation type abbreviations are defined in Table 4.

Correlograms of potential driving factors at the 1 × 1 cell level
Except for four slope aspects (the east, south, west, and north slopes) with comparatively lower Moran's I (Moran's I < 0.200), all potential driving factors showed a significantly positive spatial autocorrelation (P < 0.001) with a comparatively high Moran's I (ranging from 0.338 to 0.793) at the 1 × 1 cell level, which also decreased gradually with distance (Figure 4a-d).As seen in Figure 5, the cluster and significance maps of LISA for all damaged vegetation at the 1 × 1 aggregation level indicate that the high-high zone occupied 4146 out of 4854 grid cells with <0.05 significance levels.This means that a clustering of similar high coverage of damaged vegetation occurred in the Wenchuan earthquake-affected area, causing a highly fragmented wildlife habitat at range-wide scales and the loss of a global biodiversity hotspot [49,17].4.

LISA Map of All Damaged Vegetation at the 1 × 1 Cell Level
As seen in Figure 5, the cluster and significance maps of LISA for all damaged vegetation at the 1 × 1 aggregation level indicate that the high-high zone occupied 4146 out of 4854 grid cells with <0.05 significance levels.This means that a clustering of similar high coverage of damaged vegetation occurred in the Wenchuan earthquake-affected area, causing a highly fragmented wildlife habitat at range-wide scales and the loss of a global biodiversity hotspot [17,49].

Correlograms of potential driving factors at the 1 × 1 cell level
Except for four slope aspects (the east, south, west, and north slopes) with comparatively lower Moran's I (Moran's I < 0.200), all potential driving factors showed a significantly positive spatial autocorrelation (P < 0.001) with a comparatively high Moran's I (ranging from 0.338 to 0.793) at the 1 × 1 cell level, which also decreased gradually with distance (Figure 4a-d).As seen in Figure 5, the cluster and significance maps of LISA for all damaged vegetation at the 1 × 1 aggregation level indicate that the high-high zone occupied 4146 out of 4854 grid cells with <0.05 significance levels.This means that a clustering of similar high coverage of damaged vegetation occurred in the Wenchuan earthquake-affected area, causing a highly fragmented wildlife habitat at range-wide scales and the loss of a global biodiversity hotspot [49,17].

Spatial Autocorrelation in Residuals
In this study, we used the residuals of the standard linear regression models in the first lag constructed by Overmars et al. [19] to determine the scale threshold of the aggregation levels.The residuals' spatial autocorrelation in the first lag of all aggregation levels were significant in all the destructed vegetation and ten vegetation types with Moran's I > 0.130 except for the bamboo forest vegetation type at the aggregation level 8 × 8 (Table 6).Moreover, the residuals' spatial autocorrelation weakened at aggregation levels 4 × 4 and 8 × 8, suggesting that the aggregation level 8 × 8 may be a scale threshold for spatial autocorrelation.The numbers in bold font are not significant (P > 0.001).

Results of Spatial Autoregressive Models
We used the standard linear model based on ordinary least squares, spatial lag model (SLM), and spatial error model (SEM) as spatial autoregressive models to illustrate their differences in all damaged vegetation at aggregation level 1 × 1 (Table 7).The standard linear model contains thirteen significant variables (P > 0.05) selected using a stepwise regression from a set of 19 potential drivers with parameters including the measure of fit (R 2 ), coefficient estimate, standard error, t-test value, associated probability, and log likelihood (LIK) (Table 7(1)).We used the LIK to compare the goodness-of-fit of three spatial models.From Table 7(2), we found that the spatial lag model had smaller estimated regression coefficients and standard errors, as well as a higher R 2 , using the same variables in the standard linear model.The significance of the parameters decreased as well, and the variable "Percentage of south slope (180-270 • )" was no longer significant (P < 0.05).Moreover, the spatial lag model has a higher LIK than the standard linear model, indicating a better goodness-of-fit.
Based on the results of the spatial lag model, we excluded the insignificant variable to construct the spatial error model.From Table 7(3), we found that the pseudo R 2 and the LIK of the spatial error model decreased from 0.6685 to 0.6416 and from 64,938.3 to 61,732.3 compared with the spatial lag model, indicating that excluding the insignificant variable only results in a slight reduction of the goodness-of-fit.
The spatial autocorrelation in the residuals of spatial lag and spatial error models weakened (Figure 6).Furthermore, the difference in their spatial autocorrelations was negligible.Compared to the spatial lag and spatial error models, the standard linear model had larger residuals with a positive spatial autocorrelation.For example, as seen in Table 8, more grids with a comparatively higher/lower original value of the residuals of the standard linear model are switched from positive/negative to negative/positive or less positive/negative.Therefore, the residuals of the spatial error model are considerably lower.We list some models' outputs for all the destructed vegetation and ten vegetation types (Table 9).Based on the guidance of Overmars et al. [19], which states that the significance of ρ is a criterion of whether to apply the spatial error model, 43 out of our 44 models qualified for applying the spatial model.All ρ values at all four aggregation levels were significant except for bamboo forests at the 8 × 8 aggregation level, indicating that the vegetation types may have a patch size smaller than the cell size over the 8 × 8 aggregation level.Therefore, the aggregation level 8 × 8 may be a scale threshold for spatial autocorrelation in damaged vegetation of the Wenchuan earthquake-affected area.We list some models' outputs for all the destructed vegetation and ten vegetation types (Table 9).Based on the guidance of Overmars et al. [19], which states that the significance of ρ is a criterion of whether to apply the spatial error model, 43 out of our 44 models qualified for applying the spatial model.All ρ values at all four aggregation levels were significant except for bamboo forests at the 8 × 8 aggregation level, indicating that the vegetation types may have a patch size smaller than the cell size over the 8 × 8 aggregation level.Therefore, the aggregation level 8 × 8 may be a scale threshold for spatial autocorrelation in damaged vegetation of the Wenchuan earthquake-affected area.
no information in the literature about whether the components of research objectives can be a valid indicator that reflects the characteristics of research objectives' spatial autocorrelation.To elucidate the main factors influencing the degree of vegetation damage and to find out the key vegetation type to reflect the spatial autocorrelation of destructed vegetation is the basis of regionalized countermeasures in vegetation protection and conservation in earthquake-affected areas.Thus, statistical models should be applied to explain the correlation between vegetation destruction, its components, and the driving environmental factors using methods such as spatial generalized least-squares (GLS) or autoregressive models [31].
Vegetation destruction caused by the Wenchuan earthquake reached 1249.47 km 2 in our study area, accounting for 4.76% of the area of the nine worst-hit cities and counties [1].There was a significantly positive spatial autocorrelation in all the damaged vegetation and its components.However, the Moran's I values of the surface percentage within cells of the 11 vegetation types are clearly different at all four aggregation levels.This means that vegetation types have different patterns and different spatial characteristics within the different categories of all damaged vegetation.In fact, any categorization made with our prior knowledge can cause the differences (y h − y i ) for any distance d, independent of the location where the differences are calculated [31].Therefore, each vegetation category reflecting the spatial distribution of all the damaged vegetation or not may show similar or different spatial patterns to all the damaged vegetation.To assess the importance of each vegetation category in the whole and to determine the main damaged vegetation type contributing substantially towards the design of restoration programs by demonstrating the importance of accounting for species selection and local conditions, we set up regression models between the Moran's I of ten vegetation types and those of all the destructed vegetation across the corresponding lag distances at all four aggregation levels.We found that the Moran's I of ten vegetation types had significant relationships (including exponential function, power function, and linear function) with that of all the destructed vegetation.However, unlike the other vegetation types, only shrubs had significantly positive linear relationships with all the destructed vegetation at all four aggregation levels, indicating that the Moran's I of shrubs decreases with distance following the same rule as all destructed vegetation.In other words, shrubs can represent the characteristics of the spatial structure of all damaged vegetation to a certain extent and may play an important proxy role in designing vegetation restoration plans for the Wenchuan earthquake-affected area.Similar effects to shrub vegetation type were reported in previous studies, showing that it is the main vegetation type distributed between 1800 and 3400 m a.s.l. for headwater and animal habitat conservation efforts in our study area [1,55].The other vegetation types show a significant exponential function or power function with all destructed vegetation, indicating that their spatial autocorrelation attenuate more quickly than all the destructed vegetation with increasing distance.This means that the other vegetation types cannot reflect the spatial autocorrelation of all the destructed vegetation except shrubs when the patch size is smaller than the cell size over a certain distance.This result is consistent with other studies demonstrating significant impacts of patch size on spatial structure [19].
The spatial structures of 19 potential driving factors also showed a significantly positive spatial autocorrelation (Figure 4).This means that the characteristics of damaged vegetation may be explained by driving forces, depending on the response of land use and cover to driving factors [19].Cluster and significance maps of LISA can enable us to assess the interactions between sites in close proximity to each other by comparing their values to their neighbors and identifying the local spatial structure and instability [41].In this study, LISA maps delineated damaged vegetation zones in accordance with the types of spatial autocorrelation.Therefore, spatial autoregressive models should be used to detect the relationships between damaged vegetation and driving factors.
In all the destructed vegetation and ten vegetation types at all aggregation levels (Figure 3 and Table 6), the residuals' spatial autocorrelation of the standard linear regression is less than that of the original data, indicating that the selected driving factors used in the linear regression model are partly responsible for vegetation destruction or that damaged vegetation is a response to the spatially autocorrelated selected driving factors.However, the residuals' spatial autocorrelation is still significant, indicating that the standard linear model is not enough to explain all spatial patterns.Spatial autoregressive models can be used to reduce or remove the spatial patterns of standard linear model residuals [53].Though the pseudo R 2 in spatial models cannot be used in comparison with the traditional R 2 in the standard linear model for a measure of fit, the value of the maximized log likelihood (LIK) can be used to determine the goodness-of-fit of different models [44].In this study, spatial autoregressive models have a higher LIK than the standard linear model, indicating a better goodness-of-fit (Table 7).Comparing the residuals of the spatial lag model and spatial error model with the standard linear model (Figure 6), we found that they are considerably lower when using the spatial models.Though it may seem controversial or unsatisfactory that the prediction of a variable should use neighboring values, unlike the standard linear model which only focuses on the dominant effects caused by responses to driving forces, spatial models can deal with spatial interactions that cannot be included in the standard linear model [18,31].Moreover, compared to the spatial error model, low original values result in large negative residuals and high original values result in large positive residuals in the standard linear model (Table 8).However, our study showed that the residuals of the spatial lag model and spatial error model still had significant spatial autocorrelations.This may be due to the absence of some important environmental variables which are not available at a required spatial resolution [53,56].Therefore, more comprehensive spatial modeling techniques considering the multicollinearity of environmental variables may be a strategy to improve the model fit [53,57].
There is a scale threshold in aggregation levels that causes lower and even insignificant spatial autocorrelation at the 8 × 8 aggregation level than aggregation levels 4 × 4 (Table 6).Although the omission of variables may cause a loss of information in explaining the spatial pattern of damaged vegetation, the spatial error model is recommended because it has a higher percentage of the prediction with a higher ρ value by excluding the insignificant variable (Table 7).However, the spatial error model only explains 46.11-66.29% of the prediction at the 1 × 1 aggregation level, suggesting that we might not include all relevant environmental variables [51,53,56].A possible way to improve this is to include new variables or exclude the grids that were only slightly affected by the earthquake.
All damaged vegetation showed a significantly negative relationship with the distances to the nearest river, road, and fault (Table 6).This result is consistent with previous investigations, which highlighted that the earthquake damage degree decreased with the increase of these distances [1,24].Moreover, this result shows a significantly positive relationship with comparatively steep slopes (20-30 • and >30 • ) and percentage of high seismic intensity zone, indicating that a high potential energy in steep slopes induced secondary geo-hazards to destroy vegetation growing on the slope surface.This result is also consistent with the findings of Su and Cui [58] and Xu et al. [30], who observed that serious damage was caused in the IX, X, and XI seismic intensity zones and the areas with unstable steep slopes.
Our study has great significance in the context of vegetation protection and conservation and post-earthquake recovery in China because the Wenchuan earthquake-affected area is one of the twenty-five global hotspots for biodiversity conservation defined by the Conservation International and World Wildlife Fund [59].Our findings suggest that the spatial clustering of high-value damaged vegetation zones occurred in the Wenchuan earthquake-affected area, indicating that restoration programs should give top priority to the high-high associations in high coverage of damaged vegetation.Also, our work took an integrative approach to predict the spatial distribution of damaged vegetation after earthquakes worldwide, according to the maps of horizontal peak ground acceleration, seismic intensity, river systems, slope, and pre-earthquake vegetation distribution in earthquake-occurred areas.In addition, our results suggest that shrubs can represent the characteristics of all damaged vegetation to a certain extent.Thus, we could assess promptly the vegetation loss and can choose adaptive restoration programs according to the most representative components of the whole.Furthermore, our results might help to identify the spatial structure of land use and cover using its decisive components rather than to analyze the total data.

Figure 1 .Figure 1 .
Figure 1.The distribution of damaged vegetation in nine worst-hit cities and counties and the typical CEBRS-02B image screenshots of damaged vegetation in the Wenchuan earthquake affected area.

Figure 2 .Figure 2 .
Figure 2. The digital elevation model (DEM) of the 25 m × 25 m cell size (a), aspect (b), and slope (c) in nine worst-hit cities and counties.

Forests 2018, 9 ,
x FOR PEER REVIEW 9 of 21 calculated and compared for all combinations of the four aggregation levels, driving factors, and ten vegetation types.Within a same lag, we used correlograms pairs between the Moran's I of the percentage of all damaged vegetation and ten vegetation types to create regression models to determine their relationships.

Figure 3 .
Figure 3.The relationships between the correlograms of the Moran's I of ten vegetation types and those of all the destructed vegetation at four different aggregation levels.The points in black are significant (P < 0.001).The vegetation type abbreviations are defined in Table4.

Figure 3 .
Figure 3.The relationships between the correlograms of the Moran's I of ten vegetation types and those of all the destructed vegetation at four different aggregation levels including (a) the 1 × 1 aggregation level, (b) the 2 × 2 aggregation level, (c) the 4 × 4 aggregation level, and (d) the 8 × 8 aggregation level.The points in black are significant (P < 0.001).The vegetation type abbreviations are defined in Table4.

Figure 4 .
Figure 4.The correlograms of the Moran's I in 19 potential driving factors at the aggregation level 1 × 1. a): Correlograms with Moran's I in the distances of factors and elevation; b): Correlograms with Moran's I in seismic intensity zones; c): Correlograms with Moran's I in slope classes; d): Correlograms with Moran's I in slope aspects.The points in black are significant (P > 0.001).The driving factor abbreviations are defined in Table4.
Figure 4.The correlograms of the Moran's I in 19 potential driving factors at the aggregation level 1 × 1. a): Correlograms with Moran's I in the distances of factors and elevation; b): Correlograms with Moran's I in seismic intensity zones; c): Correlograms with Moran's I in slope classes; d): Correlograms with Moran's I in slope aspects.The points in black are significant (P > 0.001).The driving factor abbreviations are defined in Table4.

Figure 5 .
Figure 5. Cluster (a) and significance (b) maps of local indicators of spatial association (LISA) for all damaged vegetation at the 1 × 1 aggregation level.

Figure 4 .
Figure 4.The correlograms of the Moran's I in 19 potential driving factors at the aggregation level 1 × 1. (a): Correlograms with Moran's I in the distances of factors and elevation; (b): Correlograms with Moran's I in seismic intensity zones; (c): Correlograms with Moran's I in slope classes; (d): Correlograms with Moran's I in slope aspects.The points in black are significant (P > 0.001).The driving factor abbreviations are defined in Table4.

Figure 4 .
Figure 4.The correlograms of the Moran's I in 19 potential driving factors at the aggregation level 1 × 1. a): Correlograms with Moran's I in the distances of factors and elevation; b): Correlograms with Moran's I in seismic intensity zones; c): Correlograms with Moran's I in slope classes; d): Correlograms with Moran's I in slope aspects.The points in black are significant (P > 0.001).The driving factor abbreviations are defined in Table4.
Figure 4.The correlograms of the Moran's I in 19 potential driving factors at the aggregation level 1 × 1. a): Correlograms with Moran's I in the distances of factors and elevation; b): Correlograms with Moran's I in seismic intensity zones; c): Correlograms with Moran's I in slope classes; d): Correlograms with Moran's I in slope aspects.The points in black are significant (P > 0.001).The driving factor abbreviations are defined in Table4.

Figure 5 .
Figure 5. Cluster (a) and significance (b) maps of local indicators of spatial association (LISA) for all damaged vegetation at the 1 × 1 aggregation level.

Figure 6 .
Figure 6.The correlograms of the Moran's I of all damaged vegetation in the residuals of three models at the aggregation level 1 × 1.The points in black are significant (P < 0.001).The points in red are insignificant (P > 0.001).OLS: the standard linear model based on the ordinary least squares; SLM: spatial lag model; and SEM: spatial error model.

Figure 6 .
Figure 6.The correlograms of the Moran's I of all damaged vegetation in the residuals of three models at the aggregation level 1 × 1.The points in black are significant (P < 0.001).The points in red are insignificant (P > 0.001).OLS: the standard linear model based on the ordinary least squares; SLM: spatial lag model; and SEM: spatial error model.

Table 1 .
The areas of different damaged vegetation types in nine worst-hit cities and counties (km 2 ): Among all the vegetation types, we defined economic forests as the plantation that produces fresh or dry fruits, edible oils, beverages, spices, industrial raw materials, and medicinal materials.

Table 3 .
The CBERS-02B CCD images used for the interpretation of geo-hazards distribution.

Table 4 .
The damaged vegetation types and potential drivers used in the statistical analysis in nine severely damaged cities and counties.

Table 6 .
Moran's I of the residuals in the first lag at different aggregation levels.

Table 7 .
The model parameters of three different models for all damaged vegetation at aggregation level 1 × 1.

Table 8 .
The grids and percentages of the residuals of the standard linear model switched to the spatial error model at aggregation level 1 × 1.

Table 8 .
The grids and percentages of the residuals of the standard linear model switched to the spatial error model at aggregation level 1 × 1.

Table 9 .
A summary of the spatial error model for all vegetation types at four aggregation levels.