2.3.1. Calculation of ESV

The ecosystem services classification method proposed by Costanza et al. cannot be directly applied to the calculation of ecosystem service value of China. This paper refers to the equivalent factor method and revised equivalent factor table based on Chinese ecosystem services [39], and modifies it by using grain yield correction and biomass correction methods in combination with the actual production capacity in Hubei province. With an average yield of 1 hm2, 1/7 of the annual economic value of natural grain yields to cropland is a standard ecosystem service value equivalent factor [37].

The specific correction process is as follows: Select early rice, medium rice, late rice, japonica rice, wheat and corn as the main food crops. By consulting Hubei Statistical Yearbook in 2016 and the National Compilation of Information on Costs and Returns of Agricultural Products, the grain crops to yield, sown area and grain price and output value of Hubei province in 2015 are obtained. According to Formula (1), the equivalent value of a single ecosystem in Hubei province is 2516.20 CNY/hm2. The ecosystem services value coefficients are shown below (Table 1).

$$E\_n = \frac{1}{7} \sum\_{i=1}^{n} \frac{q\_i p\_i}{M} \tag{1}$$

where *En* is the economic value (USD/hm2) of 1 hm2 ecosystem service value equivalent factor, *i* is the type of food crops in the study area; *n* is the total number of major food crops in the study area; *pi* is the average price of the crop *i* in the study area. *qi* is the average yield per unit area of type *i* food crop in the study area, and M is the total sown area of all food crops in the study area.



On this basis, according to the area of each land use type, the ecosystem services value of Hubei province can be calculated as follows [38,56]:

$$ESV = \sum A\_i \times V \mathbf{C}\_i \tag{2}$$

$$ESV\_f = \sum \left( A\_i \times V C\_{if} \right) \tag{3}$$

where *ESV* is the service value of the ecosystem in the study area, *ESVf* represents the value of the service function for species *f*, and *Ai* is the area of type *i* land use type. *VCi* is the ecosystem services value coefficient of type *i* land use type, and *VCi f* is the value of the *f*th ecological service function of the *i*th land use type.

#### 2.3.2. Land Use Transfer Matrix

Land use transfer matrix is the application of Markov model in land use change. It can show the transfer direction and quantitative change of different land use types and reveal the evolution process of land use pattern. It comes from the quantitative description of system state and state transfer in system analysis. In Formula (4), *S*(*t*) is the system state at time *t*, *S*(*t*+1) is the system state at time *t* + 1, *Pij* is the transition probability matrix in the state, and the calculation formula is as Formula (5) [64,65]:

$$S\_{(t,t+1)} = P\_{\vec{i}\vec{j}} \times S\_{(t)} \tag{4}$$

$$P\_{ij} = \begin{bmatrix} P\_{11} & P\_{12} & \dots & \dots & P\_{1n} \\ & P\_{21} & P\_{22} & \dots & \dots & P\_{2n} \\ & \vdots & \vdots & \vdots & & \vdots \\ & & P\_{n1} & P\_{n2} & \dots & P\_{nn} \end{bmatrix} \tag{5}$$

where *P* is Markov probability matrix, *Pij* represents the probability of transferring from the current state *i* to another state *j* in the next time period. The low conversion probability is close to 0 and the high conversion probability is close to 1.

In this research, after reclassifying the land use type grid data with ArcGIS10.2 software, import it into ENVI software, create ENVI standard classification format, edit the header file, and obtain the land use transfer matrix every 5 years through the Change Detection Statistics tool.

#### 2.3.3. Spatial Autocorrelation Analysis

Spatial autocorrelation analysis aims to reveal the correlation and difference between the spatial distribution of a certain attribute and its neighboring regions. If the spatial correlation is positive, it indicates that the spatial distribution of the attribute has agglomeration effect. If the spatial correlation is negative, it indicates that the spatial distribution of the attribute value is significantly different. Spatial autocorrelation is divided into global spatial autocorrelation and local spatial autocorrelation. This study uses the global autocorrelation model to determine the spatial distribution pattern of ecosystem service value. The Global Moran's *I* index is usually used for this calculation, and the formula is as follows [66,67]:

$$I = \frac{\sum\_{i=1}^{n} \sum\_{j=1}^{n} \mathcal{W}\_{ij} (\mathbf{x}\_i - \overline{\mathbf{x}}) (\mathbf{x}\_j - \overline{\mathbf{x}})}{S^2 \sum\_{i=1}^{n} \sum\_{j=1}^{n} \mathcal{W}\_{ij}}, \ I \in [-1, 1] \tag{6}$$

where *I* is Moran index. *xi* and *xj* are the observed values of the *i*th evaluation unit and the *j*th ESV evaluation unit, respectively. *Wij* is the spatial weight matrix between unit *i* and unit *j*. *x* is the mean of the observed values. n is the sample size, that is, the total number of ESV evaluation units at a certain scale in the study area. The Global Moran's *I* index is between −1 to 1. When the value is closer to −1, it indicates that the greater the difference between evaluation units, the stronger the negative correlation, and the more discrete the spatial distribution of ESV. When the Global Moran's *I* index value is closer to 1, it indicates that the attribute difference between evaluation units is smaller, the positive correlation is stronger, and the spatial distribution of ESV is more concentrated. When the Global Moran's *I* index value is close to 0, it indicates that there is no correlation between evaluation units, and the spatial distribution of ESV is random.

Hot spot analysis is a common tool to identify the spatial distribution of cold spots and hot spots. Using the Getis−Ord *G*\* index to analyze the local correlation between cold spots and hot spots can further study the local performance of the change of ecosystem service value in space. Hot spot analysis can determine whether there are high−value clusters (hot spots) and low−value clusters (cold spots) in the ecosystem service value of Hubei Province, as well as the spatial clustering positions of high−value and low−value regions [68]. By calculating the score Z (standard deviation) and probability P between each patch, the spatial location where the high−value elements and low−value elements gather [33].

$$\mathbf{G}\_i^\* = \frac{\sum\_{j=1}^n \omega\_{ij} \mathbf{x}\_j - \overline{\mathbf{x}} \sum\_{j=1}^n \omega\_{ij}}{S \times \sqrt{\frac{\left[n \sum\_{j=1}^n \omega\_{ij}^2 - \sum\_{j=1}^n \omega\_{ij}\right]^2}{n-1}}} \tag{7}$$

$$
\overline{\mathfrak{x}} = \frac{\sum\_{j=1}^{n} \mathfrak{x}\_{j}}{n} \tag{8}
$$

$$S = \sqrt{\frac{\sum\_{j=1}^{n} \mathfrak{x}\_{j}^{2}}{n-1}} - \overline{\mathfrak{x}}^{2} \tag{9}$$

*xj* is plaque attribute value, and *ωij* is the spatial weight matrix between plaque *i* and plaque *j*. *n* is the total number of plaques; *G*∗ *<sup>i</sup>* is the Z score, and the Z score is a measure of statistical significance. The higher or the lower the value of *G*∗ *<sup>i</sup>* , the tighter the accumulation of hot spots or cold spots. When the ESV value is much larger than the adjacent area, a statistically significant hot spot is formed, which is called the hot spot area, that is, the area with high ESV change value. When the ESV value is much smaller than the adjacent area, a statistically significant cold spot will be formed, which is called the cold spot area, that is, the area with low ESV change value.

#### 2.3.4. The Gravity Model

The gravity model is derived from physical concepts [69]. Later, it was widely used in trade [70], energy [71], transportation [72] and other fields. This research adopts the center of gravity shift model, which can intuitively reveal the spatial evolution characteristics of the center of gravity of ESV.

If the study area consists of n units (n represents administrative units), and (*xi*,*yi*) is the geometric coordinate of the *i*th unit (*I* = 1, 2, 3,..., *n*), then the barycenter coordinate of ESV is (*x*, *y*), it can be expressed as [73]:

$$\overline{\varpi} = \frac{\sum\_{i=1}^{n} m\_i \varkappa\_i}{\sum\_{i=1}^{n} m\_i} \tag{10}$$

$$\overline{y} = \frac{\sum\_{i=1}^{n} m\_i y\_i}{\sum\_{i=1}^{n} m\_i} \tag{11}$$

The movement direction of the center of gravity can be expressed as:

$$\theta = \left[\frac{k \times \pi}{2} + \tan^{-1}\left(\frac{\overline{y\_{t\_2}} - \overline{y\_{t\_1}}}{\overline{x\_{t\_2}} - \overline{x\_{t\_1}}}\right)\right] \times \frac{180^{\circ}}{\pi} \tag{12}$$

The movement distance of the center of gravity is calculated by the following formula:

$$D = \sqrt{\left(\overline{\mathbf{x}\_2} - \overline{\mathbf{x}\_1}\right)^2 + \left(\overline{\mathcal{Y}\_2} - \overline{\mathcal{Y}\_1}\right)^2} \tag{13}$$

where *θ* represents the deflection angle of the center of gravity of ESV, (*xt*1,*yt*<sup>1</sup> ) and (*xt*<sup>2</sup> , *yt*<sup>2</sup> ) respectively represent the coordinates of the center of gravity at the beginning and end. *t*1, *t*<sup>2</sup> represents the start time and end time, and k is the adjustment coefficient. *D* represents the moving distance of the center of gravity.

#### **3. Results**
