EN中文

Journal of Geographical Sciences >

2022 , Vol. 32 >Issue 6: 1076 - 1102

DOI: https://doi.org/10.1007/s11442-022-1986-5

Research article

Spatial non-stationary characteristics between grass yield and its influencing factors in the Ningxia temperate grasslands based on a mixed geographically weighted regression model

  • SONG Xiaolong , 1 ,
  • MI Nan , 1, * ,
  • MI Wenbao 1, 2 ,
  • LI Longtang 2
Expand
  • 1.College of Agriculture, Ningxia University, Yinchuan 750021, China
  • 2.School of Geography and Planning, Ningxia University, Yinchuan 750021, China
* Mi Nan (1987‒), PhD and Associate Professor, E-mail:

Song Xiaolong (1991‒), PhD Candidate, specialized in remote sensing monitoring of grassland resources and research on forage-livestock balance. E-mail:

Received date: 2021-11-17

  Accepted date: 2022-03-09

  Online published: 2022-08-25

Supported by

Ningxia Key R&D Project(2018BEB04007)

Ningxia Colleges and Universities First-Class Discipline Construction (Grass Science Discipline) Project(NXYLXK2017A01)

Fold

Abstract

Spatial models are effective in obtaining local details on grassland biomass, and their accuracy has important practical significance for the stable management of grasses and livestock. To this end, the present study utilized measured quadrat data of grass yield across different regions in the main growing season of temperate grasslands in Ningxia of China (August 2020), combined with hydrometeorology, elevation, net primary productivity (NPP), and other auxiliary data over the same period. Accordingly, non-stationary characteristics of the spatial scale, and the effects of influencing factors on grass yield were analyzed using a mixed geographically weighted regression (MGWR) model. The results showed that the model was suitable for correlation analysis. The spatial scale of ratio resident-area index (PRI) was the largest, followed by the digital elevation model, NPP, distance from gully, distance from river, average July rainfall, and daily temperature range; whereas the spatial scales of night light, distance from roads, and relative humidity (RH) were the most limited. All influencing factors maintained positive and negative effects on grass yield, save for the strictly negative effect of RH. The regression results revealed a multiscale differential spatial response regularity of different influencing factors on grass yield. Regression parameters revealed that the results of Ordinary least squares (OLS) (Adjusted R2 = 0.642) and geographically weighted regression (GWR) (Adjusted R2 = 0.797) models were worse than those of MGWR (Adjusted R2 = 0.889) models. Based on the results of the RMSE and radius index, the simulation effect also was MGWR > GWR > OLS models. Ultimately, the MGWR model held the strongest prediction performance (R2 = 0.8306). Spatially, the grass yield was high in the south and west, and low in the north and east of the study area. The results of this study provide a new technical support for rapid and accurate estimation of grassland yield to dynamically adjust grazing decision in the semi-arid loess hilly region.

Key words： grass yield; spatial non-stationary; mixed geographically weighted regression model; temperate grassland; Ningxia

Cite this article

SONG Xiaolong , MI Nan , MI Wenbao , LI Longtang . Spatial non-stationary characteristics between grass yield and its influencing factors in the Ningxia temperate grasslands based on a mixed geographically weighted regression model[J]. Journal of Geographical Sciences, 2022 , 32(6) : 1076 -1102 . DOI: 10.1007/s11442-022-1986-5

1 Introduction

Grass yields of grasslands are fundamental for maintaining the longevity of these ecosystems, and a direct indicator of grassland ecological health. As an important aboveground material basis of grassland, it is an ideal place for grazing livestock. Due to its good palatability, natural grassland forage has high nutritional value for improving the production performance of herbivorous livestock. Further, yields have important implications for analyses of grass-livestock balances within natural grasslands across the globe. Indeed, biomass estimations of natural grasslands have long been a research hotspot corresponding fields of grassland management (Niu et al., 2003; Xie et al., 2009; Wright, 2010; Vasilis et al., 2018; Cao et al., 2019). An accurate understanding of the spatiotemporal distribution of grassland biomass also provides an important scientific basis for determining the livestock carrying capacity, and establishing sustainable animal husbandry production (Huang et al., 2021).
Ningxia is located in the ecotone of agriculture and animal husbandry in the inland of Northwest China and is a temperate continental arid and semi-arid climate area. Ningxia is also a key area for the construction of “two screens and three belts” ecological security system in China. The grassland area, which belongs to the temperate grassland region, is 244.33×104 ha (Ji et al., 2020). Temperate grassland is distributed in the loess hilly area, south of Ningxia, mainly in Yuanzhou, which is located in the central and western parts of the Loess Plateau. Yuanzhou District is a semi-arid, loess hilly region, with a large grassland area mainly consisting of temperate grasslands, including temperate meadow and typical temperate grasslands varying by region; however, desertification commenced across the grasslands of Yuanzhou as early as 2000, and large-scale disorderly reclamation and planting led to the destruction of ~466.67 km2 of natural grassland. Further destructive behaviors, such as overgrazing, Rodent hazards, poor grassland management and indiscriminate harvesting of Chinese medicinal materials, have further led to the degradation of grassland production capacity. Approximately 31% of the total grassland area in Yuanzhou District has been damaged by reclamation, followed by “grazing prohibition and house feeding,” which has become the primary animal husbandry mode for all local farmers; whereas measures for “returning farmland to forest and grassland” are also being implemented gradually (Yu et al., 2009; Dai et al., 2012; Bai et al., 2015). The implementation of such grassland management measures for over 20 years has greatly improved the local ecology (Wang et al., 2016); however, not all grasslands have been developed or utilized. Additionally, there has been withdrawal of cultivated land, shrinking of the main source of livestock diet, and a persistent competition between people and livestock for food. Accordingly, the necessity of dynamically adjusting grassland management policies to alleviate the pressure of farmers' livestock feeding on forage is an important debate; therefore, it is particularly important to reevaluate the biomass of these temperate grasslands for informing any policy adjustments. Presently, the grass and livestock industry represents one of the three largest industries in Yuanzhou, and has gradually formed a development pattern of grass planting and livestock rearing. There are 2458 large beef cattle breeding households with more than 10 heads, accounting for 38% of large-scale beef cattle breeding.
Traditional grass yield estimates employ the sample method for site data measurement. Although this method produces more accurate first-hand data of grass biomass through an abundance of sample quadrat measurements, it is time- and labor-intensive with a low estimate accuracy; moreover, the estimation method of grass yield per unit area cannot meet the practical needs of unified management across various grassland types, and it is difficult to popularize and apply on a macro scale (Cao et al., 2018; Xin et al., 2020). The recent extensive application of remote sensing technology has become an obvious solution to addressing previous limitations when estimating grass yield (Liu et al., 2018). Compared with traditional methods, it is advantageous on a spatiotemporal scale. An ever-increasing number of relevant research results and remote sensing-based models of grass yield, including statistical (Ren et al., 1998; Roy et al., 2002) and comprehensive models (Tucker et al., 1985; Anderson et al., 1993), have been put forth (He et al., 2015); however, according to the first regularity of geography, grass yield is related to adjacent environmental variables, and varies by geographical location. Since the British geographer Fotheringham posited a geographically weighted regression (GWR) analysis method, the spatial relationship between grass yield and influencing factors have been assessed (Brunsdon et al., 1996; Qin et al., 2006; Kashki et al., 2021).
In recent years, novel GWR models have provided a new remote sensing-based approach for estimating grass yields. For example, You et al. (2014) introduced the concept of GWR modeling into remote sensing estimates of grass yield for the Three-River-Source National Park. Elsewhere, Liu et al. (2019) estimated and verified the grass yield in the Xilin Gol League based on the improved GWR multi-factor model (Liu, 2019); whereas Li et al. (2016) estimated the grass yield of Qinghai Province by constructing a GWR model, combining it with ground-measured data, producing higher levels of accuracy than when using a unified weighted regression model. At present, however, most remote sensing models of grass yield are based on ground-measured sample data and constructed using various linear functions or other mathematical relationships to address the limited spatial heterogeneity considered in traditional linear regression models; although, these models do not consider the impact of model accuracy on a variety of environmental impact factors over different scales, they also do not include the spatial heterogeneity scale differences of these influencing factors, thus resulting in large deviations of estimates (Fotheringham et al., 2010; Wei et al., 2012; Zeng et al., 2016; Otheringham et al., 2017; Chao et al., 2018a; Harris et al., 2018b; Oshan et al., 2019; Li et al., 2021). Although the GWR model has been widely used among different scientific fields, there are few reports on the analysis of grass yield distributions, or its spatial non-stationary characteristics using a mixed GWR model (MGWR), which allows for the inclusion of unique bandwidths of different influencing factors, thus improving upon traditional GWR models (Yu et al., 2020), and providing a novel perspective for estimating grass yields.
Accordingly, this study focuses on Yuanzhou District, an area of warm grassland in the southern mountainous area of Ningxia Hui Autonomous Region, China. At the same time, the Liupan Mountain in Yuanzhou contains a large area of temperate meadow and temperate grasslands, which are an important part of Ningxia temperate grassland. Livestock farmers heavily depend on the loess hilly area represented by the Liupan and Yunwu mountains, which are studied for the remote sensing estimation of grass yield of temperate grassland. Therefore, measured sample point data of these grassland were collected by establishing experimental quadrats for carrying out multiscale research and accuracy verification on the grass yield of temperate grassland in combination with concurrent hydrological data, a digital elevation model (DEM), estimates of net primary productivity (NPP), and other auxiliary data. To this end, an MGWR model was used to decipher the spatial non-stationary characteristics and inform the management of local grass and forage-livestock balance dynamic policies.

2 Material and methods

2.1 Study area and quadrat setting

2.1.1 Study area

Yuanzhou District is located in the northern part of Guyuan city, in the mountainous area of southern Ningxia along the northwest edge of the Loess Plateau in China and is an important ecological barrier in the region. It is a typical semi-agricultural and semi-pastoral area, containing primarily temperate grasslands (Wang et al., 2012). Yuanzhou is the largest temperate grassland area in Ningxia. The grasslands, comprising the largest area of temperate grassland in the montane areas of southern Ningxia, can be predominantly categorized into seven types: cold Artemisia, Miscanthus, oxtail Artemisia, iron Artemisia, Baili herb, tiger stick, Miscanthus thyme. Located between 105°58ʹ-106°32ʹE (50 km) and 35°46ʹ-36°38ʹN (~100 km), it has a total area of 2739.01 km2, accounting for 4.51% of Ningxia. Located in the transition zone between the central plains agricultural area and the frontier grassland area, the town has long been an important stop along the Silk Road and maintains an extensive history and culture (Figure 1). Yuanzhou has two important nature reserves on the Loess Plateau, namely the Liupan Mountain and Yunwu Mountain. The Yunwu Mountain National Nature Reserve in this study area is the first Grassland Nature Reserve established in the Loess Plateau. Liupan Mountain, known as the “Green Island” of the Loess Plateau, is an important forest ecosystem type of nature reserve in Northwest China (Han et al., 2020). The Yunwu Mountain in Yuanzhou contains the largest reserved area of temperate grassland in the semi-arid area of the Loess Plateau, and represents the unique temperate grassland natural ecosystem of the Loess Plateau. The Yunwu Mountain is the natural ecological “background” of the Loess Plateau, the reservoir of biological resources, and a natural “treasure house” for studying the occurrence, development, and evolution law of temperate grassland ecosystem in the semi-arid area of the Loess Plateau. Yuanzhou District governs 11 townships namely Sanying, Tanshan, Huangduobao, Pengbao, Touying, Zhaike, Guanting, Zhonghe, Hechuan, Zhangyi, and Kaicheng, 149 administrative villages, and three urban sub-district offices: Beiyuan, Nanguan, and Guyan Streets. Its climate and hydrological environment are deeply affected by the northwest wind circulation because of its location in the middle latitude zone, and finally forms a temperate semi-arid continental climate type under the atmospheric mass influence of the Qinghai Tibet Plateau. The average annual temperature is 6.8℃, with an average annual rainfall between 300‒550 mm, and a large interannual variation in precipitation rate. Qingshui River, the largest and longest tributary flowing into the Yellow River in Ningxia, originates at the foot of Liupan Mountain in Yuanzhou and is also the largest river in Guyuan city. Due to the elevation of Liupan Mountain, the Qingshui River system is formed by the Jinghe and Weihe river systems. Further, the Ruhe River here belongs to the Jinghe River basin; whereas the Zhangyi River belongs to the Weihe River Basin.
Figure 1 Location and sample distribution of the study area (Yuanzhou District, Guyuan, Ningxia, China)

2.1.2 Quadrat sample locations

August typically comprises the primary growing season in temperate grasslands. Under the guidance of the field investigation method of grassland resources, and in accordance with the sample plot setting principle of Ren et al. (1998), a total of 149 experimental sample plots were arranged with the administrative village in August 2020 as the basic sampling units of grassland biomass statistics to sample across various administrative boundaries, grassland types, landforms, typical landscape characteristics, and accessibility of herbivorous livestock in Yuanzhou District (Figure 1). The quadrat setting covered all grassland areas in Yuanzhou, and field positioning was carried out using a GPS locator. A 1 m2 plot was mowed to the ground for sampling, all herbaceous plants in the plot, and recorded along with location (village and geographic coordinates), altitude, and fresh weight of each sample plot. This sampling was repeated three times using the diagonal method, for a total of 447 measured quadrats used to obtain the global grass yield through Kriging interpolation. In addition, using the previous research results for reference, downscaling of the quadrat data was performed to match the 30-meter resolution of grass yield inversion results (Arshad et al., 2021). In this study, 2698 random sample points were used as observation values for modeling and accuracy verification.

2.2 Data processing and analysis

2.2.1 Variable optimization

There are numerous influential factors on grass yield, and three requirements must be considered when selecting variables for multiple linear regression modeling: first, only those influencing factors with a significant impact on the distribution of grass yield should be considered; second, the spatial heterogeneity of these influencing factors must be accounted for; third, there must not be multicollinearity among the influencing factors selected. Here, SPPS v. 26.0 was used to analyze the correlation between all possible influencing factors—ratio resident-area index, temperature data, Air humidity data, NPP, topographic data, night light—and grass yield in the ground-measured samples. Correlation analyses were used to statistically analyze the dependency between variables, and the Pearson correlation coefficient (r) was used as an evaluation index to determine linear correlation. The correlation calculation results are shown in Figure 2. This correlation thermodynamic diagram is made by R software (v. 3.6.3).
Figure 2 Pearson correlation analysis of influencing factors on grass yield
Among the 12 influencing factors, only slope and aspect did not pass the significance test (p < 0.05); whereas all other 10 were significant at the p < 0.01 level. The highest Pearson coefficient was the DEM (0.695), followed by the relative humidity (RH; 0.596).

2.2.2 Variable determination

Based on the above Pearson coefficient calculation results, excluding the uncorrelated slope and aspect factors to grass yield, and considering the delayed impact of climate factors on grass yield (Liu et al., 2019), 10 factors were selected for the influencing factors: average July rainfall, RH, net primary productivity of vegetation (NPP), ratio resident-area index (PRI), elevation (DEM), daily temperature range, distance from gully, and distance from road (Table 1). As the constants in classical GWR and MGWR modeling represent the influence of different locations on the dependent variable when other independent variables are fixed, they can capture the influence of other external factors to a certain extent.
Table 1 Influencing factors and descriptions included in the analysis of grass yield
Variable Abbreviation Unit Variable description
Intercept Intercept g Intercept term of the model
Average rainfall in July AJR mm Average July rainfall
Relative humidity RH % Percentage of water vapor pressure in the air vs. saturated water vapor pressure at the same temperature
NPP NPP gC·m‒2·a‒1 Net primary production capacity of vegetation, refers to the total amount of organic matter accumulated by photosynthesis in unit area and unit time of green plants, minus the remaining part after autotrophic respiration
Ratio resident-area index PRI - Proportion of impervious surface in the surface area per unit area: $P R I=\frac{B L U E}{N I R}$, where blue and NIR are the pixel reflectance values of blue and near infrared wavelengths, respectively
Elevation DEM m Altitude of sample point
Daily temperature range DTR ℃ Difference between maximum and minimum daily temperatures
Distance from gully DS m Distance from sample point to valley
Distance from road DP m Distance from sample point to road
Distance from river DR m Distance from sample point to river
Night light NL - Night light distribution in the study area
In order to quickly realize the estimation of grass yield of temperate grassland in semi-arid loess hilly region of Ningxia, so as to guide local rational grazing, and at the same time, it is necessary to meet the unity of quadrat setting and spatial data of influencing factors. All data select the cross-sectional data in August 2020 as the basic data of model operation. The August 2020 data used for variable optimization included rainfall amount, elevation, RH, impervious surface index, and city light brightness (http://www.cma.gov.cn/2011qxfw/2011qsjgx/). As one of the leading climatic factors of grassland growth, rainfall has a certain delay to plant growth, hence the month of rainfall is selected as July (Liu et al., 2019). Climatological station observation data of RH, daily temperature range, and annual rainfall in the study area were acquired from China National Meteorological Data Center; elevation and remote sensing images came from the scientific data center geospatial data cloud of the computer network information center at the Chinese Academy of Sciences (http://www.gscloud.cn/sources/); NPP data of vegetation for 2020 were collected from NASA (https://search.earthdata.nasa.gov/search); spatial distribution of rivers and roads were recorded by Open Street Map (OSM) foundation (http://download.geofabrik.de/asia/china.html); night lighting data was provided by the Institute of Aerospace Information Innovation at the Chinese Academy of Sciences (http://satsee.radi.ac.cn/cfimage/nightlight/). Kriging interpolation was required for meteorological station data, the DEM required further surface analysis to obtain the slope and aspect; spatial distribution of the valley was obtained through hydrological analyses; the impervious surface index of remote sensing imagery was obtained using band math following initial preprocessing. Specifically, ENVI software was used to preprocess Landsat OLI remote sensing data of the United States (atmospheric correction, geometric correction and orthography correction) to extract urban land information. All data for the 10 influencing factors in Yuanzhou can be seen in Figure 3.
Figure 3 Spatial distribution of influencing factors across Yuanzhou District (NPP: Net primary production; PRI: Ratio resident-area index; DEM: Elevation; AJR: Average rainfall in July; NL: Night light; DP: Distance from road; DS: Distance from gully; DR: Distance from river; DTR: Daily temperature range; RH: Relative humidity)

2.2.3 Multicollinearity test

Considering the likelihood of multicollinearity between the influencing factors selected, a local variance expansion factor (VIF) was selected for testing this feature, to limit any deviation of estimation results caused by their interaction when estimating the grass yield of temperate grasslands (Eq. (1)):
$V I F=\frac{1}{1-R_{i}^{2}}$
where $R_{i}^{2}$ is the square of the determination coefficient. Thus, VIF is the reciprocal of tolerance, where the greater the VIF, the lower the tolerance between variables, and the greater the problem of collinearity. It is generally accepted that VIF should not be ≥5, or it can be relaxed to ≤10; otherwise, multiple collinearities are present among the independent variables (i.e., an independent variable can be expressed by the linear expression of one or several other independent variables; Marquardt et al., 1970). All influencing factors VIF were ≤5, and their tolerance values were small (Table 2); thus, it was concluded that no serious multicollinearity phenomenon was present among model’s influencing factors.
Table 2 Multicollinearity diagnosis results among influencing factors
Variable Tolerance VIF
AJR 0.244 4.092
RH 0.248 4.029
NPP 0.459 2.180
PRI 0.755 1.325
DEM 0.445 2.249
DTR 0.248 4.025
DS 0.792 1.262
DP 0.804 1.243
DR 0.754 1.327
NL 0.877 1.140

2.3 Methods

2.3.1 OLS model

Ordinary least squares (OLS) is a linear regression model where the weight between each observation point is equal. The OLS model is used to study the relationship between dependent and explanatory variables (Golub, 1980; Bergen et al., 2013; Jelinek et al., 2019), the starting point of all spatial regression analyses. Assuming that the linear regression relationship meets the global spatial stationary condition, the basic criterion is to find the best parameters by minimizing the residual sum of squares (RSS) between the predicted and observed values (Eq. (2)):
$y_{i}=\beta_{0}+\sum_{k} \beta_{k} x_{i k}+\varepsilon_{i}$
where yi is the value of the dependent variable at point i, β0 is the intercept, xik is the value of the k-th explanatory variable at point i, βk is the slope or regression coefficient of the k-th explanatory variable, and εi is the residual.

2.3.2 GWR and MGWR regression models

The GWR model is a geostatistical method. It was first put forth by Fotheringham et al (1998), where the geographical location of data is incorporated into the regression parameters based on the traditional global regression model and considering the spatial weight of adjacent points when estimating a local parameter (Brunsdon et al., 1996; Fotheringham et al., 1996; Fotheringham et al., 1998). Its expression is calculated according to Eq. (3):
$y_{i}=\sum_{i=0}^{n} \sum_{j=0}^{m} \beta_{j}\left(u_{i}, v_{i}\right) x_{i j}+\varepsilon_{i}$
where yi is the dependent variable, xij is the j-th independent variable of observation point i, βj(ui, vi) is the regression coefficient of the j-th independent variable in position (ui, vi), and εi is the random error term.
The GWR model is a special case of a global regression model; whereas the parameters in OLS are constant in space, those of GWR vary spatially, and the parameters of each observation point may be different, more accurately reflecting the spatial heterogeneity. The GWR expression indicates that each calibration point has its own regression model, and the corresponding parameter estimation of each calibration point is more affected by the data of nearer observation points (Qin et al., 2021). Notably, the observation data near each calibration point of the GWR model usually maintain similar properties (i.e., are spatially autocorrelated), and the selection of an appropriate bandwidth range is particularly important.
OLS, GWR, and MGWR all belong to the linear distribution regression method; however, because there is only one bandwidth in GWR, the local relationship in each model is limited to the same spatial scale and does not consider the conditional relationship between each independent and dependent variable on different spatial scales. Alternatively, the MGWR method allows the conditional relationship between predictor and response variables to change on different spatial scales (Geniaux et al., 2008; Wang et al., 2014; Fotheringham et al., 2019); that is, each variable changes on the surface of the parameter to allow the representation of different bandwidths (Zeng et al., 2016). MGWR minimizes over-fitting and reduces bias and collinearity in parameter estimation by eliminating the constraints of all relationships that vary on the same spatial scale; therefore, when GWR is being used to study spatial heterogeneity, the MGWR is the most ideal local model specification (Wei et al., 2012a). The parameter calibration process of MGWR commonly uses a back-fitting algorithm to calibrate the generalized additive model (Lianfa et al., 2019; Yu et al., 2019), defined according to Eq. (4):
$y_{i}=\sum_{i=0}^{n} \sum_{j=0}^{m} \beta_{b w j}\left(u_{i}, v_{i}\right) x_{i j}+\varepsilon_{i}$
where βbwj is the regression coefficient corrected by the effective bandwidth of the j-th independent variable.
The MGWR model is primarily established through the process of spatial weight function setting, bandwidth selection, model testing, and visual expression of coefficients.
(1) Space weight function setting
There are many types of spatial weight functions, among which bi_square and Gaussian kernel functions are the most widely used. The distance threshold kernel function sets all the weight values of observation points other than the distance threshold (bandwidth b) to 0. The Gaussian kernel function is continuous and uses all observation point data. The bi_square kernel function can be considered as a combination of these two kernel functions. Therefore, bi_square kernel function was selected as the spatial weight function in view of this characteristic. Within the bandwidth b of the calibration point i, the weight of each observation point is calculated by the continuous monotone decreasing function bi_square. The equation is as follows:
$\omega_{i j}=\left\{\begin{array}{ll}{\left[1-\left(\frac{d_{i j}}{b}\right)^{2}\right]^{2},} & d_{i j}<b \\0, & d_{i j} \geqslant b\end{array}\right.$
Similar to the Gaussian kernel function, the weight values of all observation points are within the range of [0,1] when the bandwidth b is determined. When the distance dij=0, the weight value of the corresponding observation point takes the maximum value 1. When dij increases, the weight value of the corresponding observation point gradually decreases; when dij increases to the bandwidth b, the weight value approaches 0.
(2) Bandwidth selection
The cross-validation method (CV) and Akaike Information Criterion (AIC) are two widely used bandwidth determination methods. Specifically, CV is an effective method for the criterion of least squares sum of squares when the limiting bandwidth value is 0 (Cleveland et al., 1979). The equation corresponding to this method is as follows:
$C V=\sum_{i=1}^{n}\left[y_{i}-\widehat{y_{\neq i}}(b)\right]^{2}$
where $\widehat{y_{\neq i}}(b)$ is the predicted value (excluding the calibration point) when the bandwidth is equal to b. This method can ensure that in the case of very small bandwidth b, the model is calibrated on the sample close to the calibration point, rather than on the sample with only the calibration point.
Akaike proposed a model selection criterion to measure the goodness of statistical model fitting in 1974, which is called AIC. AIC provides a standard to balance the complexity of regression model and the goodness of fitting data (Akaike et al., 1974). This criterion is one of the commonly used indicators in forecasting models using trend estimation. Its mathematical expression is:
AIC=2k-2ln(L)
where L is the likelihood function and k is the number of parameters of the regression model. The larger the likelihood function is, the better the model fit will be, and the larger the number of parameters k is, the higher the model complexity will be. In general, when the k of the model increases (that is, the complexity of the model increases), L is likely to also increase, making AIC smaller. However, when K is too large (that is, the model is too complex), the growth rate of the function L slows down, resulting in the increase of AIC; thus, the model is too complex and can easily cause over-fitting. The objective of model optimization is to find the model with the lowest AIC value. Using AIC can not only improve the applicability of the model, but also introduce penalty terms so that the model parameters are as small as possible, which helps reduce the possibility of over-fitting.
(3) Model test
Some scholars improved AIC and applied it to kernel function bandwidth selection in GWR analysis (Brunsdon et al., 2002). This criterion is called AICc, which is the test index of spatial non-stationarity of commonly used MGWR model. The AICc is calculated using the following equation:
$A I C c=2 n \log _{e}(\sigma)+\operatorname{niog}_{e}(\pi)+n\left(\frac{n+\operatorname{tr}(S)}{n-2-\operatorname{tr}(S)}\right)$
where σ is the sample standard deviation, n is the sample content, tr(S) is the trace of matrix S, which is the cap matrix of the model.
The Bayesian Information Criterion (BIC) proposed by Schwarz in 1978 is similar to AIC and can also be used to select the optimal bandwidth (Schwarz et al., 1978). The BIC was applied to GWR kernel bandwidth optimization selection (Nakaya et al., 2001). Similar to AIC, BIC also introduces penalty terms related to the number of model parameters; however, the penalty degree of BIC is greater than that of AIC. When the sample size is too large, the model can be effectively prevented from being too complex. The formula is as follows:
$B I C=k \ln (n)-2 \ln (L)$
where L is the likelihood function, n is the number of samples, and k is the number of model parameters; kln(n) is the penalty term, and the larger the sample size, the larger the penalty term. Under the condition of constant sample number, the model with fewer parameters is optimal.
(4) Visual representation of coefficients
Arc GIS ordinary Kriging interpolation method was used to visualize the coefficient values of spatial sample points obtained by model regression. This method can predict the spatial law within the whole research range based on the finite spatial sample points. Specifically, this model includes the positioning relationship between sample points and spatial structure characteristics into the analysis, which is the optimal linear unbiased method to predict the unknown from the known. The formula is as follows:
$R\left(x_{0}\right)=\sum_{1}^{n} \lambda_{i} R\left(X_{i}\right)$
where R(x0) is the value of the unknown sample point, R(Xi) is the value of known points near the unknown point, λi is the weight of the first sample point to other unknown points, and n is the number of known points.

2.3.3 Model evaluation

Here, mean error (ME), mean absolute error (MAE), and root mean square error (RMSE) were used to evaluate the accuracy of the model (Qiu et al., 2016), and calculated according to Eqs. (11-13):
$M E=\frac{1}{n} \sum_{i=1}^{n}\left[Z\left(x_{i}\right)-Z^{*}\left(x_{i}\right)\right]$
$M A E=\frac{1}{n} \sum_{i=1}^{n}\left[\left|Z\left(x_{i}\right)-Z^{*}\left(x_{i}\right)\right|\right]$
$R M S E=\sqrt{\frac{1}{n} \sum_{i=1}^{n}\left[Z\left(x_{i}\right)-Z^{*}\left(x_{i}\right)\right]}$
where n is the number of observation points used for verification, Z(xi) is the measured value, and Z*(xi) is the modeled value. Theoretically, the closer the MAE and RMSE are to 0, the better the prediction effect of the grass yield model. In some prior analyses, it has been difficult for the three indices to reach a single optimal solution; therefore, in addition to the above indices, the radius index (Eq. 14) proposed by Yang et al. (2019) was introduced to help identify the optimal model:
$\text { Radius }=100 \sqrt{\left(M E_{S}\right)^{2}+\left(M A E_{S}\right)^{2}+\left(R M S E_{S}\right)^{2}}$
where Radius is the radius-index, and MES, MAES, and RMSES are the respective values standardized by the Z-score method. The lower the radius-index, the better the model effect.

3 Results

3.1 Comparative analysis of model parameters

The calculation result parameters of the OLS global, GWR, and MGWR models are shown in Table 3. The goodness-of-fit R2 of the GWR model was as high as 0.801 compared with OLS model (0.643), and the AICC of this model was much greater than 3000 compared to the OLS model (ΔAICC = 457.844), the adjusted R2 increased from 0.642 to 0.797, further supporting that the fit of the GWR model was significantly higher than that of the OLS model. Moreover, the analysis of variance of the residuals of the two models showed that the RSS of the GWR model decreased by 722.856. When comparing the parameters of GWR and MGWR, it was found that the goodness-of-fit R2 of the MGWR model was higher (0.891), and the AICC value was significantly lower (2951.878), revealing that the MGWR results were superior to those of the GWR model. Further, the sum of squares of the residuals of the MGWR model was 23.701, an order of magnitude lower than the GWR model (236.960), indicating that more accurate regression results were obtained by using fewer influencing factor parameters. The adjusted R2 increased from 0.797 to 0.889 for the MGWR model, further supporting the superior fit compared to that of the GWR model. Overall, the significantly better performance of the MGWR model was found to be sufficient for explaining the spatial distribution of grass yield.
Table 3 Comparison of statistical parameters of different linear regression models: ordinary least squares (OLS), geographically weighted regression (GWR), and mixed GWR (MGWR)
Parameter OLS GWR MGWR
Residual sum of squares 959.816 236.960 23.701
-2 log-likelihood 3697.432 3138.206 2558.904
Classic AIC 3731.432 3270.442 2159.397
AICc 3731.660 3273.816 2951.878
BIC/MDL 3831.737 3660.555 2331.959
CV 5863.375 4370.185 -
R2 0.643 0.801 0.891
Adjusted R2 0.642 0.797 0.889
Local R2 have intuitively expressed the local model performance and mapping its distribution reveals that the performances of the GWR and MGWR models were superior to OLS (Figures 4d and 4e). The local R2 value of the GWR model was between 0.55-0.92, whereas that of the MGWR model was between 0.55 and 0.98; however, the spatial distribution of local R2 indicates that the MGWR model predicted a larger area more accurately, with less area of low-value local R2. Figures 4a-4c shows the local residuals of the OLS model (-2.59 to 3.89), the GWR model (-0.546 to 0.76), and the MGWR model (-0.18 to 0.13). The staggered structure of OLS residuals (Figure 4c) is readily apparent, whereas the high-low staggered structure in the residual distributions of the GWR and MGWR models (Figures 4a and 4b, respectively) is less obvious; however, that of the GWR model is the more of the two owing to the fixed bandwidth of the GWR model. Accordingly, the unique bandwidths of different independent variables in the MGWR model led to lower spatial distribution characteristics of residuals, higher prediction accuracy, as well as decreased noise and errors caused by a single-scale model.
Figure 4 Fitting results of different linear regression models
Table 4 presents the estimates of the OLS model parameters, which resulted in the following expression for grass yield of temperate grasslands (Eq. (15)):
Grassyield=366.649+(20.505×AJR)-(14.883×RH)-(7.321×NPP)-(6.363×PRI)+(74.194×DEM)-(40.112×DTR)-(17.548×DS)+(3.870×DP)+(9.427×DR)-(7.397×NL)
Table 4 Estimates of OLS model parameters
Variable Coefficient T-test Significance (p)
Intercept 366.649 8.482 0.000
AJR 20.505 3.867 0.000
RH -14.883 -2.593 0.010
NPP -7.321 -2.305 0.021
PRI -6.363 -2.580 0.010
DEM 74.194 27.200 0.000
DTR -40.112 -7.649 0.000
DS -17.548 -5.498 0.000
DP 3.870 1.490 0.004
DR 9.427 4.379 0.000
NL -7.397 -1.181 0.004
No sample sizes exceeded the reasonable limit, all data were valid, and all influencing factors were significant (p < 0.05). AJR, DEM, DP, and DR were significantly positively correlated with grass yield; whereas RH, NPP, PRI, DTR, DS, and NL were significantly negatively correlated with grass yield. Although the OLS model does not consider local spatial non-stationarity, the parameter estimate results show that it can adequately explain the effects of influencing factors (Table 5). The R2 of the OLS model was 0.643, indicating that the fitting degree of the model, although notable, was not optimal (Table 4). As the spatial information of the OLS model was not completely extracted, it is incapable of adequately explaining the spatial heterogeneity between grass yield and the influencing factors.
Table 5 Estimates of GWR model parameters
Variable Min Max Lwr Quartile Median Upr Quartile
Intercept 664.366 165.834 373.563 451.646
AJR -49.684 92.181 -12.087 6.980 35.114
RH -39.418 144.282 -11.053 8.416 27.694
NPP -6.232 32.852 10.403 15.679 23.962
PRI -7.563 2.873 -3.637 -2.698 -1.007
DEM -13.540 99.800 -3.165 33.316 59.446
DTR -95.718 54.511 -24.421 -3.836 12.769
DS -32.859 -0.367 -20.809 -10.443 -4.266
DP -38.467 41.975 -7.914 0.182 8.727
DR -56.344 59.046 -18.597 -2.029 18.442
NL -53.406 56.678 -10.603 -1.480 10.847
The optimal bandwidth of the GWR model was determined at 1,218, according to the minimum AICC value, and the variable coefficients were obtained (Table 5). Each influencing factor was inconsistent with the grass yield in space, and each parameter value varied within a certain range, thus objectively verifying the existence of spatial non-stationarity in grass yield observation points for temperate grasslands. The parameter estimates of influencing factors (such as AJR, DEM, DP, and DR) were positive under the OLS model, but this directionality varied in the GWR model, indicating the spatial dependence of the influencing factors on the distribution of grass yield. Comparing the 10 variables, the regression coefficients of five indexes—AJR, RH, NPP, DEM, and DP—were mostly positive; whereas the regression coefficients of PRI, DTR, DS, DR, and NL tended to be negative.
As shown in Figure 4 and Table 3, MGWR model performance was superior OLS and GWR for multiple parameters. The explanatory ability of introducing spatial effects with the GWR model reached 0.797, 15.50% higher than that of the OLS model (R2=0.642); however, the explanatory ability of the MGWR model provided an additional increase of 9.2% (to 88.9%). Notably, the explanatory power of the OLS, GWR, and MGWR models improved gradually. Additionally, the AICC value of the MGWR model was 2159.397, the minimum among the three models, further suggesting the superiority of the MGWR model to explain the spatial differentiation characteristics and estimation results of grass yield. The regression coefficient values of the different influencing factors varied greatly in space (Table 6). Standard deviation values were generally small, indicating that the correlation between the influencing factors and grass yield deviated slightly.
Table 6 Estimates of MGWR model parameters
Variable Mean STD Min Median Max
Intercept 0.827 1.170 -1.119 0.817 2.834
DEM -0.040 0.165 -0.340 0.001 0.188
NL 0.456 2.361 -4.255 0.002 8.176
PRI 0.003 0.003 -0.001 0.003 0.012
NPP 0.077 0.112 -0.213 0.042 0.585
DS -0.011 0.054 -0.245 -0.013 0.237
DP -0.073 0.156 -0.415 -0.054 0.301
DR 0.180 0.300 -0.757 0.232 0.859
DTR -1.227 1.316 -4.145 -0.715 0.561
RH -3.013 1.329 -4.979 -2.750 -1.274
AJR 1.687 1.313 -0.267 1.498 3.938

3.2 Spatial heterogeneity of influencing factors

In the practical discussion of the relationship between grass yield and environmental impact factors, both the GWR and MGWR models consider the characteristics of spatial non-stationarity, as well as the problem of spatial scale; however, the MGWR model incorporates different spatial scales for each influencing factor. While spatially predicting grass yield, MGWR also locally estimates the regression parameters of different influencing factors at each spatial location. To reduce the influence of the number of influencing factors on the relative size of the regression coefficient in the modeling process and consider the impact of influencing factors under different spatial scales on the grass yield, the local regression coefficient standardized by the MGWR model was used to describe the spatial non-stationary characteristics after comparing the different degrees of influence of each factor across the locations (Figure 5). Overall, from the regression coefficients of the indices examined, it was shown that: AJR had a positive effect on grass yield; DEM has a positive effect on alpine grassland, but a negative effect on river valley and platform area; the global effect of RH on grass yield was negative, and NL had a strong negative effect in most areas, whereas the effect of PRI on grass yield was limited.
Figure 5 Spatial patterns of regression coefficients for influencing factors of grass yield based on the mixed geographically weighted regression (MGWR) model (a. AJR: Average rainfall in July; b. PRI: Ratio resident-area index; c. DEM: Elevation; d. RH: Relative humidity; e. NL: Night light; f. DP: Distance from road; g. DS: Distance from gully; h. DR: Distance from river; i. DTR: Daily temperature range; j. NPP: Net primary production)
The MGWR model makes full use of different bandwidth scales for each variable to create a more accurate regression analysis of grass yield in these temperate grasslands. Applying Kriging interpolation to the regression coefficient surface can approximate the trends of different influencing factors at higher spatial densities (Figure 5). The value of the influencing factor coefficient of the MGWR model aligned with the true values, and thus reflected the relationships between different variables and grass yield. Specifically, Figure 5d shows the regression coefficient for AJR of the MGWR model, ranging from -0.28 to 3.94, with high distribution patterns in the south, and low in the north. The spatial distribution of the positive and negative values was divided in half across the study area, with a clear boundary between the two, thus indicating the varying effect of July rainfall on grass yield across Yuanzhou District. Figure 5b shows the regression coefficient of PRI for the MGWR model,
which ranged from -0.001 to 0.012, with positive values in the northwest, and negative values in the central and southern regions, as well as in most other areas. Figure 5c shows the regression coefficient of the DEM for the MGWR model, ranging from -0.34 to 0.188, and maintaining a spatial pattern of high in the south, and low in the north; thus, the higher altitudes of Liupan Mountain in the south, and Yunwu Mountain in the east tend to correlate to greater grass yields, while the higher the altitudes of the Qingshui River valley plain tend to have the opposite effect. Figure 5j shows the regression coefficient of RH for the MGWR model, ranging from -4.98 to -1.25, where all negative values indicate that the greater the RH, the lower the grass yield, especially in the south of the town of Touying. The negative correlation effect of RH on grass yield was the largest observed among all variables. Figure 5e presents the regression coefficient of NL for the MGWR model, ranging from -4.25 to 8.17, with a generalized spatial pattern of high in the east, low in the west, and maintaining a notably large negative value area. Although NL was positively correlated with the Yunwu Mountain area, the stronger the night light in other areas, the lower the observed grass yield. The positive correlation effect of night light on grass yield was the largest observed among all variables. Figure 5f shows the regression coefficient of DP for the MGWR model, with a coefficient value ranging from -0.4 to 0.27, and displaying mixed positive and negative geographic patterns. The results showed that the road network between the towns of Sanying, Zhangyi, Guanting, and Yuanzhou District was developed, and the grass yield decreased as the distance from the road network increased; whereas the inverse was true with other montane roads. This is consistent with remote sensing interpretation and reality. Figure 5g shows the regression coefficient of the DS for the MGWR model, with a coefficient value ranging from -0.222 to 0.217, and generally presenting a spatial pattern of medium high and low values in the north and south, respectively. In the central Qingshui River Basin, the grass yield increased with increase in proximity to the valley, while the opposite pattern was observed in the Liupan and Yunwu Mountain Areas of the loess hilly area. Figure 5h shows the regression coefficient of DR for the MGWR model, with a coefficient value ranging from -0.756 to 0.857. Generally, the northwest of Kaicheng, Zhaike, and Huangduobao Towns were negative, while all other areas were positive. Figure 5i shows the regression coefficient of DTR for the MGWR model, with a coefficient value ranging from -4.145 to 0.56, and showing a general pattern of high values in the east, and negative in the west. The positive area accounts for the majority of space, and the daily temperature range in the southeast was the largest. Further, the daily temperature range in the northwest was the smallest, the corresponding correlation coefficient was also the smallest, and the degree of negative correlation was the strongest. Figure 5a shows the regression coefficient of NPP for the MGWR model, with coefficient values ranging from -0.2 to 0.513, and most areas maintaining negative values, with patterns of high in the central region, low at north and south ends. The higher NPP values in the core area of Yunwu Mountain and the western mountainous area of Yuanzhou District, the greater was the grass yield; whereas in other areas, especially in the Qingshui River Basin, higher NPP values were associated with lower yields, with the lowest NPP coefficient values were recorded in the town of Zhangyi.
To intuitively understand the spatial scale of different models for each variable, the optimal bandwidths for each of the influencing factors in the GWR and MGWR models were analyzed and compared (Table 7). The two models maintained fundamentally different characteristics, where a constant bandwidth of the GWR model was 1218, and each bandwidth of the MGWR model was specific to the influencing factor, indicating their varied scales of spatial non-stationary characteristics; whereas the GWR model did not reflect the effects of these different spatial scales. Indeed, the ability of the MGWR model to accurately reflect the differential spatial scale of each influencing factors, while the GWR model only reflects the average value of the influencing factors, leads to a significantly improved accuracy effect. The bandwidth of the GWR model was consistent for all influencing factors (1218); whereas the spatial scales of different influencing factors varied greatly through the operation of the MGWR model, with PRI maintaining the largest value. The large spatial scales of DEM and NPP are marked by small spatial heterogeneity, and the relatively smaller spatial scales of all other variables (<50) is indicative of the significant spatial heterogeneity present in Yuanzhou District. In the MGWR model results, all regression coefficients of the 10 influencing factors were significant (P < 0.05).
Table 7 Bandwidth comparison between classical GWR and MGWR models
Variable MGWR GWR
Intercept 43 1218
DEM 150 1218
NL 43 1218
PRI 1348 1218
NPP 113 1218
DS 49 1218
DP 43 1218
DR 50 1218
DTR 46 1218
RH 43 1218
AJR 46 1218

3.3 Spatial predictions of grass yield

According to the above analyses and model parameters, grass yield distribution maps for these temperate grasslands were drawn for the OLS, GWR, and MGWR models examined (Figure 6), and revealing similar spatial trends among each method. The grass yield estimates of the GWR and MGWR models (Figures 6b and 6c, respectively) were smoother primarily because the incorporation of influencing factors on the spatial non-stationary characteristics of grass yield. Conversely, when compared with the MGWR model, the spatial scale effect of various influencing factors on the degree of influence of dependent variables in the OLS model (Figure 6a) was not well described, resulting in a relatively poor prediction ability.
Figure 6 Prediction results of grass yield for temperate grasslands according to each linear regression model
For the MGWR model, the estimated grass yield was generally higher in the south and west, lower in the north and east, and maintained the highest levels of overall estimation accuracy. The highest modeled grass yield estimates were in the south of Kaicheng, and a small part of Zhangyi (993-1310 g·m‒2) notably higher than that estimated by the GWR model. The second highest yields were seen in the smaller, central areas of Zhangyi and Kaicheng (798-992 g·m‒2), also higher than similar regional estimates from GWR. Notably, the north of Kaicheng and the west of Zhonghe were low local yield areas (631-797 g·m‒2); however, it can be seen from the figure that the grass yield in the southwest of Zhangyi was the same as that in the surrounding areas of the urban construction area of Yuanzhou District, as there are many aquaculture industries, large cultivated land areas, and relatively flat terrains located here, and the resulting yield is quite different from that of the surrounding montane areas. Additionally, most areas of Touying, some areas of Sanying, and a few areas of Huangduobao in the Qingshui River Basin also maintained lower grass yields (500-630 g·m‒2), and the yields across most areas of Hechuan, Guanting, Zhaike, and Tanshan were generally the lowest observed (296-499 g·m‒2); The predicted grass yields in smaller montane regions of these towns were higher than that in the surrounding areas; whereas those in Touying to north of Yuanzhou District, were also very low, consistent with the ground observations.

3.4 Model accuracy verification

To further verify the advantages of the MGWR model in the modeled estimates of grass yields in these temperate grasslands, model accuracy was assessed. Initially, scatterplots of the observed and predicted values for the OLS, GWR, and MGWR models were created (Figure 7).
Figure 7 Scatterplots of the observed vs. predicted grass yields
The OLS model’s fitness accuracy was relatively high (R2 = 0.7816), with relatively compact result in areas with lower outputs, and maintaining a clear linear relationship (Figure 7a). According to the OLS model accuracy test parameters (Table 8), the RMSE was 104.945; however, in areas with larger grass yields, the results were overgeneralized, and held a unique relationship; thus, the OLS model is suitable for remote sensing estimates of grass yield in areas <800 g·m‒2.
Table 8 Statistical parameters for the different linear regression models analyzed
Model ME MAE RMSE Radius
OLS Model -7.0112 75.7841 104.9458 184.0735
GWR Model -4.6232 70.5916 97.6234 156.5453
MGWR Model -3.9770 65.2953 92.6180 39.0543
The GWR modeled parameter results showed a significantly improved prediction of grass yield over the OLS model, with an increase in R2 by 0.0287, and a significantly reduced RMSE to 97.6234 (Figure 7b and Table 8). Accordingly, it can be seen from Figure 6b that the prediction accuracy of the grass yield using the GWR model was significantly improved, its scatterplot was more compact, and the linear trend was stronger. There was no over-dispersion, as with the OLS model, in areas with higher yields; thus, the GWR model can be applied to yield superior remote sensing estimates across different grass yield areas.
According to the statistics of the MGWR model parameters (Table 6), and the relevant test parameters of the MGWR model (Table 8), the remote sensing estimates of grass yield by the MGWR model held clear spatial non-stationary characteristics. Different bandwidth scales explained the impact of different influencing factors on grass yield, and the correlated R2 reached 0.8306, 0.0203 higher than that of the GWR model. Similarly, the RMSE wassignificantly reduced by 5.0054. The MGWR model highlights non-grassland areas, such as buildings and water bodies, and the resulting predicted grass yields were more consistent with observations, presenting a more compact linear trend, especially across the areas with larger grass yields (Figure 7c). Notably, because the MGWR model considered spatial non-stationary and multiscale advantages, produced superior predictions across different yield ranges. Overall, the MGWR was more suitable than either the OLS or GWR for remote sensing estimation of grass yields in temperate grasslands.

4 Discussion

4.1 The performance of MGWR model considering multiscale spatial non-stationarity is superior

In this study, different model methods were used to estimate grass yield, and MGWR model was selected as the method for estimating grass yield of temperate grassland in semi-arid loess hilly region. At the same time, the spatial non-stationarity characteristics between different influencing factors and grass yield were analyzed, and the effects of different spatial scales were explored. To explore the prediction accuracy of the OLS more comprehensively, GWR, and MGWR models, previous research results were summarized, and statistical indices were introduced to evaluate the models’ fitting effect (Mei et al., 2004; Kumar et al., 2013; Fotheringham et al., 2017; Mansour et al., 2020; Arshad et al., 2021; Chen et al., 2021; Liu et al., 2021; Yang et al., 2021). Conventionally, MAE and RMSE are two of the most commonly used indicators for measuring model accuracy (Yang et al., 2019). MAE represents the average value of the absolute error; whereas RMSE is commonly used to reflect dataset dispersion and is the square root of the average differences between predicted and observed values, squared.
The MAE values of the OLS, GWR, and MGWR models were decreased; whereas the RMSE values were also decreased, respectively, indicating the superior fitting effect of the MGWR model (Table 8). The significantly lower radius value of the MGWR model also indicated that the MGWR >>> GWR > OLS when explaining the relationship between influencing factors and grass yield. For the information criterion indices, where smaller values represent superior model effects, the simulated results of the OLS model were AICC (3731.660) and AIC (3731.432), those of the GWR model were AICC (3273.816) and AIC (3270.442), and those of the MGWR model were AICC (2951.878) and AIC (2159.397; Table 3). Accordingly, it was clear that the simulated effects of the MGWR and GWR models were superior to the OLS models. The fitting degrees of the three models were 0.643, 0.801, and 0.891 for OLS, GWR, and MGWR, respectively, while the adjusted fitting degrees were 0.642, 0.797, and 0.889. The goodness-of-fit was consistent with the adjusted goodness-of-fit results, as well as the radius-index performance results (Figure 7); therefore, the MGWR model more effectively calculated the weight values of different variables through the kernel function to limit interference factors and noise. Further, the different bandwidths provided to explain the spatial influence and distribution of influencing factors on the grass yield of temperate grasslands, more effectively revealed the scale of impact for each variable; whereas in the GWR model, these scaled effects were masked due to the constant bandwidth.

4.2 The spatial distribution of grass yield is the result of the joint action of climate factors and human activities

The introduction of mixed geographically weighted regression into remote sensing estimation of grass yield is an attempt to improve the accuracy of remote sensing estimation of grass yield, and can expand the types of remote sensing estimation models of grass yield. Although the variables selected in this study are highly correlated with grass yield (P < 0.01), and have a multiscale impact. However, in the traditional sense, the factors affecting grass yield are various vegetation indexes, such as normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), soil adjusted vegetation index (SAVI), modified soil adjusted vegetation index (MSAVI) and optimal soil adjusted vegetation index (OSAVI) (Xu et al., 2007; Luo et al., 2015; Luo et al., 2010). Previous studies have introduced meteorological factors such as temperature and rainfall into the estimation model as parameters for grass yield estimation (Dong et al., 2018; Wei et al., 2012b). However, other studies have shown that the influencing factors of grass yield are not only affected by climate change and meteorological factors, but also disturbed by many human factors, such as urbanization, road traffic, and light (Dong et al., 2018; Tian et al., 2021). In view of this, the data of PRI, DP and NL are introduced in this study. The reality is that the above influencing factors are rarely considered in previous studies. Therefore, this study is also an attempt to introduce the remote sensing estimation data of grass yield, which integrates meteorological factors and human factors.

4.3 The prospect of future research on the spatial difference method of mountain meteorological elements, research scope, grassland type division and so on

Climate change affects the species composition, growth rate and species accumulation process of grassland ecosystem (Mackey, 1994; Qian et al., 2010; Li et al., 2013). Therefore, temperature and precipitation are indispensable key climatic factors in the process of grassland biomass growth and accumulation (Qiu et al., 2009; Xu et al., 2018), and their accuracy is very key to the study of grassland biomass accumulation. As an important means of spatial estimation of grassland biomass, remote sensing technology can estimate regional grassland biomass from different spatial scales (Sun et al., 2013). However, when estimating grass yield, MGWR model needs to consider many key climate factors affecting grassland biomass accumulation, and climate factors need to match remote sensing spatial data through various spatial interpolation methods. Due to the different accuracy of meteorological data obtained by different spatial interpolation methods, it may affect the grass yield estimation of MGWR model, and then affect the accuracy of the estimation results. Existing studies show that different spatial interpolation methods have advantages and disadvantages, which need to be optimized in combination with the sensitivity of different research areas, especially for areas with large fluctuations of meteorological elements in complex mountainous climate environment (Liu, 2019; Zhang, 2019b). In this study, the methods used assume smooth continuity in space, without considering other environmental factors such as mountain vertical differentiation (Sun et al., 2008; Wang, 2008; Zhan Shi, 2014). Although the area of high-altitude area in the study area is small, there may still be error results of low accuracy of meteorological data interpolation results (Liu, 2012). What are the effects of different interpolation methods on the accuracy of data, Optimizing the best spatial interpolation method of meteorological elements for estimating grass yield of temperate grassland, so that the data interpolation results reflect better physical significance, which are major scientific problems that must be further studied in the future.
The sample plots in this study were located the temperate grasslands of the semi-arid loess hilly region in Ningxia, and although the prediction accuracy of grass yield in the involved regions was high, the sample size of this study is quite limited, and the spatial scale was restricted to the temperate grasslands. The temperate grassland in the selected study area is typical in the semi-arid loess hilly area; whereas the grassland types in this region of China are complex and diverse. Therefore, it is necessary to expand the scope of the study area in future research, and include a wider array of different grassland areas, increase the number of grass yield sample points observations, and further verify the model prediction in the loess hilly area across a large spatial range in an effort to provide greater theoretical and technical support for the popularization of the MGWR model in grass yield prediction, and the balanced management of grassland and forage-livestock within natural grasslands. Future research should promote the remote sensing estimation of grass yield in a larger research scope. Using the same technology will expand the application in the same type of semi-arid loess hilly area, to guide the monitoring of grassland resources, especially the dynamic adjustment of animal husbandry policy in the future.

5 Conclusions

The MGWR model was generally superior to the OLS and GWR models in terms of the fitting effect. Its goodness-of-fit (R2) was ≤0.891, and the adjusted R2 was ≤0.889. The AIC value and the sum of squares of residuals were significantly lower for the MGWR model compared to those of the other two models, indicating that the model is more suitable for the correlation analyses between grassland yield and influencing factors in the semi-arid loess hilly region of China, and is capable of overcoming limitations with regression weighted models used in traditional grass yield estimates.
The spatial scales varied among the influencing factors of the MGWR model. The results further showed that PRI maintained the largest spatial scale on grass yield of these temperate grasslands, followed by DEM > NPP > DR > DS > DTR > AJR > DP > RH > NL (with the spatial scales of night light, distance from road, and RH being the most limited). In the estimation of grass yield, the spatial scale bandwidths of the above influencing factors were 1348, 150, 113, 50, 49, 46, 46, 43, 43, and 43, respectively.
According to the statistical results of MGWR model parameters, and the spatial distribution of the regression coefficients of influencing factors, the impact of influencing factor-related indicators was spatially heterogeneous, with different spatial scales on grass yield of temperate grasslands in semi-arid loess hilly regions. Among them, RH had a negative effect on grass yield; whereas all other influencing factors had both positive and negative effects on grass dependent on the location. Thus, the multiscale differential spatial response regularity of different influencing factors on grass yield was highlighted, producing more realistic and reliable spatial analysis results than those of the GWR model.
According to the estimated results of temperate grassland yield using different regression models, the study area presented a spatial distribution pattern of high in the south and west, low in the north and east. The remote sensing estimates of grass yield by the MGWR model had obvious spatial non-stationary characteristics, and different bandwidth scales more accurately explained the impact of each influencing factor, as indicated by its larger R2 value, increased by 0.0203, and reduced RMSE by 5.0054 compared to the GWR model; thus, a superior spatial prediction performance of grass yield was realized, and the varying effects associated with the change in geographic position of each influencing factor according to their different spatial scales was obtained. This effect change has different spatial scales for different influencing factors.

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
Akaike H T, 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19(6): 716-723.

DOI

[2]
Anderson G L, Hanson J D, Haas R H, 1993. Evaluating Landsat thematic mapper derived vegetation indices for estimating above-ground biomass on semiarid rangelands. Remote Sensing of Environment, 45(2): 165-175.

DOI

[3]
Arshad A, Zhang W, Zhang Z, et al., 2021. Reconstructing high-resolution gridded precipitation data using an improved downscaling approach over the high altitude mountain regions of Upper Indus Basin (UIB). Science of The Total Environment, 784(2021): 147140.

DOI

[4]
Bai X, Wen Z, An S, et al., 2015. Evaluating sustainability of cropland use in Yuanzhou County of the Loess Plateau, China using an emergy-based ecological footprint. Plos One, 10(3): 1-10.

[5]
Bergen S, Sheppard L, Sampson P D et al., 2013. A national prediction model for PM2.5 component exposures and measurement error-corrected health effect inference. Environmental Health Perspectives, 121(9): 1017-1025.

DOI PMID

[6]
Brunsdon C, Fotheringham A S, Charlton M E, 1996. Geographically weighted regression: A method for exploring spatial nonstationarity. Geographical Analysis, 28: 281-298.

DOI

[7]
Brunsdon C, Fotheringham A S, Charlton M E, 2002. Geographically weighted summary statistics: A framework for localised exploratory data analysis. Computers Environment & Urban Systems, 26(6): 501-524.

[8]
Cao X, Sun L, Zhao Z, et al., 2018. Application of MODIS remote sensing products in the estimation of grass yield in Sanjiang Source Area. Remote Sensing for Land & Resources, 30(4): 115-124. (in Chinese)

[9]
Cao Y N, Wu J S, Zhang X Z, et al., 2019. Dynamic forage-livestock balance analysis in alpine grasslands on the northern Tibetan Plateau. Journal of Environmental Management, 238(15): 352-359.

DOI

[10]
Chao L, Zhang K, Li Z, et al., 2018. Geographically weighted regression based methods for merging satellite and gauge precipitation. Journal of Hydrology, 558: 275-289.

DOI

[11]
Chen Y, Zhu M K, Zhou Q, et al., 2021. Research on spatiotemporal differentiation and influence mechanism of urban resilience in China based on MGWR Model. International Journal of Environmental Research and Public Health, 18(3): 1056.

DOI

[12]
Cleveland W S, 1979. Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association, 74(368): 829-836.

DOI

[13]
Dai Y L, Tang W, Wang S, et al., 2012. Monitoring report on natural grassland vegetation change after returning grazing land to grassland in southern mountainous area of Ningxia. Ningxia Journal of Agriculture and Forestry Science and Technology, 3(11): 118-119, 122. (in Chinese)

[14]
Dong R L, Jia W, Yu G H, 2018. Research progress on prediction models of grass yield and livestock carrying capacity of grassland. Acta Agrestia Sinica, 26(5): 1043-1051. (in Chinese)

[15]
Fotheringham A S, Brunsdon C, 2010. Local forms of spatial analysis. Geographical Analysis, 31(4): 340-358.

DOI

[16]
Fotheringham A S, Charlton M, Brunsdon C, 1996. The geography of parameter space: An investigation of spatial non-stationarity. Geographical Information Systems, 10(5): 605-627.

[17]
Fotheringham A S, Charlton M E, Brunsdon C, 1998. Geographically weighted regression: A natural evolution of the expansion method for spatial data analysis. Environment and Planning A: Economy and Space, 30(11): 1905-1927.

DOI

[18]
Fotheringham A S, Yang W B, Kang W, 2017. Multiscale geographically weighted regression (MGWR). Annals of the American Association of Geographer, 107(6): 1247-1265.

DOI

[19]
Fotheringham A S, Yue H, Li Z, 2019. Examining the influences of air quality in China’s cities using multi scale geographically weighted-regression. Transactions in GIS, 23(1): 1444-1464.

DOI

[20]
Geniaux G, Napoléone C, 2008. Semi-parametric tools for spatial hedonic models: An introduction to mixed geographically weighted regression and geoadditive models. Hedonic Methods in Housing Markets, 101-127.

[21]
Golub G H, Loan C F, 1980. An analysis of the total least squares problem. SIAM Journal on Numerical Analysis, 17(6): 883-893.

DOI

[22]
Han X, Wang P, Wang J et al., 2020. Evaluation of human-environment system vulnerability for sustainable development in the Liupan mountainous region of Ningxia, China. Environmental Development, 34(5734): 100525.

DOI

[23]
He B, Li X, Quan X et al., 2015. Estimating the aboveground dry biomass of grass by assimilation of retrieved LAI into a crop growth model. IEEE Journal of Selected Topics in Applied Earth Observations & Remote Sensing, 8(2): 550-561.

[24]
Huang L, Ning J, Zhu P, et al., 2021. The conservation patterns of grassland ecosystem in response to the forage-livestock balance in North China. Journal of Geographical Sciences, 31(4): 518-534.

DOI

[25]
Jelinek A, Zalud L, Jilek T, 2019. Fast total least squares vectorization. Journal of Real-Time Image Processing, 16(2): 1-17.

DOI

[26]
Kashki A, Karami M, Zandi R, et al., 2021. Evaluation of the effect of geographical parameters on the formation of the land surface temperature by applying OLS and GWR, A case study Shiraz City, Iran. Urban Climate, 37(4): 100832.

DOI

[27]
Kumar S, Lal R, Liu D, et al., 2013. Estimating the spatial distribution of organic carbon density for the soils of Ohio, USA. Journal of Geographical Sciences, 23(2): 280-296.

DOI

[28]
Li N, Li T, Wang Y et al., 2016. Monitoring grassland yield in Qinghai Province based on Landsat8 remote sensing image. Journal of Qinghai University, 34(5): 63-68. (in Chinese)

[29]
Li W, Dong F, Ji Z, 2021. Research on coordination level and influencing factors spatial heterogeneity of China’s urban CO2 emissions. Sustainable Cities and Society, 75(1): 103323.

DOI

[30]
Li X Z, Han G D, Guo C Y, 2013. Impacts of climate change on dominant pasture growing season in central Inner Mongolia. Acta Ecologica Sinica, 33(13): 4146-4155. (in Chinese)

DOI

[31]
Liu H, Huang B, Gao S, et al., 2021. Impacts of the evolving urban development on intra-urban surface thermal environment: Evidence from 323 Chinese cities. Science of The Total Environment, 771(1): 144810.

DOI

[32]
Liu H X, 2019. Spatio-temporal analysis of grassland productivity and research on remote sensing estimation and prediction of grass yield in Inner Mongolia[D]. Shandong: Shandong University of Science and Technology. (in Chinese)

[33]
Liu J, Huang X, He X, et al., 2018. Estimation of grassland yield and carrying capacity in Qinghai Province based on MODIS data. Pratacultural Science, 35(10): 2520-2529. (in Chinese)

[34]
Liu K, 2012. Study on the mechanism of loess landslide induced by earthquake in Yuanzhou District, Guyuan city[D]. Xi’an: Chang’an University. (in Chinese)

[35]
Luo C, Yu X, Liu C, et al., 2015. Estimation of grass yield in large region on geographically weighted regression model. ISPRS: International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XL-7/W3(7): 9-13.

[36]
Luo L, Wang Z M, Ren C Y, et al., 2010. Models for estimation of grassland production and spatial inversion based on MODIS data in Songnen Plain. Transactions of the Chinese Society of Agricultural Engineering, 26(5): 182-187. (in Chinese)

[37]
Mackey Brendan G, 1994. Predicting the potential distribution of rain-forest structural characteristics. Journal of Vegetation Science, (4): 43-54.

[38]
Mansour S, Al K A, Al Said A, et al., 2020. Sociodemographic determinants of COVID-19 incidence rates in Oman: Geospatial modelling using multiscale geographically weighted regression (MGWR). Sustainable Cities and Society, 65(2): 102627.

DOI

[39]
Marquardt D W, 1970. Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation. Technometrics, 12(3): 591-612.

[40]
Mei C L, He S Y, Fang K T, 2004. A note on the mixed geographically weighted regression model. Journal of Regional Science, 44(1): 143-157.

DOI

[41]
Nakaya T, 2001. Local spatial interaction modelling based on the geographically weighted regression approach. GeoJournal, 53: 347-358.

DOI

[42]
Niu Z, Ni S, 2003. Study on models for monitoring of grassland biomass around Qinghai Lake assisted by remote sensing. Acta Geographica Sinica, 58(5): 695-702. (in Chinese)

[43]
Oshan T M, Li Z, Kang W et al., 2019. MGWR: A python implementation of multiscale geographically weighted regression for investigating process spatial heterogeneity and scale. International Journal of Geo-information, 8(6): 269.

[44]
Qian S, Fu Y, Pan F F, 2010. Climate change trend and grassland vegetation response in growing season in Sanjiangyuan area. Scientia Sinica Terrae, 40(10): 1439-1445. (in Chinese)

DOI

[45]
Qin N, Ren N, Huang N et al., 2021. Application of geographically weighted regression model in the estimation of surface air temperature lapse rate. Journal of Geographical Sciences, 31(3): 389-402.

DOI

[46]
Qin W, Wang J, 2006. Exploring spatial relationship non-stationary based on GWR and GIS. Proceedings of SPIE: The International Society for Optical Engineering, 6420: 64201X-8.

[47]
Qiu L, Kai W, Long W et al., 2016. A comparative assessment of the influences of human impacts on soil Cd concentrations based on stepwise linear regression, classification and regression tree, and random forest models. Plos One, 11(3): e0151131.

[48]
Qiu X F, Qiu Y P, Zeng Y, 2009. Distributed modeling of monthly mean air temperature of rugged terrain of Chongqing. Advances in Earth Science, 24(6): 621-628. (in Chinese)

[49]
Ren J Z, 1998. Research Methods of Grassland Science. Beijing: China Agricultural Press, 1-16. (in Chinese)

[50]
Roy D P, Borak J S, Devadiga S, et al., 2002. The MODIS land product quality assessment approach. Remote Sensing of Environment, 83(1/2): 62-76.

DOI

[51]
Schwarz G, 1978. Estimating the dimension of a model. The Annals of Statistics, 6(2): 461-464.

[52]
Sun C M, Liu T, Tian T, et al., 2013. Remote sensing estimation and application of grassland NPP based on MODIS data in southern China. Acta Prataculturae Sinica, 22(5): 11-17. (in Chinese)

[53]
Sun R H, Chen L D, 2008. Study on vertical differentiation and digital identification of the landscape in mountain areas. Proceedings of the 5th Chinese Young Ecologists Symposium, 326-331. (in Chinese)

[54]
Tian J, Xiong J N, Zhang Y C, et al., 2021. Quantitative assessment of the effects of climate change and human activities on grassland NPP in Altay Prefecture. Journal of Resources and Ecology, 12(6): 743-756.

[55]
Tucker C J, Vanpraet C L, Sharman M J, et al., 1985. Satellite remote sensing of total herbaceous biomass production in the Senegalese Sahel. Remote Sensing of Environment, 17(3): 233-249.

DOI

[56]
Vasilis D, Hauke H, Fabian L, et al., 2018. Predicting methane yield by linear regression models: A validation study for grassland biomass Bioresource Technology. Bioresource Technology, 265: 372-379.

DOI

[57]
Wang C, Zhen L, Du B, 2016. Assessment of the impact of China’s Sloping Land Conservation Program on regional development in a typical hilly region of the loess plateau: A case study in Guyuan. Environmental Development, 21: 66-76.

DOI

[58]
Wang K B, Li J P, Shangguan Z P, 2012. Biomass components and environmental controls in Ningxia grasslands. Journal of Integrative Agriculture, 11(12): 2079-2087.

DOI

[59]
Wang S J, Fang C L, Ma H T, et al., 2014. Spatial differences and multi-mechanism of carbon footprint based on GWR model in provincial China. Journal of Geographical Sciences, 24(4): 612-630.

DOI

[60]
Wang X, 2008. Study on biomass and species diversity on environmental gradient of alpine grassland[D]. Urumqi: Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences. (in Chinese)

[61]
Wei C H, Qi F, 2012. On the estimation and testing of mixed geographically weighted regression models. Economic Modelling, 29(6): 2615-2620.

DOI

[62]
Wei Y X, Wang L W, Shi Y C, et al., 2012. Net primary productivity of grassland resources monitoring based on remote sensing data in Qinghai Province. Scientia Geographica Sinica, 32(5): 621-627. (in Chinese)

[63]
Wright I A, 2010. Field and laboratory methods for grassland and animal production research. Grass & Forage Science, 57(2): 189-189.

[64]
Xie Y, Sha Z, Yu M, et al., 2009. A comparison of two models with Landsat data for estimating above ground grassland biomass in Inner Mongolia, China. Ecological Modelling, 220(15): 1810-1818.

DOI

[65]
Xu B, Yang X C, Tao W G, et al., 2007. Remote sensing monitoring upon the grass production in China. Acta Ecologica Sinica, 27(2): 405-413. (in Chinese)

DOI

[66]
Xu X, Xu Y, Sun Q Q, et al., 2018. Comparison study on meteorological spatial interpolation approaches in Kangdian region of China. Journal of Central China Normal University, 52(1): 122-129. (in Chinese)

[67]
Yang S H, Liu F, Song X D, et al., 2019. Mapping topsoil electrical conductivity by a mixed geographically weighted regression kriging: A case study in the Heihe River Basin, Northwest China. Ecological Indicators, 102(C): 252-264.

[68]
Yang Z, Xia J, Zou L, et al., 2021. Efficiency and driving force assessment of an integrated urban water use and wastewater treatment system: Evidence from spatial panel data of the urban agglomeration on the middle reaches of the Yangtze River. Science of The Total Environment, 805(2): 150232.

DOI

[69]
You H Y, Luo C F, Liu Z J, et al., 2014. Study on the method of grass yield model in the Source Region of Three Rivers with multivariate data. IOP Conference Series Earth and Environmental Science, 17(1): 012031.

[70]
Yu H, Fotheringham A S, Li Z, et al., 2019. Inference in multiscale geographically weighted regression. Geographical Analysis, 52(1): 1-20.

DOI

[71]
Yu H, Fotheringham A S, Li Z, et al., 2020. On the measurement of bias in geographically weighted regression models Spatial Stats. Spatial Statistics, 38: 100453.

DOI

[72]
Yu H J, Lu X M, Yue X X, et al., 2009. Dynamic change analysis and evaluation of forest resources in Yuanzhou District of Guyuan City. Modern Agricultural Science and Technology, (14): 181-183. (in Chinese)

[73]
Zeng C, Yang L, Zhu A X, 2016. Mapping soil organic matter concentration at different scales using a mixed geographically weighted regression method. Geoderma, 281: 69-82.

DOI

[74]
Zhang Q, 2019. Temporal and spatial patterns and dynamic changes of grass yield in Hulunbeier grassland[D]. Hohhot: Inner Mongolia University. (in Chinese)

Options
Outlines

/

〈 〉