key: cord-0972748-vbjcsage
authors: Zhou, Chunshan; Zhang, Rongrong; Ning, Xiaoju; Zheng, Zhicheng
title: Spatial-Temporal Characteristics in Grain Production and Its Influencing Factors in the Huang-Huai-Hai Plain from 1995 to 2018
date: 2020-12-09
journal: Int J Environ Res Public Health
DOI: 10.3390/ijerph17249193
sha: 32a832962710933e99a8da18e553ad5f531b81c2
doc_id: 972748
cord_uid: vbjcsage

The Huang-Huai-Hai Plain is the major crop-producing region in China. Based on the climate and socio-economic data from 1995 to 2018, we analyzed the spatial–temporal characteristics in grain production and its influencing factors by using exploratory spatial data analysis, a gravity center model, a spatial panel data model, and a geographically weighted regression model. The results indicated the following: (1) The grain production of eastern and southern areas was higher, while that of western and northern areas was lower; (2) The grain production center in the Huang-Huai-Hai Plain shifted from the southeast to northwest in Tai’an, and was distributed stably at the border between Jining and Tai’an; (3) The global spatial autocorrelation experienced a changing process of “decline–growth–decline”, and the area of hot and cold spots was gradually reduced and stabilized, which indicated that the polarization of grain production in local areas gradually weakened and the spatial difference gradually decreased in the Huang-Huai-Hai Plain; (4) The impact of socio-economic factors has been continuously enhanced while the role of climate factors in grain production has been gradually weakened. The ratio of the effective irrigated area, the amount of fertilizer applied per unit sown area, and the average per capita annual income of rural residents were conducive to the increase in grain production in the Huang-Huai-Hai Plain; however, the effect of the annual precipitation on grain production has become weaker. More importantly, the association between the three factors and grain production was found to be spatially heterogeneous at the local geographic level.

China is an important food-producing country in the world, as well as a large food consumer. China's food self-sufficiency rate has reached more than 95% [1] . Although the current grain supply and demand in China maintains a balance of total quantity and a surplus in harvest, the small per capita arable land area, low mechanization level, and family-based business units determine that the current arable land has an extremely limited potential for increasing grain production [2, 3] . From 1850-1900 to 2006-2015, the mean land surface air temperature increased by 1.53 • C [4] . Climate change has already affected food security due to warming, changing precipitation patterns, and the greater frequency of some extreme events [4] . Under the dual constraints of climate change and the process of urbanization, unpredictable meteorological disasters, limited arable land resources, huge population pressure, and the diversified consumer demand of residents directly generated strong demands for stable grain production [5, 6] . Additionally, the outbreak of COVID-19 in 2019 led to labor shortages 

The Huang-Huai-Hai Plain (HHH), the second largest plain in China, is located at 32 • -40 • N and 114 • -121 • E, with a land area of about 4 × 10 5 km 2 , spanning seven provinces and cities of Beijing, Tianjin, Hebei, Shandong, Henan, Anhui, and Jiangsu ( Figure 1 ). The HHH belongs to the warm temperate semi-humid monsoon climate zone and is one of the most sensitive areas to climate change in China. The winter climate in this area is typically dry and cold, spring is dry, with less rain and much evaporation, and summer is characterized by high temperatures and heavy rainfall, including high intense rainfall that often leads to summer floods [24] . As an important agricultural region, this area has a long history of farming and is an important grain production base for food security in China, with its sown area of 20.4% of the nation's farmland and 23.6% of the whole nation's grain yield [25] .

Here, wheat and maize occupy a larger proportion in the structure of grain production. The annual double-cropping system of winter wheat and summer maize is the most popular planting pattern. In addition, the yield of wheat and maize accounts for approximately 61% and 31% of the total national output, respectively [26] . Therefore, it is necessary to determine the pattern change rules and influencing factors of grain production in upgrading HHH grain production and ensuring national food security. 

The Huang-Huai-Hai Plain (HHH), the second largest plain in China, is located at 32°-40° N and 114°-121° E, with a land area of about 4 × 10 5 km 2 , spanning seven provinces and cities of Beijing, Tianjin, Hebei, Shandong, Henan, Anhui, and Jiangsu ( Figure 1 ). The HHH belongs to the warm temperate semi-humid monsoon climate zone and is one of the most sensitive areas to climate change in China. The winter climate in this area is typically dry and cold, spring is dry, with less rain and much evaporation, and summer is characterized by high temperatures and heavy rainfall, including high intense rainfall that often leads to summer floods [24] . As an important agricultural region, this area has a long history of farming and is an important grain production base for food security in China, with its sown area of 20.4% of the nation's farmland and 23.6% of the whole nation's grain yield [25] .

Here, wheat and maize occupy a larger proportion in the structure of grain production. The annual double-cropping system of winter wheat and summer maize is the most popular planting pattern. In addition, the yield of wheat and maize accounts for approximately 61% and 31% of the total national output, respectively [26] . Therefore, it is necessary to determine the pattern change rules and influencing factors of grain production in upgrading HHH grain production and ensuring national food security. 

The historical socio-economic data used in this paper were collected from the Statistical Yearbook published by the National Bureau of Statistics of China, which includes annual data on the grain yield per unit area, the sown area of grain crops, the amount of agriculture fertilizer application, the amount of effective irrigation area, the amount of pesticide application, the amount of mechanical power, and the per capita annual income of rural residents in the HHH from 1995 to 2018. According to China's statistics, grain production includes corn, rice, soybeans, wheat, potatoes, and sweet potatoes.

The historical climate data were collected from the Chinese meteorological data hub (https://data.cma.cn), including the annual precipitation, annual temperature, and annual sunshine duration of 123 meteorological stations in the HHH from 1995 to 2018. The meteorological data need 

The historical socio-economic data used in this paper were collected from the Statistical Yearbook published by the National Bureau of Statistics of China, which includes annual data on the grain yield per unit area, the sown area of grain crops, the amount of agriculture fertilizer application, the amount of effective irrigation area, the amount of pesticide application, the amount of mechanical power, and the per capita annual income of rural residents in the HHH from 1995 to 2018. According to China's statistics, grain production includes corn, rice, soybeans, wheat, potatoes, and sweet potatoes.

The historical climate data were collected from the Chinese meteorological data hub (https:

//data.cma.cn), including the annual precipitation, annual temperature, and annual sunshine duration of 123 meteorological stations in the HHH from 1995 to 2018. The meteorological data need to be processed by SQL-Server (Microsoft Corporation, Redmond, WA, USA) and ArcGIS10.2 software (Environmental Systems Research Institute Inc, Redlands, CA, USA), which can analyze the average annual temperature, annual precipitation, and average annual sunshine hours of each city in different years. Grain production is influenced by a variety of natural and socio-economic factors. In this paper, based on the previous literature [27] [28] [29] , the grain yield per unit sown area (kg/km 2 ) was taken as the dependent variable, and climate and socio-economic factors were taken as the independent variables. Among them, the climate indexes included the annual average temperature ( • C), the annual average precipitation (mm), and the annual average sunshine duration (h), and the social-economic factors included the proportion of effective irrigation area (%), amount of fertilizer application per unit sown area (t/km 2 ), amount of pesticide application per unit sown area (t/km 2 ), amount of mechanical power per sown area (kw/km 2 ), and per capita annual income of rural residents (RMB). Table 1 summarizes the data used in this study. 

Exploratory spatial data analysis is a collection of techniques for describing and visualizing spatial distributions, determining atypical locations or spatial outliers, discovering spatial associations, clusters, or hot spots, and to infer spatial characteristics or other forms of space heterogeneity [30] . In general, global and local spatial autocorrelation (or hot spots analysis) is often used to explore the spatial characteristics of observations [31] .

Global spatial autocorrelation is used to test the spatial correlation of the observations of spatial units within the study area [32] , and is mainly measured by the Global Moran's I, which was first proposed by Moran [33] . The Moran's I can be calculated using Equation (1):

where I represents Moran's I, n is the number of spatial units (in this study, n = 59), x i and x j are the observations of spatial units I and j, respectively, x is the average value of observations of spatial units, and wij is the spatial weight matrix, where w ij = 1 if spatial units I and j share a common border and w ij = 0 otherwise. The values of Global Moran's I range from −1 to 1. If I < 0, it means there is a negative spatial correlation in the space, while if I > 0, it means there is a positive spatial correlation, an d if I = 0, it means there is no spatial correlation.

The significance of Moran's I is usually measured by Z statistics using Equation (2):

where E(I) and Var(I) are the expected value and variance of Moran's I, respectively.

The Getis-Ord G i * is commonly used for hot spot analysis, which can identify clustering relationships at different spatial locations. Compared with the local spatial autocorrelation, the Getis-Ord G i * is more sensitive to the identification of cold and hot spots, and can fully reflect the high or low value distribution relationship between a certain geographic element and other surrounding elements [34] . The formula is [34] [35] [36] [37] :

where x = 1 n n j=1

x j , n is the number of spatial units (in this study, n = 59), and w ij is the spatial weight matrix, where w ij = 1 if spatial units I and j share a common border and w ij = 0 otherwise. The significance of G i * is usually measured by Z statistics using Equation (4):

where E(G i * ) and Var(G i * ) are the expected value and variance of G i * , respectively. If Z(Gi * ) is significantly positive, it indicates that the observations around the spatial unit i are relatively high (higher than the average), and are high-value clusters in the space, belonging to hot spots; on the contrary, if Z(G i * ) is significantly negative, it indicates that the observations around the spatial unit i are relatively low (lower than the mean), and are low-value clusters in the space, belonging to cold spots. The larger (or smaller) the Z(G i * ) is, the more intense the clustering of high (or low) values. A Z(G i * ) near zero indicates no apparent spatial clustering.

The gravity center model is used to measure the overall distribution of a certain attribute in a region. It can provide a concise and accurate feature of the distribution of the attribute in the space, and can indicate the general trend and central location of its distribution. We assumed that a large region (such as an administrative region) consists of several subregions, and so the gravity center of grain production in the region could be calculated by the grain production and geographic coordinates of each sub-region. The formula is [38] :

In Equation (2), X i and Y i represent the geographic coordinates of the ith subregion. M i represents the grain yield per unit sown area of the subregion. X and Y represent the gravity center of grain production in a large region. Using formula (3), the moving distance of the gravity center in grain production can be obtained, which can reflect the evolution of the gravity center of a property in a region;

In Equation (3), D ij is the gravity center movement distance (km) of grain production from j to i years. (X i , Y i ) and (X j , Y j ) are the gravity center coordinates of grain production in the i and j years. R is typically 111.111, which represents the coefficient of the spherical longitude and latitude coordinates converted to plane distance.

When the data have spatial autocorrelation effects, the residuals are no longer independent of each other; thus, it is not appropriate to use the ordinary least square regression (OLS) model. Instead, the spatial lag model (SLM) or the spatial error model (SEM) should be used for analysis. The formula is [39] [40] [41] :

SLM : y = ρw ij y + xβ + µ

SEM :

In Equations (4) and (5), y is the dependent variable, x is the explanatory variables, W ij is the space weight matrix, ρ is the spatial hysteresis parameter, β is the parameter vector, µ is the random interference term, ε is the regression residual vector, and λ is the autoregression parameter.

Geographically weighted regression models are superior to traditional regression models such as ordinary least squares (OLS). The geographical weighted regression (GWR) model can fully consider the spatial characteristics of each influencing factor, and more accurately show the spatial relationship between independent and dependent variables [42] . The form of a GWR model is as follows:

In Equation (6), Y i represents the grain production in region i, β 0(ui,vi) represents a constant, β λ(ui,vi) represents the regression coefficient, (u i ,v i ) represents the geographic location of the cities i, X iλ represents the parameter value of the λ independent variable of city i, and ε represents the random error. The optimal bandwidth distance can be obtained automatically in GWR4.0 corrected by finite correction of the Akaike Information Criterion (AICc). The smaller the AICc value, the higher the goodness of fit of the model will be [43] . 

We classified the grain yield data of each urban unit of sown area into four types according to 2000-3500, 3500-5000, 5000-6500, and 6500-8000 kg/km 2 using ArcGIS10.2 software. For the convenience of analysis, four time sections of 1995, 2005, 2015, and 2018 were selected for study, for which the spatial distribution characteristics of the grain production pattern were discussed ( Figure 3 ).

It can be seen from Figure 3 that the spatial variation in grain yield per unit of sown area in each city is very significant, and the overall grain yield shows an increasing trend. Specifically, the number of cities in the HHH where the grain yield per unit of sown area remained within the range of 2000-3500 kg/km 2 continued to decrease, and was mainly distributed in Zhangjiakou, Cangzhou, Zhengzhou, Luoyang, Nanyang, Sanmenxia, and Bozhou in 1995, Zhangjiakou and Sanmenxia in 2005, and only in Zhangjiakou in 2015. In 2018, the number of cities of this type decreased to zero.

The number of cities within the range of grain yield per unit of sown area of 3500 to 5000 kg/km 2 also showed a downward trend, and the type area shrank from 30 cities in 1995 to 5 cities in 2018, and presented a layout trend of agglomeration to dispersion in spatial distribution.

The number of cities in the range of 5000-6500 kg/km 2 showed a rising-falling-rising trend. Among them, the number of cities in this type of area increased from 19 to 35 in 1995-2005, reduced from 35 to 28 in 2005-2015, and increased from 28 to 40 in 2015-2018. The spatial distribution of this type overall showed a tendency varying between scattered and clustered development, and appeared roughly from east to west, and from south to north.

The trend of unit of area sown to grain production maintained in the 6500-8000 kg/km 2 interval of urban change was different from the above three kinds. In 1995, the number of cities within this range was 3, which increased to 5 in 2005 and 24 in 2015, and reduced to 14 in 2018, which reflects a characteristic of a sharp increase followed by a slow decline. However, they were still mainly distributed in the central and southern parts of the HHH.

In terms of spatial distribution, the areas with high grain yield per unit of sown area in the Huang-Huia-Hai Plain were mainly distributed in the east and south, while the areas with low grain yield per unit of sown area were mainly distributed in the west and north of the HHH, which indicated that the grain production capacity of the eastern and southern regions of the HHH was higher, while that of the northern and western regions was lower. 

We classified the grain yield data of each urban unit of sown area into four types according to 2000-3500, 3500-5000, 5000-6500, and 6500-8000 kg/km 2 using ArcGIS10.2 software. For the convenience of analysis, four time sections of 1995, 2005, 2015, and 2018 were selected for study, for which the spatial distribution characteristics of the grain production pattern were discussed ( Figure 3 ).

It can be seen from Figure 3 that the spatial variation in grain yield per unit of sown area in each city is very significant, and the overall grain yield shows an increasing trend. Specifically, the number of cities in the HHH where the grain yield per unit of sown area remained within the range of 2000-3500 kg/km 2 continued to decrease, and was mainly distributed in Zhangjiakou, Cangzhou, Zhengzhou, Luoyang, Nanyang, Sanmenxia, and Bozhou in 1995, Zhangjiakou and Sanmenxia in 2005, and only in Zhangjiakou in 2015. In 2018, the number of cities of this type decreased to zero.

The number of cities within the range of grain yield per unit of sown area of 3500 to 5000 kg/km 2 also showed a downward trend, and the type area shrank from 30 cities in 1995 to 5 cities in 2018, and presented a layout trend of agglomeration to dispersion in spatial distribution.

The number of cities in the range of 5000-6500 kg/km 2 showed a rising-falling-rising trend. Among them, the number of cities in this type of area increased from 19 to 35 in 1995-2005, reduced from 35 to 28 in 2005-2015, and increased from 28 to 40 in 2015-2018. The spatial distribution of this type overall showed a tendency varying between scattered and clustered development, and appeared roughly from east to west, and from south to north.

The trend of unit of area sown to grain production maintained in the 6500-8000 kg/km 2 interval of urban change was different from the above three kinds. In 1995, the number of cities within this range was 3, which increased to 5 in 2005 and 24 in 2015, and reduced to 14 in 2018, which reflects a characteristic of a sharp increase followed by a slow decline. However, they were still mainly distributed in the central and southern parts of the HHH.

In terms of spatial distribution, the areas with high grain yield per unit of sown area in the Huang-Huia-Hai Plain were mainly distributed in the east and south, while the areas with low grain yield per unit of sown area were mainly distributed in the west and north of the HHH, which indicated that the grain production capacity of the eastern and southern regions of the HHH was higher, while that of the northern and western regions was lower. 

We classified the grain yield data of each urban unit of sown area into four types according to 2000-3500, 3500-5000, 5000-6500, and 6500-8000 kg/km 2 using ArcGIS10.2 software. For the convenience of analysis, four time sections of 1995, 2005, 2015, and 2018 were selected for study, for which the spatial distribution characteristics of the grain production pattern were discussed ( Figure 3 ).

It can be seen from Figure 3 that the spatial variation in grain yield per unit of sown area in each city is very significant, and the overall grain yield shows an increasing trend. Specifically, the number of cities in the HHH where the grain yield per unit of sown area remained within the range of 2000-3500 kg/km 2 continued to decrease, and was mainly distributed in Zhangjiakou, Cangzhou, Zhengzhou, Luoyang, Nanyang, Sanmenxia, and Bozhou in 1995, Zhangjiakou and Sanmenxia in 2005, and only in Zhangjiakou in 2015. In 2018, the number of cities of this type decreased to zero.

The number of cities within the range of grain yield per unit of sown area of 3500 to 5000 kg/km 2 also showed a downward trend, and the type area shrank from 30 cities in 1995 to 5 cities in 2018, and presented a layout trend of agglomeration to dispersion in spatial distribution.

The number of cities in the range of 5000-6500 kg/km 2 showed a rising-falling-rising trend. Among them, the number of cities in this type of area increased from 19 to 35 in 1995-2005, reduced from 35 to 28 in 2005-2015, and increased from 28 to 40 in 2015-2018. The spatial distribution of this type overall showed a tendency varying between scattered and clustered development, and appeared roughly from east to west, and from south to north.

The trend of unit of area sown to grain production maintained in the 6500-8000 kg/km 2 interval of urban change was different from the above three kinds. In 1995, the number of cities within this range was 3, which increased to 5 in 2005 and 24 in 2015, and reduced to 14 in 2018, which reflects a characteristic of a sharp increase followed by a slow decline. However, they were still mainly distributed in the central and southern parts of the HHH.

In terms of spatial distribution, the areas with high grain yield per unit of sown area in the Huang-Huia-Hai Plain were mainly distributed in the east and south, while the areas with low grain yield per unit of sown area were mainly distributed in the west and north of the HHH, which indicated that the grain production capacity of the eastern and southern regions of the HHH was higher, while that of the northern and western regions was lower. Figure 3 reflects the static distribution pattern of grain production in the HHH in 1995, 2005, 2015, and 2018, but does not reflect the dynamic change trend. Therefore, we used Equation (2) to calculate the barycenter of grain production in the HHH from 1995 to 2018 (Figure 4) , and then used Equation (3) to calculate the barycenter movement distance ( Table 2) to analyze the dynamic change characteristics of the grain production pattern.

According to Figure 4 , the grain production center in the HHH generally shifted from southeast to northwest in Tai'an, and gradually stabilized in the central area of the HHH with the passage of time. Specifically, from 1995 to 1997, the grain production center of HHH was distributed in Jining city, which moved first to the northwest and then to the southwest. From 1998 to 2000, the movement direction of the barycenter in grain production remained highly stable; that is, it continued to move in the northwest direction. The barycenter of grain production began to enter Tai'an City. From 2015 to 2018, the grain production center of gravity showed characteristics of moving from southeast to northwest, but it was still stable in the territory of Tai'an City, that is, the southeast of Figure 3 reflects the static distribution pattern of grain production in the HHH in 1995, 2005, 2015, and 2018, but does not reflect the dynamic change trend. Therefore, we used Equation (2) to calculate the barycenter of grain production in the HHH from 1995 to 2018 (Figure 4) , and then used Equation (3) to calculate the barycenter movement distance ( Table 2) to analyze the dynamic change characteristics of the grain production pattern.

According to Figure 4 , the grain production center in the HHH generally shifted from southeast to northwest in Tai'an, and gradually stabilized in the central area of the HHH with the passage of time. Specifically, from 1995 to 1997, the grain production center of HHH was distributed in Jining city, which moved first to the northwest and then to the southwest. From 1998 to 2000, the movement direction of the barycenter in grain production remained highly stable; that is, it continued to move in the northwest direction. The barycenter of grain production began to enter Tai'an City. From 2015 to 2018, the grain production center of gravity showed characteristics of moving from southeast to northwest, but it was still stable in the territory of Tai'an City, that is, the southeast of the HHH. The dynamic shift in the grain production center in the HHH indicates that the regional grain production capacity has the characteristics of non-stationarity in time and non-equilibrium in space at the same time, and the shift in the grain production center from southeast to northwest indicates that the grain production capacities in the western and northern parts of the HHH were continuously enhanced. In addition, the grain production center in the HHH tended to be stable over time, and it was concentrated in the border area between Jining and Tai'an. the HHH. The dynamic shift in the grain production center in the HHH indicates that the regional grain production capacity has the characteristics of non-stationarity in time and non-equilibrium in space at the same time, and the shift in the grain production center from southeast to northwest indicates that the grain production capacities in the western and northern parts of the HHH were continuously enhanced. In addition, the grain production center in the HHH tended to be stable over time, and it was concentrated in the border area between Jining and Tai'an. Table 2 shows that the barycenter of grain production in the HHH as a whole moved from the southeast to northwest from 1995 to 2005, with a distance of 17.7 km. From 2005 to 2018, the barycenter of grain production moved to the northwest with a distance of 9.4 km, which was significantly smaller than that from 1995 to 2005, which confirmed that the barycenter of grain production in the HHH showed good time-stability characteristics over time; however, this could not cover up the disequilibrium in the spatial characteristics of grain production. On the whole, from 1995 to 2018, the center of gravity of grain production moved 26.1 km to the northwest. The center was still stable in the territory of Tai'an city, that is, to the southeast of the HHH. Table 2 shows that the barycenter of grain production in the HHH as a whole moved from the southeast to northwest from 1995 to 2005, with a distance of 17.7 km. From 2005 to 2018, the barycenter of grain production moved to the northwest with a distance of 9.4 km, which was significantly smaller than that from 1995 to 2005, which confirmed that the barycenter of grain production in the HHH showed good time-stability characteristics over time; however, this could not cover up the disequilibrium in the spatial characteristics of grain production. On the whole, from 1995 to 2018, the center of gravity of grain production moved 26.1 km to the northwest. The center was still stable in the territory of Tai'an city, that is, to the southeast of the HHH. 

Based on the grain yield data per unit sown area, the Moran's I value, Z statistic, and P-value were calculated using Geoda software, and the spatial correlation characteristics of grain production are shown in Table 3 .

It can be seen from Table 3 that the Moran's I values from 1995 to 2018 were all greater than 0 and significant at the threshold level of 5%, indicating that the grain yield per unit sown area in the HHH was not randomly distributed but positively correlated. This indicates that the grain production layout showed strong spatial clustering characteristics. From 1995 to 1997, the Moran's I value decreased from 0.4114 to 0.1718. This indicates that the agglomeration and development trend in grain production in the HHH weakened during this period. From 1997 to 2002, the Moran's I value showed a fluctuating trend, ranging from 0.1718 to 0.3356. In this period, the grain production experienced a process of both agglomeration and dispersion development, but the change trend was small. From 2002 to 2006, the Moran's I value showed "up-down" fluctuation characteristics twice, and the change trend was more intense. This indicates that the grain production pattern of the HHH was greatly changed during this period. From 2007 to 2011, the increase in Moran's I value indicated that the grain production distribution in the HHH was in an increasingly concentrated state in this period. From 2011 to 2012, the Moran's I value changed from 0.2897 to 0.2401, a decrease by 0.0496 units, indicating that the clustering characteristics of grain production layout weakened during this period. The Moran's I value showed a continuous rising trend, indicating that the clustering characteristics of grain production distribution were enhanced from 2012 to 2014. The value of Moran's I changed from 0.3727 to 0.2315 in 2014-2018, indicating that the agglomeration characteristics of grain production distribution during this period weakened. On the whole, the Moran's I value experienced a change process of "down-up-down" in this period; however, the agglomeration and distribution characteristics of grain production did not change.

The evaluation of the global spatial correlation feature has the defect of ignoring the instability of local spatial processes. Therefore, the local spatial correlation characteristics of grain production in the HHH can be analyzed by observing the 

First, SPSS software was used to eliminate the collinearity of variables, and finally four indexes were extracted, namely, the effective irrigation area (EIA), amount of fertilizer application per unit of sown area (AFA), per capita annual income of rural residents (PCI), and annual average precipitation (PRE). In addition, 3.3 proves that the grain production pattern in the HHH has dependent characteristics; therefore, the influence of the spatial spillover effect cannot be ignored. Secondly, the GeoDa software was used to obtain the parameter estimation results of the spatial lag model, spatial error model, and OLS model. Finally, the statistical values of LMLAG and R-LMLAG in the OLS results were significantly higher than those of LMERR and R-LMERR at the 10% level; as such, it was appropriate to use the spatial lag model to explore the key factors affecting grain production (Table 4) .

According to Table 4 , the coefficients of EIA were 0.505, 0.415, 0.532, and 0.588 in 1995, 2005, 2015, and 2018, respectively, and were effective at the 1% threshold level. This indicates that there According to Figure 5 , the numbers of hot spots, sub-hot spots, sub-cold spots and cold spots were 13, 21, 21, and 4, respectively, in 1995, and the cold spots were mainly distributed in Sanmenxia, Luoyang, and Nanyang. The sub-cold spots were distributed in Beijing, Zhangjiakou, Chengde, Baoding, Jiaozuo, Xinxiang, Zhengzhou, Kaifeng, Zhoukou, Zhumadian, Xinyang, Fuyang, Bengbu, Huainan, and Bozhou. Sub-hot spots were concentrated in Tianjin, Qinhuangdao, Langfang, Cangzhou, Hengshui, Xingtai, Shijiazhuang, Dezhou, Binzhou, Dongying, Heze, Anyang, Hebi, and Puyang. Hot spots were distributed in the eastern part of the HHH. In particular, Yantai, Qingdao, Weifang, Jinan, Rizhao, Linyi, Zaozhuang, Jining, Lianyungang and Xuzhou were the most typical ones.

In 2005, the numbers of cities in the four categories were 16, 14, 16, and 13 . Compared with 1995, the number of cities in hot spots and cold spots increased, indicating that the agglomeration of grain production increased from 1995 to 2005. The cold spots were concentrated in the north and southwest of the HHH. The sub-cold spots were distributed in Qinhuangdao, Tangshan, Tianjin, Langfang, Baoding, Xinyang, Fuyang, Zhoukou, Xuchang, Zhengzhou, and Jiaozuo. The sub-hot spots continued to shrink, and were mainly distributed in Shijiazhuang, Cangzhou, and Hengshui. The spatial distribution of hot spots is clear, and the main distribution was in the middle and eastern part of the HHH. In 2015, the numbers of hot spots, sub-hot spots, sub-cold spots, and cold spots were 13, 20, 14, and 12 cities, respectively.

Compared with 2005, the number of cities in hot spots and cold spots decreased by 3 and 1, respectively, indicating that the local spatial clustering characteristics of grain production in the HHH were weakened from 2005 to 2015. In 2018, the hot spots, sub-hot spots, sub-cold spots, and cold spots included 10, 17, 22, and 10 cities, respectively. Compared with 2015, the number of cities in hot spots and cold spots decreased by three and two, and the number of cities in sub-hot spots increased by eight. This indicates that the local spatial agglomeration characteristics of grain production in the Huang-Huai-Hai Plain gradually weakened from 2005. In general, hot spots spread from the east of the HHH to the southeast and the central region, and sub-hot spots were mainly distributed in the central region. The cold spots and sub-cold spots were mainly distributed in the north and south regions.

First, SPSS software was used to eliminate the collinearity of variables, and finally four indexes were extracted, namely, the effective irrigation area (EIA), amount of fertilizer application per unit of sown area (AFA), per capita annual income of rural residents (PCI), and annual average precipitation (PRE). In addition, 3.3 proves that the grain production pattern in the HHH has dependent characteristics; therefore, the influence of the spatial spillover effect cannot be ignored. Secondly, the GeoDa software was used to obtain the parameter estimation results of the spatial lag model, spatial error model, and OLS model. Finally, the statistical values of LMLAG and R-LMLAG in the OLS results were significantly higher than those of LMERR and R-LMERR at the 10% level; as such, it was appropriate to use the spatial lag model to explore the key factors affecting grain production (Table 4) .

According to Table 4 , the coefficients of EIA were 0.505, 0.415, 0.532, and 0.588 in 1995, 2005, 2015, and 2018, respectively, and were effective at the 1% threshold level. This indicates that there was a positive correlation between the EIA and grain yield; that is, an increase in the EIA will bring about an increase in the grain yield. At the same time, the coefficient of EIA on the whole was on the rise, indicating that it has a stronger positive effect on grain yield. The continuous improvement of irrigation and water conservancy facilities in the HHH is the reason for this phenomenon.

The coefficients of fertilizer application per unit of sown area in 1995 and 2005 were 0.008 and 0.208, respectively, which were both significant at the 1% threshold level, and the coefficients in 2015 and 2018 were 0.124 and 0.023, respectively, and were significant at the 5% threshold level. This shows that the positive effect of chemical fertilizer application per unit of sown area on grain production experienced a changing process of decline after rising first, reflecting that the increase in chemical fertilizer application per sown area from 1995 to 2005 significantly increased the grain production in the HHH.

From 2005 to 2018, the boosting effect on grain production was alleviated. The reason for this phenomenon may be that the increase in the use of chemical fertilizers in the low-level agricultural development stage has a significant effect on increasing grain production. As the amount of fertilizer input remains at a high level, the effect of increasing the amount of fertilizer input in the future is limited as regards improving the grain yield. The coefficients of rural residents' per capita annual income in the four time sections were 0.405, 0.170, 0.124, and 0.050, which passed the significance tests at the critical value levels of 1%, 5%, 10%, and 10%, respectively.

The increase in per capita annual income continued to weaken the boost to grain production. The reason is that agriculture has been the main source of income for Chinese farmers to maintain their livelihoods for a long time. As farmers' income levels increase, they have more funds to purchase agricultural machinery, fertilizers, pesticides, seeds, and other production materials, thereby achieving the goal of increasing grain yield. With the rapid advancement of urbanization, farmers' income channels have become increasingly diversified, which has greatly reduced their dependence on agriculture, and the tendency of farmers to invest in non-agriculture has become more obvious. Therefore, to some extent, the increase in PCI has an inhibitory effect on grain output.

The coefficient of PRE was 0.508 in 1995 and passed the significance test of the 1% critical value level, indicating that the increase in PRE in that year played a promoting role in grain production. The coefficients of PRE were 0.016 and 0.148 in 2005 and 2015, both of which failed to pass the 10% significance test, indicating that the PRE increase did not have an obvious positive promotion effect on grain production. The underlying reason may be that with the continuous improvement of irrigation facilities, the impact of changes in PRE on grain production became weaker. *** means significant at the 1% threshold level; ** means significant at the 5% threshold level; * means significant at the 10% threshold level.

Based on the above research findings, the EIA, AFA, and PCI had a significant impact on the grain production in the HHH. However, SLM cannot explain the specific degree and scope of the impact of these three factors in space, and so a GWR model was explored to investigate the spatial difference in the influence of these three factors. The three factors-EIA, AFA, and PCI-were put into the model according to the criteria of AICc minimization.

The spatial distribution of Local R-squared values derived from the GWR model is displayed in Figure 6 . The GWR results show that three factors explain 81.1% of the variance in grain production. Geographic variations in these factors describe a difference in the combined statistical influence of the three variables on grain production across cities in the HHH, from a very weak relationship (near 0.30) to a strong relationship (>0.80). We found that 50.94% of cities maintained local R-squared values of more than 50%. The predictive power of the model shows characteristics of increasing from east to west. The local R-squared map suggests that the predictive power of the analysis was greatest in relation to the northern part (Chengde, Baoding, and Zhangjiakou) and the southwestern part (Jiaozuo, Jiyuan, and Sanmenxia). The lower R-squared values demonstrate a poorer regression fit in the eastern parts of the HHH, such as Weihai, Yantai, and Qingdao.

To explore the strength of the influence of each of the three factors on grain production, we created maps for each factor, which represent the geographic distribution of their regression coefficient values across HHH, according to the results of the GWR modeling. The mapped regression coefficients are divided into five classifications through Natural Breaks. In Figure 7 , white is used to indicate cities without data, and gradient shading is used to show cities with a significant relationship between variables and grain production. (Jiaozuo, Jiyuan, and Sanmenxia). The lower R-squared values demonstrate a poorer regression fit in the eastern parts of the HHH, such as Weihai, Yantai, and Qingdao. To explore the strength of the influence of each of the three factors on grain production, we created maps for each factor, which represent the geographic distribution of their regression coefficient values across HHH, according to the results of the GWR modeling. The mapped regression coefficients are divided into five classifications through Natural Breaks. In Figure 7 , white is used to indicate cities without data, and gradient shading is used to show cities with a significant relationship between variables and grain production.

The proportion of the effective irrigated area had a positive impact on grain production in the HHH, and its impact intensity presents a spatial distribution characteristic of "low west and high east". A total of 75.47% of the cities showed a significant positive relationship between the effective irrigate area and production, mainly in Hebei and Henan provinces. The reason may be that these two provinces are mostly inland, with relatively drier climate and less rainfall; therefore, irrigation is mainly used to supply the water requirement of crops.

AFA had a positive effect on 94.34% of the cities in the HHH, which can significantly increase the food production of 39.62% of the cities, mainly in the central, northern, and southern regions of the HHH, such as Zhangjiakou, Puyang, and Fuyang. AFA had a negative effect on 5.66% of the cities, mainly in the eastern part of the HHH, such as Qingdao and Weifang. The impact of AFA on grain production is limited. In particular, with the long-term investment of chemical fertilizers by Chinese farmers on cultivated land, the impact of various chemical fertilizers applied by farmers on soil fertility has become nearly saturated. The majority of the chemical fertilizers play a role in maintaining soil fertility after application, so even if the input of chemical fertilizers is increased, the positive effect on crop yield is not clear enough.

The per capita income of rural residents had a positive impact on 86.79% of the cities in the HHH, of which only 9.43% passed the significance test, mainly in the northeastern part of the HHH, such as Chengde, Qinhuangdao, Tangshan, Langfang, and Cangzhou. The per capita income of rural The proportion of the effective irrigated area had a positive impact on grain production in the HHH, and its impact intensity presents a spatial distribution characteristic of "low west and high east". A total of 75.47% of the cities showed a significant positive relationship between the effective irrigate area and production, mainly in Hebei and Henan provinces. The reason may be that these two provinces are mostly inland, with relatively drier climate and less rainfall; therefore, irrigation is mainly used to supply the water requirement of crops.

AFA had a positive effect on 94.34% of the cities in the HHH, which can significantly increase the food production of 39.62% of the cities, mainly in the central, northern, and southern regions of the HHH, such as Zhangjiakou, Puyang, and Fuyang. AFA had a negative effect on 5.66% of the cities, mainly in the eastern part of the HHH, such as Qingdao and Weifang. The impact of AFA on grain production is limited. In particular, with the long-term investment of chemical fertilizers by Chinese farmers on cultivated land, the impact of various chemical fertilizers applied by farmers on soil fertility has become nearly saturated. The majority of the chemical fertilizers play a role in maintaining soil fertility after application, so even if the input of chemical fertilizers is increased, the positive effect on crop yield is not clear enough.

The per capita income of rural residents had a positive impact on 86.79% of the cities in the HHH, of which only 9.43% passed the significance test, mainly in the northeastern part of the HHH, such as Chengde, Qinhuangdao, Tangshan, Langfang, and Cangzhou. The per capita income of rural residents had a negative impact on 13.21% of cities, and was concentrated in the southwestern region of HHH, that is, the western and southern regions of Henan Province, such as Sanmenxia, Nanyang, and Xinyang. As rural residents flow into developed urban areas, the non-agricultural income they earn from moving into cities has gradually become an important source of livelihood, and their dependence and emphasis on agriculture has gradually decreased; abandonment can even occur, which causes a huge negative impact on grain production. of HHH, that is, the western and southern regions of Henan Province, such as Sanmenxia, Nanyang, and Xinyang. As rural residents flow into developed urban areas, the non-agricultural income they earn from moving into cities has gradually become an important source of livelihood, and their dependence and emphasis on agriculture has gradually decreased; abandonment can even occur, which causes a huge negative impact on grain production. 

Even our results suggested that the effect of precipitation on grain production became weaker, we should also be aware of the relationship between rainfall and groundwater; that is, rainfall becomes groundwater through infiltration, providing sufficient water for irrigating crops. With the improvement in irrigation facilities, irrigation has gradually become an indispensable and important means for stable grain production. This may also be the reason why the direct impact of PRE on grain production is gradually weakening and why the EIA is gradually increasing. Future research should focus on related research in this area. The reasons for the insignificant effect of temperature and sunshine duration may be attributed to two points: first, compared with precipitation, the heat resources of the HHH can well meet the growth needs of winter wheat and summer corn, and the yield is less affected by temperature and sunshine duration. At the same time, an increase in temperature and sunshine duration may increase evaporation, thereby offsetting the effect of precipitation. It is also possible to attribute this effect to precipitation rather than temperature and sunshine duration. Second, it may be that the crops in the research area are not subdivided. Different crops have inconsistent requirements for temperature and sunshine duration, which may weaken the effects of temperature and sunshine duration. This will also be the focus of future research. China's food self-sufficiency rate has reached more than 95%. The concept of "with grain in the hand, the heart is unharried" has been made reality. However, in the face of the impact of global warming and COVID-19, as well as the requirements of new urbanization, as a country with one of the largest populations in the world, ensuring the stability of China's food supply is not only related to national security, but is also related to the stability of the world.

Therefore, based on the results of this paper, the following policy suggestions are put forward to increase grain production in the HHH, an important grain production base for food security in China to maintain China's food security. The suggestions include the following: to increase the construction investment for basic farmland infrastructures, such as irrigation facilities; to cultivate and promote good varieties and treatments, and implement soil testing and formula fertilization; to standardize the rural land market; to promote the transfer of rural land in an orderly manner; and to realize the large-scale management of cultivated land. 

Even our results suggested that the effect of precipitation on grain production became weaker, we should also be aware of the relationship between rainfall and groundwater; that is, rainfall becomes groundwater through infiltration, providing sufficient water for irrigating crops. With the improvement in irrigation facilities, irrigation has gradually become an indispensable and important means for stable grain production. This may also be the reason why the direct impact of PRE on grain production is gradually weakening and why the EIA is gradually increasing. Future research should focus on related research in this area. The reasons for the insignificant effect of temperature and sunshine duration may be attributed to two points: first, compared with precipitation, the heat resources of the HHH can well meet the growth needs of winter wheat and summer corn, and the yield is less affected by temperature and sunshine duration. At the same time, an increase in temperature and sunshine duration may increase evaporation, thereby offsetting the effect of precipitation. It is also possible to attribute this effect to precipitation rather than temperature and sunshine duration. Second, it may be that the crops in the research area are not subdivided. Different crops have inconsistent requirements for temperature and sunshine duration, which may weaken the effects of temperature and sunshine duration. This will also be the focus of future research. China's food self-sufficiency rate has reached more than 95%. The concept of "with grain in the hand, the heart is unharried" has been made reality. However, in the face of the impact of global warming and COVID-19, as well as the requirements of new urbanization, as a country with one of the largest populations in the world, ensuring the stability of China's food supply is not only related to national security, but is also related to the stability of the world.

Therefore, based on the results of this paper, the following policy suggestions are put forward to increase grain production in the HHH, an important grain production base for food security in China to maintain China's food security. The suggestions include the following: to increase the construction investment for basic farmland infrastructures, such as irrigation facilities; to cultivate and promote good varieties and treatments, and implement soil testing and formula fertilization; to standardize the rural land market; to promote the transfer of rural land in an orderly manner; and to realize the large-scale management of cultivated land.

However, the impact of relevant agricultural policies issued by the country also needs to be considered in the future. From 1995 to 2018, the overall increase in grain production in HHH and the gradual reduction in spatial differences well reflects the background and national policies of the country in different periods. Since 1995, with the advancement of agricultural technology, the agricultural development of the HHH has been weakened by natural conditions, and the grain production has increased significantly [13] . At the same time, due to the popularization of fertilizers, pesticides, and irrigation, the gap in grain production in various regions in HHH is also narrowing.

Due to rapid urbanization, the conversion of fertile irrigated land to non-agricultural land seems to pose a potential threat to the food security of the HHH, and even to the whole of China [10] . Moreover, farmers are gradually migrating to cities in search of higher incomes due to the urban-rural development gap [44] , which has caused the arable land in the HHH to be abandoned, resulting in the expansion of spatial differences in grain production in different regions in the HHH, and also threatening national food security. In order to ensure national food security, the central government proposed the construction of "high-standard basic farmland projects" and "agricultural modernization" to promote the large-scale, intensive, and modernized management of arable land in the HHH, thereby increasing the grain production of the HHH, and promoting the development of sustainable agriculture.

At the same time, we should not ignore the resistance, created by resource endowments, to sustainable agriculture in the HHH. In particular, the precipitation cannot meet the water demand of the crops [8, 45] , and water shortage is one of the major factors threatening the high and stable production of wheat [20] . The water consumption greatly exceeds the precipitation, and groundwater must be extracted to make up for the deficiency so as to maintain high yields [21] . Some areas in the HHH even appear salinized due to unreasonable irrigation [10] , which poses huge challenges for the minimization of environmental impacts and the development of sustainable agriculture [46, 47] . Therefore, the government must not only build irrigation facilities, but more importantly, must promote water-saving irrigation technologies, improve water resource utilization efficiency, implement drip irrigation and sprinkler irrigation [46] , etc., so as to achieve sustainable agricultural production in the HHH. At the same time, in the long run, in the face of the encroachment of arable land in the promotion of urbanization, the amount of arable land in the future will also face severe challenges [10] . Promoting intensive agricultural production and improving the level of intensive use of agricultural production are also of great significance to ensure future food production [48] . Moreover, these actions can also enhance the ability to respond to natural disasters in the future. However, a sustainable food and agriculture system is one which is environmentally sound, economically viable, socially responsible, nonexploitative, and which serves as the foundation for future generations [49] [50] [51] . With the long-term development of intensive agriculture production in the Huang-Huai-Hai Plain, agricultural practices ranging from the development of irrigation projects to the use of agrichemicals have often had negative environmental impacts, such as wildlife kills, pesticide residues in drinking water, soil erosion, groundwater depletion, and salinization [52] . Substituting environmentally sound inputs for those which are damaging is an important step in addressing these problems [49] . In view of this, the Ministry of agriculture of China started the construction of the Key Laboratory of Agricultural Environment in the Huang-Huai-Hai Plain, with the objectives of scientific research, environmental monitoring, detection analysis and technical services, in 2012, aiming to carry out research on regional agricultural pollution prevention by means of agricultural non-point source pollution prevention, the environmental protection of producing areas, and the development and application of environment-friendly inputs. However, for the farmers who are the main body of agricultural production, whether these agricultural technology inputs will increase agricultural production costs and reduce agricultural income will be an important factor affecting the promotion of agricultural technology and the development of sustainable agriculture. Therefore, whether the economic, social and environmental benefits generated by agricultural production in the Huang-Huai-Hai Plain under the influence of agriculture technology can achieve a balance will be the focus of our future research.

In this paper, exploratory spatial data analysis, the gravity center model, and the spatial lag model were used to explore the spatial-temporal variation and influencing factors of the grain production pattern in the HHH from 1995 to 2018. The main conclusions were drawn as follows:

The grain production pattern in the HHH has the characteristics of being non-equilibrated in space and non-stationary in time. The spatial non-equilibrium is reflected in the shift of the grain production center from the southeast to the northwest of Tai'an city. The high-level areas of grain production capacity were mainly distributed in the east and south, while the low-level areas were distributed in the west and north. The non-stationarity of time is reflected in the rising trend in the grain production capacity and the weakening of the non-stationarity of time in the grain production center over time;

The global and local spatial agglomeration characteristics of grain production in the HHH were significant. The global spatial correlation characteristics underwent a "decrease-growth-decrease" change process, and the local spatial correlation characteristics demonstrated a concentrated distribution. Specifically, the hot spots were mainly distributed in the central and eastern regions of the HHH, and the cold spots were distributed in the north and southwest. The global and local spatial autocorrelation characteristics showed that the polarization of grain production in local areas has gradually weakened and the spatial difference has gradually decreased in the HHH, which indicates that its agricultural production has gradually shifted in the direction of sustainable development;

The impact of social-economic factors on grain production was constantly strengthened and the influence of climate factors on grain production was gradually weakened. EIA, AFA, and PCI helped to increase the grain yield per unit of sown area in the HHH; however, the effect of the PRE on grain production became weaker as time went on. We adopted the GWR model to prove that the EIA, AFA, and PCI had clear spatial heterogeneity in the intensity and direction of the local scale. The results showed that the EIA had a larger impact on grain production in the HHH compared with other factors, with the percentage of significance at 75.47%.

Study on grain self-sufficiency rate in China: An analysis of grain, cereal grain and edible grain

Food demand and the food security challenge with rapid economic growth in the emerging economies of India and China

Poverty-Food Security Nexus: Evidence from a Survey of Urban Slum Dwellers in Kolkata

Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems

Spatial shifts in grain production increases in China and implications for food security

The grain potential of cultivated lands in Mainland China in 2004

Spatial and temporal differentiation of COVID-19 epidemic spread in mainland China and its influencing factors

Temperature and Precipitation Suitability Evaluation for the Winter Wheat and Summer Maize Cropping System in the Huang-Huai-Hai Plain of China

Water consumption in summer maize and winter wheat cropping system based on SEBAL

Changes in quantity and quality of cropland and the implications for grain production in the Huang-Huai-Hai Plain of China. Food Secur

Strip rotary tillage with subsoiling increases winter wheat yield by alleviating leaf senescence and increasing grain filling

Summer maize growth under different precipitation years in the Huang-Huai-Hai Plain of China

Influencing factor decomposition of grain production at county level in Huang-Huai-Hai region based on LMDI

A Study on the Driving Factors of Food Production in Huang-Huai-Hai Plain Based on Path Analysis

Analysis on spatial spillover effect and influence factors of grain yield per hectare at county level in Huang-Huai-Hai region

Coupling Relationship Between Land Use Transitions and Grain Yield in the

The relationship between grain output and fertilizer input in wheat-corn cropping area of the Huang-Huai-Hai plain

Wheat yield benefited from increases in minimum temperature in the Huang-Huai-Hai Plain of China in the past three decades

Contributions of climatic and crop varietal changes to crop production in the North China Plain, since 1980s

The impacts of climate change on wheat yield in the Huang-Huai-Hai Plain of China using DSSAT-CERES-Wheat model under different climate scenarios

ETc of winter wheat across HuangHuai-Hai Plain in the future with the modifed DSSAT model

Regional Differences of Winter Wheat Phenophase and Grain Yields Response to Global Warming in the Huang-Huai-Hai Plain in China Since 1980s

Impact of thermal time shift on wheat phenology and yield under warming climate in the

Spatio-temporal variations in precipitation on the Huang-Huai-Hai Plain from 1963 to 2012

Projection of extreme precipitation in the context of climate change in Huang-Huai-Hai region

Agricultural Vulnerability and Adaptation to Climatic Changes in Malaysia: Review on Paddy Sector

An Integrated Analysis of Impact Factors in Determining China's Future Grain Production

The spatio-temporal change of agricultural hydrothermal conditions in China from 1951 to 2010

Analysis of sensitivity of main grain crops yield to climate change since 1980 in Henan Province

Spatial Data Analysis in the Social and Environmental Sciences

Spatial-temporal patterns of PM2.5 concentrations for 338 Chinese cities

Methods and Models

The interpretation of statistical maps

Local spatial autocorrelation statistics: Distributional issues and an application

Evolution of Patterns in the Ratio of Gender at Birth in Henan province

Spatial-Temporal Evolution of PM 2.5 Concentration and its Socioeconomic Influence Factors in Chinese Cities in 2014-2017

Spatial differentiation of daily travel carbon emissions in small-and medium-sized cities: An empirical study in Kaifeng

A Center-of-Gravity-Based Approach to Estimate Slow Power and Frequency Variations

Applied Spatial Econometrics: Raising the Bar

Evaluation of Parameter Estimation Methods for Fitting Spatial Regression Models

Hedonic Housing Prices In Paris: An Unbalanced Spatial Lag Pseudo-Panel Model With Nested Random Effects

A Computer Movie Simulating Urban Growth in the Detroit Region

Identifying the socioeconomic determinants of population exposure to particulate matter (PM 2.5 ) in China using geographically weighted regression modeling

What Factors Influence Rural-To-Urban Migrant Peasants to Rent out Their Household Farmland? Evidence from China's Pearl River Delta

Water-yield relations and water-use efficiency of winter wheat in the North China Plain

Water management and crop production for food security in China: A review

Agricultural sustainability and intensive production practices

Governing farmland conversion: Comparing China with the Netherlands and Germany

Expanding the Definition of Sustainable Agriculture. Sustain

What is Sustainable Agriculture? A Systematic Review

Agricultural sustainability: Concepts, principles and evidence

Sustainable Intensification in Agriculture: Premises and Policies

Author Contributions: C.Z. and X.N. developed the main ideas of the study, gathered the data, performed the model construction and estimation, and wrote the manuscript. R.Z. helped to collect and process data. Z.Z. participated in revising the manuscript and proofreading the article. All authors have read and agreed to the published version of the manuscript. 

The authors declare no conflict of interest.

Publisher's Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.