*2.2. Methods*

To investigate the spatiotemporal complexity of temperature and its driving factors, the correlation dimension method and the Geogdetector method were used. It can be seen from Figure 3, we first showed the spatiotemporal pattern of temperature; then, we analyzed the complexity of temperature on the daily, seasonal, and annual scales by using the Correlation Dimension (CD) method;finally, the individual contribution rates and interactional contribution rates of driving factors to the temperature slope (TS) by using Geogdetector method were detected.

**Figure 3.** The framework of this study.

#### 2.2.1. Trend Analysis Method

Trend analysis is the most studied and most popular quantitative forecasting method by far. It is based on a known historical data to fit a curve, so that this curve can reflect the growth trend of things themselves, and then to predict the future according to this growth trend curve. Commonly used trend models include linear trend models, polynomial trend models, linear trend models, log trend models, power function trend models, exponential trend models, and so on [27]. In this study, we use the linear trend method to analyze the change trend of the time series:

$$y(t) = at + b,\tag{1}$$

where *y* represents the time series, *t* represents the time, *a* represents the linear slope, and *b* represents the intercept.

If *a* > 0, it indicates that the time series is increasing, if *a* = 0, it means that the time series is not changing, and, if *a* < 0, it indicates that the time series is decreasing. The size of *a* indicates the degree of change in time series.

#### 2.2.2. Kriging Interpolation Method

Temperature is a regionalized variable, which is changing with the variation of space position. In order to analyze the distribution of the plum rainfall in different years, Kriging interpolation is employed. Kriging interpolation (or space local estimation) is named by D. G. Krige, who is a mining engineer in South Africa, and it is an optimal interpolation method [28]. The original data of the regional variables and the structural characteristics of the variance function is used to estimate the value of non-sampling points unbiasedly and optimally [28]. In general, Kriging interpolation contains several types, namely Ordinary Kriging, Universal Kriging, Co-Kriging, and so on. Ordinary Kriging is shown below:

Assume that *<sup>Z</sup>*(*x*) is a regionalized variable that satisfies two-stage stationary hypotheses and intrinsic hypothesis. *m* is mathematical expectation, with covariance function and variance function all existing at the same time. The relation between them is indicated below:

$$E[Z(\mathbf{x})] = \mathfrak{m},\tag{2}$$

$$\mathbf{C}(h) = \mathbb{E}[Z(\mathbf{x})Z(\mathbf{x} + h)] - m^2,\tag{3}$$

$$\gamma(h) = \frac{1}{2}E[(Z(\mathbf{x}) - Z(\mathbf{x} + h))^2].\tag{4}$$

Assuming that there are *no* measured points in the neighborhood of *x*, namely *x*1, *x*2, ... , *xn*, for which the sample value is *<sup>Z</sup>*(*xi*)(*<sup>i</sup>* = 1, 2, 3, ... , *<sup>n</sup>*), the formula can be defined as follows:

$$Z^\*(\mathbf{x}) = \sum\_{i=1}^n \lambda\_i Z(\mathbf{x}\_i),\tag{5}$$

where λ*i* is a weight coefficient that presents the contribution degree of the observed values of *<sup>Z</sup>*(*xi*) to estimate the values of *<sup>Z</sup>*<sup>∗</sup>(*x*). Two points need to be noticed about this formula: on the one hand, the estimated value of *<sup>Z</sup>*<sup>∗</sup>(*x*) must be unbiased, namely the mathematical expectation of the deviation is zero; on the other hand, it must be optimal, namely the difference between the estimated value and the actual value is the smallest.

#### 2.2.3. Correlation Dimension (CD)

Since the appearance of fractal theory, fractal dimension has been welcomed by scholars as one of the quantitative indicators to describe whether the dynamic system has chaotic characteristics. There are different types of fractal dimensions, such as topological dimension, Hausdorff dimension, information dimension, and correlation dimension. As for the correlation dimension, Grassberger and Procaccia [29] proposed an analysis method for experimental time series data in 1983, which is to obtain the fractal dimension through the relationship between the integral *C (r)* and the distance *r* on the reconstructed phase space through the univariate time series. This method is called a G–P method. Because it is particularly suitable for experimental observation data and the algorithm is simple and easy to implement, it has been widely used. In this study, we use this method when calculating the correlation dimension.

The correlation dimension (CD) is usually applied to analyze time series and determine if it exhibits a chaotic dynamic characteristic [30,31]. Considering {*<sup>x</sup>*1, *x*2, *x*3, ... , *xi*, ... }, the equal interval time series of daily temperature, and the first *m* data are extracted, and they determine the first point in the *m*-dimensional space, which is denoted as *X*1. Then, remove *x*1, and take m data *x*2, *x*3, ... , *xm*+1, and the second point is composed of this set of data in *m*-dimensional space, which is recorded as *X*2. According to this, a series of phase points *X*1, *X*2, ... , *XN* are formed. Given the number *r*, and check how many point pairs (*Xi, Xj*) distance is less than *r*, and the ratio of the number that point pairs distance is less than *r* to the total number of point pairs *N* is denoted as *C(r)* [17]. It can be expressed as the following formula:

$$\mathbf{C}(r) = \frac{1}{N^2} \sum\_{\substack{i,j=1 \\ i \neq j}}^N \mathbf{i}\_{i,j} \neq \mathbf{1} \quad \theta(r - |\mathbf{X}\_i - \mathbf{X}\_j|), \tag{6}$$

where θ(x) is the Heaviside function, which is defined as:

$$\theta(x) = \begin{cases} 1 & \text{ $x > 0$ } \\ 0 & \text{ $x < 0$ } \end{cases} \tag{7}$$

If *r* is too large, the distance of all point pairs will not exceed it. In addition, this *r* cannot measure the correlation between phase points. In addition, appropriate reduce *r*, the following formula may exist:

$$C(r) \ltimes r^d.\tag{8}$$

If this relationship exists, *d* is a dimension called the correlation dimension, denoted as *D2*:

$$D\_2 = \lim\_{r \to 0} \frac{\ln C(r)}{\ln r}.\tag{9}$$

The limit here mainly represents a direction in which *r* is reduced, and it is not mean that the *r* must be close to 0. In the scale transformation of the actual system, there are scale restrictions in the magnitude of both directions. Exceeding this limit is beyond the featureless scale area.

Figure 4 shows the results of ln(*r*) versus ln(*r*), and the correlation dimension (*d*) versus embedding dimension (*m*) used the measured data of temperature in this paper. It is apparent that the correlation dimension, *D*2, is given by the slope coefficient of ln(*r*) versus ln*<sup>r</sup>*. According to (ln*<sup>r</sup>*, ln*C*(*r*)), *D*2 can be obtained by the least squares method using a log–log grid (as shown in Figure 4a).

To detect the chaotic behavior of the system, the correlation dimension has to be plotted as a function of the embedding dimension (as shown in Figure 4b).

The MATLAB 2014a software (Manufacturer, City, US State abbrev. if applicable, Country) was used to calculate Correlation Dimension. First, we calculated the time series of daily average temperature, monthly average temperature, and annual average temperature of each meteorological station during the period of 1980–2012, and then calculated the CD value of each station on different time scales through programming using MATLAB software (The MathWorks, Natick, MA, USA).

**Figure 4.** (**a**) a plot of ln(*r*) versus ln(*r*) and (**b**) the correlation dimension (*d*) versus embedding dimension (*m*).
