Assessment of Ecosystem Services: Spatio-Temporal Analysis and the Spatial Response of Influencing Factors in Hainan Province

Abstract

The impact of human activities on ecosystems is receiving increasing attention because their mechanisms of action are complex; the spatial response of ecosystem service drivers still needs to be explored further. This study evaluated three ecosystem services—water yield, soil conservation, and carbon storage—in Hainan Province from 2000 to 2020; we analyzed the spatial and temporal changes of the ecosystem services, and the spatial heterogeneity of the influencing factors. The results were as follows: (1) The average water yield, soil conservation and carbon storage of Hainan Province from 2000 to 2020 were 42.36 billion,$8.01\times {10}^8 \mathrm{t}$ and 1.52 ×${10}^7$ t, respectively. Overall, the ecosystem services were relatively weak at lower elevations. (2) There were obvious hot spots and cold spots in the water yield and soil conservation, and the hot spot distribution of carbon storage was not obvious. (3) There were differences between the ecosystem services for different land use types; trade-off relationships only appeared between unused land and ecosystem services. (4) The precipitation, normalized difference vegetation index and elevation factors had great impacts on the ecosystem services. Most of the human activity factors showed a significant nonlinear enhancement effect during their interaction. Population and elevation had obvious spatial differentiation effects on the water yield and carbon storage services.

1. Introduction

‘Ecosystem services’ refer to the direct or indirect benefits that human beings receive from ecosystems to meet their own needs [1]. As the basis of sustainable development, the concept of ecosystem services has attracted extensive attention from scholars [2]. The Millennium Ecosystem Assessment showed that about 60% of global ecosystem services exhibited a trend of degradation [3]. With the continuous advancement of urbanization, human activities are becoming increasingly widespread; under the background of sustainable development, it is particularly important to evaluate the service functions and influencing factors of ecosystems, so as to guide the government when allocating ecological engineering.

In order to evaluate ecosystem services, researchers first need to accurately quantify the services provided by ecosystems. There are three evaluation methods: the value evaluation method [4], the material quality evaluation method [5,6,7], and the energy evaluation method [8]. The value evaluation method is greatly affected by subjective factors and fails to realize spatialization, while the energy evaluation method is not often applied because of its complex calculation process [9]. With the in-depth understanding of ecosystem service theory and the development of information technology, composite “modeling and spatialization” models have been developed on the basis of the material quality evaluation method [10]. These composite models include the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) Model, the Aries Model [11], and the MIME Model [12]. Compared to other composite models, the InVEST model is more mature and widely used. The InVEST model has the characteristics of simple data acquisition, quantitative evaluation, and spatial results [13]. Therefore, this study selected the InVEST model to evaluate the ecosystem services of Hainan Province. This model can overcome the defect that the value evaluation method cannot explain the process mechanism, and its results are more scientific. Scholars have used the InVEST model to study the Yeerqiang River Basin [14], the Yellow River Basin [15], the Sanjiangyuan Basin [16], the Maine area [17], and the Shanxi coalfield [18]. This research has focused on the water yield [15,19], carbon storage [20], soil conservation [21], habitat quality [22], and biodiversity [23]. To sum up, scholars have mainly focused on arid and semi-arid watershed areas, while relatively little research has been conducted on humid areas. Whether the achievements made in arid areas can be applied to humid areas remains to be tested; there are relatively few studies on island ecosystem services. Therefore, this study focused on Hainan Province, a typical humid area. The InVEST model was used to study ecosystem services in order to refine case studies in wet areas. At the same time, it can provide reference for island ecosystem services in the same latitude of other countries.

Ecosystem services and their changes are affected by many driving factors [24]. The influencing factors of ecosystem services will vary with different research scales. For example, on the global scale, climate change is the main controlling factor leading to the change of ecosystem services [25]. On the regional scale, natural factors such as topography, vegetation cover, and soil conditions will affect ecosystem services [26]. Human factors also affect ecosystem services, such as land use change, policy change, population change, and human activities [27]. Exploring ecosystem service responses to various driving factors is critical to achieving ecosystem sustainability. Therefore, the question of how to analyze the spatial heterogeneity of ecosystem services has attracted the attention of many scholars. A previous study used correlation analysis to investigate the natural social driving factors of ecosystem services [28]. In another study, the analytic hierarchy process was used to reveal the impact of land use and socioeconomic factors on ecosystem services [29]. In addition, the relationship between the climate, human activities, and runoff was analyzed using the GAM model in a previous work [30]. To sum up, most scholars conduct relevant research on the influencing factors of ecosystem services from the overall perspective, and rarely analyze the spatial response of driving factors to ecosystem services from the perspective of different spatial scales. In order to address this knowledge gap, this study used geographic detectors to identify the main driving factors of ecosystem services, and then used multi-scale geographically weighted regression (MGWR) to analyze the impacts of the main factors on ecosystem services at different spatial scales, so as to reveal the spatial responses of ecosystem driving factors in Hainan Province.

Hainan Province established an international tourism island in 2009, and the Chinese central government proposed to establish the China (Hainan) Pilot Free Trade Zone in Hainan in 2018. The increase of tourism and the rapid economic development caused by the construction of the free trade port are major characteristics of Hainan Province [31]. However, once the island’s ecosystem is damaged, it becomes difficult to repair the ecological environment, and how to coordinate the growing economy and maintain the environmental quality of Hainan Province is an urgent problem. Based on the actual situation of Hainan Province, this paper used the InVEST model to evaluate its ecosystem services, and the Moran index and G* to analyze the spatial distribution characteristics of ecosystem services. We examined the trade-off relationship between ecosystem services of different land use types through correlation analysis, and explored the spatial responses of the ecosystem service driving factors based on the geographical detectors and the MGWR model. The objectives of this paper were (1) to analyze the spatio-temporal changes of ecosystem services in Hainan Province from 2000 to 2020, (2) to identify the location of clustering in the space of ecosystem services and conduct hot spot analysis, (3) to analyze the relationships between the ecosystem services of various land use types, and (4) to analyze the spatial response mechanisms of factors affecting ecosystem services. Achieving these objectives could reveal new directions for exploring spatial response mechanisms, and could provide a scientific reference for Hainan province in order to construct a national ecological civilization pilot zone.

2. Materials and Methods

2.1. Study Area

Hainan Province (18.80°–20.10° N, 108.37°–111.03° E) is located at the southernmost part of China, with a land area of about 35,300 ${\mathrm{km}}^2$. It is dominated by plains and platforms, exhibiting the geomorphic characteristics of high elevations in central Hainan Island and low elevations at the coast, with an average altitude of 120 m, as shown in Figure 1. Hainan Island has a typical tropical monsoon marine climate, with an annual precipitation of 1000–2600 mm [32]. It has an obvious rainy season and a dry season, and the annual average temperature is 23.8–26.2 °C. There are three major rivers on Hainan Island: Nandu River, Changhua River, and Wanquan River. Together, their drainage areas account for about 47% of the total of Hainan Island. The forest coverage rate of Hainan Province is as high as 62.1%, at the forefront of the country. With 1944 km of coastline, Hainan Island is rich in tourism resources. In order to promote the construction of a national ecological civilization pilot area, the red line area of ecological protection will not be less than 27.3% of the land area of the whole province.

2.2. Data Sources

The data used this study are shown in detail in

Table 1. Data sources.

Type	Resolution	Sources
Landuse	30 m	Resource science and satellite center (http://www.resdc.cn/, accessed on 6 May 2021)
DEM	30 m	Geospatial data cloud (http://www.gscloud.cn/, accessed on 6 May 2021)
Precipitation, soil data	1 km	National Earth System Science Data Center (http://www.geodata.cn/, accessed on 6 May 2021)
Potential evapotranspiration data	1 km	NASA (https://modis.gsfc.nasa.gov/, accessed on 6 May 2021)
Hainan statistical yearbook	/	Hainan Provincial People’s Government (https://www.hainan.gov.cn/, accessed on 6 May 2021)

. The annual potential evapotranspiration data were extrapolated from the eight-day synthetic product data. After preprocessing, all of the data were uniformly projected and resampled to a resolution of 1 km.

2.3. Research Methods

2.3.1. Water Yield

The InVEST model is a comprehensive evaluation model of ecosystem services and trade-offs. This model is based on the principle of water balance, combined with climatic factors and soil physical and chemical properties in order to calculate water yield. The water yield is the precipitation of each pixel minus the water after actual evaporation. The model parameters are set with reference to relevant research [33], and the results are verified with reference to the data of Hainan Water Conservancy Yearbook. The calculation formula is as follows:

$$Y_{xj}=\left(1-\frac{AE{T}_{xj}}{P_x}\right)\cdot P_x$$

(1)

where $Y_{xj}$ is the water yield of the x-th pixel of j ecosystem type (mm), $AE{T}_{xj}$ is the annual actual evapotranspiration of the x-th pixel of j ecosystem type (mm), and $P_x$ is the annual average precipitation of the x-th pixel (mm).

$$\frac{AE{T}_{xj}}{P_x}=1+\frac{PE{T}_{xj}}{P_x}-{\left[1+{\left(\frac{PE{T}_{xj}}{P_x}\right)}^{\omega }\right]}^{\frac{1}{\omega }},$$

(2)

where $PE{T}_{xj}$ is the annual potential evapotranspiration of the x-th pixel of j ecosystem type (mm), and $\omega$ is the non-physical parameters of natural climate soil properties.

$$\omega =Z\frac{AWC\left(x\right)}{P\left(x\right)}+1.25$$

(3)

where $AWC\left(x\right) \mathrm{is} \mathrm{the}$available water content of vegetation (mm); its value is determined by the soil texture and effective soil depth. Z is the empirical constant. These parameters were derived from Table 1.

2.3.2. Soil Conservation

The calculation principle of soil conservation for the InVEST model is based on the revised universal soil loss equation. Its calculation is divided into two steps: one is the sediment accumulation of confluence, and the other is the amount of erosion reduced after surface interception and conservation measures [34]. The calculation process is as follows:

$$\mathrm{SEDRET}=R\times K\times LS\times \left(1-C\times P\right)+SEDR$$

(4)

where SEDRET is the amount of soil retention$\left(t\right)$, and R is the rainfall erosivity factor ($\mathrm{MJ}\cdot \mathrm{mm}{(\mathrm{ha}\cdot \mathrm{h})}^{-1}$), which is calculated by using the Wischmeier formula [35] based on the annual precipitation. K is the soil erodibility factor ($\mathrm{t}\cdot \mathrm{ha}\cdot \mathrm{h}{(\mathrm{MJ}\cdot \mathrm{ha}\cdot \mathrm{mm})}^{-1}$), which is calculated by the erosion–productivity impact calculator model in ArcGIS; LS is the slope length factor, and is dimensionless; C is the vegetation cover and crop management factor, and is dimensionless, with a range of 0–1. P is factors of the water and soil conservation factor, with a range of 0–1, and $SEDR$ is the sediment retention (t).

$$K=\left\{0.2+0.3exp\left[-0.0256\times \mathrm{sand}\times \left(1-\frac{\mathrm{silt}}{100}\right)\right]\times {\left(\frac{ \mathrm{silt} }{ \mathrm{clay} + \mathrm{silt} }\right)}^{0.3}\times \left[1-\frac{0.25C}{C+\mathit{exp}\left(3.72-2.95C\right)}\right]\times \left[1-\frac{0.7sn}{sn+\mathit{exp}\left(22.9sn-5.51\right)}\right]\right\}$$

(5)

$$sn=1-\frac{\mathrm{sand}}{100}$$

(6)

where sand, silt, clay and C are the percentage contents of soil sand, silt, clay and organic carbon, respectively (%).

In the process of soil erosion, intercepted by the attachments on the underlying surface, not all eroded soil can be transported to the catchment area by surface runoff; that is, the amount of erosion is reduced after surface interception and maintenance measures are taken. The calculation formula is as follows:

$$SEDR=S{E}_x{\displaystyle \sum }_{y=1}^{x-1}USL{E}_y{\displaystyle \prod }_{z=y+1}^{x-1}\left(1-S{E}_z\right)$$

(7)

where $S{E}_x$is the sediment retention of the x-th pixel, $USL{E}_y$ is the amount of actual sediment of the y-th pixel (t), and $S{E}_z$ is the sediment retention on the upslope of the z-th pixel.

2.3.3. Carbon Storage

Carbon storage is an important aspect of terrestrial ecosystems; it provides services and plays an important role in the global carbon balance and global climate regulation. Carbon storage is calculated by adding the carbon reserves of different periods and different land types in the five carbon pools [20]. The carbon storage in terrestrial ecosystem is generally divided into four basic carbon pools, which are aboveground carbon(all living vegetation biomass on the soil), belowground carbon (the living root system of plants), soil carbon (mineral soil organic carbon and organic soil carbon) and dead organic carbon (the amount of carbon stored in litter, fallen trees and dead trees) [9]. As the fifth carbon pool represents specific products and will not enter the carbon in the atmosphere, the fifth carbon pool is not considered in this paper. The calculation formula is as follows:

$$C=C_{\mathrm{above}}+C_{\mathrm{below}}+C_{\mathrm{soil}}+C_{\mathrm{dead}}$$

(8)

where C is the total carbon storage for ecosystems$, C_{\mathrm{above} }$ is the aboveground carbon storage$, C_{\mathrm{below} }$ is the belowground carbon storage, $C_{\mathrm{soil} }$ is the soil carbon storage, and $C_{\mathrm{dead} }$ is the carbon storage of dead organic matter; the unit is (t).

2.3.4. Identification of Ecosystem Service Hot spots

The cold spot or hot spot areas of the ecosystem service water yield, soil conservation and carbon storage functions are given by the degree of aggregation in spatial distribution, and the ecosystem service in hot spots has a high value [36]. The Z value of Getis-Ord G* statistics helps us to identify the differences between cold spots and hot spots with different confidence levels of 99%, 95%, and 90% [37]. In this study, the spatial autocorrelation Moran’s I index is used to analyze whether there is spatial aggregation of the water yield, soil conservation and carbon storage, and then the hot spot analysis is used to analyze the hot spots of ecosystem services in Hainan Province. The calculation formula is as follows:

$$G_i^{*}=\frac{{\sum }_j^n{W}_{ij}{x}_j}{{\sum }_j^n{x}_j}$$

(9)

$$Z\left(G_i^{*}\right)=\frac{G_i^{*}-E\left(G_i^{*}\right)}{\sqrt{VAR\left(G_i^{*}\right)}}=\frac{{\sum }_j{W}_{ij}-\overline{x}{\sum }_j^n{W}_{ij}}{\sqrt{\frac{n{\sum }_j^n{W}_{ij}^2-{\left({\sum }_j^n{W}_{ij}\right)}^2}{n-1}}}$$

(10)

$$\overline{x}=\frac{{\sum }_j^n{x}_j}{n}$$

(11)

where $W_{ij}$ is the weight matrix used to represent patches i and j in space, $x_j$ is the attribute value of j, n is the total number of patches, $\overline{x}$ is the average of $x_j$, and the range of $Z\left(G_i^{*}\right)$ is −2.58 and 2.58.

2.3.5. Ecosystem Service Trade-Off and Synergy Analysis

Correlation analysis is usually used to study the correlation between variables [38]. This study used this method to study the relationships between the water yield, soil conservation, and carbon storage of different land use types. The magnitude of the value indicates the strength of the correlation. A positive value indicates that one ecosystem service would be enhanced with the enhancement of another service, exhibiting a collaborative relationship. The calculation formula is as follows:

$$R_{ab}=\frac{{\sum }_{i=1}^n\left(ES{1}_i-\overline{ES1}\right)\left(ES{2}_i-\overline{ES2}\right)}{\sqrt{{\sum }_{i=1}^n{\left(ES{1}_i-\overline{ES1}\right)}^2{\sum }_{i=1}^n{\left(ES{2}_i-\overline{ES2}\right)}^2}}$$

(12)

where$R_{ab}$is the correlation coefficient of services for two ecosystems, $ES{1}_i$ and $ES{2}_i$ represent two ecosystems value of the i-th pixel, $\overline{ES1} \mathrm{and} \overline{ES2}$ are the mean values of the two ecosystems, and the number of samples is n.

2.3.6. Driving Force Analysis Based on Geo-Detectors

Geographical phenomena have spatial heterogeneity. Geographical detectors can reveal the driving forces behind geographical phenomena. There is a spatial similarity between a dependent variable and its independent variable [39]. The driving mechanism of ecosystem service change is comprehensive, which is the result of the joint action of natural factors and economic and social factors [28]. The principle of factor detection is as follows:

$$\begin{array}{c}q=1-\frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}\end{array}$$

(13)

where $N_h$ and N are the number of layers (h) and units in the whole area; ${\sigma }_h^2$ and ${\sigma }^2$ are the variance of layer h and the whole region, respectively; q is the factor interpretation degree; and the value range is [0, 1]. The larger the value of q, the more obvious the influence of the independent variable on the dependent variable.

Interaction detection refers to the interpretation strength of dependent variables under the joint action of two factors. The interaction relationships were divided into five categories (

Table 2. Classification of the interaction relationships.

Judgment Basis	Interaction Type
$q\left(\mathrm{X}1{\displaystyle \cap }\mathrm{X}2\right)<Min\left(q\left(\mathrm{X}1\right),q\left(\mathrm{X}2\right)\right)$	Nonlinear weakening
${\left(q\left(\mathrm{X}1\right),q\left(\mathrm{X}2\right)\right)}_{min}$$<q\left(\mathrm{X}1{\displaystyle \cap }\mathrm{X}2\right)$$<{\left(q\left(\mathrm{X}1\right),q\left(\mathrm{X}2\right)\right)}_{Max}$	Single-factor nonlinear weakening
$q\left(\mathrm{X}1{\displaystyle \cap }\mathrm{X}2\right)>Max\left(q\left(\mathrm{X}1\right),q\left(\mathrm{X}2\right)\right)$	Two-factor enhancement
$q\left(\mathrm{X}1{\displaystyle \cap }\mathrm{X}2\right)=q\left(\mathrm{X}1\right)+q\left(\mathrm{X}2\right)$	Two-factor independence
$\left(\mathrm{X}1{\displaystyle \cap }X2\right)>q\left(\mathrm{X}1\right)+q\left(\mathrm{X}2\right)$	Nonlinear enhancement

$q\left(\mathrm{X}1{\displaystyle \cap }\mathrm{X}2\right)$ is the interaction between X1 and X2. Nonlinear weakening means that when X1 and X2 act together; the combined interpretation degree of the two is less than that of a single X1/X2.

).

This study used the “factor detection” and “interactive detection” modules of the model to detect the impact of factors on ecosystem services. Precipitation and vegetation coverage have good protective effects on the water yield, soil conservation, and carbon storage, and the impact of human activities on the ecosystem cannot be ignored. As Hainan Province is the largest special economic zone in China and has experienced the rapid development of tourism in recent years, combined with the actual situation of Hainan Province, this paper selects precipitation (X1), normalized difference vegetation index (NDVI) (X2), net primary productivity (NPP) (X3), the digital elevation model (DEM) (X4), and slope (X5) as natural condition factors. Cultivated land area (X6), population (X7), gross domestic product (GDP) per capita (X8), land use degree (X9), urbanization rate (X10), and overnight tourist population (X11) are selected as social and economic factors.

2.3.7. Spatial Response Analysis Based on MGWR

MGWR (multi-scale geographically weighted regression) is a software application developed by Stewart Fotheringham. It is used to calibrate the geographically weighted regression (GWR) model, which can be used to explore the spatial relationship between dependent variables and independent variables [40]. MGWR adopts the backward fitting method for the spatial scale of variables, such that each selected variable presents different scale characteristics. This is more in line with the spatial heterogeneity of geographical processes. The smaller the bandwidth of the variable, the smaller the spatial influence scale and the stronger the spatial heterogeneity [41]. The calculation formula is as follows:

$$Ln{y}_i={\beta }_{bw0}\left(u_i,v_i\right)+\sum {\beta }_{bwk}\left(u_i,v_i\right)x_{ik}+{\epsilon }_i$$

(14)

where ${\beta }_{bw0}(u_i, v_i)$ is the $bw0$ regression intercept under the bandwidth, ${\beta }_{bwk}(u_i, v_i)$ is the optimal bandwidth $bwk$ regression coefficient of the k-th variable, $x_{ik}$is the i-th value of the k-th variable$, \mathrm{and} {\epsilon }_i$ is a random error term.

3. Results

3.1. Spatio-Temporal Changes in Ecosystem Services

Obvious spatiotemporal heterogeneity was found in the three ecosystem services in Hainan Province. From the perspective of quantity, the water yields in 2000, 2010, and 2020 were 56.13 billion m³, 42.06 billion m³, and 28.91 billion m³, showing a trend of continuous decline (

Table 3. Temporal changes in the ecosystem services during 2000–2020.

	Water Yield (Billion m3)	Soil Conservation (Million t)	Carbon Storage (Million t)
2000	56.13	16.1	80.8
2010	42.06	18.0	80.4
2020	28.91	11.5	79.2
2000–2020 (change)	−27.22	−4.7	−1.6

). The water yield is affected by the landform, climate and hydrology. It is mainly affected by precipitation and potential evaporation; the precipitation levels in Hainan Province in 2000, 2010, and 2020 were 2153 mm, 2138 mm, and 1987 mm, respectively, and exhibited a continuous downward trend. The potential evaporation is 1117 mm, 1778 mm and 1718 mm, respectively. There is heavy precipitation and low potential evaporation, which leads to the difference in the temporal change of the water yield. From the perspective of space, the water yield in the east was higher and gradually decreased to the west, in which the cities of Qionghai, Wenchang, and Wanning had the higher water yields (Figure 2A). Due to the large amount of precipitation and relatively weak evapotranspiration in the eastern region, the precipitation in the western region is relatively low. The land use intensity and spatial distribution characteristics of soil attributes also lead to the strong spatial heterogeneity of the water yield. Soil conservation in Hainan Province was $1.61\times {10}^7 \mathrm{t}$, $1.80\times {10}^7 \mathrm{t},$ and$1.15\times {10}^7 \mathrm{t}$ in 2000, 2010, and 2020, respectively, showing a trend of first rising and then declining (Table 3). In 2002, Hainan Province implemented the Grain for Green Program, giving priority to the return of farmland to forests in the soil erosion area in the middle of Hainan Province and desertification areas in the west of Hainan Province. Therefore, soil conservation showed an upward trend. In 2018, the Chinese central government proposed to establish the China (Hainan) Pilot Free Trade Zone in Hainan, and a large amount of construction land expansion occupied the original forest land and cultivated land, resulting in a downward trend in soil conservation. The areas with high soil conservation were mainly located in the central part of Hainan Province, such as Wuzhishan city and Qiongzhong city (Figure 2B). The regional precipitation was more abundant and the rainfall erosion was higher than in the northwest, but due to the higher vegetation coverage in the area, the soil conservation capacity was strong. The soil conservation areas were distributed in the northern part of Hainan Province, and the changes in soil conservation from 2000 to 2020 were not obvious. The carbon storage of Hainan Province was $8.08\times {10}^8 \mathrm{t},$ $8.04\times {10}^8 \mathrm{t}$, and$7.92\times {10}^{8 }\mathrm{t}$ in 2000, 2010, and 2020 respectively, showing a decreasing trend (Table 3), but the overall rate was relatively slow. High carbon storage values were found in the interior of Hainan Island, such as in Wuzhishan city and Qiongzhong city (Figure 2C), while the amount of carbon storage in the surrounding areas was relatively low. The reason for this finding is that more arable land and forests are located in the interior of Hainan Island, significantly enhancing the carbon sequestration, while the periphery contains more cultivated land, building land and water areas. The aboveground carbon density of cultivated land is low, and the carbon density of water areas is 0 where the amount of carbon storage is relatively low.

3.2. Analysis of the Ecosystem Service Hot Spots

Before the analysis of the cold and hot spots, ArcGIS spatial autocorrelation analysis was used to judge whether ecosystem services have aggregation characteristics in space. In 2020, the Moran index of ecosystem services in the study area was greater than 0 (

Table 4. Spatial autocorrelation index of the ecosystem services in 2020.

	Moran’s I	Z	P
WY	1.03	72.29	0.000
SC	0.265	49.28	0.000
CS	0.163	98.94	0.000

Moran’s I is Moran’s index, I < 0 is negatively correlated, I > 0 is positively correlated, I = 0 is not correlated. The Z value is the multiple of the standard deviation; if Z > 2 58, it is classified as an extremely significant hot spot area; if Z < −2 58, it is classified as an extremely significant cold spot area; if −1 ≤ Z ≤ 1, it is classified as insignificant area. The smaller the p value, the less likely the observed spatial pattern is to produce a random process.

), indicating that the ecosystem has spatial correlation; the Z values were greater than 2.58 and the p values were less than 0.01, indicating that the spatial distribution of the ecosystem services in the study area has aggregation characteristics. Therefore, the hot spot analysis of ecosystem services can be carried out. The hot spots of the water yield service are mainly distributed in the north of Hainan Province (Figure 3), accounting for about 20.7%. These areas have more precipitation and less potential evaporation. The cold spots are mainly in the southwest of Hainan Province, such as Dongfang city. This is similar to the spatial distribution of the water yield. The soil conservation service in Hainan Province is mainly made up of cold spots, accounting for about 97%. Hot spots are mainly distributed in Qiongzhong city and other areas with high altitudes. This area is mainly forest land, with high vegetation coverage and few human disturbance activities; the soil conservation capacity was strong. The cold and hot spots served by carbon storage showed complex characteristics.

3.3. Ecological System Services Trade-Offs/Synergies Analysis

Water yield, soil conservation, and carbon storage showed different correlations in different land use types (

Table 5. Relationships between the ecosystem services under different land use types.

	WY-SC	WY-CS	SC-CS
Cultivated land	/	0.37	0.37
Forest land	/	0.36	0.36
Grassland	/	0.38	0.39
Water	−0.085	/	/
Construction land	0.13	0.44	0.44
Unused land	0.2	−0.15	−0.11

The “/” indicates that the representatives are not related to the two factors. WY-SC: the service relationship between water yield and soil conservation; WY-CS: the service relationship between water yield and carbon storage; SC-CS: the service relationship between soil conservation and carbon storage.

). The service relationship between the water yield and carbon storage (WY-CS) and the service relationship between the soil conservation and carbon storage (SC-CS) of grassland, cultivated land, construction land, and forest land passed the significance level test of 0.01, with all of them exhibiting values greater than 0.3. That is, water yield and soil conservation service capabilities will increase with the increase of carbon storage. The total amount of ecosystem services of construction land was low, but due to the ecological protection during the development of Hainan Island, including ecological engineering projects such as artificial wetlands and parks, water yield and carbon storage had a significant correlation on the construction land. The correlations in the WY-CS and SC-CS in unused land were −0.15 and −0.11, respectively, revealing a trade-off relationship; that is, increases in water yield service capabilities reduced the soil conservation. Regional vegetation with a high water yield generally exhibited strong carbon storage. Therefore, WY-CS had a strong synergistic relationship. The WY-CS and SC-CS in the water did not pass the significance test. As water can be stored directly without carbon fixation, there was no obvious correlation between these ecosystem services.

3.4. Ecological System Services Driving Factor Space Response Analysis

Various ecosystem services in Hainan Province showed obvious spatial heterogeneity. This differentiation was commonly affected by natural factors and human activities. Because the influence mechanism of the ecosystem services in different years varied little [42], this article only explained the driving factors in 2020. In order to avoid the distortion of geographical detectors caused by multi collinearity, this paper uses stepwise regression analysis to analyze the multi collinearity relationship of variables. Using the variance expansion factor (VIF) to judge the collinearity between variables, when 0 < VIF < 10, there is no multicollinearity; when 10 ≤ VIF < 100, there is strong multicollinearity; and when VIF ≥ 100, there is severe multicollinearity. The results showed that there was no obvious multicollinearity among the factors (

Table 6. Results of the multiple co-curve test of influencing factors.

	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10	X11
VIF	1.361	1.437	1.883	2.186	1.778	1.417	1.178	3.273	2.367	2.568	1.589

).

3.4.1. Single-Factor Attribution of Ecosystem Services

On the whole, the precipitation (X1) and cultivated land area (X6) explained the ecosystem services with a large degree of interpretation, at 0.593 and 0.375, respectively (Figure 4). It was found that the ecosystem services were not only affected by natural factors but were also affected by human activities.

From the perspective of a single ecosystem service, the amount of regional precipitation directly affected the water yield. Therefore, precipitation had the most impact on the water yield service (0.555). When there was little change in precipitation, the slope affected the water yield capacity. The greater the slope, the poorer the water yield effect. Therefore, the contribution rate of the slope to water yield services was low (X5, 0.072). The degree of land use represented the development level of a region, and its factor interpretation of the water yield services was second only to precipitation (X9, 0.337). Therefore, when considering the construction of a new urban area, planners should also focus on protecting the environment. With the continuous development of Hainan’s tourism industry, the interpretation of the water yield services by the number of overnight tourists was 0.266 (X11). This finding indicates that the government needs to develop ecological tourism according to local conditions during the future development process. The DEM interpretation of the soil conservation service was the greatest (X4, 0.114). The areas with strong soil conservation functions were mainly concentrated in areas with high altitudes and high slope changes in the central part of Hainan Province. In general, the contribution rates of various factors to the soil were relatively low, indicating that the influencing factors of soil conservation services were more complicated. The contribution of the NDVI (X2) to carbon storage was 0.262, and areas with high NDVI values were generally areas with high vegetation coverage, so the amount of carbon was high. The per capita GDP (X8) had a secondary impact on carbon reserves of 0.146, indicating that in the process of rapid economic development, population increases lead to increases in construction land and decreases in the vegetation area. The relationships between ecosystem services should be considered in planning and development.

3.4.2. Ecological System Service Interaction Detection and Analysis

Interactive detection can explain whether a two-factor combination will enhance or weaken the interpretation of variables. The results showed that the interpretation of the interaction between the driving force factors in Hainan Province was greater than the degree of interpretation of any single factor (Figure 5), and the combined relationships had two combination methods: nonlinear enhancement and dual-factor enhancement.

Overall, the spatial pattern of ecosystem services was constrained by both natural and socioeconomic factors. The NPP (X3) ∩ human activities (the interaction between NPP and human activity factors) showed significant nonlinear enhancement relationships with the three ecosystem services (the interaction of two factors exceeded the explanatory power of the sum of the two single factors). Precipitation (X1) ∩ human activity also showed nonlinear relationships with soil conservation and carbon storage, which meant that human activities could increase the impact of precipitation on ecosystem services.

From the perspective of various services, the interaction between the precipitation of the water yield service and other factors exhibited a good interpretation of the water yield service in Hainan Province, with all of the values exceeding 0.56. Among them, the value of the precipitation ∩ NDVI (X2) was as high as 0.58, and it showed a nonlinear enhancement relationship, indicating that at the same precipitation level, the NDVI would have a greater impact on the spatial pattern of the water yield service. The NPP ∩ urbanization rate factor was 0.4, indicating that the water yield service in Hainan Province is affected by both human activities and NPP. In the soil conservation service, the soil conservation service explained by the interaction between precipitation and NDVI had a value of 0.13. However, because these single factors were not strong enough to explain the soil conservation, their interaction was not significant. In the carbon storage service, the interaction between NDVI and other factors had a greater interpretation of Hainan carbon storage service, with all of the values being about 0.3, of which NDVI ∩ NPP reached the maximum at 0.34. NPP ∩ human activity factors all showed nonlinear enhancement relationships, indicating that carbon fixation services were sensitive to human activity. The urbanization rate was significantly increased by several times compared with its single-factor interpretation. This indicated that the impact of urbanization on the water yield was significantly magnified under the comprehensive effect of other factors.

3.4.3. Spatial Difference Analysis of the Main Controlling Factors of Ecosystem Services

According to the previous driving force analysis results, because the driving mechanism of the soil conservation service was relatively complex, this article does not discuss its spatial differences. This paper used the MGWR model to analyze the spatial and geographic relationship between the main controlling factors (precipitation, population, elevation, and NDVI) of the water yield and carbon storage services. The data points with insignificant (p < 0.05) spatial and geographical relationships for each driving factor were excluded. The spatial correlation coefficients of each driving factor and the spatial response of ecosystem services were obtained.

On the whole, the spatial response relationships between the DEM, population, and precipitation and water yield were more obvious (Figure 6), while the response relationship of the GDP was relatively weak, and its coefficient ranged from 0.036 to 0.045, exhibiting a decreasing trend. The elevation had a significant negative correlation in the north of Hainan Province, indicating that the lower the elevation, the lower the water yield. In central Hainan Island, the opposite result was found. The reason may be that Wuzhishan city, Qiongzhong city, and Baoting city are located at higher elevations, and are less affected by human activities. Therefore, the water yield effect is obvious in the areas with higher altitudes, and the results were reversed in northern Hainan Province, especially in areas with lower elevations. Areas such as Haikou city have a high intensity of human activities, thus affecting the service function of the water yield. The effect of precipitation on the spatial heterogeneity of the water yield in Hainan Province showed a decreasing trend from southeast to northwest, indicating that the precipitation increased by 1 mm and the average water yield depth increased by 0.284 mm.

Overall, the correlation coefficient between the effects of the DEM and population size on carbon storage services was relatively high (Figure 7), with a clear circular structure. The response relationship of the GDP was weak; its coefficient ranged from −0.081 to 0, and only the southeastern part of Hainan Province passed the significance test of 0.05. Spatially, the impact of GDP on carbon storage increased from southwest to northeast. The coefficients of the DEM ranged from 0.086 to 0.839, and its spatial response was similar to that of GDP. The spatial response of the population to carbon storage was high in the southern regions. The reason may be that the rapid economic development and large population in the northern and southern regions have led to the increase of construction land and cultivated land, thereby significantly impacting carbon storage. The central area of the island has a higher elevation, a relatively weak land use change trend, and a relatively weak spatial response to carbon storage.

4. Discussion

4.1. InVEST Model Feasibility Analysis

The simulation of the water yield, soil conservation, and carbon storage after the localization of the InVEST model was the basis of this study. The results were similar to the results obtained by the Hainan Provincial Water Resources Bulletin and the existing research results reported in Hainan Province [43,44,45], indicating that the model achieved good results in the simulation of ecosystem services. The simulation results of the water yield included both the surface water and groundwater, so the results were compared to the total water resources data in the Hainan Water Conservancy Yearbook [46]. The total water resources of Hainan Province in 2020 were 26.36 billion m³, and the simulation performed in this study yielded a total of 28.914 billion m³. In 2002 it was reported that the soil conservation in Hainan Province was $1.48\times {10}^7$ tons [43], while the simulation result obtained in the present study indicated that the soil conservation in 2000 was $1.61\times {10}^7$ t. In summary, the model has high accuracy, and can be used to analyze the spatial heterogeneity and influencing factors of ecosystem services.

4.2. Spatial Heterogeneity of the Ecosystem Services

From 2000 to 2020, precipitation showed a trend, and potential evapotranspiration showed a trend of rising first and then decreasing, mainly because they were significantly affected by precipitation and potential evapotranspiration [47]. Wuzhishan city and Baozhong county in Hainan Province had relatively high precipitation, but due to the high forest coverage in this area, the annual average actual evapotranspiration was also relatively high, such that the water yield in these areas was not the highest. In the past 20 years, the water yield of construction land has remained at a consistently high level. This may be due to the impervious subsurface layer created during the urbanization process, which increases the surface water yield and reduces the infiltration, such that the water yield capacity is relatively high. This result is similar to that of Li [48].

Geographical environmental factors have a significant impact on the soil conservation service, and the soil conservation capacity increases first and then decreases with the increase of altitude. This is similar to the results of Cheng [29]. At the same time, soil conservation is also affected by human activities. When Hainan Province implemented the Grain for Green Program, soil conservation showed an upward trend. The Chinese central government proposed the establishment of the China (Hainan) Pilot Free Trade Zone in Hainan, resulting in a downward trend in soil conservation. The results are consistent with the Hainan statistical yearbook [46].

The change of the land use type is the main direct cause of the change of carbon reserves in Hainan Province, which is similar to the results of Liu [7]. We can effectively improve carbon sequestration by converting cultivated land and unused land into forest land.

4.3. Ecosystem Service Trade-Offs and Synergies

The relationship between carbon storage and soil conservation is mainly synergistic, and there are differences in spatial patterns. The vegetation coverage level in areas with high carbon reserves is high, which inhibits soil erosion. At this time, there is a synergistic relationship between carbon sequestration and soil conservation services, which is similar to the research results of Wang [38]. In the future, it will be necessary to strengthen the research on the scale characteristics and regional differences of watershed ecosystem service trade-offs, build a grid-township-county watershed multi-level evaluation unit, and clarify the ecosystem service capacity and trade-offs on various scales.

4.4. Influencing Factors of Ecosystem Services

In this study, the geographic detectors were used to analyze the contribution of each driving factor to the ecosystem services, and the MGWR model was used to analyze the spatial responses of the effects of the main controlling factors at different spatial scales. By innovating on the traditional use of regression models and the global analysis of driving factors, this method can provide a new technology for the analysis of the driving factors of ecosystem services in Hainan Province. Natural and human activity factors such as precipitation, the NDVI, and land use intensity all play important roles in the spatial heterogeneity of various ecosystem services in Hainan Province [26], and the results of these articles are similar to those of this study. The differences in the single-factor driving force of water storage at different scales reveals the scale effect of water storage changes. Therefore, the selection of scales needs to be considered fully when optimizing regional water resource allocation [49]. Land use change will lead to changes in ecosystem types, which will lead to changes in ecosystem services. However, changes in ecosystem services are not necessarily caused by changes in land use types, but may also be caused by changes in ecosystem quality [50]. The spatial differentiation of the effects of the population and DEM on the water yield and carbon storage was relatively obvious, which was consistent with the results of a previous study on the spatial responses of the effects of human activities and the natural environment on ecosystem services [51]. Rainfall had a great impact on ecosystem services in Hainan Province, and the explanatory degrees of the land use intensity and population number factors on ecosystem services all exceeded 0.35. While rainfall is difficult to control due to the macroscopic effects of the global climate, land use intensity is a controllable management factor. Ecosystem services can be therefore be improved by changing human activities. The explanatory power of each factor increased significantly after the interaction with human activity factors. This indicates that the impact of human activity factors on ecosystem services cannot be ignored. In order to coordinate the relationship between the economy and the environment, this paper put forward the following suggestions for Hainan Province: (1) Hainan Province should adhere to the principle of “ecological priority and green development”, such as strictly implementing the ecological protection red line and preserving ecosystems such as woodland, wetland, grassland wetland, and grassland; (2) in the process of promoting urbanization, Hainan Province should strengthen the construction of ecological cities, such as through the construction of artificial lakes, parks, and green belts; (3) Hainan Province should optimize agricultural production methods, and should transform a single water-intensive agricultural model into a modern agricultural model with less water consumption, such as sustainable agriculture.

4.5. Limitations

First, when using the InVEST model to evaluate ecosystem services in this study, the biological coefficients and soil physicochemical properties referred to the existing calculation methods and related literature, such as soil depth, soil and water conservation measure factors, and carbon density data. This introduced errors into the results.

Second, when using the geographic detection model, variables must be classified. This paper selected the natural breakpoint method to classify these variables, which may have introduced some errors.

Third, the spatial trade-off and spatial intensity of ecosystem services need to be explored further. Therefore, future work should analyze the ways in which different data discretization methods affect the experimental results. There may be differences in the impact mechanisms of ecosystems at different geographic scales, so it is necessary to study the driving forces of ecosystem services from different temporal and spatial scales.

5. Conclusions

Based on the InVEST model, three ecosystem services—water yield, soil conservation, and carbon storage—in Hainan Province were simulated; then, the temporal and spatial changes and hot spots of the ecosystem services were analyzed, and the relationship between the ecosystem services of different land use types was analyzed. Finally, the driving factors of the ecosystem services were analyzed using the geographic detectors combined with the MGWR model to identify the spatial differences of the dominant factors. The results were as follows:

(1)

From 2000 to 2020, the average annual water yield and carbon storage in Hainan Province were 42.36 billion m³ and $8.01\times {10}^8 \mathrm{t}$, respectively, and both services exhibited a downward trend; the multi-year soil conservation amount was $1.52\times {10}^7 \mathrm{t}$.

(2)

The hot spots of the water yield service are mainly distributed in the north of Hainan Province; the soil conservation service in Hainan Province is mainly cold spots. Cold and hot spots served by carbon storage showed complex characteristics.

(3)

There was a trade-off relationship between the CS-WY and CS-SC in unused land, the other ecosystem services of different land use types were synergistic.

(4)

Precipitation and land use intensity have strong explanatory power for the water yield service, and the NDVI and elevation factors are the main controlling factors for soil conservation and carbon storage. The interaction among driving force factors in Hainan Province was greater than the explanatory degree of any single factor, and most of the human activity factors showed obvious nonlinear enhancement during the interaction.