Comparison of Machine Learning Methods for Predicting Soil Total Nitrogen Content Using Landsat-8, Sentinel-1, and Sentinel-2 Images

: Soil total nitrogen (STN) is a crucial component of the ecosystem’s nitrogen pool, and accurate prediction of STN content is essential for understanding global nitrogen cycling processes. This study utilized the measured STN content of 126 sample points and 40 extracted remote sensing variables to predict the STN content and map its spatial distribution in the northeastern coastal region of Hebei Province, China, employing the random forest (RF), gradient boosting machine (GBM), and extreme gradient boosting (XGBoost) methods. The purpose was to compare the ability of remote sensing images (Landsat-8, Sentinel-1, and Sentinel-2) with different machine learning methods for predicting STN content. The research results show the following: (1) The three machine learning methods accurately predicted the STN content and the optimal model provided by the XGBoost method, with an R 2 of 0.627, RMSE of 0.127 g · kg − 1 , and MAE of 0.092 g · kg − 1 . (2) The combination of optical and synthetic aperture radar (SAR) images improved prediction accuracy, with the R 2 improving by 45.5%. (3) The importance of optical images is higher than that of SAR images in the RF, GBM, and XGBoost methods, with optical images accounting for 87%, 76%, and 77% importance, respectively. (4) The spatial distribution of STN content predicted by the three methods is similar. Higher STN contents are distributed in the northern part of the study area, while lower STN contents are distributed in coastal areas. The results of this study can be very useful for inventories of soil nitrogen and provide data support and method references for revealing nitrogen cycling.


Introduction
Nitrogen is one of the key essential nutrients for plant growth and development [1]. Low levels of nitrogen can negatively impact plant growth, while excessive nitrogen levels can lead to reduced ecosystem productivity and environmental pollution [2,3]. Soil is a crucial nitrogen pool in terrestrial ecosystems, playing a fundamental role in the global nitrogen cycle [4]. However, the rapid development of the economy has caused changes in land use types, particularly the conversion of natural ecosystems to artificial ones [5,6], boosting (XGBoost) [26], and generalized boosted machines (GBM) [27], have emerged to help explain nonlinear relationships. However, choosing the best modeling method for a given region has always been a challenge for soil property mapping [28].
The primary objective of this study was to map the STN content in the coastal wetlands of northeast Hebei, China, using Landsat-8, Sentinel-1, and Sentinel-2 images, and to evaluate the effectiveness of different remote sensing sensors. Landsat-8, Sentinel-1, and Sentinel-2 images were obtained for generating predictors (multispectral remote sensing bands, remote sensing indices, and backscatter coefficients). We utilized RF, GBM, and XGBoost methods to compare the prediction accuracy of different combinations of these predictor variables in predicting STN content. Furthermore, we assessed the potential of different remote sensing sensors, such as Landsat-8, Sentinel-1, Sentinel-2, and different combinations of sensors, for mapping STN content. We then investigated the importance of the generated predictor variables. Finally, we plotted the spatial distribution of STN content in the study area based on the optimal models. This study helps to explore the most suitable remote sensing imagery and machine learning methods for predicting STN content in coastal areas. The map of STN content provides a better understanding of land resources, helps assess land suitability for different uses, and aids in land planning and decision-making.

Study Area
The study area is located in the northeast of Hebei, China (38.92 • -40.32 • N, 118.14 • -119.85 • E) ( Figure 1) and covers an approximate area of 7387 km 2 . It is a typical plateau continental climate, with a mean annual temperature and mean annual precipitation of 12.4 • C and 1086.6 mm, respectively. The region is mainly composed of plains, with mountains in the north and elevations ranging from 0 to 1091 m, with an average elevation of 43 m. The eastern and southern parts of the study area are adjacent to the sea, located in a transitional zone between land and sea. It possesses abundant wetland resources and is an ecosystem with distinct environmental features. STN content in this area is influenced by both terrestrial and marine factors. Additionally, this area is an important region within the Bohai economic circle with rapid economic development. The land use in this region is complex and changes rapidly, and there is a large uncertainty in STN content.

Satellite Imagery and Processing
The remote sensing data used for modeling included Landsat-8, Sentinel-1, and Sentinel-2 images; the specific parameter information is shown in Table 1. Three Landsat-8 images covering the study area were downloaded from Geospatial Data Cloud

Satellite Imagery and Processing
The remote sensing data used for modeling included Landsat-8, Sentinel-1, and Sentinel-2 images; the specific parameter information is shown in Table 1. Three Landsat-8 images covering the study area were downloaded from Geospatial Data Cloud (https: //www.gscloud.cn/, accessed on 10 November 2022). Radiometric calibration and atmospheric correction were performed on Landsat-8 images. Two Sentinel-1 images (single-look complex (SLC) products) covering the study area were downloaded from the Google Earth Engine platform and are in IW (interferometric wide swath) mode. Three Sentinel-2 images were downloaded from the European space agency (ESA) website (https://www.esa.int/, accessed on 12 November 2022) as Level-2A products, which were already atmospherically corrected with the Sen2Cor processor and PlanetDEM digital elevation model. These images were then mosaiced and clipped to obtain an optical image covering the study area. A total of 40 remote sensing variables were extracted from all the remote sensing images, including 13 derived from Landsat-8, 22 from Sentinel-2, and 5 from Sentinel-1 (Table 2). Spectral bands 1 to 7 (from 0.43 to 2.29 µm) of Landsat-8 images, 10 spectral bands from the Sentinel-2 images (B2, B3, B4, B5, B6, B7, B8, B8A, B11, and B12), and 2 polarization modes from the Sentinel-1 images (Vertical-vertical polarization and Verticalhorizontal polarization) were utilized. In addition, we calculated vegetation indices from Landsat-8 and Sentinel-2 (normalized difference vegetation index (NDVI), ratio vegetation index (RVI), and difference vegetation index (DVI)), which were reported to be strongly correlated with STN content [29]. The bare soil index (BSI) has a negative correlation with STN content, indicating that a higher degree of land surface bareness corresponds to a lower STN content. BSI can be used as a remote sensing method to monitor land surface bareness and can be combined with soil sampling data to analyze the spatial distribution of STN content [30]. The normalized difference built-up index (NDBI) was selected to reflect the impact of human buildings on STN content, as there are urban and village areas in the study area [31]. The normalized difference water index (NDWI) was initially used for water body monitoring [32], and later studies have used it to predict STN content, which has a strong positive correlation with STN content [33]. Sentinel-2 images have three red edge bands between the visible and near-infrared bands, which are more sensitive to monitoring plant photosynthesis. STN content is closely related to vegetation growth status and type, and thus the red edge bands can be used to estimate STN content [34]. This study used three red edge bands to calculate six commonly used red edge indices to predict STN content. All covariates were resampled to have a similar scale and the same cell size of 30 × 30 m [35,36]. Table 2. List of all remote sensing variables in the study for STN prediction.

Soil Sampling and Analysis
A total of 126 soil samples (0-30 cm) were randomly collected within the study area in 2020 (Figure 1), with a straight-line interval of sampling points for approximately 5 km. The geographic coordinates, vegetation types, land uses, and soil types were duly recorded at each sampling site. Three soil samples were collected and thoroughly mixed at each sampling point to form a composite sample for determining the STN content at that sampling point. All soil samples were air-dried for three weeks, subsequently crushed, and sifted through a 2 mm sieve. STN content was measured using the Kjeldahl method [37].
A descriptive statistical analysis of the target STN content was performed. The statistical properties of the measured STN content at the sampling sites are presented in Table 3. The measured STN content is defined as moderately variable (with a coefficient of variation (CV) value of 59.86%), ranging from 0.052 to 2.396 g·kg −1 , with an average of 0.745 g·kg −1 . The standard deviation (SD) of the STN content was 0.446 g·kg −1 .

Predictive Models
The RF, GBM, and XGBoost methods are currently the most commonly used treebased machine learning methods. Models based on these three methods were implemented through the "train" function in the "caret" package in R-4.2.3-software, and the model parameters were optimized using the grid search method. The final modeling used the parameter combination that resulted in the minimum prediction error.

Random Forest
RF is a commonly used machine learning algorithm. It is a model composed of a random collection of independently trained decision trees [38]. The training data for each decision tree is obtained by random sampling with replacement from the original dataset, and the final model's prediction is the average of all decision tree results ( Figure 2). The advantages of the RF method have the ability to (1) handle nonlinear relationships between multiple predictors, (2) identify and correct overfitting problems, thereby improving prediction accuracy, (3) handle high-dimensional data and automatically deal with missing and outlier values, and (4) output the importance of each predictor to the model's prediction, which further helps understand the influencing factors of the soil properties [39].

Gradient Boosting Machine
The GBM method is also comprises decision trees similar to RF. However, the method is a weighted iterative method for generating decision trees, so the trees i GBM model can be non-independent [40]. The GBM method first generates a decision using the original dataset, then calculates the prediction error, and adjusts the sa weights based on the prediction error. When generating the next round of decision the GBM method prioritizes the sampling of samples with larger prediction errors t hance the model's ability to fit these difficult-to-predict samples. The steps are as fo [41]:

Gradient Boosting Machine
The GBM method is also comprises decision trees similar to RF. However, the GBM method is a weighted iterative method for generating decision trees, so the trees in the GBM model can be non-independent [40]. The GBM method first generates a decision tree using the original dataset, then calculates the prediction error, and adjusts the sample weights based on the prediction error. When generating the next round of decision trees, the GBM Remote Sens. 2023, 15, 2907 7 of 22 method prioritizes the sampling of samples with larger prediction errors to enhance the model's ability to fit these difficult-to-predict samples. The steps are as follows [41]: Step 1: Initialize the model with a constant value: where F 0 (x) is the function initially assumed by GBM, γ is an initial constant, n is the total number of samples, i is the index of the sample, and l(y, F(x)) is the loss function.
Step 2: Looping: m = 1 to M (where m represents the iteration number and M is the predetermined number of iterations, i.e., the number of trees). a.
Compute residuals: where r im represents the residual of the i sample in the m iteration.
b. Fit a decision tree h m (x) to the residuals.
c. Compute multiplier γ m : d. Update the model: Repeat Step 2 iteratively for M times and output F M (x). The GBM method calculates the contribution of each feature to the loss function and then weighs the features according to their contribution. This can increase the model's attention to important features, reduce attention to unimportant features, and improve the accuracy and generalization ability of the models.

Extreme Gradient Boosting
The XGBoost method is an optimization of the GBM method, which introduces regularization to prevent overfitting and improve model generalization performance on top of the original GBM method [42]. The regularization term is added to the loss function, and the new loss function becomes where Ω( f m ) is the regularization term for the m iteration. It also incorporates a novel algorithm for splitting nodes that can speed up training and improve model accuracy [43].

Recursive Feature Elimination
Some of the remote sensing variables may not provide useful information for predicting the target STN content, as they may be redundant or highly correlated. It is necessary to select the subset of features that can best represent the characteristics of the soil to improve the prediction accuracy of the models and to reduce computation and data storage costs. Recursive feature elimination (RFE) is a commonly used feature selection method that can be used to determine which remote sensing variables are most important in building STN content prediction models [44]. In this study, we used the "rfFuncs" method to sort the model as an argument, making 40 iterations the preset number of features, which decreased from 40 to 1. RFE removed the least important feature and retrained the model using the remaining features at each iteration. The trained model was then used to predict the validation set, and the root mean square error (RMSE) was calculated. The number of features corresponding to the smallest RMSE was selected, and the feature variables were outputted. RFE is performed using the "rfe" function of the R software.

Model Validation
We constructed STN content models using three different machine learning methods and various combinations of predictor variables. The combinations of the different factors are shown in Table 4. Model I, Model II, and Model III used Landsat-8-, Sentinel-1-, and Sentinel-2-derived predictors, respectively, to predict STN content. Model IV and Model V were combinations of Landsat-8-derived predictors with Sentinel-1-and Sentinel-2-derived predictors, respectively, and Model VI used a combination of Sentinel-1-and Sentinel-2derived predictors. Model VII included all predictor variables. Figure 3 shows an overview of the flowchart for STN content mapping using these experimental models. For the predictive performance of these models, we used a 10-fold cross-validation method [45]. For the 10-fold cross-validation, we randomly divided the observed dataset into 10 groups [46]. In each of the 10 folds, one group was designated as the test dataset and the other nine groups were used as the training set [47]. Three validation criteria were calculated to evaluate the performance of the model: the RMSE, the mean absolute error (MAE), and the coefficient of determination (R 2 ). These validation criteria are calculated from the following [25]: where n represents the number of samples,ŷ l and y i represent the predicted and observed values at site i, respectively, and y represent the mean of observed values.
where represents the number of samples, and represent the predicted and observed values at site , respectively, and represent the mean of observed values.

Model Evaluation and Comparison
The performances of the RF, GBM, and XGBoost methods based on different combinations of predicting STN content are shown in Table 5. The different methods and variable combinations significantly affected the modeling performance. For the RF and GBM methods, Model I (R 2 = 0.446 vs. R 2 = 0.410, respectively), Model II (R 2 = 0.411 vs. R 2 = 0.391, respectively), and Model III (R 2 = 0.409 vs. R 2 = 0.394, respectively) were better predicted by the RF, indicating that the RF method is better than GBM in predicting STN content using single-type remote sensing data. However, the GBM method performed better than RF in Model IV (R 2 = 0.479 vs. R 2 = 0.459, respectively), Model V (R 2 = 0.496 vs. R 2 = 0.463, respectively), Model VI (R 2 = 0.488 vs. R 2 = 0.457, respectively), and Model VII (R 2 = 0.533 vs. R 2 = 0.475, respectively), indicating that the GBM method is more suitable than RF for predicting STN content using multiple source remote sensing data. Among the three machine learning methods, whether using single or multiple data types as prediction variables, the XGBoost method has the highest prediction accuracy for STN content. When comparing different combinations of predictive variables, Model I (R 2 = 0.446 and R 2 = 0.410 for RF and GBM methods, respectively) performed better than Model II (R 2 = 0.411 and R 2 = 0.391 for the RF and GBM methods, respectively) and Model III (R 2 = 0.409 and R 2 = 0.394 for the RF and GBM methods, respectively), indicating that Landsat-8 images have a better predictive capability than the Sentinel-1 and Sentinel-2 images when modeling using the RF and GBM methods; Model II and Model III have similar prediction levels, indicating that Sentinel-1 and Sentinel-2 images have similar predictive capabilities. For the models that were established using the XGBoost method, the R 2 of Model I and Model III is 15.5% (from 0.431 to 0.498) and 21.6% (from 0.431 to 0.524) higher than that of Model II, respectively, indicating that optical imagery performs better than SAR imagery. When compared with single-type remote sensing data, the combination of different types of data can improve prediction accuracy, and the addition of each type of data has a different degree of improvement in model accuracy. For example, when Sentinel-1 and Sentinel-2 predictors were added to Model I to form Model IV and Model V, the highest R 2 increased by 16.83% and 20.98%, respectively. When Sentinel-2 predictors were added to Model II to form Model VI, the highest R 2 increased by 37.5%. This indicates that the added data contain valuable information that is different from the original data. This improvement is similar across all three machine learning methods.
The models built using all prediction factors have the highest prediction accuracy. Model VII, which combines three types of remote sensing data as prediction factors, showed the most significant improvement. For example, when compared to Model II and constructed from a single data source, Model VII, based on the XGBoost method, improved the R 2 value from 0.431 to 0.627, an increase of 45.48% compared to Model IV, which combines two types of data. The R 2 value of Model VII (0.627) is 15.05% higher than that of Model IV (0.545). Model VII is based on three methods (R 2 = 0.475, R 2 = 0.533, and R 2 = 0.627 for the RF, GBM, and XGBoost methods, respectively), which can explain the variation in STN content of 47.5%, 53.3%, and 62.7%, respectively. From the distribution of the measured value and the predicted value scatterplot (Figure 4), the STN content predicted by the XGBoost method is closer to the measured value, and the fitted straight line between the measured value and the predicted value is closer to the 1:1 line, followed by GBM, and RF is the worst.

Relative Importance of Variables
Model VII was built based on three machine learning methods using a combination of all predictive variables. The predictive variables sorted by relative importance are shown in Figure 5 (percentages were used to enhance comparability). Variables with less than 1% importance are not displayed in the graph, as they may be due to chance. Variables are ranked roughly the same in importance. For example, six out of the top nine important variables were duplicated in the GBM and XGBoost methods, namely VV, L_B5, S_NDVIRE2, VH, S_NDVIRE3, and S_B5. Five duplicated variables were found in the RF and XGBoost methods, namely L_NDWI, S_NDVIRE2, L_B5, VH, and VV. L_NDWI were the most significant explanatory variables in both models, accounting for 14.24% and 21.44% of the relative importance in predicting STN content. Four variables were duplicated among the top nine most important variables in all three methods, namely, VV, L_B5, S_NDVIRE2, and VH.
Model VII, which was constructed using the RF method, showed that Landsat-8 imagery (relative importance of 63%) is the main explanatory variable for STN content, followed by Sentinel-2 (24%) and Sentinel-1 (13%). Similarly, the XGBoost method-established Model VII also indicates that Landsat-8 has the highest relative importance (44%), followed by Sentinel-2 (33%), and Sentinel-1 has the lowest importance (23%). This suggests that Landsat-8 imagery has a stronger explanatory power for STN content than Sentinel series imagery for both the RF and XGBoost methods. However, GBM method-established Model VII, Landsat-8, and Sentinel-2 have similar explanatory powers. For all the models established using the RF, GBM, and XGBoost methods, the relative importance

Relative Importance of Variables
Model VII was built based on three machine learning methods using a combination of all predictive variables. The predictive variables sorted by relative importance are shown in Figure 5 (percentages were used to enhance comparability). Variables with less than 1% importance are not displayed in the graph, as they may be due to chance. Variables are ranked roughly the same in importance. For example, six out of the top nine important variables were duplicated in the GBM and XGBoost methods, namely VV, L_B5, S_NDVIRE2, VH, S_NDVIRE3, and S_B5. Five duplicated variables were found in the RF and XGBoost methods, namely L_NDWI, S_NDVIRE2, L_B5, VH, and VV. L_NDWI were the most significant explanatory variables in both models, accounting for 14.24% and 21.44% of the relative importance in predicting STN content. Four variables were duplicated among the top nine most important variables in all three methods, namely, VV, L_B5, S_NDVIRE2, and VH.
Model VII, which was constructed using the RF method, showed that Landsat-8 imagery (relative importance of 63%) is the main explanatory variable for STN content, followed by Sentinel-2 (24%) and Sentinel-1 (13%). Similarly, the XGBoost method-established Model VII also indicates that Landsat-8 has the highest relative importance (44%), followed by Sentinel-2 (33%), and Sentinel-1 has the lowest importance (23%). This suggests that Landsat-8 imagery has a stronger explanatory power for STN content than Sentinel series imagery for both the RF and XGBoost methods. However, GBM method-established Model VII, Landsat-8, and Sentinel-2 have similar explanatory powers. For all the models established using the RF, GBM, and XGBoost methods, the relative importance of Sentinel-1 is the lowest, at 13%, 24%, and 23%, respectively. This indicates that optical imagery is more helpful for predicting STN content. The same rules were observed in others, from Model I to Model VI (Figures A1-A3).

Spatial Distribution Pattern of STN Content
Based on the RF, GBM, and XGBoost methods, the established Model VII was selected to predict the STN content in the entire study area, and a spatial distribution map of the STN content in the study area was drawn ( Figure 6). The spatial patterns of STN content predicted by the three methods are similar, and a strong spatial heterogeneity for STN content was observed on all of the distribution maps, with higher STN content in the northern part of the study area and lower content in coastal areas. Based on the statistical analysis, the predicted STN contents from different models show similarities. For instance, in coastal areas, the majority of the pixels have STN contents ranging from 0.4-0.5 g•kg −1 . As we move from coastal wetlands to farmland areas to mountainous areas, the peak of the distribution curve shifts towards higher STN contents, indicating that the STN content in inland areas is significantly higher than that in coastal areas.
The three methods predict STN content in the study area, and the descriptive statistics are shown in Table 6. The SD values for STN content predicted by the RF model are lower than that predicted by GBM and XGBoost, indicating the robustness of the RF model is the highest. The predicted average STN content of each model is higher than the actual value.

Method
Minimum (g•kg −1 ) Maximum (g•kg −1 ) Mean (g•kg −1 ) SD (g•kg −1 ) Figure 5. The relative importance of variables used for the STN content prediction in Model VII based on RF, GBM, and XGBoost methods.

Spatial Distribution Pattern of STN Content
Based on the RF, GBM, and XGBoost methods, the established Model VII was selected to predict the STN content in the entire study area, and a spatial distribution map of the STN content in the study area was drawn ( Figure 6). The spatial patterns of STN content predicted by the three methods are similar, and a strong spatial heterogeneity for STN content was observed on all of the distribution maps, with higher STN content in the northern part of the study area and lower content in coastal areas. Based on the statistical analysis, the predicted STN contents from different models show similarities. For instance, in coastal areas, the majority of the pixels have STN contents ranging from 0.4-0.5 g·kg −1 . As we move from coastal wetlands to farmland areas to mountainous areas, the peak of the distribution curve shifts towards higher STN contents, indicating that the STN content in inland areas is significantly higher than that in coastal areas.  The three methods predict STN content in the study area, and the descriptive statistics are shown in Table 6. The SD values for STN content predicted by the RF model are lower than that predicted by GBM and XGBoost, indicating the robustness of the RF model is the highest. The predicted average STN content of each model is higher than the actual value.

Accuracy and Influencing Factors of STN Content Prediction Models
The results demonstrate that the prediction methods, different types of data, and different combinations of data significantly influence the STN content predictive accuracy. The study did not find that the RF method-established models consistently outperformed GBM in predicting STN content using different variable combinations. Therefore, it is necessary to calibrate and evaluate competitive prediction models based on specific experimental datasets under different model combinations. The XGBoost method outperforms the RF and GBM methods in terms of prediction accuracy, a finding supported by Tien Dat Pham's research [27]. However, Zhang et al. used RF, GBM, and XGBoost to study STN content in tobacco planting areas and found that GBM performed the best, followed by RF, with XGBoost performing the worst [48]. Some scholars have found that the three methods perform similarly [49]. This discrepancy may be due to differences in STN sample quantity as well as the types of remote sensing variables used. There is no consistent conclusion about the model performance of the RF, GBM, and XGBoost methods. It seems that no single machine learning method is most suitable for all ecosystems, so it is important to choose different methods based on different regions and remote sensing variables. As the RF method calculates an average value of the output values from multiple trees as the model's prediction result, it is not sensitive to outliers [50]. This means that the RF method ignores the effect of extremely high or low STN content values on the prediction of STN content in the study area, resulting in a small predicted range for STN content in the entire region. GBM and XGBoost are both iterative models, with each model's prediction based on the residuals of the previous model. The models are sensitive to outliers, as a large outlier may affect the residuals of each model and result in a wider predicted range of STN content [51]. Based on measured soil data, the STN content ranged from 0.052 to 2.396 g·kg −1 . Among the three machine learning algorithms used in this study, the XGBoost method was found to be more accurate, supporting the results that XGBoost has better accuracy.
Our research findings demonstrate the crucial importance of three types of remote sensing imagery, namely Landsat-8, Sentinel-1, and Sentinel-2, in predicting STN content. The accuracy of the model based on the derived variables extracted from different remote sensing images is different. Although both Landsat-8 and Sentinel-2 are optical remote sensing images, the information contained in the images is different due to differences in their center wavelength, bandwidth, and overlapping bands [15], and the difference in image acquisition time also leads to differences in the information contained in the images [52], which can lead to different prediction capabilities for STN content. The models (Model II) based on Sentinel-1-derived variables had lower accuracy compared to the models (Model I and Model III) based on optical image-derived variables. This suggests that the predictive ability of optical images is superior to SAR images in the study area, which is consistent with previous research [53]. However, the Sentinel-1 data helped to improve the accuracy of the models, and the study found that when the predictors extracted from SAR images were added, the model accuracy improved, indicating that Sentinel-1 imagery contains useful information beyond Landsat-8 and Sentinel-2 [15]. There is a study that found the inclusion of Sentinel-1 imagery improves the model accuracy, contributing 9% and 7% to the RF and BRT models, respectively, supporting the results of this experiment [53]. The inclusion of different sensor data in the model significantly improved its accuracy, indicating that Landsat 8, Sentinel-1, and Sentinel-2 images contain different valuable information. Previous studies on predicting soil properties mainly used a single sensor, such as Landsat [54] and Sentinel-2 [35,45], without considering the feasibility of radar sensors. In this study, the better prediction accuracy obtained from the combination of optical and SAR images demonstrates the usefulness of SAR data in predicting STN content. The combination of optical and radar sensors has great potential in predicting soil properties [55].
Among all the prediction models, Model VII, based on the XGBoost method, has the highest prediction accuracy and can explain 62.7% of the variability in STN content. Our prediction model has achieved higher accuracy compared to other scholars' predictions of STN content. For example, Wadoux et al. established an RF model using French LUCAS data, which can only explain 20% of the variability in STN content [54]. ZHOU et al. used data from the Second National Land Survey from 1979 to 1985 to construct XGBoost, RF, and weighted model averaging methods to predict STN content across China, with R 2 values of 0.34, 0.38, and 0.41, respectively [56]. Although public soil datasets have wide coverage and rich attribute information for the sampling points, they have poor timelines due to their long sampling time and cannot be matched with recent remote sensing images, thus, they can only be used to invert soil properties during the sampling period. When compared to using public datasets, field sampling provides controllable data for specific experiments, with the spatial and temporal consistency between the samples and remote sensing images being the most important advantage, which had a certain effect on improving the accuracy of the model.

Relative Importance of Variables
In this study, among the RF, GBM, and XGBoost method-established optimal models, the prediction factors provided by Landsat-8 data accounted for 63%, 37%, and 44% of the importance of all variables, respectively. The prediction factors provided by Sentinel-2 data accounted for 24%, 39%, and 33% of the importance of all variables, respectively. The prediction factors provided by Sentinel-1 data accounted for 13%, 24%, and 23% of the importance of all variables, respectively. It can be seen that optical images are the most important in explaining the variability of STN content. Among the two optical images, Landsat-8 has greater importance than Sentinel-2 in all but the GBM prediction results, where their importance is similar. This suggests that Landsat-8 data have a greater impact on the study area than Sentinel-2. Additionally, when compared to optical images, SAR images have lower importance in predicting STN content. This result is supported by the study of Zhou et al. [53].
The spectral bands, the remote sensing indices of optical imagery, and the backscatter coefficients of radar imagery are extracted through remote sensing images to help explain the spatial variation in STN content in the soil-vegetation system. They can capture the relationship between soil properties and vegetation to reflect soil information to some extent. Remote sensing images represent a valuable dataset for explaining spatial changes in the soil in natural vegetation areas [57]. In the RF, GBM, and XGBoost prediction models, remote sensing indices accounted for 69%, 45%, and 55% of the model's contribution, respectively. Remote sensing indices contribute more to the models than band reflectance. Specifically, in the RF and XGBoost models, remote sensing indices have a higher contribution rate than band reflectance, indicating that remote sensing indices can better characterize soil information and have higher values for the prediction of STN content. In addition, band reflectance is more important than remote sensing indices in GBM models, indicating that the role of band reflectance is relatively more important in the GBM model. In the RF and XGBoost models, the importance of NDWI is highest, and soil moisture promotes the accumulation of ecosystem STN content [58]. Xu's research also proves that NDWI strongly correlates with STN content [33]. The sampling sites in this study are mostly located in intertidal flats, paddy fields, and marsh land cover types, where soil moisture is high, which is expected to result in the high importance of NDWI. Remote sensing images are more sensitive to vegetation and are indirectly sensitive to soil properties; the data were collected from September to October, when vegetation was flourishing, and the calculated vegetation indices were relatively high. The importance of vegetation indices is high in this study, among which the red edge index S_NDVIRE2 has a large contribution to the models. The red edge indices are recognized as the most suitable remote sensing indices for reflecting vegetation growth [59], which means that they can estimate soil properties better through the vegetation medium. This view is also supported in this study.

Spatial Distribution of STN Content
The spatial distribution of STN content predicted in this study is similar to that of the 0-20 cm STN content data set in China by Zhou et al. [56]. The spatial distribution of STN content predicted by the three modeling methods is also similar, further indicating that the results of this study are in line with reality. High levels of STN content were mainly distributed in the northern mountainous areas that were covered with dense vegetation. Correspondingly, low levels of STN content were mainly found in the southern and central regions that were dominated by high levels of human activity areas such as coastal and urban areas, indicating that areas with dense tree cover are more conducive to the accumulation of STN content, which is consistent with the results of Zhou et al. [60]. Cropland STN content is second only to the northern mountainous areas, which may be due to the fertilization of farmland, which causes some nitrogen elements to infiltrate into the soil [15,61]. In addition, a large amount of agricultural waste, residue, and feces will be produced during the agricultural production process, and these materials will degrade into organic matter and release nitrogen elements, further increasing the STN content [62]. In terms of the prediction results, the spatial distribution of STN based on the XGBoost method shows that the STN content ranges from 0 to 2.01 g·kg −1 . These distribution ranges are consistent with those presented in the STN maps produced by Xu et al. [33] and Li et al. [63].
However, there are some uncertainties in our study. Firstly, time-series images provide more information compared to single-time remote sensing images, reduce uncertainty and improve model accuracy, and future studies might consider extracting time-series remote sensing variables for predicting STN content [18]. This study resampled the spatial resolution of all the remote sensing variables to 30 m without considering the impact of spatial resolution on modeling and inversion accuracy. However, different spatial resolutions can lead to different mixes of land features, thereby affecting the accuracy of the modeling and inversion. Therefore, it is crucial to determine an appropriate spatial resolution for predicting soil properties. In the next step of the research, we will pay attention to considering the impact of spatial resolution on soil property prediction and explore suitable strategies for selecting spatial resolution to enhance the accuracy and reliability of the model.

Conclusions
This study combined three commonly used remote sensing images, including two multispectral images (Landsat-8 and Sentinel-2) and one SAR image (Sentinel-1), and used three decision tree-based machine learning methods (RF, GBM, and XGBoost) to predict the STN content in a coastal area. The spatial distribution of the STN content was mapped. Our conclusions are summarized as follows:

•
The application of SAR and optical images proved useful for predicting STN content, and their combination showed enhanced model accuracy. The RF, GBM, and XGBoost methods demonstrated maximum improvements of 16%, 36%, and 45%, respectively; • The XGBoost method had higher accuracy than the RF and GBM methods. The optimal model was built using the XGBoost method, with an R 2 of 0.627, RMSE of 0.127 g·kg −1 , and an MAE of 0.092 g·kg −1 ; • Optical imagery is more helpful than SAR imagery in predicting STN content. In the models established by the RF and XGBoost methods, Landsat-8 had the highest relative importance (63% and 44%, respectively), followed by Sentinel-2 (24% and 33%, respectively). In the model established by the GBM method, the importance of Landsat-8 and Sentinel-2 was similar but higher than that of Sentinel-1; • The STN content predicted by the three models has a certain degree of similarity for spatial distribution. The predicted range of STN content is from 0 to 2.01 g·kg −1 . These maps showed significant spatial variability. The STN content is high in the densely forested areas in the north and low in the paddy wetlands in the southeast.