A New Method for Estimating Soil Fertility Using Extreme Gradient Boosting and a Backpropagation Neural Network

Abstract

Soil fertility affects crop yield and quality. A quick, accurate evaluation of soil fertility is crucial for agricultural production. Few satellite image-based evaluation studies have quantified soil fertility during the crop growth period. Therefore, this study proposes a new approach to the quantitative evaluation of soil fertility. Firstly, the optimal crop spectral variables were selected using the integration of an extreme gradient boosting (XGBoost) algorithm with variance inflation factor (VIF). Then, based on the optimal crop spectral variables where the red-edge indices were introduced for the first time, the estimation models were developed using the backpropagation neural network (BPNN) algorithm to assess soil fertility. The model was finally adopted to map the soil fertility using Sentinel-2 imagery. This study was performed in the Conghua District of Guangzhou, Guangdong Province, China. The results of our research are as follows: (1) five crop spectral variables (inverted red-edge chlorophyll index (IRECI), chlorophyll vegetation index (CVI), normalized green-red difference index (NGRDI), red-edge position (REP), and triangular greenness index (TGI)) were the optimal variables. (2) The BPNN model established with optimal variables provided reliable estimates of soil fertility, with the determination coefficient (R²) of 0.66 and a root mean square error (RMSE) of 0.17. A nonlinear relation was found between soil fertility and the optimal crop spectral variables. (3) The BPNN model provides the potential for soil fertility mapping using Sentinel-2 images, with an R² of 0.62 and an RMSE of 0.09 for the measured and estimated results. This study suggests that the proposed method is suitable for the estimation of soil fertility in paddy fields.

1. Introduction

Soil fertility refers to the capacity for soil to offer different nutrients for crop growth, significantly affecting crop yield [1,2]. Scientific and reasonable evaluation of soil fertility can provide a reference for land-use planning and fertilizer prescriptions, guiding agricultural production [3]. Soil fertility is typically evaluated using soil sampling and laboratory chemical analysis to determine the chemical properties (soil pH, soil organic matter (SOM), available phosphorus (AP), available potassium (AK), and total nitrogen (TN)). It proves to be time-consuming and expensive for deriving spatially explicit estimates across a large study area [4,5]. In contrast, remote sensing techniques can be used to evaluate soil fertility using spectral variables and obtain spatially explicit estimates of soil fertility relatively quickly.

Current research methods for soil fertility estimation in arable land using remote sensing can be divided into soil spectrum-based and crop spectral index-based methods. The soil spectrum-based method evaluates the correlation between soil spectral indices and soil fertility. Rossel et al. [6] evaluated the soil fertility of sugarcane using a decision tree algorithm based on visible/near-infrared (vis-NIR) soil spectra and terrain attributes, which showed that the method could improve the efficiency of soil fertility evaluation. Munnaf et al. [7] developed a method for assessing soil fertility indices based on online vis-NIR soil spectroscopy. The results demonstrated that the method can be effective in assessing soil fertility. Yang et al. [8] evaluated soil fertility of paddy fields in southern China using vis–NIR soil spectral indices and partial least-squares regression. According to the results, vis–NIR spectroscopy improved the efficiency of estimating soil fertility. Wang et al. [9] explored the association of a comprehensive soil fertility index with soil spectral curves for agricultural soils in different states under optimal observation conditions. The results indicated that soil spectral indices were suitable for estimating soil fertility. However, arable land in southern China has few bare soil areas; therefore, it is difficult to obtain soil fertility using soil spectral indicators.

The crop spectral index-based method evaluates the correlation between spectral vegetation indices and soil fertility. Zeeshan et al. [10] found that a higher normalized difference vegetation index (NDVI) value corresponds to better soil fertility. Wang et al. [11] set thresholds for the NDVI based on cotton growth and performed density partitioning to obtain information on different levels of soil fertility in cotton fields. The research results showed that the NDVI could be used to estimate soil fertility. Duan et al. [12] used the mean and the coefficient of variation (CV) of the NDVI of arable land for three consecutive years to determine the magnitude and stability of soil fertility, respectively. A larger mean value and a smaller CV indicated higher soil fertility of the arable land. In these studies, soil fertility was estimated qualitatively. Moreover, only vis–NIR spectral indices were used to estimate soil fertility, and other spectral indices variables (e.g., red-edge) were not considered; thus, the accuracy of soil fertility estimation was not very high.

The contribution of this work was to design a novel method for the quantitative estimation of soil fertility in paddy fields using crop spectral variables. In this method, the machine learning algorithms (extreme gradient boosting (XGBoost) and backpropagation neural network (BPNN)) were used for the inaugural time to quantitatively assess soil fertility. For the first time, the red-edge index contributes to the soil fertility evaluation model, a method which provides the potential to quantitatively estimate soil fertility at both the soil sample level and regional scale. The study was performed in the Conghua District of Guangzhou, Guangdong Province, China using Sentinel-2 data.

2. Materials and Methods

In this study, a novel method for quantitative estimation of soil fertility was developed. The method combines the XGBoost algorithm for variable selection with the BPNN algorithm for soil fertility estimation. Figure 1 is the flow chart for the method. The methodological framework involves data collection, data pre-processing, determination of optimal crop spectral indicators, construction of the soil fertility estimation model, and the spatial mapping of soil fertility.

2.1. Study Area

The study area is located in the Conghua District of Guangzhou, Guangdong Province, China (113°17′E–114°04′E, 23°22′N–23°56′N). According to the statistical yearbook of the Conghua District, the annual average temperature was 22.3 °C, the annual rainfall was 1297.5 mm, and the annual sunshine hours were 1976.8 h in 2021. The main soil types in the Conghua District are red loam, yellow loam, lateritic red soil, and rice soil. The total area of arable land in Conghua is 205 km²; most of it is located in the northwestern area of the Conghua District. Arable land planted with rice was chosen because rice accounts for the highest proportion of cropland in the study area.

2.2. Data Sources and Preprocessing

2.2.1. Soil Samples

A total of 150 sample points were obtained in the rice-growing area of Conghua District in September 2017 using stratified random sampling and considering different land units, soil types, land-use patterns, and agricultural facility construction levels. The 150 points fell into three groups: 90 samples (yellow plots of Figure 2) were introduced for model training, 30 samples (black plots of Figure 2) for evaluating model accuracy, and 30 samples (red plots of Figure 2) for validating mapping accuracy. At each sample point, five soil sub-samples of the topsoil (0–20 cm) were collected in an X shape, mixed, and used as the soil sample of this point. Basic information, such as the crop type, was recorded during sample collection, and the latitude and longitude of the sample points were recorded by a Global Positioning System (GPS) receiver. The samples underwent air drying at room temperature and were milled and screened via a 100-mesh sieve (0.15 mm) to observe soil properties. We used the regulation of classification of paddy soil fertility and fertilizer technology (DB43/T 2087-2021) and other information [5] to obtain the soil fertility index (SFI) of the 150 sample points using Equation (1):

$$SFI={\displaystyle \sum }{W}_i\times N_i$$

(1)

where $W_i$ and $N_i$ represent the weight coefficient and membership degree of the ith indicator (soil pH, SOM, AP, AK, and TN). We followed the measurement procedures and methods described in Lu [13].

Table 1. Descriptive statistics of the soil properties.

Soil Properties	Min	Max	Mean	SD	Skewness	Kurtosis	CV (%)
pH	4.90	8.20	5.84	0.44	1.23	4.95	7.52
SOM	6.42	68.90	23.72	8.88	1.10	3.43	37.43
TN	0.37	2.14	0.86	0.42	1.22	0.83	47.95
AP	6.80	140.8	43.89	24.64	1.14	1.35	56.15
AK	2.00	235.00	74.30	49.19	1.06	0.54	66.21

shows the descriptive statistics for soil properties. $W_i$ was obtained from correlation analysis between the indicators, and $N_i$ was obtained from an affiliation function of the ith indicator. The boxplot of the SFI of the sample points is presented in Figure 3.

2.2.2. Satellite Image Data and Preprocessing

According to the rice key growth stages in the study area, 6 Sentinel-2 images (Figure 4) ranging from 17 September 2017 to 11 November 2017 were acquired. The band information (13 bands) and the spatial resolutions are listed in

Table 2. Band information of the Sentinel-2 images.

Band	Description	CW (nm)	SR (m)	Band	Description	CW (nm)	SR (m)
B1	Coastal aerosol	443	60	B8	NIR-1	842	10
B2	Blue	490	10	B8A	NIR-2	865	20
B3	Green	560	10	B9	Water vapor	945	60
B4	Red	665	10	B10	SMIR-Cirrus	1375	60
B5	Red edge-1	705	20	B11	SMIR-1	1610	20
B6	Red edge-2	740	20	B12	SMIR-2	2190	20
B7	Red edge-3	783	20

Note: Central wavelength, CW; Spatial resolution, SR.

. The images were acquired from the Google Earth Engine platform. They had been preprocessed, including radiometric, geometric, and atmospheric corrections and orthorectification.

2.3. Methods

2.3.1. Determination of Optimal Crop Spectral Variables for Estimating SFI

(1) Acquisition of crop spectral variables

The characteristics related to crop growth and the fertility of arable land, which is closely related to vegetation indices [14]. Twenty-seven crop spectral variables were calculated on the Google Earth Engine Platform (https://earthengine.google.com/ (accessed on 12 April 2022)) using the Sentinel-2 images. The details of the crop spectral variables are listed in

Table 3. Crop spectral variables for estimating SFI.

Vegetation Index	Formulation in Sentinel-2	References	Vegetation Index	Formulation in Sentinel-2	References
NDVI	(B8 − B4)/(B8 + B4)	Haboudane et al. [15]	MCARI	((B5 − B4) − 0.2 × (B5 − B3)) × (B5/B4)	Daughtry et al. [16]
MTCI	(B6 − B5)/(B5 − B4)	Dash et al. [17]	MCARI1	1.2 × (2.5 × (B8 − B4) − 1.3 × (B8 − B3))	Haboudane et al. [15]
MGRVI	((B3)2 − (B4)2)/((B3)2 + (B4)2)	Bendig et al. [18]	MCARI2	1.5 × (2.5 × (B8 − B4) − 1.3 × (B8 − B3)/((2.0 × B8 + 1)2) − (6.0 × B8 – 5 × ((B4)0.5)) − 0.5)0.5
REP	705 + 35 × ((((B7 + B4)/2) − B5)/(B6 − B5))	Frampton et al. [19]	MTVI1	1.2 × (1.2 × (B8 − B3) − 2.5 × (B4 − B3))
IRECI	(B7 − B4)/(B5/B6)		MTVI2	1.5 × (1.2 × (B8 − B3) − 2.5 × (B4 − B3)/((2.0 × B8 + 1)2) − (6.0 × B8 − 5 × ((B4)0.5)) − 0.5)0.5
RVI	B8/B4	Birth et al. [20]	NDREI	(B8 − B5)/(B8 + B5)	Gitelson et al. [21]
EVI	2.5 × ((B8 − B4)/(B8 + 6 × B4 − 7.5 × B2 + 1))	Huete [22]	NGRDI	(B3 − B4)/(B3 + B4)	Tucker [23]
DVI	B8 − B4	Jordan [24]	NIRv	((B8 − B4)/(B8 + B4)) × B8	Badgley et al. [25]
SAVI	(B8 − B4) × 1.5/(B8 + B4 + 0.5)	Huete [26]	OSAVI	(B8 − B4)/(B8 + B4 + 0.16)	Rondeaux et al. [27]
MASVI	0.5 × (2 × B8 + 1 − ((2 × B8 + 1)2 − 8 × (B8 − B4))0.5)	Qi et al. [28]	SELI	(B8A − B5)/(B8A + B5)	Pasqualotto et al. [29]
CIG	(B8/B3) − 1	Anatoly et al. [30]	TCARI	3 × ((B5 − B4) − 0.2 × (B5 − B3) × (B5/B4))	Haboudane et al. [31]
CIRE	(B8/B5) − 1		TCI	1.2 × (B5 − B3) − 1.5 × (B4 − B3) × (B5/B4)0.5
CVI	(B8 × B4)/((B3)2)	Meng et al. [32]	TGI	−0.5 × (190 × (B4 − B3) − 120 × (B4 − B2))	Hunt et al. [33]
TVI	0.5 × (120 × (B8 − B3) − 200 × (B4 − B3))	Broge et al. [34]

.

(2) Determining the optimal crop spectral variables

Feature selection is a key step in regression analysis to improve prediction accuracy and reduce redundant indicators. Compared to other screening algorithms (e.g., random forest, deep learning), related research [35,36] showed that the extreme gradient boosting (XGBoost) algorithm offers characteristics such as interpretability, computationally efficient, and being less prone to over-fitting under small sample size condition. Because of the relatively small sample size in this study (n = 150), the XGBoost algorithm was employed for selecting optimal crop spectral variables to estimate SFI. The XGBoost is a machine learning algorithm based on a decision-tree ensemble and gradient boosting framework. It gives importance scores for each feature (FI) in each iteration of the training process, so as to indicate the importance of each feature to the training of the model. The FI is directly used as a basis for feature selection [37]. The specific screening steps are as follows [38]:

(1)

A classification model is built on the basis of all the features.

(2)

Based on the information from the generated model process, the FI is obtained and ranked in descending order. FI is calculated as follows [38]:

$$\mathrm{IG}\left(\mathrm{T},\mathrm{F}\right)=\mathrm{H}\left(\mathrm{T}\right)-\mathrm{H}\left(\mathrm{T}|\mathrm{F}\right)=-{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{J}}{\mathrm{p}}_{\mathrm{i}}{\log }_2{\mathrm{p}}_{\mathrm{i}}-{\displaystyle \sum }_{\mathrm{F}}\mathrm{p}\left(\mathrm{F}\right)*{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{J}}\mathrm{p}\left(\mathrm{i}|\mathrm{F}\right){\log }_2\mathrm{p}\left(\mathrm{i}|\mathrm{F}\right),$$

(2)

where $\mathrm{H}\left(\mathrm{T}\right)$ and $\mathrm{H}\left(\mathrm{T}|\mathrm{F}\right)$ denote the entropy of parent and child nodes on the basis of the F-feature segmentation, separately; ${\mathrm{p}}_{\mathrm{i}}$ represents the score of the labeled samples at the node.

(3)

A subset of features is generated by selecting a number of features with the highest FI values.

(4)

Classification experiments were performed on the subset of features to examine their classification ability.

(5)

Repeat steps (3) and (4) until all features have been selected.

(6)

Check the classification for all subsets and choose the optimal subset of features (namely, the subset having relatively high area under the curve values and fewer features).

2.3.2. Model Construction

Two algorithms (multiple linear regression (MLR) and a backpropagation neural network (BPNN)) were employed for determining the association of optimal crop spectral variables with the SFI, and accuracy assessments were carried out. We present a brief summary of each algorithm.

(1) Multiple linear regression model

MLR refers to a linear regression model that uses multiple explanatory variables to depict the linear correlation between independent and dependent variables. The model can describe the degree of influence of a variable on the soil properties. MLR has been widely used for predicting soil properties [39,40,41]. In this study, MLR was performed in SPSS software. The definition of MLR is as follows [40]:

$$\mathrm{Y}={\displaystyle \sum }_{i=1}^n{\beta }_i{\check{X}}_i+a$$

(3)

where ${\check{X}}_i$ denotes the ith optimal crop spectral variable, ${\beta }_i$ represents the regression coefficient of the ith variable, and $a$ represents the intercept.

(2) Backpropagation Neural Network model

The BPNN model is a multilayer feed-forward network that uses fine-tuning of the weights according to the error rate of the former epoch. It comprises an input layer, an output layer, and several hidden layers (Figure 5). The model uses a gradient descent algorithm and backpropagation algorithm to iteratively adjust the weights and biases of the network. The training ends when the predicted value is as close as possible to the actual value. The learning process comprises forward propagation of the input signal and backward propagation of the error [42,43,44].

(1) Forward propagation

As for neural networks, forward propagation requires the calculation of both neuron input and output values. The output value ($o_j$) was written as:

$$O_j=f_i{\displaystyle \sum }({\omega }_{ji}{O}_i+{\theta }_j)$$

(4)

where $o_k$ is the output layer information (each of the SFI); $f_i$ means the transfer function of the hidden layer to the output layer, where the Purelin function is chosen as $f_i$ by the current research [42]; ${\omega }_{kj}$ suggests the weight of the hidden layer to the output layer; ${\theta }_j$ represents the threshold value in the output layer.

When the output value ($o_k$) of the hidden layer was transferred to the output layer, the $o_k$ was written as:

$$O_k=f_j{\displaystyle \sum }({\omega }_{kj}{O}_j+{\theta }_k)$$

(5)

where $o_i$ means the input layer information (crop spectral variables); $o_j$ refers to the hidden layer information; ${\omega }_{kj}$ denotes the weight of the input layer to the hidden layer; $f_j$ signifies the transfer function of the input layer to the hidden layer, where the trainlm function is selected by the current research [42]; ${\theta }_k$ suggests the threshold value in the hidden layer.

(2) Error back propagation

In cases where the predicted value is different significantly from the measured value, the difference can be transferred to the error in the backpropagation process. The backpropagation process utilizes the Levenberg–Marquardt algorithm for modifying connection weights from the output layer to the input layer to decrease the mean squared error (MSE).

$$MSE=\frac{1}{N}{\displaystyle \sum }{\left(O-O_k\right)}^2$$

(6)

where $o$ and $O_k$ represent the measured and predicted SFI, respectively; $N$ means the number of training samples.

In this study, the number of neuron nodes of the hidden layer (H) is decided using the empirical formula [43]:

$$H=2n+1$$

(7)

where $n$ refers to the number of input units.

2.3.3. Accuracy Metrics

The coefficient of determination (R²), concordance correlation coefficient (CCC), root mean square error (RMSE), and the ratio of performance to interquartile range (RPIQ) were employed for assessing the performance of the SFI estimation models using the training and validation set. The metrics are expressed as follows [45,46]:

$${\mathrm{R}}^2=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\overline{\mathrm{y}}\right)}^2}$$

(8)

$$\mathrm{CCC}=\frac{2\frac{1}{n}{\sum }_{i=1}^n\left(y_i-\overline{\mathrm{y}}\right)\left({\hat{y}}_i-\widehat{y}\right)}{\frac{1}{n}{\sum }_{I=1}^n{\left(y_i-\overline{\mathrm{y}}\right)}^2+\frac{1}{n}{\sum }_{i=1}^n{\left({\widehat{y}}_i-\widehat{y}\right)}^2+{\left(\overline{\mathrm{y}}-\widehat{y}\right)}^2}$$

(9)

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{y}}_{\mathrm{i}}-\widehat{{\mathrm{y}}_{\mathrm{i}}}\right)}^2}{\mathrm{n}}}$$

(10)

$$\mathrm{RPIQ}=\frac{\mathrm{IQ}}{\mathrm{RMSE}}$$

(11)

where ${\mathrm{y}}_{\mathrm{i}}$ is the measured SFI and $\widehat{{\mathrm{y}}_{\mathrm{i}}}$ is the estimated SFI of the ith sample point, n represents the number of samples, and $\overline{y}$ and $\widehat{y}$ denote the average value of observations and estimations, respectively. IQ signifies the interquartile range (IQ = Q3 − Q1) of the observed values. Q1 and Q3 denote the first and third quartile, respectively.

3. Results

3.1. Optimal Crop Spectral Variables for Estimating SFI

The spatial distribution of the crop spectral variables in the study area was obtained using the Google Earth Engine platform and the calculation formula (Table 3). The results are shown in Figure 6.

Numerous experiments showed the prediction error in the XGBoost algorithm stabilized with the shrinkage feature weight (eta), the maximum depth (max_depth), and the number of iterations (nround) being 0.4, 10, and 150, separately. The results showed that 6 crop spectral variables (inverted red-edge chlorophyll index (IRECI), chlorophyll vegetation index (CVI), normalized green-red difference index (NGRDI), red-edge position (REP), triangular greenness index (TGI), and optimized soil-adjusted vegetation index (OSAVI)) derived optimal initial results (Figure 7). Subsequently, the variance inflation factor (VIF) was employed for eliminating the collinearity between these characteristic variables with the screening criteria VIF < 10 [47]. Finally, five crop spectral variables (IRECI, CVI, NGRDI, REP, and TGI) were determined as the optimal crop spectral variables for estimating the SFI.

3.2. Model Construction and Accuracy Evaluation

The MLR and BPNN models were adopted for determining the association of optimal crop spectral variables (IRECI, CVI, NGRDI, REP, and TGI) and the SFI using 90 training samples. For the BPNN algorithm, referred to relevant literature [48,49,50] and through numerous experiments, the number of neuron nodes of the hidden layer was eventually set to 11, the number of iterations was set to 5000, and the learning rate and learning objective were set to 0.01. Figure 8 indicates that the BPNN model provides more accurate estimates than the MLR because the values in the scatter plot approach the 1:1 line.

Table 4. The accuracy assessment results for the soil fertility index.

Model	Data Set	R2	RMSE	CCC	RPIQ
MLR	training	0.03	0.26	0.17	0.76
	validation	0.02	0.28	0.02	0.72
BPNN	training	0.84	0.06	0.92	3.60
	validation	0.66	0.17	0.81	1.16

presents the accuracy assessment findings of SFI based on 30 validation samples. The RMSE values of the BPNN are smaller than those of the MLR, whereas the R² values are larger, indicating that the BPNN model is optimal. Thus, the BPNN model is used to estimate the SFI based on the five spectral variables selected by XGBoost.

3.3. Soil Fertility Index Map

The map of the SFI in the Conghua District obtained from the BPNN model is presented in Figure 9. The SFI value is generally concentrated within 0.20–0.60. The soil fertility is lower in the west but higher in the northeast of the study area.

Figure 10 shows the measured and estimated SFI for 30 sample plots (red dots in Figure 2). The accuracy metrics (R² of 0.62, RMSE of 0.09) show that the proposed model provides potential for mapping the SFI in the Conghua District.

4. Discussion

This study determined the quantitative relationship between crop spectral variables and the SFI. The results indicated that the proposed method has great potential to evaluate soil fertility using crop spectral variables.

4.1. Comparison with Other Similar Studies

Previous studies [6,7,8,9] have focused primarily on descriptions of the relationship between soil spectral indicators and soil fertility. Due to the complexity of soil components, the response spectrum of soil fertility may be disturbed, leading to difficulties in identification. In addition, regionally, in southern China there is more vegetation cover and fewer bare soil areas, which also makes it difficult to use soil spectra for soil fertility monitoring. Thus, some scholars used vegetation spectra to assess soil fertility. However, these studies [10,11,12,51,52] lacked quantitative information on soil fertility and their reliability could not be verified by ground-truth data. In this study, we selected five soil fertility indicators (pH, SOM, TN, AP, and AK) based on previous studies [5,53] and calculated the SFI using a fuzzy approach. The relationship model between crop spectral variables and SFI was established, and quantitative evaluation results of soil fertility were obtained. This quantitative method for SFI estimation provided reasonable accuracy at the sample point level (R² = 0.66) and the regional scale (R² = 0.62). The BPNN algorithm had higher estimation accuracy than MLR, suggesting the marked nonlinear relationship between crop spectral variables and the SFI.

Current studies that used crop spectral indices to estimate soil fertility [10,11,12,51,52] focused primarily on indices in the vis-NIR spectral range and did not consider other spectral indices (e.g., red-edge bands). References [54,55] showed that red-edge bands provided abundant information not contained in the red, green, and short-wave infrared bands. The bands can be utilized for identifying and monitoring the chlorophyll content, phenological growth status, health status of vegetation, and heavy metal pollution and can also reflect soil fertility [55]. Red-edge indices (e.g., IRECI, NDREI, and REP) were evaluated in this study to evaluate their ability to estimate the SFI. The XGBoost algorithm selected five crop spectral variables (IRECI, CVI, NGRDI, REP, and TGI) as the optimal variables for estimating the SFI. The results showed that five crop spectral variables (IRECI, CVI, NGRDI, REP, and TGI) could explain 66% of the variance in soil fertility.

4.2. Prospects for Future Studies

This study considered only paddy fields. Future research may consider various arable land types (such as dry land and irrigated land) to improve the generalizability of the model. In addition, we evaluated soil fertility using only spectral vegetation indices and did not consider other crop growth indicators (such as gross primary productivity, net primary productivity, or leaf area index). Additional crop growth indicators should be introduced to future research for improving SFI estimation accuracy.

Notably, the BPNN model was applied to SFI estimation, while the large uncertain-ties of weights and threshold may potentially affect the accuracy of the model. In future research, more parameter optimization algorithms (e.g., whale optimization, particle swarm optimization) should be introduced to promote efficiency and stability in the BPNN model. Finally, some mixed pixels existed in the images, although the spatial resolution was 10 m. It is unclear as to whether the association of spectral variables with SFI holds true for mixed pixels. Thus, further research is required to support our conclusions.

5. Conclusions

The research here proposes a new approach to the quantitative estimation of soil fertility using crop spectral variables during the crop growth period. The method combines the XGBoost algorithm for variable selection with the BPNN algorithm for soil fertility estimation. The five optimal crop spectral variables (IRECI, CVI, NGRDI, REP, and TGI) were screened, which was the first time soil fertility using red-edge indices was assessed. Based on the five optimal crop spectral variables, BPNN algorithm was used to construct the model to realize the quantitative estimation of soil fertility. The research result showed that the proposed method is reliable with the R² of 0.62 and RMSE of 0.09 at the regional scale. To the best of our knowledge, this research is the first to provide an efficient solution to quantify soil fertility based on crop spectral variables.