Agricultural non-point source pollution (ANSP) has been a long-term research hotspot in academia. In recent years, with the rapid development of agricultural economy and the increase of food crop production for 12 successive years, China has become a large user of pesticides and fertilizers among the developing countries in the world (He, 2020). The consequences of inappropriate behaviors in agricultural production, such as leaching of pesticides and chemical fertilizers, excretion of livestock and aquaculture industry, and the disordered discharge of rural domestic garbage, not only gradually eroded the agricultural ecological environment, but also had a huge impact on food security (Wittman et al., 2017). This suggests the importance of effective ANSP prevention in the rapid growth of agricultural economy for China’s long-term development. Therefore, more research efforts are being devoted to the ANSP patterns and its reduction strategies.

However, it is a challenge to measure ANSP and understand the mechanism of its emission due to the complexity in its formation, the uncertainty of emission direction and the difficulty of ANSP monitoring and control (Zhang et al., 2004; Collins et al., 2016; Sidemo-Holm et al., 2018). In the past decade, scholars have paid intensive attention to the calculation of ANSP. In the related research, pollution source analysis is the focus of pollution calculation. Emergy and economic analyses were used to evaluate the agricultural production in Tongxiang region, including traditional rice monoculture, integrated farming, and non-grain production systems (Su et al., 2020). An integrated assessment system including the modified DRASTIC model (standing for Depth to the water table, Recharge, Aquifer media, Soil media, Topography) was developed to evaluate the pollution caused by agricultural system (Wu et al., 2020). Energy dispersive X–ray fluorescence spectroscopy was used to measure the concentration levels of heavy metals in soil samples in Turkey, which showed that Cr, Ni, Zn, As, and Pb levels were highly influenced by agricultural practices (Baltas et al., 2020). The SWAT (Soil and Water Assessment Tool) model was used to link between pollution input in the upstream watershed and pollutant load response at the watershed outlet to identify critical source areas at the lake basin scale (Shang et al., 2012). Equivalent standard pollution loads method was used to estimate the amount of total nitrogen and phosphorus loss from agricultural sources in Shandong Province (Gao et al., 2010). These studies have provided effective reference for ANSP measurement. However, due to late start of pollution research in China and lack of extensive and long-term monitoring data, previous research is mostly concentrated on basin, scales, fixed types of pollution, on-site investigations and regional trials in the measurement process (Hu et al., 2015;Ma et al., 2015; Lu et al., 2017). Such research results can not significantly reflect the complexity of pollution, resulting in deviations in ANSP explanation at the national or regional level and thus lack of proof for the formulation and implementation of policies related to pollution prevention and control at this level. The inventory analysis method is applicable to a large research area. For on hand, the basic data of the inventory analysis method comes from public statistics. For another hand, the key coefficients of analysis methods are mainly determined by the comprehensive research of Department of Environmental Science and Engineering of Tsinghua University and the Provincial data of National First Pollutant Source Census Agricultural Source Coefficient Manual. Specifically, the inventory analysis method measures the pollution emission unit by unit through establishing the response relationship between agricultural activities and pollution discharges (Ongley et al., 2010). When compared with other methods, the inventory analysis method can more comprehensively reflect the discharge of various types of ANSP, suggesting its potential application value (Lai et al., 2004).

Currently, research attention has been focused on the inventory analysis methods, such as quantitative analysis of ANSP in terms of spatial distribution, emission efficiency, impact mechanism, and ANSP prevention and control strategies, as well as the relationship of ANSP with the level of economic agglomeration and the internal structure of agriculture (Ge and Zhou, 2011; Rao et al., 2011; Pan and Ying, 2013; Qiu et al., 2018). However, few reports are available about the use of the inventory analysis method to measure China’s ANSP emission and further explore its temporal and spatial changes. Most previous analyses of China’s ANSP at the spatial level paid little attention to its long-term variation law or its convergence characteristics, which, however, are important to reflect the polarization, imbalance and spatial distribution of economic indicators. Furthermore, they failed to consider the obvious clustering development characteristics of China’s agricultural economy and non-point source pollution at the dimensional level. With the development of new economic geography and spatial analysis, the spatial econometrics, which focuses on the potential spatial connection among regional units, could solve the problem inaccurate estimation of marginal effects by traditional models and improve the measurement of evolution and judgment of the influencing mechanism of the ANSP (Lesage and Pace, 2007). Therefore, the long-term variation law and the spatial convergence characteristic of the ANSP should be explored to further reveal the actual situation and the influencing mechanism of the ANSP.

In view of this, this paper uses the inventory analysis method to measure the ANSP emission intensity of China and its provincial units from 1999 to 2017, and analyzes its spatial convergence characteristics. Furthermore, a spatial econometric model is constructed to study the factors influencing ANSP. Finally, the evolution law of ANSP in China is explored. The aim is to accurately evaluate the level of ANSP emission intensity of China, clarify its impact mechanism, and provide the basis for the formulation and optimization of ANSP prevention strategies in China.

The intensity of ANSP can reflect the degree of agricultural pollution accumulation per unit land area, suggesting that it can be used to represent the level of regional ANSP (Li et al., 2017a). The research has shown that ANSP can be calculated by inventory analysis method using the main pollution-producing units shown in Table 1 (Chen et al., 2006). As is shown in Fig. 1, the pollution-producing units use precipitation as the carrier, enter the water network through the process of surface runoff and underground infiltration, and cause agricultural non-point source pollution (Wang, 2009).

Figure 1.  Accounting framework of the agricultural non-point source pollution

The ANSP emission intensity can be calculated by the following formula:

where y represents the intensity of ANSP (kg/ha); E, the total emission of ANSP (kg); Area, the area of agricultural land, (ha). q, the number of pollution-producing units; EL_m, the pollution-producing unit of different categories; m, the pollutant production unit; S_m, the pollution intensity coefficient of unit m; η_m, the coefficient of regional resource utilization and management of unit m; C_m, the emission coefficient of unit m after considering regional resource utilization and management factors. Among them, the pollution intensity and emission coefficient are determined by the comprehensive effects of regional environment, resource utilization and management characteristics, according to the comprehensive ‘National First Pollution Source Survey Agricultural Source Coefficient Manual’ and correlation coefficients of studies by Liang and Ma, etc. (Liang, 2009; Ma, 2013). PE_m is the product of EL_m and S_m, which indicates the potential pollution of each pollution-producing unit when the external factors are ignored.

Spatial dependence test is one of the preconditions for spatial analysis of ANSP. It can be divided into two categories: global spatial dependence test and local spatial dependence test.

Specifically, the global spatial dependence test is used to analyze the spatial dependence of the overall regional system. In this study, we use the Moran’s I index for testing as shown by the following formula (Anselin, 1988):

where I, is the Moran’s I of the intensity of ANSP; y_i and y_j are the intensity of ANSP of the adjacent regions i and j, respectively; $\overline y $, the average of intensity of ANSP; W_ij, the spatial weight matrix; n, the statistics of the study region. The value of I ranges from –1 to 1. When I is close to 1, y exhibits a positive overflow effect in space; when I is close to 0, y is spatially independent of each other; when I is close to –1, y is spatially adjacent to each other. Since the choice of spatial weight matrix is exogenous, this paper investigates the spatial dependence of ANSP emission intensity in China by combining Rook Contiguity weights (r), K-Nearest Neighbors weight (k) and Euclidean Distance weight (d) according to the research of Anselin (1988). Rook Contiguity weights used only common boundaries to define neighbors, and r₁, r₂, r₃ are used to investigate different order weights.

The significance level of Moran’s I can be tested by Z in the following formula:

Where E(I) and VAR(I) are the expected values and variances of Moran’s I, respectively. When |Z| > 1.96, the spatial characteristics are verified to be significant.

Furthermore, the Moran scatter plot of local spatial dependence test was used to analyze whether there is a ‘club’ spatial association of ANSP. The research unit located in the first quadrant of the scatter plot represents the situation of (high, high). The second quadrant represents the situation of (low, high), the third quadrant represents the situation of (low, low), and the fourth quadrant represents the situation of (high, low) (Herrerias et al., 2013).

Convergence analysis can be divided into two parts: α convergence and β convergence. The former refers to the convergence of horizontal quantities, indicating that the intensity gap of ANSP emissions in various regions has narrowed over time, while the latter refers to the convergence of growth rate. When the growth rate is higher in the low pollution area than in the high pollution area at the initial period and tends to present a common equilibrium state, the absolute β convergence exists; when the regional characteristics are considered, and the pollution in a region tends to be in its own stable state, the conditional β convergence exists. Only when the β convergence exists, high-value areas and low-value areas tend to be in a common steady state. Therefore, β convergence is the premise of α convergence. In β-convergence analysis, the benchmark for absolute β-convergence comparison is based on the pollution growth rate in other regions, while the comparison criterion for conditional β convergence is the equilibrium state of pollution growth rate in their own region. This paper focuses on the analysis of absolute β convergence, with a purpose to investigate the pollution growth trends in different regions.

The absolute β convergence can be measured by the classical or spatial measurement methods developed based on the neoclassical growth model of Barro (Barro and Sala-i-Martin, 1992). Later, scholars found the defects in the classical measurement methods at the technical level. For example, the macroeconomic research failed to consider the interaction of economic activities in different regions, thus neglecting the role of spatial elements in convergence and leading to bias in estimation results. In this research, a spatial measurement method is used to investigate the convergence mechanism of ANSP intensity in China, followed by comparison of its results with those of classical measurement methods (standard convergence analysis, hereafter).

According to Bernard (Bernard and Jones, 1996), the standard convergence analysis model can be established as follows:

where, i indicates the province, municipality or autonomous region, u represents the error term, and $\Delta {\rm{ln}}\left({{y_{iT}}} \right)$ is the growth rate of the ANSP emission intensity of the province i during the T period after taking the logarithm. The convergence coefficient α₀ and α₁ can be used for calculation of convergence rate κ and the convergence half life cycle to the steady-state value τ. Among them, the convergence speed is $\kappa = - \left({1/T} \right){\rm{ln}}\left({1 + {\alpha _1}} \right)$, and if α₁ < 0 and passes the significance test, the absolute convergence exists. The time span T in this paper is set to 1 year, and the half life cycle of convergence to the steady state value is $\tau = \left({{\rm{ln}}1/2} \right)/{\rm{ln}}\left({1 - \kappa } \right)$.

However, the standard convergence model does not consider the spatial effect and cannot reveal the effect of economic activities in adjacent areas on the convergence mechanism of ANSP intensity. According to the theory of space economy, dependence and difference between adjacent regions should be considered in economic activities and they can be evaluated using the Spatial Lag Model (SLM) and Spatial Error Model (SEM) (Lesage and Pace, 2007). Correspondingly, the spatial measurement of absolute β convergence can be set as follows:

where, W is a spatial weight matrix; ν and ε, error vectors; ρ, the spatial lag coefficient; λ, the spatial error coefficient, y, T, and α share the same definition of formula (4). $\varepsilon :N(0,{\sigma ^2}{I_n})$ indicates that ε is considered to obey the Gaussian distribution. When ρ = 0 and λ = 0, the model is set to the standard convergence model; when ρ ≠ 0 and λ = 0, the model is set to SLM; when ρ = 0 and λ ≠ 0, the model is set to SEM; If ρ and λ ≠ 0 at the same time, the space factor should be considered. In the model test process, the Moran’s I value of the model residual is a preliminary criterion for discriminating the model type. If the residual of the spatial regression model is significantly auto-correlated, the model cannot well reflect the spatial correlation, and LM test should be further used to analyze whether the spatial lag term or the spatial error term exists. If the LM test results pass the 5% significant level, the sample should be analyzed by using the spatial measurement model.

(1) Spatial measurement model

The spatial effect model was used to analyze the spatial effects of various factors on the ANSP intensity. Spatial interactions can generally be characterized in three forms: Spatial Durbin Model (SDM), SLM and SEM. In this study, the sample is first assumed to fit the SDM model as follows:

where Y is the ANSP intensity; WY, the hysteresis term of Y; X, the explanatory variable of Y; n, the number of research units; ε, the random error term; W, the spatial distance weight matrix (Lesage and Pace, 2007); β, ρ and θ, the spatial regression coefficients, respectively.

The validity of model setting is evaluated by using the LM test to determine whether a spatial lag term and a spatial error term exist. If the LM-test result passes the test at a significance level of 5%, then Wald-test and LR-Test are used to determine if SDM is applicable. When $ {H}_{0}{}^{1}: \theta =0$, formula (6) is simplified to SLM; when $ {H}_{0}{}^{2}: \theta +\rho \theta =0$, formula (6) is simplified to SEM; if the null hypothesis is rejected, SDM estimation needs to be further performed.

Due to the introduction of the spatial distance weight matrix, the regression coefficients in the model cannot directly reflect the influence of various factors on the ANSP intensity. Thus, the total effect, direct effect and indirect effect of each explanatory variable need to be further calculated by using the formulas rewritten from formula (6) according to the research of Lesage (Lesage et al., 2007).

(2) Variable description

According to relevant research, the intensity of ANSP can be affected by factors such as agricultural economic scale, agricultural internal structure, precipitation and runoff, rural infrastructure investment, agricultural research and development (R&D) investment, and labor-force quality (Ge and Zhou, 2011; Liu et al., 2013; Lu et al., 2018). In this study, these previously verified influencing factors are introduced into the model as explanatory variables for empirical research of the spatial effects of factors on ANSP and these influencing factors are presented in Table 2.

The research area of this paper is China’s 31 provinces (municipalities, autonomous regions, excluding Hong Kong, Macau and Taiwan). The study period is from 1999 to 2017. The data used in this paper are mainly obtained from China Statistical Yearbook (NBS, 2000–2019a), China Statistics Yearbook on Environment (NBS and MEE, 1998–2018), China Population and Employment Statistics Yearbook (NBS, 2000–2019b), China Rural Statistical Yearbook (NBS, 2000–2019c), National Agricultural Science and Technology Statistics Collection^①, and statistics of various provinces, municipalities and autonomous regions in China; precipitation data and hydrological data are derived from the National Climate Center of China National Meteorological Administration and the Data Center of the Ministry of Water Resources. The price data are converted to comparable prices in 1999.

Based on the panel data of China’s 31 provinces, municipalities and autonomous regions (excluding Hong Kong, Macau and Taiwan) from 1999 to 2017, this paper uses the inventory analysis method to measure the ANSP emission intensity and the results are shown in Table 3.

In Table 3, it is shown that: 1) from 1999 to 2017, the ANSP intensity exhibits an upward trend, with its arithmetic mean increased from 134.55 kg/ha in 1999 to 139.74 kg/ha in 2017. 2) The ANSP emission intensity varies in different provinces and cities in each period. From 1999 to 2007, the ANSP emission intensity shows a decrease in most provinces and cities, followed by an increase from 2007 to 2013 and a decrease from 2013 to 2017. 3) The gap in ANSP intensity has been narrowing over time in various provinces and cities. The difference between the maximum and minimum values was dropped from 527.97 kg/hm² in 1999 to 429.26 kg/ha in 2017, but the ANSP intensity was still high in each period. The distribution of high and low values is relatively stable. The high values are mainly concentrated in Shanghai, Tianjin, and Jiangsu while the low values are mainly distributed in Inner Mongolia, Qinghai, and Xinjiang, indicating the potential connection of ANSP with economic-geographic spatial distribution. From the perspective of the temporal and spatial changes, the upward trend of the ANSP and the different values of spatial distribution pattern are basically consistent with the previous studies (Chen et al., 2006; Liang, 2009; Qiu et al., 2018).

The global spatial autocorrelation results are shown in Table 4 and three tendencies are obtained as follows: 1) China’s ANSP intensity shows a significant positive spatial dependence. Most of the Moran’s I values are positive and have passed the 5% significance test, indicating that the ANSP intensity is positively correlated between adjacent areas, and changes of this intensity in neighboring provinces will affect the ANSP in the study province, leading to changes in emission intensity. 2) The spatial dependence of China’s ANSP intensity shows a gradual increasing trend in the time dimension. Regardless of the adjacency relationship, the Moran’s I index basically shows a trend of ‘decline-rise-fall’. 3) The spatial dependence of China’s ANSP intensity is gradually weakened with spatial expansion with the highest spatial positive correlation for r₁. Under the ‘adjacent neighbor relationship’, the spatial dependence is significant for r₂. However, the degree of dependence is significantly weakened and r₃ shows no statistically significant spatial dependence.

Fig. 2 shows the local indicators of spatial association (LISA) cluster for China’s ANSP intensity of China in 1999 and 2017. In 1999, 10 provinces and autonomous regions were shown to have passed the significant level test, including 3 (high, high), namely Jiangsu, Zhejiang and Anhui, and 7 (low, low), namely Inner Mongolia, Ningxia, Qinghai, Xinjiang, Tibet, Sichuan, and Gansu. In 2017, the local spatial pattern changed slightly. Among the provinces and autonomous regions that passed the significance level test, Anhui was excluded from the list (high, high) and Shaanxi was included in the list (low, low) compared in 1999. Collectively, two stable ‘clubs’ were formed in China for ANSP intensity, the low-intensity club in the northwest region and the high-intensity club in the eastern coast.

Figure 2.  Local Indicators of Spatial Association (LISA) cluster for China’s agricultural non-point source pollution intensity in 1999 (a) and 2017 (b); data excluding Hong Kong, Macau and Taiwan

Based on the above analysis, the overall spatial and local spatial dependence of China’s ANSP intensity is roughly the same, with the main form being the local ‘club’ spillover. But, how does this spillover influence the interpretation of the intensity changes in agricultural non-point source pollution? This question was further investigated by performing spatial metrology analysis.

The spatial measurement model was used to analyze the convergence intensity of ANSP in China from 1999 to 2017, and the results of standard convergence model analysis were also listed for comparison. According to Anselin’s research (Anselin, 1988), the residuals of the model without considering spatial factors were tested by using the Moran’s I index and the optimal model was determined by LM test to obtain the simplified formula (6). In Table 5, it can be seen that the Moran’s I index of the regression residual of the ANSP intensity is significant in the standard convergence model, indicating that the spatial factors should be considered in convergence analysis. In the spatial error model (SEM), the error coefficient λ is not significant, and the residual Moran’s I index passes the spatial correlation test, implying that SEM does not conform to the basic assumptions of the spatial model. Besides, the estimated coefficients of SLM are all significant, and the ρ value is not 0 under the significant test, indicating that the spatial correlation of the error term u is reasonable in the assumption model. Meanwhile, the Moran’s I index of SLM residual is not significant, indicating that the residual term of the spatial lag model does not have a spatial correlation and SLM can reflect the authenticity of the spatial convergence mechanism of ANSP. Furthermore, R² is higher in SLM than in other models, even the Robust LM values of each model did not pass the significance test, and the LM value was highly significant, which indicate that the spatial correlation of the residuals may be derived from the hysteresis correlation. Therefore, the SLM is defined as the optimal model setting.

Table 5 shows the absolute β convergence results, which can reflect the spatial lag correlation characteristics of ANSP in China. It can be seen that the growth rate of regional ANSP intensity is not only affected by the emission level of its initial period, but also by the initial levels of surrounding areas. The α₁ value is negative and significant, indicating an absolute convergence for the ANSP emission intensity under the effect of spatial hysteresis. At the same time, the spatial analysis revealed a convergence rate of 0.0064, indicating an increase of 0.64% versus 0.0060 from the ordinary convergence analysis. The convergence half life cycle is shortened from 115 yr to 108 yr, indicating the spatial factor plays an obvious role in the change of ANSP intensity. Additionally, a comparison with the convergence cycle of China’s regional economic growth found that the semi-life cycle of ANSP convergence is longer than that of the economic growth, which was assumed to be 87 yr by Pan (2010). Collectively, in China, the average rate is lower in shortening the gap between low-emission areas and high-emission areas in ANSP intensity than in that between underdeveloped areas and developed areas in economic development.

Due to significant geospatial characteristics of China’s agricultural development, the adjacent provinces in each region have strong convergence at agricultural internal structure and economic development level. The interaction of labor, capital and other factors in space, the demonstration role of the agricultural production input behavior and the competition behavior between local governments have also caused the complexity and diversity in the factors affecting pollution emission. Since no single factor can fully explain the spatial characteristics of ANSP emission intensity in the real economic society, spatial analysis has been used to study the influencing factors of pollution in the present study. For instance, the spatial dependent items of independent variables are incorporated into the spatial econometric model, and the multi-collinearity problem of variables no longer becomes a misplaced factor in the model, thus overcoming the invalidity or bias of the traditional econometric model (Li et al, 2010).

Firstly, the Hausman Test rejected the assumption of no difference between fixed and random effects, and random-effects model was fitted. As shown in Table 6, the LM test rejects the assumption of no spatial lag terms and spatial error terms in the explanatory variables at the 1% significance level, indicating that a spatial measurement model should be established. Furthermore, both the Wald test and the LR test of the model have passed the test of 1% significance level, indicating that the SDM setting is suitable for the sample data. The comparison of R² indicated that the fitting effect of SDM is improved versus that of OLS. In terms of variable fitness, most variables are highly significant except for the scale of agricultural economy, lagging precipitation and runoff, lagging agricultural infrastructure investment, and lagging agricultural R&D investment. The spatial autoregressive coefficient ρ is significantly positive, indicating the spatial correlation of these factors with ANSP exists in 31 provinces, municipalities and autonomous regions in China. Therefore, the external influencing factors on ANSP cannot be ignored in the research.

Furthermore, direct and indirect effects of various factors on the intensity of ANSP discharge were analyzed and the calculation results are shown in Table 7. In the estimation of the direct effect, most variables have passed the test at the 5% level except for the variable of agricultural R&D input. In the estimation of the indirect effect, the agricultural population size and the labor-force quality have both passed the test at the 5% level.

The coefficients for direct, indirect and total effects of the agricultural economy scale are 0.1322, –0.0064 and 0.1259, respectively. Consistent with the general assumptions, the direct effect is a significant reflection of the ANSP growth caused by the expansion of agricultural economy. The indirect effect is negative but not significant, indicating that the agricultural economic development is competitive between regions. With the expansion of agricultural economy in the study region, a “siphon effect” can be easily produced in neighboring regions, which may promote the flow of agricultural resources, and reduce agricultural chemical inputs in adjacent regions, but this negative effect is not obvious.

The direct, indirect and total effect coefficients of the internal structure of agriculture are –0.0594, –0.0271 and –0.0865 respectively. The direct effect is negative and significant, indicating a negative correlation between the increase in the proportion of planting industry and its ANSP intensity in the region. Generally, with a reduction in the proportion of ANSP-intensive industries, the emission intensity will be relatively reduced. According to the agronomic characteristics in the agricultural industry, an increase in the proportion of planting industry and a decline in the proportion of aquaculture suggest a lower pollution intensity. The negative indirect effect indicates that an appropriate adjustment in the agriculture internal structure would promote ANSP reduction in neighboring regions, for one region can learn from other regions in pollution reduction mode and effectively enhance the pollution reduction effect in agricultural restructuring.

The coefficients for direct, indirect and total effects of precipitation and runoff are significantly positive, but the total effect is not significant, indicating that the increase in precipitation will not only increase the ANSP emission intensity in the study region, but also produce the same effects in adjacent areas. Precipitation and runoff is one of the important natural factors affecting ANSP. An increase of precipitation will inevitably lead to more serious soil erosion and intensify the nitrogen and phosphorus pollution. These impacts also tend to spread spatially due to the fluidity of the water, but the infrastructure construction for farmland water conservancy will have a certain role in promoting the intensity of pollution mitigation. Therefore, the positive effects of precipitation and runoff cannot fully reflect the overall impact of this variable on changes in non-point source pollution emissions.

The direct, indirect and total effect coefficients of rural infrastructure investment are –0.0298, 0.0076 and –0.0222, respectively. The direct effect of infrastructure investment is negative and significant, indicating that it has increased the capacity of agricultural production and pollution control, thereby promoting the decline of ANSP intensity in the study region. The indirect effect is positive but not significant, indicating a positive correlation between the increase of rural infrastructure investment in the study region and the intensity of ANSP in neighboring regions. The rural infrastructure investment at the provincial level has a competitive effect, leading to an increase in its scale. The regional agricultural production capacity has induced the accumulation of economic factors from the surrounding areas to the study region, giving rise to a certain negative impact on the agricultural infrastructure investment and an indirect increase in the ANSP intensity in the surrounding areas, but this effect is not strong.

The effects of agricultural R&D investment are not significant, indicating its limited impact on the regional ANSP intensity. In fact, the effect of new technology in environmental protection was neutral on the intensity of agricultural pollution. Though higher agricultural technology and pollution prevention capabilities in economically developed regions, their high-consumption and high-emission growth methods have not been completely changed. This is one of the reasons that eastern areas mostly top the ANSP list in China.

The quality of labor has a significant and negative impact on ANSP, and both direct and indirect effects are significantly negative. The negative direct effect indicates that the improvement of labor-force quality in the study region has promoted the decline of its ANSP intensity. Higher-quality agricultural workers have a stronger sense of acquiring new knowledge, which contribute to improving agricultural production reduction methods and reducing the environmental damage behavior in agricultural production processes. The indirect effect of labor-force quality is negative and higher than the direct effect, indicating ANSP could be reduced through improving labor-force quality, and good demonstration effects could be formed among regions. The study region can learn from the successful experience in neighboring regions to improve the quality of agricultural labor and achieve the agricultural development benefits from reducing ANSP.

Based on the results from this research, the main conclusions are summarized as follows.

(1) The ANSP intensity in China is generally on the rise, and the spatial difference is obvious. The high emission intensity is concentrated in the eastern and central regions and relatively low in the western region.

(2) Spatial autocorrelation analysis reveals a significant spatial agglomeration of ANSP intensity from 1999 to 2017 as well as a gradual increase in the degree of agglomeration. The ANSP intensity is highly consistent with the economic-geographical distribution pattern, which shows a higher level in eastern China versus a lower level in western China.

(3) The convergence analysis of ANSP intensity is applicable to the spatial lag model. Compared with the convergence analysis in tradition model, the convergence rate was accelerated after considering the spatial factors, but slower than that of regional economic growth.

(4) The increase in infrastructure investment, labor-force quality and planting industry is conducive to the reduction of ANSP intensity in the region, while the expansion of agricultural economy scale and precipitation runoff has a positive impact on ANSP. The direct effect of agricultural R&D inputs is not significant. Meanwhile, improving the quality of labor force can also significantly reduce the ANSP intensity in adjacent areas, but the impact of precipitation runoff on the ANSP of adjacent areas is significant and positive.

This paper explores the intensity of agricultural non-point source pollution (ANSP) and its spatial convergence at national and provincial levels in China from 1999 to 2017, and reveals spatial factors affecting the ANSP by constructing a spatial econometric model. In an effort to reduce the ANSP intensity in China, three suggestions are proposed as follows:

(1) Agricultural pollution emissions have been showing a varying upward trend in most provinces of China since 1999, indicating that there are still some inefficiencies and environmental losses in China’s agricultural development, and there is a great space to reduce agricultural pollution. Therefore, in the further long-term development, we should further reduce the use of pesticides and chemical fertilizers, strengthen the prevention and control of agricultural pollution, accelerate the transformation of traditional agriculture to modern agriculture, and promote the sustainable development of agriculture.

(2) The convergence rate of ANSP in China is lower than that of regional economic growth, due to its convergence acceleration induced by spatial spillover factors. Therefore, the transfer and diffusion of modern agricultural technology should be accelerated to accelerate the rate of convergence of ANSP intensity though technology spillovers. In particular, high-ANSP intensity areas should further strengthen the education and training of workers and gradually promote the overall decline in the intensity of agricultural pollution.

(3) The influence of various factors on ANSP has obvious spatial correlation characteristics. When formulating agricultural emission reduction policies, the spatial interaction between various influencing factors should not be ignored. In order to reduce the ANSP intensity in the region, it is important to optimize the agricultural industrial structure, upgrade the level of human capital and increase the infrastructure investment. On this basis, the regional agricultural economic development should be rationally designed to reduce competitive effect driving by the increasing of agricultural investment. At the same time, all regions should further enhance exchanges and cooperation and strengthen the demonstration effect of ANSP reduction measures among regions.

Though the research results has further enhanced by exploring the spatial convergence characteristic and the influencing mechanism of the ANSP, we are cautious in claiming any causal effect between the influencing factors and the ANSP. Further research should be conducted with county data and survey data of farmers from the major agricultural producing areas.