Tree-Structured Parzan Estimator–Machine Learning–Ordinary Kriging: An Integration Method for Soil Ammonia Spatial Prediction in the Typical Cropland of Chinese Yellow River Delta with Sentinel-2 Remote Sensing Image and Air Quality Data

Abstract

Spatial prediction of soil ammonia (NH₃) plays an important role in monitoring climate warming and soil ecological health. However, traditional machine learning (ML) models do not consider optimal parameter selection and spatial autocorrelation. Here, we present an integration method (tree-structured Parzen estimator–machine learning–ordinary kriging (TPE–ML–OK)) to predict spatial variability of soil NH₃ from Sentinel-2 remote sensing image and air quality data. In TPE–ML–OK, we designed the TPE search algorithm, which encourages gradient boosting decision tree (GBDT), random forest (RF), and extreme gradient boosting (XGB) models to pay more attention to the optimal hyperparameters’ high-possibility range, and then the residual ordinary kriging model is used to further improve the prediction accuracy of soil NH₃ flux. We found a weak linear correlation between soil NH₃ flux and environmental variables using scatter matrix correlation analysis. The optimal hyperparameters from the TPE search algorithm existed in the densest iteration region, and the TPE–XGB–OK method exhibited the highest predicted accuracy (R² = 85.97%) for soil NH₃ flux in comparison with other models. The spatial mapping results based on TPE–ML–OK methods showed that the high fluxes of soil NH₃ were concentrated in the central and northeast areas, which may be influenced by rivers or soil water. The analysis result of the SHapley additive explanation (SHAP) algorithm found that the variables with the highest contribution to soil NH₃ were O₃, SO₂, PM₁₀, CO, and NDWI. The above results demonstrate the powerful linear–nonlinear interpretation ability between soil NH₃ and environmental variables using the integration method, which can reduce the impact on agricultural nitrogen deposition and regional air quality.

1. Introduction

As an alkaline gas, ammonia (NH₃) can effectively reduce the occurrence of acid rain [1,2]. NH₃ can be used to synthesize fertilizers, which greatly improves agricultural productivity and meets growing food demand [3]. However, the large-scale volatilization of NH₃ into the atmosphere leads to air pollution and nitrogen deposition, which leads to the eutrophication of water bodies and the destruction of biodiversity [4,5]. As an indirect source of nitrogen dioxide (N₂O) emissions, NH₃ has an important impact on global warming. It is estimated that the potential of carbon dioxide (CO₂) to cause global warming is only about 1/265 times that of N₂O in one hundred years [6,7], and global NH₃ emission is projected to increase rapidly to 132 Tg by 2100 [8], with the major emission source of NH₃ being agriculture. The nitrification of NH₃ will lead to soil acidification and eutrophication [9]. NH₃ is highly dependent on environmental variables compared to other gases due to its volatile characteristics [9,10,11]. Therefore, there is an urgent need for rapid and effective methods to predict NH₃ emissions from the soil with environmental variables to provide data support for alleviating NH₃ emissions or studying NH₃ from soil-influencing factors.

Traditional soil NH₃ measurement requires a lot of sampling work, which consumes significant manpower resources and causes cost waste. Additionally, laboratory determination requires using a concentrated sulfuric acid (H₂SO₄) and boric acid (H₃BO₃) solution. The time period for obtaining soil NH₃ content is long, and it can easily cause secondary environmental pollution. Spatial prediction is an important method that can be used to realize the continuous distribution of NH₃ in farmland soil in comparison with the high cost and difficulty of a soil NH₃ content survey. It is essential to develop a high-performance prediction model with efficient environment variables for spatial prediction of soil NH₃ flux, which is of profound significance to alleviating global warming and improving air quality.

Air quality (such as PM_2.5, NO₂, and SO₂) has been shown to closely interact with soil ammonia. Because of the important role of NH₃ in aerosol nucleation and haze [12], the gradual decrease in the formation of ${\mathrm{NH}}_4^{+}$ aerosols will lead to a significant increase in the residence time of NH₃ in the atmosphere [13], while excessive NH₃ stays in the aerosol phase, which may reduce the removal of the condensed phase [14]. On the contrary, the reaction of SO₂ and NH₃ with active radicals and NO_x may also affect the formation of O₃ and secondary aerosol intermediates [15,16]. Gu et al. (2022) found that the reduction in acidic precursors sulfur dioxide and nitrogen oxides in recent years has also reduced the chemical sink of NH₃ in the atmosphere [13]. This is because gaseous NO₂ is a reaction intermediate that reacts with NH₃; NH₃ has a good adsorption effect on NO₂. The concentration of SO₂, NO₂, CO, and other precursors is significantly higher, which may impact the production of ${\mathrm{NH}}_4^{+}$ in the environment [14].

A large number of methods have been widely used in the field of soil gas flux prediction. Geostatistical interpolation is a basic method that can estimate the target of any coordinate with zero deflection and minimum variance [17,18]. For example, Roberto et al. collected 50 soil CO₂ flux samples using the dynamic concentration method and interpolated the discrete points to obtain the spatial variability result of CO₂ flux in the survey area based on the ordinary kriging (OK) model [19]. However, since only the spatial coordinates have no information about other auxiliary variables, the semivariogram constructed by the kriging method depends on the measured point pairs, resulting in low spatial prediction accuracy of soil gas flux [20].

With the development of artificial intelligence, the machine learning (ML) method is often used to construct the relationship between environmental variables and targets, and is applied extensively in the prediction of farmland soil gas flux [20,21]. The commonly used models include support vector machine (SVM), artificial neural network (ANN), random forest (RF), and gradient-boosted regression (GBR). For example, Abbasi et al. used the RF model with multi-source input variables (i.e., air temperature, soil organic matter, soil moisture, soil total carbon, solar radiation, etc.) to predict CO₂ fluxes from the soil in the inorganic fertilizer environment in combination with field measured CO₂ fluxes [22]. Morad et al. used the RF model, multivariate adaptive spline curve (MDSC), and general linear model (GLM) to predict soil NO, CO₂, and nitrogen oxide fluxes under no-tillage and conventional tillage, respectively, in the case of temperature, soil moisture, and crop straw rate as environmental variables and found that the prediction results of RF model were the best [21]. Daniel et al., used quantile regression forests to predict the distribution of soil–atmosphere CH₄ and CO₂ fluxes in forest watersheds in different seasons by inputting DEM-derived topographic attributes [23]. The ML method is an efficient and robust substitution compared with traditional soil gas measurement methods [22,24].

However, the existing ML models suffer from two shortcomings in explaining the relationship between environmental variables and soil gases. On the one hand, the prediction performance of ML models is affected by hyperparameters. Previous studies have shown that hyperparameters can significantly affect the performance and accuracy of prediction models. Additionally, these hyperparameters also affect each other. It is difficult to capture the global optimal hyperparameters of the prediction model by using artificial trial-and-error experimentation and default hyperparameters. The automatic parameter adjustment method with high efficiency and robustness is an important auxiliary method to realize automatic ML [19]. For example, Yan et al., used a grid search (GS) algorithm to optimize and find the hyperparameters of the RF model, light gradient boosting machine (LGBM), and adaptive boosting models to adjust the parameter values so that the model performed best [25]. However, the GS algorithm is very time-consuming because it traverses all parameters in an isolated manner [26] and cannot meet the requirements of identifying global optimality [27]. Zhu et al., used an efficient Bayesian optimization method to find the optimal combination of hyperparameters of the extreme gradient boosting (XGB) model with higher prediction accuracy than GS [26]. On the other hand, the interpretation of the relationship between environmental variables and soil gas by ML is monotonous. The ML model has certain limitations on the data set and is suitable for nonlinear data sets [21], ignoring the spatial autocorrelation between adjacent observation data [28].

A new prediction model, the integration kriging model, has been proven to be more accurate in predicting soil properties in many current studies [29,30,31,32,33,34,35]. For example, Guo et al. found that the integration method has advantages in predicting the spatial variability of soil organic matter in the case of comparing the methods of RF, stepwise linear regression (SLR), and the combination of RF and residual kriging by inputting topographic attributes, such as geological units, climate factors, and vegetation indices [28]. Based on 73 years of NO₂ observation data from 14 monitoring stations, Chen et al. predicted the spatial–temporal characteristics of NO₂ concentrations in Taiwan using the NO₂ data interpolated kriging method and the land use regression model with local non-traditional geographical predictors [36].

However, whether the integration of the TPE-based ML model and geostatistical model can significantly improve the spatial prediction accuracy and explain the relationship between environmental variables and soil NH₃ flux is still unclear (i.e., nonlinear and linear coupling characteristics). Therefore, to overcome the limitations of the hyperparameter adjustment and spatial autocorrelation neglect of ML modeling, the objective of this study is to demonstrate the increased predicted model performance using TPE–ML–OK-based models with Sentinel-2 remote sensing image and air quality data, which involves establishing ML models, determining the optimal hyperparameters using the TPE algorithm, building a residual ordinary kriging model, and revealing the relationships between key variables and soil NH₃ flux via explainable analysis. The environment variables (i.e., Sentinel-2 remote sensing image and air quality data) related to soil NH₃ were collected first. Then, the tree-structured Parzen estimator (TPE) search was applied to find the optimal hyperparameters of the gradient boosting decision tree (GBDT), random forest (RF), and extreme gradient boosting (XGB) models, and the residual ordinary kriging model was used to further predict and mapping the soil NH₃ flux. Finally, the relationships between environmental variables and soil NH₃ were evaluated using the SHapley additive explanation (SHAP) method.

2. Materials and Methods

Figure 1 shows the flowchart in this study and includes three phases: (i) collection and preprocessing of sampling and environmental variable data (vegetation index and texture features from Sentinel-2 L2A image, interpolation analysis of air quality stations data, and soil NH₃ sampling); (ii) construction of integration method (machine learning model, TPE hyperparameter optimization, and residual ordinary kriging); and (iii) accuracy evaluation of prediction results, soil mapping using TPE–ML–OK models, and importance analysis of environmental variable using the SHAP method.

2.1. Overview of Study Area

The study area is situated in the Yellow River Delta (YRD), East China (36°55′–38°10′N, 118°07′–119°10′E). It has a temperate continental monsoon climate, and the terrain inclines from southwest to northeast along the Yellow River (Figure 2). Soil types are divided into five categories: cinnamon soil, Shajiang black soil, fluvo-aquic soil, saline soil, and paddy soil. Based on the statistical yearbook of Dongying in 2020, the average temperature in the study area is 14.1 °C, the annual precipitation is 665.3 mm, and the annual sunshine hours are 2998.5 h. The study area mainly contains oil, natural gas, oil shale, and many other mineral resources, which is the main production area of the Shengli Oilfield. However, the soil pollution caused by oilfield development in this area indirectly affects the survival rate of soil microorganisms and, thus, the absorption, transformation, and utilization of NH₃ by some soil microorganisms, which changes the NH₃ content in the soil. Spatial prediction of soil NH₃ is of great significance in improving microclimate change, air pollution, and nitrogen loss from farmland in the study area.

2.2. Soil Sample Collection and Analysis

In the farmlands of the study area, 134 sampling points were uniformly and randomly arranged from 22 September 2022 to 29 September 2022. The NH₃ volatilization was measured per day at 8:00–10:00 a.m. and 3:00–5:00 p.m. The concentration of NH₃ was determined using the continuous airflow enclosure method, and a cylindrical plexiglass chamber (20 cm in diameter and 30 cm in height) was installed in the soil of each sample plot. We inserted the air chamber into the soil of 5 cm and connected the air in the glass chamber from the air outlet to the pump-suction gas-phase speedometer (FZ2800-2, Guangzhou, China) through the pipeline. The measurement interval of each sample point was 10 min, the measurement was carried out in three parallels, and the initial and final NH₃ content (unit ppm) was recorded on the spot. The initial NH₃ content was converted into emission flux (unit: kg N ha⁻¹ d⁻¹) according to Equation (1) [37,38].

$$AE=c\times v\times 14\times 1440\times {10}^6\times {10}^{-9}/(mv\times 0.0177)$$

(1)

AE represents the volatilization flux of NH₃ in farmland soil, c represents the content of soil NH₃ (m mol mol⁻¹) from the pump-suction gas-phase speedometer, v represents the air velocity (0.25 L min⁻¹), 14 represents the molar mass of NH₃ (g mol⁻¹), 1440 represents the conversion from min⁻¹ to d⁻¹ (1 d 1440 min), 10⁶ represents the conversion from m⁻² to ha⁻¹, 10⁻⁹ represents the conversion from mg N into kg N, mv represents the molar volume of the NH₃, and 0.0177 is the covered area of the air chamber (m²).

2.3. Environmental Variables

2.3.1. Remote Sensing Data

Based on the Sentinel-2 L2A remote sensing image (23 September 2022), the remote sensing index was calculated and obtained. The data comes from the European Space Agency (https://scihub.copernicus.eu/, accessed on 12 November 2022). The L2A image data were preprocessed (i.e., atmospheric correction), and four vegetation index factors were calculated using band 2, band 3, band 4, band 5, band 6, and band 8, including visible atmospheric resistant index (VARI), soil-adjusted vegetation index (SAVI_red), normalized difference water index (NDWI), and plant senescence reflectance index (PSRI) (

Table 1. Indexes derived from Sentinel-2 L2A image.

Index	Formula	Reference
VARI	$(band3-band4)/(band3-band4+band2)$	[39]
SAVIred	$(band8-band5)\times 1.5/(band8+band5+0.5)$	[40]
NDWI	$(band3-band8)/(band3+band8)$	[41]
PSRI	$(band4-band3)/band6$	[42]

); the resolution of all vegetation variables was 10 m × 10 m.

Texture features can reflect additional information, such as hue changes and attributes of remote sensing images, and have been widely used in the extraction or classification of ground object information [43,44,45]. For example, many studies have used the GLCM method to calculate the relative frequency between the pixel brightness values of remote sensing images to obtain spatial information gain, which effectively improves the interpretation and classification accuracy of ground objects [46,47,48]. Based on the 11 bands of the preprocessed Sentinel-2 L2A image, the texture features were also calculated and analyzed using filter/co-occurrence measures in ENVI 5.3 software. A total of 88 texture features were obtained, including eight transformations: entropy (ENT), correlation (CORR), mean (MEAN), homogeneity (HOM), contrast (CTRA), dissimilarity (DIS), angular second moment (SECM), and variance (VAR). In order to simplify the variable dimension, the first three principal components were obtained using principal component analysis (PCA) in ENVI 5.3 software, and the total variance contribution rate that can be explained is 65.83%. The features with high load values of PC1, PC2, and PC3 were mainly concentrated in band 6–band 9; band 1–band 4; and band 5, band 10, and band 11 after texture transformation (Figure 3). These principal components (PC1–PC3) representing image texture features can reflect soil physical information such as spatial information (i.e., pattern, shape, size, etc.), geometric morphology, texture, and coverage on the soil surface, and have important implications for the basic soil conditions of NH₃ emissions.

2.3.2. Air Quality

The air quality in the study area indirectly interfered with soil microbial interaction, affecting soil NH₃ emissions. In order to be consistent with the sampling time of soil NH₃, the selected time of air quality data was the monthly average of September 2022. Six air quality pollutants (CO, NO₂, O₃, PM_2.5, PM₁₀, and SO₂) were obtained as environmental variables for spatial prediction of NH₃. The data were obtained from the real-time air quality release system (http://218.58.213.53:8081/dyfb_air/fb_web/dongying.html, accessed on 10 November 2022) of Dongying Environmental Protection Bureau, including 40 air quality monitoring stations. The air quality data of all stations were processed using interpolation analysis with a resolution of 10 m × 10 m in ArcGIS 10.7 software.

2.4. Machine Learning Model

2.4.1. Random Forest

Random forest (RF) is a tutorial integration algorithm with a decision tree as a base estimator [49]. In the modeling process of RF, different subsets are randomly sampled from the provided dataset and then used to build multiple decision trees. According to bagging rules, the results of several decision trees were synthesized to predict soil NH₃ fluxes, that is, the average values of the results output by several base estimators. The learning ability of the RF model is strong, and it can overcome the over-fitting problem caused by a highly dimensional variable input. The RF model has been widely used in soil, ecology, and other fields [50,51,52]. The core hyperparameters of the RF model that need to be optimized are shown in

Table 2. Information description of core hyperparameters of RF model.

Hyperparameters	Type	Range	Explanation
n_estimators	int	(50, 80)	Number of trees
max_depth	int	(10, 17)	The maximum depth of a tree
min_impurity_decrease	float	(0, 0.1)	Minimum impurity in the node split
min_samples_split	int	(2, 8)	The minimum samples required for node re-split
max_features	int, float, string	(0, 64, “log2”, sqrt”, “auto”)	The number of features needed to find the best segmentation

. For example, max_depth and max_features have a significant impact on the prediction accuracy and splitting depth of the basic estimator.

2.4.2. GBDT

Gradient boosting decision tree (GBDT) is one of the ensemble ML models that uses classification and regression tree (CART) as the basic model [53]. Different from the bagging rule of the RF model, the GBDT model meets the basic process of the boosting rule, calculates the loss function based on the results of the previous base estimator, and uses the loss function to affect the next base estimator (Equation (2)) [53].

$$f_n(x)=f_{n-1}(x)+T_n(x)$$

(2)

$f_n(x)$ and $f_{n-1}(x)$ represent the n-th (n ≥ 1) and n − 1 weak learners, respectively. $T_n(x)$ represents a new learner based on the residual of the n − 1 weak learner. The current i-th learner $f_i(x)$ in the GBDT model is affected by all previously trained weak learners [49,54,55] and finally integrated multiple base estimators and output soil NH₃ prediction values. Because GBDT is optimized in the base estimator, loss function, fitting residual, and random sampling, it has a stronger learning ability than most boosting algorithms. GBDT is one of the most stable machine learning algorithms in actual scenarios. The core hyperparameters and information of the GBDT model are shown in

Table 3. Information description of core hyperparameters of GBDT model.

Hyperparameters	Type	Range	Explanation
n_estimators	int	(90, 120)	Number of trees
learning_rate	float	(0.2, 0.3)	The learning speed
subsample	float	(0.6, 0.72)	The proportion of subsampling
max_depth	int	(5, 10)	The maximum depth of a tree
max_features	int, float, string	(2, 16, “log2”, “sqrt”, “auto”)	The number of features needed to find the best segmentation
min_impurity_decrease	float	(2, 4)	The amount of information gained to consider when splitting nodes

.

2.4.3. XGB

The eXtreme gradient boosting (XGB) is a new machine-learning model proposed by Chen et al. in 2016. It has been optimized based on the GBDT algorithm, making XGB have better results in dealing with regression and classification problems [51,52]. As an advanced boosting algorithm, XGB builds the next weak learner based on the results of the previous weak learner in each iteration [53]. The advantage of the XGB algorithm is that the objective function with a regular term controls overfitting [54]. In addition, similar to the RF model, XGB also supports column sampling, which extracts some data for training during the iteration process so as to control overfitting by reducing the amount of calculation. Compared with the first-order Taylor expansion of GBDT, XGB adopts the second-order Taylor expansion, which can approximate the real loss function more accurately. Equation (3) shows the objective function ($O^{\left(t\right)}$) [51].

$$O^{(t)}\cong {\displaystyle \sum {}_{i=1}^n}[L(y_i,{\widehat{y}}_i^{(k-1)})+g_i{f}_k(x_i)+\frac{1}{2}{h}_i{f}_k^2(x_i)]+\mathrm{\Omega }(f_k)$$

(3)

where $y_i$ and ${\widehat{y}}_i^{}$ represent the i-th measured value and the predicted value at step t, respectively, L represents the loss function, Ω represents the regular term, $x_i$ and $f_k$ represent the i-th input variables and the k-th tree, respectively, and $g_i$ and $h_i$ represent the first and second derivatives of the L, respectively. Then, the final objective function (Equation (4)) at the j-th leaf $I_j=\left\{i|q(x_i)\right\}=j$ is as follows:

$$O^{(t)}={\displaystyle \sum {}_{j=1}^T}[({\displaystyle \sum {}_{i\in I_j}{g}_i){\omega }_j+\frac{1}{2}({\displaystyle \sum {}_{{}_{i\in I_j}}{h}_i+\lambda })}{\omega }_j^2]+\gamma T,$$

(4)

where T represents the total leaf nodes of the XGB model, ${\omega }_j$ and $\gamma$ represent the scores on the j-th leaf and complexity of leaf nodes, respectively, and $\lambda$ is a compromise parameter that scales the penalty. The core hyperparameters and information of the XGB model are shown in

Table 4. Information description of core hyperparameters of XGB model.

Hyperparameters	Type	Range	Explanation
subsample	float	(0.85, 1)	Construct the sampling rate of each tree to the sample
num_round	int	(30, 50)	The number of trees
eta	float	(0.2, 0.3)	Learning rate
lambda	float	(0, 4)	Regularization section for processing XGB
min_child_weight	float	(0.1, 2.8)	The sum of weights of the minimum leaf node sample
colsample_bytree	float	(0.92, 1)	The proportion of features used in training out of all features
colsample_bynode	float	(0.95, 1)	Sub-sampling rate of columns split per node
max_depth	int	(1, 4)	Maximum depth of a tree

.

2.4.4. TPE Search Algorithm

Not only is the internal structure of the above prediction model complex, but many hyperparameters also need to be adjusted. Moreover, the hyperparameters have a great influence on the prediction performance of the model [55,56]. Hyperparameters usually refer to the parameters shared between different modules, functions, classes, or objects in the model. These parameters can be used by multiple components to achieve data transfer, interaction, and sharing of information. The traditional manual adjustment will take a lot of time and cost, and the effect is not good. Studies have shown that compared with the GS algorithm, the Bayesian optimization (BO) algorithm cannot easily plunge the local optimum and is more efficient. Bayesian optimization is considered the most advanced superparameter optimization framework at present and is widely used in various fields of AutoML [57]. It uses the frequency of the minimum value to judge. The frequency reflects the probability of the minimum value to some extent. The higher the frequency, the greater the probability that the function has a minimum value. Tree-structured Parzen estimator (TPE) is one of the most advanced BO algorithms based on tree-structured Parzen density estimation. First, the TPE algorithm defines two probability density functions (Equations (5) and (6)) based on $y^{*}$ (i.e., the value at the quantile r in the y-set).

$$r=p(y<y^{*})$$

(5)

$$p(x | y)=\{\begin{array}{l}l(x)\hspace{1em}if\hspace{1em}y<y^{*}\\ g(x)\hspace{1em}if\hspace{1em}y\ge y^{*}\end{array}$$

(6)

Then, the TPE optimization algorithm selects EI (expected improvement) as the evaluation criterion (Equation (7)):

$$E{I}_{y^{*}}(x)={\displaystyle {\int }_{-\infty }^{y^{*}}(y^{*}-y)p(y | x)dy={\displaystyle {\int }_{-\infty }^{y^{*}}(y^{*}-y)\frac{p(x | y)p(y)}{p(x)}}}dy$$

(7)

Finally, we can obtain Equation (8):

$$E{I}_{y^{*}}(x)=\frac{r{y}^{*}l(x)-l(x){\displaystyle {\int }_{-\infty }^{y^{*}}p(y)dy}}{rl(x)-(1-r)g(x)}\propto {\left(r+\frac{g(x)}{l(x)}(1-r)\right)}^{-1}$$

(8)

In Equation (8), we can see a direct correlation between the value of EI and the value of $g(x)/l(x)$. In order to improve EI, we make the probability of $l(x)$ as large as possible and the probability of $g(x)$ as small as possible so as to select the most suitable x value. In each iteration, the algorithm returns $x^{*}$ with the maximum EI value, and the returned $x^{*}$ will participate in the next iteration so as to repeatedly find the best hyperparameter combination of the soil NH₃ prediction model. In this study, packages (i.e., xgboost, sklearn.model_selection, sklearn.ensemble, hyperopt, etc.) are imported into the Python 3.9 programming environment to build the above algorithms (i.e., RF, GBDT, XGB, and TPE).

2.5. The Integration Method

The integration method was used to further improve the prediction accuracy of soil NH₃ flux. First, the residual value $r_{TPE-ML}(x_i)$ between the true value $z\left(x_i\right)$ of soil NH₃ and the predicted value $z_{TPE-ML}(x_i)$ of the TPE-optimized ML method (i.e., TPE-GBDT, TPE-RF, and TPE-XGB) was calculated (Equation (9)).

$$r_{TPE-ML}(x_i)=z(x_i)-z_{TPE-ML}(x_i)$$

(9)

Then, the ordinary kriging (OK) interpolation algorithm was used in ArcGIS 10.7 software to obtain the residual spatial continuity distribution of the TPE-optimized machine learning model. The formula of the OK model (Equations (10) and (11)) is as follows [58]:

$$\widehat{\gamma }(h)=\frac{1}{2N(h)}{\displaystyle \sum _{i=1}^{N(h)}{[z(x_i)-z(x_i+h)]}^2}$$

(10)

$${\widehat{z}}_{OK}(x_i)={\displaystyle \sum _{i=1}^n{\lambda }_{i}z(x_i)}$$

(11)

where $\widehat{\gamma }(h)$ represents the experimental semivariogram in the kriging interpolation process, $N(h)$ represents the number of sampling points separated by h (where h is the lag distance between $z(x_i)$ and $z(x_i+h)$, and ${\lambda }_i$ is the optimal weight. ArcGIS 10.7 software was used to realize the spatial interpolation calculation of the OK model. Finally, the spatial prediction results of soil NH₃ flux in the hybrid geostatistical model were obtained by adding the spatial prediction value of the TPE-optimized machine learning model to the residual kriging value (Equation (12)) [59].

$$z_{TPE-ML-OK}(x_i)=z_{TPE-ML}(x_i)+{\displaystyle \sum _{i=1}^n{\lambda }_i{z}_{TPE-ML}(x_i)}$$

(12)

where $z_{TPE-ML-OK}(x_i)$ includes the comprehensive prediction results of tree-structured Parzen estimator–random forest–ordinary kriging (TPE–RF–OK), tree-structured Parzen estimator–gradient boosting decision tree–ordinary kriging (TPE–GBDT–OK), and tree-structured Parzen estimator–extreme gradient boosting–ordinary kriging (TPE–XGB–OK), representing the optimal spatial estimation of soil NH₃ flux.

2.6. Verification

In this study, the datasets were randomly divided into 70% for the training set and 30% for the test set to build a prediction method. For accuracy evaluation, root mean square error (RMSE) (Equation (13)) and coefficient of determination (R²) (Equation (14)) were used. $x_i$ and ${\widehat{x}}_i$ are, respectively, the real value of soil NH₃ flux and the prediction value of the super parameter optimization model, and ${\overline{x}}_i$ and ${\overline{\widehat{x}}}_i$ are, respectively, the average of the real value of soil NH₃ at all sample points and the average of the estimated value of the model.

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{({\widehat{x}}_i-x_i)}^2}}$$

(13)

$${\mathrm{R}}^2=\frac{{[{\displaystyle {\sum }_{i=1}^m(x_i-{\overline{x}}_i)({\widehat{x}}_i-{\overline{\widehat{x}}}_i)}]}^2}{{\displaystyle {\sum }_{i=1}^m{(x_i-{\overline{x}}_i)}^2{\displaystyle {\sum }_{i=1}^m{({\widehat{x}}_i-{\overline{\widehat{x}}}_i)}^2}}}$$

(14)

Because the ensemble learning models, deep learning models, and other models are extremely complex, abstract, and difficult to understand, the importance score of each environmental variable was output and visualized using the SHapley additive explanation (SHAP) analysis so that we can determine which environmental variables play a leading role in affecting soil NH₃ flux. As a unified interpretation framework, SHAP can calculate the importance score of each feature in the data so as to explain the model. In this study, the SHAP package was imported into the Python 3.9 programming environment to explain the predicted method and visualize the importance of different environment variables.

3. Results and Analysis

3.1. Statistical Analysis

Figure 4 shows the normalized numerical distribution of NH₃ flux in soil and auxiliary variables. The soil NH₃ fluxes are concentrated in the median area, and the data structure is normally distributed. PSRI, PC1, and PC2 also showed significant normal distribution. However, the values of other variables showed skewed or bimodal distribution. There were significantly nonlinear characteristics among different environmental variables. The linear correlation between similar variables was strong, the positive linear correlation between vegetation index factors (i.e., SAVI_red and NDWI, VARI and PSRI, and SAVI_red and PSRI) was high (r > 0.2), the negative linear correlation between soil texture factors (PC3 and PC1) was strong (r = 0.39), and there was a high correlation between air quality variables. Obviously, there was a positive linear correlation between soil NH₃ flux and O₃ (r = 0.21) and a weak linear correlation with other variables. The above correlation analysis results show complex fitting characteristics between soil NH₃ flux and environmental variables.

3.2. Spatial Prediction of Soil NH₃ Using TPE–ML–OK Method

3.2.1. Evaluation of Hyperparameter Optimization Process

Figure 5, Figure 6 and Figure 7 show the sampling point iteration trend of GBDT, RF, and XGB models in the TPE optimization process, respectively. Obviously, the sampling points in the TPE optimization process have an aggregation trend, and there is a region with the densest distribution of sampling points in each hyperparameter iteration process. For example, some hyperparameters (min_impuity_decrease of the GBDT model and min_child_weight of the XGB model) increase with the increase in the value, and the distribution of sampling points change from sparse to dense. Other hyperparameters (n_estimators of the RF model and colsample_bynode of the XGB model) show the trend of sampling points from dense to sparse (that is, it is difficult to reach the optimal value, even if the hyperparameter value is increased). These sampling-intensive areas are significantly different from the sparse sampling areas, and the optimal parameters often appear near the sampling-intensive areas. The lowest RMSE values obtained by the RF, GBDT, and XGB models in the TPE optimization process are relatively close, which are 11.96 kg N ha⁻¹ d⁻¹, 11.77 kg N ha⁻¹ d⁻¹, and 11.70 kg N ha⁻¹ d⁻¹, respectively, indicating that the TPE algorithm has strong generalization in the hyperparameter optimization of different models. TPE always selects the next set of sampling points according to the results of previous sampling points [60]. After multiple iterations, the probability that the optimal hyperparameters appear in the dense sampling area is much higher than in the sparse area.

3.2.2. Prediction Accuracy of TPE–ML Model for Soil NH₃

The scatter and residual distribution of the measured and predicted values of soil NH₃ concentration and the regression accuracy of the model are shown in Figure 8 and

Table 5. The accuracy evaluation of prediction models for soil NH₃ flux.

Model	RMSE (kg N ha−1 d−1)	R2 (%)	Model	RMSE (kg N ha−1 d−1)	R2 (%)
TPE–GBDT	10.40	65.90%	TPE–GBDT–OK	8.28	75.48%
TPE–RF	8.85	71.78%	TPE–RF–OK	7.11	80.92%
TPE–XGB	8.13	74.22%	TPE–XGB–OK	6.42	85.97%

, respectively. The spatial prediction accuracies of the TPE-optimized machine learning model for soil NH₃ fluxes are TPE–XGB (RMSE = 8.13 kg N ha⁻¹ d⁻¹ and R² = 74.22%) > TPE-RF (RMSE = 8.85 kg N ha⁻¹ d⁻¹ and R² = 71.78%) > TPE–GBDT (RMSE = 10.40 kg N ha⁻¹ d⁻¹ and R² = 65.90%) (Table 5). The TPE–XGB model had the strongest ability to fit the nonlinear relationship between soil NH₃ flux and environmental variables in comparison with the other models. The residuals of the training sets of the three models satisfy the normal distribution, indicating a better-fitting performance of the models. The scatter distribution of the test set of the TPE–GBDT model is dispersive, and the prediction accuracy is significantly affected by the anomalous samples. In addition, the fitting curves of the residuals of the test sets of the three models show left skewness (i.e., the predicted residuals are concentrated in the negative region), indicating that the model predictions are underestimated.

3.2.3. The Semi-Variation Analysis of Residual Values from TPE–GBDT, TPE–RF, and TPE–XGB Models

Furthermore, the spatial semi-variation characteristics of the predicted residuals of the three models were analyzed in ArcGIS 10.7 software. The nugget effect is an important indicator used for measuring the spatial variability of the sample points, usually divided into low spatial dependence (NE ≤ 25%) and moderate spatial dependence (25% ≤ NE ≤ 75%) [61]. The NE of the soil NH₃ flux of residual sample points for OK exhibited moderate spatial dependence, making it appropriate for kriging interpolation (

Table 6. The semi-variation analysis of residuals of soil NH₃ using TPE–GBDT, TPE–RF, and TPE–XGB models.

Variogram	Function	Nugget (Co)	Sill (Co + C)	Nugget/Sill [Co/(Co + C)] (%)
Residuals of TPE–GBDT	Exponential	120.57	163.07	73.94
Residuals of TPE–RF	K-Bessed	122.10	178.82	68.28
Residuals of TPE–XGB	Exponential	119.80	161.17	74.33

). In the process of semivariogram analysis, the residual sampling points showed a nonuniform distribution, and the semivariogram value (γ) tended to be stable as the step size increased. The residual sample points have the highest density (i.e., the strongest spatial autocorrelation) between the distance (h) from 0.1 to 0.6 (Figure 9).

3.2.4. Spatial Prediction and Mapping Using the TPE–ML–OK Method

As shown in Table 5, the combined prediction accuracy of the TPE–ML–OK method was higher than that of the monotonic machine learning model, and the three integration models characterized more than 75% of the spatial variability (R²) of soil NH₃ fluxes. Moreover, the test set scattering points of the three integration models were more compactly distributed near the diagonal, and their residuals were closer to normal distribution and concentrated in the low-value region (Figure 8). The TPE–XGB–OK model exhibited the highest prediction accuracy (i.e., RMSE = 6.42 kg N ha⁻¹ d⁻¹, and R² = 85.97%) in comparison with the other models, which indicated that the TPE–XGB–OK model could explain the linear and nonlinear relationships between soil NH₃ and environmental variables. The above results exhibited that the ML models assisted by the TPE algorithm have a more uniform distribution of residuals with fewer extremes, and the smoothing effect of residual kriging interpolation promotes the higher prediction accuracy of the TPE–ML–OK method.

Figure 10 shows soil NH₃ fluxes spatially mapped using three integration prediction models. The predicted range of the TPE–XGB–OK model (0.7–82.36 kg N ha⁻¹ d⁻¹) was closest to the range of the true values (1.31–81.89 kg N ha⁻¹ d⁻¹), whereas the mapping results of the TPE-GBDT-OK model appeared to feature a significant underestimation. The spatial mapping results of the three hybrid models were similar, with the TPE–XGB–OK model predicting a more natural spatial excess of soil NH₃ fluxes. Soil NH₃ fluxes in the study area were smaller in the south as a whole, and the areas with larger concentrations were mainly distributed in the central and northeastern of the study area. Especially near the Yellow River estuary, the high soil NH₃ fluxes may be influenced by river transport effects or soil moisture.

4. Discussion

4.1. Importance Analysis of Environmental Variables on Soil NH₃

SHAP analysis is used to explore the contribution trajectory of environmental variables in the prediction process based on the TPE–XGB model. As shown in Figure 11, among the three types of variables of vegetation, air quality, and soil texture, the most significant impact on the spatial variability of soil NH₃ flux is air quality (i.e., cumulative contribution ratio is 60.52%), followed by vegetation (i.e., cumulative contribution ratio is 30.79%) (Figure 11a). Among all the environmental variables, O₃ had the highest contribution to the prediction accuracy of soil NH₃ (26.29%), and the cumulative importance contribution rate of the first five environmental variables exceeded 70%. The mean SHAP values of other environmental variables such as VARI, PC3, SAVI_red, PC1, and PM_2.5 were all below 1, which had a low importance on the prediction accuracy of soil NH₃. The characterization of the driving force of environmental variables on soil NH₃ flux at each sample point is shown in Figure 11b. The SHAP value of the samples with low O₃ value is negative (i.e., the closer the color is to blue, the smaller the feature value is), and the samples with high O₃ value have positive SHAP values (i.e., the closer the color is to red, the larger the feature value is). That is, O₃ shows a significant positive driving effect on soil NH₃ (i.e., an increase in O₃ content in the air will promote the emission of soil NH₃). Similarly, NDWI is also positively correlated with soil NH₃. Inversely, the SHAP values corresponding to the samples with higher SO₂, PM₁₀, and PSRI values are negative. This result indicates that these variables had a negative driving effect on soil NH₃.

4.2. The Mechanism of Influence of Environmental Variables on Soil NH₃

Air quality is an important driving factor affecting NH₃ fluxes from farmland soil in the study area. According to the SHAP analysis, the positive driving effect of O₃ is the most significant. The mechanism of the response of O₃ on soil N emissions is that it changes soil microbial activity through plant-mediated processes, which may have an important impact on net gas fluxes, because it does not penetrate the soil and has a direct impact on the components of the soil ecosystem [62]. For example, elevated O₃ concentrations reduce net photosynthesis [63,64] and dry matter accumulation [65] by oxidatively damaging plant cell membranes and chloroplasts [63,64,66,67,68,69], which in turn alters the soil subsurface processes of plant root respiration and microbial activity [41,70,71,72,73], ultimately affecting the nitrogen (N) and carbon (C) cycles in soil [74,75] and greenhouse gas emissions (CO₂, CH₄, N₂O, etc.) [76]. Recio et al. (2020) found that the high O₃ concentration significantly reduced the absorption of NH₃ and other emissions from soil by plants [9], thereby increasing the emission of some nitrogen gases. For example, the increase in O₃ concentration will significantly increase the peak CH₄ flux of wheat and rice [77], and even under the highest O₃ concentration treatment, the cumulative emission of soil N₂O will double [78]. With the increase in O₃ concentration, the emission of N element gas also increased, which proved that O₃ was a positive driving force for soil NH₃ in this study.

In addition to air quality variables, soil moisture (NDWI) also had a significant effect on the spatial variability of soil NH₃ flux in this study. On the one hand, soil moisture has a significant positive effect and sensitivity on soil NH₃ emission [79,80]. For example, since soil NH₃ flux is a physical process affected by the concentration of ${\mathrm{NH}}_4^{+}$-N in soil solution [81], the increase in soil moisture will promote the diffusion of NH₃ in soil, resulting in the increase of soil NH₃ emission [79,82,83,84]. On the other hand, soil moisture controls the loss of NH₃ by affecting urea-hydrolyzing microorganisms and urease activity [84,85,86]. The increase in the liquid diffusion rate in moist soil is beneficial to the absorption of matrix by microorganisms and promotes the diffusion of NH₃ in soil (and the emission of soil NH₃) [87]. This is similar to the results of this study. In the estuary area of the Yellow River in the northeast of the study area, due to the dense river network and the high moisture of the farmland soil, the soil NH₃ concentration is high. In addition, other variables such as soil salinity, soil pH, temperature, and wind speed may also affect soil NH₃ emissions [88,89,90,91].

4.3. Limitations

We developed a TPE–ML–OK method with Sentinel-2 remote sensing image and air quality data to predict the spatial distribution of NH₃ in farmland soil. However, this study has some limitations that should be considered in the future. First, the environmental variables that dominate soil NH₃ emissions under different land use types are significantly different. For example, differences in land cover and land management methods can significantly affect soil NH₃ emissions. The highest total environmental NH_x concentration in rural areas may be related to intensive agricultural activities (such as fertilization and irrigation) [92], while soil pH, mineral nitrogen content, soil temperature, and humidity in the agro-pastoral ecotone play a key role in regulating NH₃ emissions [93].

Second, the high fertilization rate was a significant environmental variable that promoted soil NH₃ emission. The accumulation of chemical fertilizers and organic fertilizers in the soil leads to an increase in the transfer of NH₃ from the liquid phase to the gas phase and promotes the volatilization of NH₃ [94,95]. For example, the high nitrogen fertilizer input will lead to an increase in NH₃ volatilization, and the more nitrogen fertilizer applied, the more nitrogen loss caused by NH₃ volatilization [96], which becomes the main source of ammonia emission increasing in the atmosphere [97]. On the contrary, the volatilization of NH₃ is limited by washing the fertilizer into the deep soil so that it is adsorbed [98].

Third, the influence of meteorological conditions on soil NH₃ emission should be considered. NH₃ volatilization loss is affected by many factors, such as temperature, precipitation, and precipitation [88,89,90,91,99]. Cold weather is not conducive to the emission of NH₃, and strong precipitation may also lead to the removal of ammonia in the soil, resulting in a low concentration of ammonia discharged from the soil [97]. With the increase in temperature, the initial emission of NH₃ will increase. When the temperature exceeds the optimum level, NH₃ volatilization may not be significantly affected by temperature changes. [95] Other soil properties also have positive or negative effects on NH₃ emissions. For example, acidic loam is also not conducive to the emission of NH₃ [97]; that is, soil pH is positively correlated with NH₃ emission, and soil pH will affect the content of ${\mathrm{NH}}_4^{+}$-N in the ammonia conversion reaction system, thus affecting the emission process of soil NH₃ [93]. The higher the soil pH, the more NH₃ volatilization [96]. The increase in clay content leads to decreased electrical conductivity, which limits the vertical transfer of nitrogen [95].

In addition, spatial interpolation is used to obtain the spatial continuity distribution of air quality stations, which is a stable and effective method in small-scale areas. However, in future studies, improving the spatial accuracy of air quality data should be considered in order to further improve its response to the spatial prediction of soil NH₃ and reduce the uncertainty of prediction. For example, in large-scale areas, the use of aerosol images (such as MODIS and Himawari) and a large number of monitoring station data, combined with machine learning models, can effectively improve the spatial prediction accuracy of air quality, thereby reducing the uncertainty of air quality variables in predicting soil NH₃ flux.

5. Conclusions

In this study, we developed an integration method to predict the spatial distribution of NH₃ flux in farmland soil. The GBDT, RF, and XGB models assisted by the TPE algorithm found the optimal hyperparameters, and the probability that the optimal hyperparameters appear in the dense sampling area of TPE iterations. Compared with a single ML model, the three integration models exhibited higher predicted performance for the spatial variability of soil NH₃ fluxes. Moreover, the TPE–XGB–OK model effectively explained the linear and nonlinear relationships between environmental variables and soil NH₃. The areas with high flux of NH₃ in farmland soil were concentrated centrally and in the northeast. Further importance analysis of environmental variables exhibited that air quality and soil moisture have a significant positive driving effect on soil NH₃. That is, O₃ promotes soil NH₃ flux by altering the soil subsurface processes of plant root respiration and microbial activity, and soil moisture also has a direct promotion effect on soil NH₃ emission. In the future, the integration method can be generalized to spatio–temporal fusion prediction, supporting global warming alleviation and high-quality agriculture development in coastal areas.