Spatiotemporal Evolution and the Influencing Factors of Tourism-Based Social-Ecological System Vulnerability in the Three Gorges Reservoir Area, China

With the rapid development of global tourism, identifying the vulnerability of tourism-based social-ecological systems (SESs) has become an important topic in sustainable development research. This study aimed to quantitatively evaluate the vulnerability, spatiotemporal evolution characteristics, and influencing factors of tourism-based SESs in the counties of the Three Gorges Reservoir Area (TGRA). A comprehensive evaluation system containing 46 indicators was constructed using a model that combines a social–economic–ecological model and a pressure–state–response model (SEE-PSR). The entropy and composite index methods were used to calculate the vulnerability values of the indicators, and Geodetector was used to explore the factors influencing system vulnerability in the whole study area. The results showed the following: (1) The mean value of the composite vulnerability index of the TGRA from 2010 to 2018 was 0.4849, indicating a moderate vulnerability state. The system vulnerability of the study area gradually decreased from moderately high to moderately low. (2) There were obvious differences in the spatiotemporal evolution of vulnerability in different counties; high and moderately high vulnerability continued to decrease, moderately low and low vulnerability increased, and moderate vulnerability showed a trend of increasing and then decreasing. Meanwhile, the relative differences in vulnerability among counties were small but gradually increasing. (3) System vulnerability was mainly caused by the social subsystem. Six factors, including the growth rate of the number of tourists and the amount of fiscal expenditure, were more likely to contribute to system vulnerability than other factors. The interaction types were mainly nonlinear enhancement types, supplemented by two-factor enhancement. This study presents an approach for evaluating the vulnerability of tourist destinations and constructing an evaluation index system. In this way, it has reference value for reducing regional vulnerability and promoting sustainable development.


Introduction
Vulnerability studies took shape in the context of natural hazard research in the 1960s [1,2]. Timmerman [3], who worked in the field of geology, first introduced the concept of vulnerability in 1981. In the 1990s, international scientific programs and agencies, such as the IHDP, IGBP, and IPCC, began to include vulnerability as an important research component [4,5]. Since then, vulnerability research has become a frequent topic in the field of global environmental change and sustainable development [6]. It has gradually been applied to many fields, such as climate change, disaster management, ecology, public health, and engineering [7]. Its connotation has gradually evolved from a singlefactor, single-dimensional concept of natural vulnerability into a composite, multifactor,

Study Area
The Yangtze River TGRA refers to the area flooded by the construction of the Three Gorges Hydropower Station, which required the relocation of many residents. With a total area of about 58,000 km 2 and a geographical range of 28°52′ N-31°72′ N, 105°82′ E-111°66′ E [26], it consists of four districts and counties in Hubei Province and 22 districts and counties in Chongqing ( Figure 1). The region has a subtropical monsoonal humid climate with four distinct seasons and average annual rainfall of 1000-1400 mm, concentrated in April-September [27]. Its average annual temperature is 17 °C-19 °C. The landscape is dominated by mountains and hills, accounting for 74% and 21.7% of the total reservoir area, respectively [28]. Soil erosion is severe in the reservoir area, and landslides, mudslides, and earthquakes occur frequently [29]. In 2018, the resident population of the TGRA was 21,025,900, with an urbanization rate of 64.42% and a per capita disposable income of USD 4119. This is slightly lower than the average per capita disposable income of USD 4266 in China. The study area is rich in tourism resources, including natural landscapes such as mountains, water, gorges, forests, and springs, as well as in cultural attractions such as traditional houses, festivals, ethnic costumes, and music in areas where ethnic minorities, such as the Shizhu and Badong, gather [30]. There are 177 A-class scenic spots, including 13 5A scenic spots, such as Dazu Rock Carvings, Three Gorges Dam, and Badong Shenlong Stream. These tourism resources have created favorable conditions for tourism development in the TGRA. In 2018, the region received 564,096,800 tourists accounting for 21% of the total number of tourists received in China and generated USD 39.58 billion, which accounted for 21% of the total tourism revenue. In recent years, tourism development in the TGRA has faced both opportunities and challenges under the national development strategies of Yangtze River Economic Belt, One Belt One Road, and All-for-One Tourism. In 2018, the region received a total of 564,096,800 tourists accounting for 21% of the total number of tourists received in China and generating a total tourism income of USD 39.58 billion, which accounted for 21% of the total tourism revenue.

Data Sources
The basic data for this study included statistical data and remote sensing data. Statistical data were mainly obtained from municipal and county-level statistical yearbooks; some data came from statistics and bulletins provided by the relevant government departments. Thirty-three indicators (e.g., natural population growth rate, urbanization rate, and tourist growth rate) were extracted from the statistical yearbooks of Chongqing and the corresponding districts and counties in Hubei Province (i.e., Chongqing Statistical Yearbook, Yichang Statistical Yearbook, Enshi Prefecture Statistical Yearbook, and China County Statistical Yearbook). Annual precipitation and underground water resources were obtained from the Chongqing Water Resources Bulletin, Enshi Prefecture Water Resources Bulletin, and Yichang Water Resources Bulletin. Data for the following were obtained by submitting "Government Information Disclosure Applications" to the relevant municipal government departments: the area of built-up areas (from the Housing and Urban-Rural Development Commissions of Chongqing, Enshi Prefecture, and Yichang), the area of planted forests (from the Forestry Bureau of Chongqing, the Forestry Bureau of Enshi Prefecture, and the Forestry and Garden Bureau of Yichang), and wastewater and waste gas emissions (from the Ecological Environment Bureaus of Chongqing, Enshi Prefecture, and Yichang). Information about A-class tourist attractions came from the following official websites: the Chongqing Culture and Tourism Development Committee (http://whlyw.cq.gov.cn/, accessed on 2 April 2021), Enshi Tujia and Miao Culture and Tourism Bureau (http://wtxgj.enshi.gov.cn/, accessed on 2 April 2021), and Yichang Culture and Tourism Bureau (http://whly.yichang.gov.cn/, accessed on 2 April 2021). Missing individual data were interpolated based on adjacent-year values. Regarding remote sensing data, Normalized Difference Vegetation Index (NDVI) data were obtained from the MODIS vegetation index product MOD13Q1 provided by NASA, with a temporal resolution of 16 d and a spatial resolution of 500 m. Since vegetation in the TGRA is the most luxuriant around August, NDVI data for August 2010, 2014, and 2018 were selected. The MOD13Q1.HDF data were concatenated using MRT (MODIS Reprojection Tool) and then projected, mosaicked, cropped, and resampled in ArcGIS.

Construction of Evaluation Index System
A combined SEE-PSR framework was adopted based on previous studies [8,31]. Combining the vulnerability concept with the actual situation of the study area, an evaluation system comprising 46 indicators was constructed based on the principles of scientific, systematic, typical, and accessible index selection (Table 1). This evaluation system aimed to comprehensively reflect the pressures faced by tourism-based SESs, the present state, and the dynamic response measures in each study period. Reflects investment in pollution control (−) 0.0204 Note: (1) "+" and "−" indicate positive and negative correlations between the indicators and system vulnerability. To a certain extent, the definition of the nature of the indicators, of both its positive and negative effects, is subjective and relative. For example, the tourism growth elasticity coefficient, which is the most weighted indicator in this study, reflects the degree to which tourism income growth rate contributes to GDP growth rate. The higher the value is, the greater the tourism economy contributes to the GDP growth in the study area. However, what is emphasized in this study is that the instability of the tourism economy will not benefit from the steady growth of the GDP in the study area. In addition, it is also emphasized in this article that there is the possibility of landslides and floods under the circumstance of precipitation in the TGRA. (2) Tourism economic density = total tourism revenue/total land area [8].
(3) Tourism growth elasticity coefficient = tourism revenue growth rate/GDP growth rate [32]. (4) The industrial structure diversification index D is calculated as D = n ∑ t=1 P t log P t , where P t is the proportion of the output of each industry, and n is the number of industries [21]. To eliminate the problem of different dimensions among the indicator data, the raw data needed to be standardized. Since positive and negative vulnerability indicators have different effects on tourism-based SESs, the range standardization method was used to deal with the positive and negative indicators, as follows: positive indicator: negative indicator: where x ij is the standardized value of the jth indicator of the ith year of a district or county, x ij is the original value of the jth indicator of the ith year of a district or county, and λ jmax and λ jmin are the maximum and minimum value of the original data of the jth indicator of all years of a district or county.

Determining the Indicator Weights
The entropy method is an objective weighting method that determines weight coefficients based on the amount of information carried by each evaluation factor. It avoids interference from subjective factors and overcomes information overlap between indicators. Therefore, this study used the entropy method to assign weights to each vulnerability index; the calculation steps and formulas are as follows: weight of the jth indicator value in ith year: indicator information entropy: information entropy redundancy: indicator weight coefficient: where m is the evaluation year, and n is the number of indicators.

Composite Vulnerability Index
The composite vulnerability index [33] was used to calculate the composite vulnerability of 26 districts and counties; it is calculated as follows: single-indicator vulnerability index: composite vulnerability index of a district and county in the ith year:

Vulnerability Classification
Since the criteria for classifying the vulnerability of tourism-based SESs are not unified, referring to the literature (e.g., on urban vulnerability [34] and ecological vulnerability [35]), this study combined the study area's vulnerability characteristics to classify the vulnerability index into five levels: S i ≤ 0.36 is low vulnerability, 0.36 < S i ≤ 0.45 is moderately low vulnerability, 0.45 < S i ≤ 0.54 is moderate vulnerability, 0.54 < S i ≤ 0.63 is moderately high vulnerability, and 0.63 < S i is high vulnerability.

Geodetector
Geodetector is a recent statistical method developed by Professor Wang et al. [36] of the Chinese Academy of Science to detect spatially stratified heterogeneity and reveal the driving forces behind it. It includes the four modules of factor detection: interaction detection, risk detection, and ecological detection. This method has been extensively utilized in the literature [37,38]. It is a software that is compiled by Excel and easy to use. Besides, it can be downloaded for free from the following website: http://www. geodetector.cn/ (accessed on 2 April 2021). The independent variables in Geodetector that are continuous need to be discretized into the type of variables for spatial partitioning. Thus, the K-means clustering algorithm was used in this study. Moreover, factor detection and interaction detection were used to analyze the intensity of influence of each influencing factor on the vulnerability of the TGRA system. The core idea of the two detectors is that the relevant factors, which has affected the vulnerability of the TGRA system, are spatially different [39]. A certain factor has significant spatial consistency with the distribution of the vulnerable TGRA system, which means that this factor is of significantly decisive significance to the spatial distribution of the vulnerability of the TGRA system.
(1) Factor detection was used to measure the influence of each influencing factor on the vulnerability of tourism-based SESs, modeled as follows: where q is the effect of the influencing factor on system vulnerability; L is the number of subregions; N h and N are the number of cells in layer h and the study area as a whole, respectively; σ 2 h and σ 2 are the variance in layer h and the study area as a whole, respectively; and q has a value range of [0, 1]-the larger the value, the stronger the effect of the influencing factor on system vulnerability, and vice versa.
(2) Interaction detection was used to detect whether the influencing factors X1 and X2 could act together to increase or decrease the effect on system vulnerability. The interactions could be divided into five categories ( Table 2). Table 2. Types of interaction between two covariates.

Graphical Representation Criterion Interaction
Sustainability 2021, 13, x FOR PEER REVIEW 10 of 19 Table 2. Types of interaction between two covariates.

Overall Vulnerability of the SESs in the Study Area
From 2010 to 2018, the mean value of the social subsystem vulnerability index Single-factor nonlinear weakening

Graphical Representation Criterion
Linear enhancement Note: ( ( 1), ( 2)) Min q X q X : gets the minimum value from ( 1), ( 2) q X q X ; ( ( 1), ( 2)) Max q X q X : gets the maximum value from ( 1), ( 2) q X q X ; ( 1)+ ( 2) q X q X : ( 1), ( 2) q X q X are summed; ( 1 2) q X X  : ( 1), ( 2) q X q X interact.  From 2010 to 2018, the mean value of the social subsystem vulnerability index showed a continuous decreasing trend (Table 3), with a decrease of 47.83%. Its maximum value was 0.2927 in Fuling in 2010; its minimum value was 0.0797 in Xingshan in 2018. The number of people receiving the subsistence security allowance and the unemployment rate in the whole study area first decreased and then increased from 2010 to 2018. Furthermore, with the development of the tourism economy, the density of tourists increased from 47,162 to 147,058 persons/km 2 , and the ratio of tourists to local residents increased from 8:1 to 30:1, placing greater pressure on the social system and the regional carrying capacity of the tourism destinations. However, with the growth of local fiscal funds from USD 15.61 billion to USD 37.08 billion; investment in transportation, education, and social security; and an increase in tertiary industry employees from 3.7947 million to 6.0284 million, pressure from internal and external sources combined with a series of response measures eventually shaped the social subsystem vulnerability in a rapid and then slowly declining process. The mean value of the economic subsystem vulnerability index also continued to decline, showing a significantly larger decline from 2010 to 2014 than from 2014 to 2018. The maximum value was 0.2362 in Shizhu in 2010, and the minimum value was 0.0916 in Wanzhou in 2014. From 2010 to 2018, the regional GDP of the whole study area increased from USD 92.93 billion to USD 255.17 billion. Fiscal revenues increased from USD 8.76 billion to USD 21.12 billion, the per capita GDP increased from USD 4590 to USD 10,108, the per capita disposable income of urban and rural residents increased by about 2.36 times, and the Engel's coefficient decreased by 17.95%. These favorable factors all helped weaken the vulnerability of the economic subsystem in the study area, as economic strength increased and living standards improved. In addition, the economic density of tourism increased from 1974 to 73,659,600 yuan/km 2 , and the proportion of total tourism revenue to GDP increased from 9.80% to 24.05%. While this highlights the important role of tourism in promoting the economic development of the TGRA, it also points to the possibility of vulnerability and volatility, especially in the period following a crisis event.

Overall
The mean value of the ecosystem vulnerability index showed a "rising and then falling" trend. The maximum value was 0.2242 in Fengdu in 2014, and the minimum was 0.0704 in Nanan in 2018. From 2010 to 2018, the total retail sales of consumer goods and total wastewater discharge increased, going up by 138.88% and 54.44%, respectively. Thus, with the socioeconomic development of the area's tourist destinations, resource consumption and pollution increased the pressure on and vulnerability of the ecosystem. However, there has also been increased awareness of energy conservation and environmental protection among governments, residents, tourists, and factories in the tourist areas. This is seen in the increasing rate of urban domestic sewage treatment (from 85.59% to 92.02%). Moreover, waste gas emissions showed a trend of rising then falling-from 1.2302 million tons in 2010 to 2.347 million tons in 2014, and then down to 1.2378 million tons in 2018. The trajectory of the good air quality rate is the opposite of that of total emissions; in fact, the good air quality rate in 2018 was 2.01% higher than in 2010.

Composite Vulnerability of SESs
The mean value of the composite vulnerability index of the SESs decreased from 2010 to 2018 by 42.37%. The evolution of its vulnerability level was in the order of moderately high vulnerability-moderate vulnerability-moderately low vulnerability. The maximum value was 0.6468 in Badong in 2010, and the minimum was 0.3219 in Xingshan in 2018. The situations for subsystems with the highest vulnerability varied from year to year. Among them, the social subsystem progressed the fastest, changing from the most vulnerable subsystem to the least vulnerable one, with the largest coefficient of variation (0.3381). Meanwhile, the ecological subsystem had the smallest coefficient of variation, indicating that it was more stable than the other subsystems. From 2014 to 2018, the economic subsystem showed a small decrease in vulnerability, thus leading it to becoming the most vulnerable subsystem in 2018. The mean value of the multiyear composite vulnerability index in the study area was 0.4849, which is moderately vulnerable.

Spatiotemporal Evolution Characteristics of the Vulnerability Types of the SESs
From 2010 to 2014, there was an overall decreasing trend in the vulnerability of the area's districts and counties ( Figure 2). Spatially, the dominant vulnerability type shifted from moderately high to moderate vulnerability, with moderately low and moderate vulnerability replacing high vulnerability. In 2010, there were 4 districts and counties with high vulnerability, 15 with moderately high vulnerability, and 7 with moderate vulnerability. In 2014, there were 4 districts and counties with moderately high vulnerability, 14 with moderate vulnerability, and 8 with moderately low vulnerability. Wulong had the largest decreases in vulnerability, and Badong, Fuling, and Dadukou decreased from high vulnerability to moderately low vulnerability as well. However, Changshou and Yubei rose from moderate vulnerability to moderately high vulnerability. Zigui, Yiling, Banan, Yuzhong, Shapingba, and Jiangbei had no differences in level, showing moderately high and moderate vulnerability. Other districts and counties decreased by 1-2 vulnerability levels.   The overall vulnerability level of the districts and counties also showed a decreasing trend from 2014 to 2018. Spatially, the dominant vulnerability type shifted from moderate to moderately low vulnerability, with low and moderately low vulnerability replacing moderately high vulnerability. There were 4 districts and counties with moderately high vulnerability, 14 with moderate vulnerability, and 8 with moderately low vulnerability in 2014. Meanwhile, there were 6 districts and counties with moderate vulnerability, 15 with moderately low vulnerability, and 5 with low vulnerability in 2018. Zigui, Xingshan, Kaizhou, Zhong, Changshou, Yubei, Banan, and Jiangjin all decreased by two levels, containing types ranging from low to moderately high vulnerability. Wanzhou, Shizhu, and Nanan rose from moderately low to moderate vulnerability. Wushan, Fengjie, Yuzhong, Shapingba, Dadukou, and Jiangbei had no differences in level, showing moderately low and moderate vulnerability. Other districts and counties decreased by one vulnerability level each. The highest vulnerability value was reduced from 0.5779 in Yubei to 0.5013 in Wanzhou, and the lowest was reduced from 0.3975 in Dadukou to 0.3219 in Xingshan.
In addition, based on the characteristics of the numerical changes in the vulnerability of each district and county from 2010 to 2018, they can be divided into gradually declining, stable, and fluctuating. First, the districts and counties belonging to the gradually declining type were Xingshan, Yiling, Badong, Wushan, Wuxi, Fengjie, Yunyang, Kaizhou, Zhong, Fuling, Fengdu, Wulong, Banan, Jiangjin, Beibei, Jiulongpo, and Dadukou. These 17 districts and counties had an average drop of 36.10%, among which Zhong and The overall vulnerability level of the districts and counties also showed a decreasing trend from 2014 to 2018. Spatially, the dominant vulnerability type shifted from moderate to moderately low vulnerability, with low and moderately low vulnerability replacing moderately high vulnerability. There were 4 districts and counties with moderately high vulnerability, 14 with moderate vulnerability, and 8 with moderately low vulnerability in 2014. Meanwhile, there were 6 districts and counties with moderate vulnerability, 15 with moderately low vulnerability, and 5 with low vulnerability in 2018. Zigui, Xingshan, Kaizhou, Zhong, Changshou, Yubei, Banan, and Jiangjin all decreased by two levels, containing types ranging from low to moderately high vulnerability. Wanzhou, Shizhu, and Nanan rose from moderately low to moderate vulnerability. Wushan, Fengjie, Yuzhong, Shapingba, Dadukou, and Jiangbei had no differences in level, showing moderately low and moderate vulnerability. Other districts and counties decreased by one vulnerability level each. The highest vulnerability value was reduced from 0.5779 in Yubei to 0.5013 in Wanzhou, and the lowest was reduced from 0.3975 in Dadukou to 0.3219 in Xingshan.
In addition, based on the characteristics of the numerical changes in the vulnerability of each district and county from 2010 to 2018, they can be divided into gradually declining, stable, and fluctuating. First, the districts and counties belonging to the gradually declining type were Xingshan, Yiling, Badong, Wushan, Wuxi, Fengjie, Yunyang, Kaizhou, Zhong, Fuling, Fengdu, Wulong, Banan, Jiangjin, Beibei, Jiulongpo, and Dadukou. These 17 districts and counties had an average drop of 36.10%, among which Zhong and Yunyang had the largest and smallest drops of 46.47% and 18.96%, respectively. Second, the districts and counties belonging to the stable type were Yuzhong, Shapingba, and Jiangbei; the coefficients of variation of their composite vulnerability indexes were 0.0233, 0.0805, and 0.0329, respectively. This shows that these three districts had better stability. Finally, six districts and counties belonged to the fluctuating type, among which Wanzhou, Shizhu, and Nanan showed a decreasing and then increasing change process, while Zigui, Changshou, and Yubei showed the opposite process.

Spatial Variation Characteristics of the Vulnerability of SESs among Counties
Calculating the mean value of the composite vulnerability of each district and county during the study period (Figure 4), it can be seen that Badong (0.5145) and Banan (0.5034) were in the top two positions while Yiling (0.4704) and Wulong (0.4660) were at the bottom. Similarly, the standard deviation and coefficient of variation were calculated to reflect the absolute and relative differences among the districts and counties during the study period ( Figure 4). As shown in Figure 5, the standard deviation and coefficient of variation of Yuzhong were the smallest, with values of 0.0109 and 0.0233, indicating the smallest dispersion and smallest variation of vulnerability. However, the standard deviation and coefficient of variation of Wulong were the largest, with values of 0.1518 and 0.3258. The coefficients of variation for the whole study area were 0.0898, 0.1052, and 0.1230 in 2010, 2014, and 2018, respectively. This shows that the relative differences in vulnerability among counties were small but gradually increasing. The mean values of the standard deviation and coefficient of variation of the composite vulnerability index for all districts and counties over the study period were 0.0956 and 0.1976, respectively. It is noteworthy that the mean value of the composite vulnerability of Wulong was the smallest among all districts and counties, but its standard deviation and coefficient of variation were the largest. This difference is consistent with the actual situation in which Wulong has risen and developed rapidly as a result of tourism. Specifically, in 2007, the Southern China Karst (with the Wulong Karst as one of its three components) applied to become a World Natural Heritage Site, and, in 2011, Wulong was reconfirmed as a 5A tourist attraction. Since then, Wulong's popularity has increased dramatically, attracting a large number of tourists and thus leading to the high vulnerability of its SESs. Later, Wulong gradually decreased to low vulnerability under the combined effects of a series of response measures.

Main Influencing Factors
Based on the q values of each influencing factor calculated by Geodetector, the top 10 indicators in 2010, 2014, and 2018 were selected as the main factors affecting the vulnerability of SESs in the counties of the TGRA in descending order ( Table 4). The top 10 influencing factors from 2010 to 2018 had numerical evolutions of the vulnerability indicators within the social, economic, and ecological subsystems of 4-3-7, 3-3-3, and 3-4-0, respectively. This indicates that social subsystem vulnerability was the main cause of SES vulnerability. In addition, the numerical evolutions of the vulnerability indicators within the pressure, state, and response layers were 4-3-4, 2-1-2, and 4-6-4, respectively, indicating that system vulnerability was caused by a lack of power in the response layer and an overload in the pressure layer. Notably, the growth rates of tourist number (V3), fiscal expenditure amount (V12), population density (V33), accommodation and catering business volume (V31), industrial structure diversification index (V32), and forest coverage rate (V39) all appeared twice, indicating that they were more likely to contribute to system vulnerability than other factors. Therefore, adopting relevant response measures will help reduce the system vulnerability of the study area, such as improving tourism reception capacity, increasing financial investment (e.g., to improve livelihoods and public service facilities), improving the star level and service quality of accommodations and catering, diversifying the industrial structure, continuing to increase the forest coverage rate, and moderately controlling the population density of tourist destinations.

Interactions between Factors
Due to the large amount of data and limited space, the overall interactions are summarized, and some data characteristics are explained. In 2010, the types of interactions between factors included 542 groups of two-factor enhancement types and 493 groups of nonlinear enhancement types, indicating that any two-factor interaction will increase the vulnerability of the study area. The top 10 groups of interaction factors, represented by green space per capita (V43) and domestic wastewater treatment rate (V45), were all nonlinear enhancement types (Table 5); their interaction values were higher than 0.8 and showed high interaction effects. When the proportion of school students to the resident population (V5), urban registered unemployment rate (V10), social security and employment expenditure to total fiscal expenditure (V14), total retail sales of consumer goods (V34), and domestic sewage treatment rate (V45) interacted with other factors, the number of groups with interaction values higher than 0.8 were six, five, five, five, and five, respectively. This indicates that they are more likely to interact than other factors. The year 2014 included 221 groups of two-factor enhancement types and 814 groups of nonlinear enhancement types. The top 10 groups of interaction factors, represented by the growth rate of the number of tourists (V3) and the ratio of the number of school students to the resident population (V5), were also nonlinear enhancement types. When the proportion of education expenditure to total fiscal expenditure (V15), domestic sewage treatment rate (V45), proportion of school students to the resident population (V5), number of rated road miles (V11), and Engel's coefficient (V21) interacted with other factors, the number of groups with interaction values higher than 0.8 were three, three, two, two, and two, respectively. The year 2018 included 400 groups of two-factor enhancement types and 635 groups of nonlinear enhancement types. The first 10 groups of interactions included three groups of two-factor enhancement types represented by the growth rate of the number of tourists (V3) and the share of social security and employment expenditures in total fiscal expenditures (V14), as well as seven groups of nonlinear enhancement types represented by the growth rate of the number of tourists (V3) and the density of tourists (V7). When the share of social security and employment expenditures in total fiscal expenditures (V14), growth rate of tourist numbers (V3), natural population growth rate (V1), and urban registered unemployment rate (V10) interacted with other factors, the number of groups with interaction values higher than 0.8 were 21, 15, 6, and 4, respectively. From 2010 to 2018, the interaction type was mainly the nonlinear enhancement type, supplemented by two-factor enhancement types, and the mean value of the interaction (0.5604, 0.4854, 0.5253) decreased and then increased. Therefore, adopting comprehensive response measures will help reduce the negative effects of increasing the interaction values. Such measures can include continuously improving the reception and service capacity of tourist destinations, enhancing their ability to prevent and cope with tourism risks, expanding the number of miles of district and county-level roads (or building direct highways to tourist attractions in districts and counties with high tourism development), increasing financial investment in education development, creating more jobs, providing vocational skills training, encouraging residents to work in the tourism service industry, controlling the natural population growth rate, increasing the incomes of local residents, and improving locals' consumption structure.

Discussion
(1) Constructing a sound evaluation index system and selecting a suitable research method are the foundations of evaluation-focused research [8,40]. Based on the vulnerability concept, previous literature, and the principles of indicator selection, this study constructed a comprehensive evaluation system consisting of 46 indicators, aiming to reflect the actual vulnerability status of the study area as much as possible. This study then calculated the vulnerability values of each system and its subsystems using the entropy method, which is an objective assignment method and composite index method. Despite the rationality of this research process, some limitations exist. These include a lack of a deep understanding of vulnerability; moreover, the rationality of the selection of indicators and methods need to be justified and tested.
Considering that some studies [41,42] have used multiple weighting methods at the same time, it is proposed that, if a certain vulnerability evaluation supports impor-tant decision-making or special needs, a multimethod, multi-indicator, or multiscale approach could be used to ensure the scientificity and feasibility of the results. (2) There are few studies on vulnerability in the TGRA. Research has mainly been conducted from a single perspective, such as farmers' livelihoods, agriculture, the drawdown zone environment, landslide hazard risk, or ecological and economic vulnerability. Ma [43], for example, evaluated the ecological vulnerability of the TGRA from 2001 to 2010. Comparing the present study's findings with that study's spatial distribution map of ecological vulnerability in 2010, we found that the two maps basically matched, albeit with some differences. The differences could be due to the fact that Ma selected more indicators with "natural attributes" while the present study selected more indicators under the influence of human society according to the SEE-PSR model. Thus, the vulnerability of Wulong, Wuxi, Fengjie, and Wushan, which are mainly mountainous areas, is higher in the present study. From 2010 to 2014, the mean system vulnerability and economic subsystem vulnerability of the TGRA showed a decreasing trend, which is consistent with Liu et al. [40]. The difference, however, is that the SEE vulnerability in Chongqing in Liu's study showed a spatial distribution feature of high in the west and low in the northeast and southeast. In the present study, however, no specific distribution feature was formed. The reason for this difference could be because Liu et al. assigned higher weights to the ecological subsystems, such that the northeastern and southeastern regions of Chongqing-which have better natural conditions, higher vegetation cover, and less human activity-had lower vulnerability while the more developed western regions had higher vulnerability. In addition, the higher comprehensiveness of the evaluation index system and the smaller relative differences in weights among the subsystems may be important reasons for the insignificant spatial distribution characteristics in this study. (3) SESs vulnerability in tourism destinations deserves more attention and in-depth study.
With the overall positive development of tourism worldwide, there are increasing flows of people, materials, capital, and information between tourism-based SESs and the external environment. This positively affects local socioeconomic development, facilitates cultural exchange, increases employment opportunities, and alleviates poverty. However, excessive tourism development and poor resource protection have exacerbated conflicts among populations, environments, resources, and development. Therefore, enhancing tourism destinations' resilience to stresses and risks while reducing the vulnerability of their systems will contribute to regional sustainable development. Given the complexity and dynamics of tourism-based SESs [8,16], it is difficult to accurately measure their vulnerability. Evaluation research on their vulnerability is still in the process of exploration, validation, and revision with regard to evaluation index system construction, research method selection, concept cognition, threshold determination, and influencing factors. Additional work in this area will therefore deepen the theory and methods of vulnerability research. (4) The Geodetector model has both advantages and disadvantages. The model is good at detecting spatial heterogeneity and revealing the driving forces behind it; thus, it was employed to analyze the factors affecting the vulnerability of SESs in the TGRA as a whole. However, its properties also hindered analyzing the factors that affected vulnerability in each district and county. Therefore, future work can explore and simulate the factors affecting the vulnerability of each district and county, the influencing mechanism of the vulnerability of tourism-based SESs, and future vulnerability scenarios in each area.

Conclusions
Adopting an SEE-PSR model, this study constructed an evaluation system comprising 46 indicators to analyze the spatiotemporal evolutionary characteristics and influencing factors of the vulnerability of SESs in the TGRA using the composite index method and Geodetector. The main conclusions are as follows: (1) From 2010 to 2018, the mean value of the composite vulnerability index of the TGRA was 0.4849, which is a moderate vulnerability state. (2) The system vulnerability of the study area gradually decreased from moderately high vulnerability to moderately low vulnerability, and the vulnerability of socioeconomic subsystems showed a continuous decreasing trend while the ecological subsystem shows a trend of rising then falling. (3) The spatiotemporal evolution of vulnerability in counties showed obvious differences. In terms of the vulnerability types of districts and counties, high vulnerability and moderately high vulnerability types kept decreasing, and moderately low vulnerability and low vulnerability types kept increasing, while moderate vulnerability showed a trend of increasing and then decreasing. In addition, changes in the system vulnerability of districts and counties could be divided into gradually decreasing, stable, and fluctuating. Spatially, the relative differences in vulnerability among counties were small but gradually increasing. (4) The factor detection results showed that system vulnerability was mainly caused by the social subsystem but also by a lack of power in the response layer and an excessive load in the pressure layer. Six factors, including the growth rate of the number of tourists and the amount of fiscal expenditure, were more likely to contribute to system vulnerability than other factors. (5) The interaction detection results showed that the interaction type was mainly the nonlinear enhancement type, supplemented by the two-factor enhancement type, and the mean value of interaction showed a trend of decreasing and then increasing.
This study focused on the vulnerability of special tourism-based social-ecological systems in the context of rapid global tourism development. It exploratively constructed a tourism-based vulnerability evaluation index system and used it for an empirical investigation. In this way, this study may enrich theoretical approaches and practical cases in vulnerability research. In addition, this study clarified the extent, spatiotemporal evolution characteristics, and main influencing factors of the vulnerability of SESs in the TGRA. This study suggests that the response measures intended to reduce regional vulnerability and enhance the ability of tourism destinations to resist pressure and risk will contribute to regional sustainable development.

Data Availability Statement:
The data used in this study are available from the corresponding author upon request. The data are not publicly available due to privacy restrictions.