Method for inverting Cu content of agricultural land soil by using multispectral data

Abstract

The invention discloses a method for inverting the Cu content of agricultural land soil by using multispectral data. The method comprises the following steps: acquiring Sentinel-2 satellite Level-2A remote sensing image data; the method comprises the following steps: constructing a decision tree by using a soil pH value, humidity and organic matter content, and dividing samples into different range groups through the decision tree; calculating the importance of the characteristic variables by adopting a continuous projection algorithm SPA and a Boruta algorithm, and finding the most suitable spectral index form; a gradient boosting tree XGBoost algorithm is adopted as an inversion model of the Cu content in the soil, and inversion models of a pH interval, a humidity interval and the organic matter content are established in combination with an Optuna framework; and classifying through the decision tree, and inputting the spectral index of the soil to be detected into the inversion model with good prediction performance to obtain a prediction value. According to the method, the accuracy of inversion of the content of the heavy metal Cu in the soil is improved by means of a comprehensive resampling strategy, a multi-dimensional spectrum transformation form, classified data, classified modeling and the like.

Description

Technical Field

[0001] The present invention relates to the technical field of soil inversion, and in particular to a method for inverting the Cu content of agricultural soil using multispectral data. Background Art

[0002] Applying satellite remote sensing data to extract spectral information from remote sensing image data for inversion of soil heavy metal content is the forefront of monitoring and early warning of soil pollutant content in agricultural land. However, the physical and chemical properties of soil, such as pH, humidity, and organic matter, will affect the spectral reflectance. Differences in soil pH, humidity, and organic matter content will cause different spectral response characteristics. Existing studies have shown that soil pH, humidity, and organic matter have a strong correlation with soil spectral reflectance. Therefore, when using these spectral reflectance data for inversion, in order to improve the accuracy of data analysis, the influence of these factors needs to be considered. Although Sentinel-2 satellite multispectral image data has the advantages of free, easy to obtain, and multiple time series, it is difficult to meet the data requirements of the high-precision quantitative inversion model of soil heavy metal content due to its low spectral resolution and relatively weak spectral response of soil heavy metals. The inversion accuracy of existing inversion methods is generally low. Summary of the invention

[0003] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a method for inverting the Cu content of agricultural soil using multispectral data, which can improve the inversion precision and accuracy.

[0004] The present invention is achieved through the following technical solutions:

[0005] A method for inverting the Cu content of agricultural soil using multispectral data comprises the following steps:

[0006] S1. Data collection and processing: Acquire Level-2A remote sensing image data processed by Sen2Cor from the Sentinel-2 satellite in the study area, process and download them through the Google Earth Engine (GEE) cloud platform, and obtain surface reflectance spectral data after filtering and masking, atmospheric correction, geometric correction and orthorectification correction; conduct field sampling and analysis of the soil in the study area, test the performance indicators of the sampled samples including Cu content, pH value, humidity and organic matter content, and form calibration data sets and validation data sets;

[0007] S2. Construct a decision tree: take soil pH as a factor and pH decision tree correlation coefficient R phAs the evaluation index, two pH value optimal nodes pH1 and pH2 are obtained, thereby obtaining three pH intervals [pHmin, pH1), [pH1, pH2) and [pH2, pHmax]; taking these three pH intervals as classification targets, a pH value decision tree is constructed, and the pH value decision tree is used to divide the samples into corresponding pH intervals according to the pH value; wherein pHmin represents the minimum pH value observed in the soil sample, pHmax represents the maximum pH value observed in the soil sample, pHmin<pH1<pH2<pHmax,

[0008] Taking soil moisture WETNESS as the factor and humidity decision tree correlation coefficient R WETNESS As the evaluation index, two optimal humidity nodes WETNESS1 and WETNESS2 are obtained, thereby obtaining three humidity intervals [WETNESSmin, WETNESS1), [WETNESS1, WETNESS2), and [WETNESS2, WETNESSmax]. Taking these three humidity intervals as classification targets, a humidity decision tree is constructed, and the humidity decision tree is used to divide the samples into corresponding humidity intervals according to humidity. Among them, WETNESSmin represents the minimum humidity value observed in the soil sample, WETNESSmax represents the maximum humidity value observed in the soil sample, WETNESSmin<WETNESS1<WETNESS2<WETNESSmax,

[0009] Taking organic matter content SOM as a factor and the organic matter decision tree correlation coefficient R SOM As the evaluation index, two optimal nodes of organic matter content, SOM1 and SOM2, were obtained, thereby obtaining three organic matter content intervals [SOMmin, SOM1), [SOM1, SOM2), and [SOM2, SOMmax]. Taking these three organic matter content intervals as classification targets, an organic matter decision tree was constructed, which was used to divide the samples into corresponding organic matter content intervals according to the organic matter content. Among them, SOMmin represents the minimum moisture value observed in the soil sample, SOMmax represents the maximum moisture value observed in the soil sample, SOMmin<SOM1<SOM2<SOMmax,

[0010]

[0011] Where Bi represents the spectral reflectance of the i-th band of the Sentinel-2 satellite, where i = 1, 2, 3, 4, 5, 6, 7, 8, 8A, 9, 11, 12, corresponding to the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 bands of the Sentinel-2 satellite; SpH、 S WET and S SOM represents the soil Cu content at the specified soil pH value, humidity and organic matter content, mg / kg; Rn represents the correlation coefficient R between the spectral reflectance of the i-th band in the n-th branch of the Sentinel-2 satellite and the soil Cu content at the specified pH value, humidity or organic matter content, where the n-th branch refers to one of the three intervals of pH value, humidity or organic matter content, n = 1, 2, 3, and the correlation coefficient X i and Y i are the i-th observation values ​​of X and Y respectively, and are the average values ​​of X and Y, respectively, where X represents the spectral reflectance, Y represents the soil Cu content at a specified soil pH, moisture or organic matter content, and k is the number of samples;

[0012] S3. Screening spectral index: By constructing two-dimensional spectral index form and three-dimensional spectral index form for the original bands of Sentinel-2, using the original bands of Sentinel-2, two-dimensional spectral index and three-dimensional spectral index, M feature variables are obtained; then the continuous projection algorithm SPA is used to screen out the top N most important bands; then the Boruta algorithm is used to further refine these N bands. The Boruta algorithm is a feature selection method based on random forests. The Boruta algorithm introduces shadow features, which are generated by disrupting the order of the original features (N bands), adding noise and randomization (these shadow features are used to simulate randomly selected features and compared with the original features). By calculating and comparing the importance scores of each original feature and its corresponding shadow feature, the truly important key features are identified; the original bands of Sentinel-2, the two-dimensional spectral index and the three-dimensional spectral index find the most suitable spectral index form through the above SPA and Boruta feature screening algorithms;

[0013] S4. Construction and training of inversion model: Taking the key features as independent variables, the gradient boosting tree XGBoost algorithm is used to construct an inversion model for Cu content in soil, and three pH intervals, three moisture intervals, and three organic matter content inversion models are established respectively; the pH value decision tree, moisture decision tree, and organic matter decision tree are used to classify the data in the data set; the calibration data set is used to train the regression model, and the validation data set is used to make predictions and calculate the root mean square error RMSE to evaluate the prediction accuracy; the optimal parameter configuration of each model is determined according to the value of the loss function;

[0014] S5. Prediction of soil Cu content: The soil to be tested is classified according to its pH value, humidity, and organic matter content through the pH value decision tree, humidity decision tree, and organic matter decision tree. The inversion model with the best model performance is selected from the inversion models corresponding to the pH value interval, the inversion model corresponding to the humidity interval, and the inversion model corresponding to the organic matter content interval, and the soil spectral index screened in the corresponding branch is input into the inversion model, and the inversion model outputs the Cu content value in the soil. That is, the influencing factor corresponding to the model with the best inversion result is the influencing factor with the least influence on the spectral reflectance. Selecting this model for prediction can improve the accuracy of the prediction. The factor interval corresponding to the model is also the factor interval to which the sample belongs. The spectral index of the soil to be tested input is the most suitable spectral index form of the inversion model found in step S3.

[0015] Furthermore, in step S4, the Optuna framework is used in model training to evaluate the performance of each model and return the corresponding evaluation index by defining the objective function, and the search space of the hyperparameters is defined in the objective function.

[0016] Furthermore, the two-dimensional spectral index includes a difference index DI, a ratio index RI and a normalized difference index NDIs, DI(Ri, Rj))=Ri-Rj,

[0017] The three-dimensional spectral index includes TBI1, TBI2, TBI3, TBI4, and TBI5, which represent the spectral index of the first, second, third, fourth, and fifth transformation forms, respectively. Where Ri, Rj and Rk are the reflectances of any three bands selected from all the bands of the Sentinel-2 satellite.

[0018] Furthermore, in step S1, the method for processing and downloading through the Google Earth Engine (GEE) cloud platform includes: defining a function for masking the cloud layer and the cumulus cloud layer in the Sentinel-2 remote sensing image in Google Earth Engine (GEE), the function using the QA60 band containing the quality assessment information of the image, screening out the image set between a period of time, excluding the images with a cloud cover exceeding 5%, obtaining the image set after filtering and masking processing and taking the average value thereof;

[0019] Furthermore, in step S1, the method for processing and downloading through the Google Earth Engine (GEE) cloud platform also includes: using cubic convolution interpolation in the Google Earth Engine (GEE) cloud platform to resample the B1 band and B5 to B12 bands in the Level-2A data to a resolution of 10m.

[0020] Furthermore, in step S2, when constructing the humidity decision tree, the surface reflectance spectral data of the Sentinel-2 satellite is transformed by tasseled cap, and the multi-spectral band data is reduced to three main components: brightness, greenness and humidity, so as to extract the humidity index T CW (Tasseled cap Wetness) formula to build a decision tree;

[0021] Humidity index T CW The formula is:

[0022] T CW =0.0649B1+0.1363B2+0.2802B3+0.3072B4+0.5288B5+0.1379B6-0.0001B7-0.0807B8-0.1389B8A

[0023] -0.0302B9-0.4064B11-0.5602B12

[0024] Among them, B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, and B12 are the reflectances of the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, and B12 bands respectively.

[0025] Furthermore, in step S2, the method for obtaining the pH value optimal nodes pH1 and pH2 is as follows: select two pH values, pHa and pHb, respectively, pHmin<pHa<pHb<pHmax, form three pH intervals [pHmin, pHa), [pHa, pHb) and [pHb, pHmax], assign samples to corresponding intervals according to their pH values, repeat the process iteratively, and ensure that the number of samples in each interval is not less than the set value (such as 70), and the evaluation index R pH When the maximum value is reached, the corresponding decision tree is the best, and the corresponding pHa is pH1, and pHb is pH2.

[0026] Furthermore, in step S2, the method for obtaining the optimal humidity nodes WETNESS1 and WETNESS2 is as follows: select two humidity values, WETNESSa and WETNESSb, WETNESSmin<WETNESSa<WETNESSb<WETNESSmax, to form three humidity intervals [WETNESSmin, WETNESSa), [WETNESSa, WETNESSb), [WETNESSb, WETNESSmax], assign samples to corresponding intervals according to their humidity, repeat the process iteratively, and ensure that the number of samples in each interval is not less than the set value (such as 70), and the evaluation index R WETNESS When the maximum value is reached, the corresponding decision tree is the best. At this time, the corresponding WETNESSa is WETNESS1, and WETNESSb is WETNESS2.

[0027] Furthermore, in step S4, the hyperparameters that need to be adjusted in the XGBoost algorithm include: num_boost_round, eta, max_depth, subsample, colsample_bytree, gamma, min_child_weight, rel_lambda, and alpha. The value ranges of each hyperparameter are set as follows: num_boost_round is [500,2000], eta is [0.01,0.3], max_depth is [3,10], subsample is [0.5,1], colsample_bytree is [0.5,1], gamma is [0,1], min_child_weight is [1,10], rel_lambda is [0,1], and alpha is [0,1].

[0028] Furthermore, the determination coefficient R 2 , RMSE and RPD were used to evaluate the accuracy of each soil Cu content regression model. Where n represents the number of soil samples, y i and represent the measured value and model predicted value of the i-th soil sample in the validation data set, respectively. Represents the measured average.

[0029] The present invention adopts a comprehensive resampling strategy to uniformly resample all bands to a spatial resolution level of 10m, thereby ensuring that more accurate spatial information can be obtained when inverting soil heavy metals, solving the spatial limitation of multispectral in small-scale soil research, and helping to improve the accuracy of inversion results, especially in specific areas or farmlands that require fine inversion; by using multispectral original spectral data and spectral transformation forms, such as two-dimensional spectral index and three-dimensional spectral index, potential band associations in multispectral data are mined, and the signal response of the spectrum to soil heavy metal Cu is enhanced, which helps to improve the sensitivity of the soil Cu content inversion model to spectral information and reduce the risk of missing important information; by establishing a decision tree with soil pH value, humidity (WETNESS) and organic matter content (SOM) as branch standards, effective classification of soil samples is achieved, and further classification modeling is performed, thereby reducing the influence of soil physical and chemical property indicators on soil Cu content inversion. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 The figure is a flow chart of a method according to an embodiment of the present invention.

[0031] Figure 2 Schematic diagram of pH interval division in an embodiment of the present invention.

[0032] Figure 3 Schematic diagram of humidity interval division in an embodiment of the present invention.

[0033] Figure 4 This is a schematic diagram of determining pH branch nodes when constructing a pH value decision tree in an embodiment of the present invention.

[0034] Figure 5 This is a schematic diagram of determining humidity branch nodes when constructing a humidity decision tree in an embodiment of the present invention. DETAILED DESCRIPTION

[0035] A method for inverting Cu content in agricultural soil using multispectral data, such as Figure 1 , including the following steps:

[0036] S1. Data collection and processing: Level-2A remote sensing image data processed by Sentinel-2 satellite and Sen2Cor in the study area were obtained, and processed and downloaded through the Google Earth Engine (GEE) cloud platform to obtain surface reflectance spectral data after filtering and masking processing, atmospheric correction, geometric correction and orthorectification.

[0037] Sentinel-2 can provide high-resolution optical images. Level-1C is the raw data containing satellite observation information of the earth's surface. Sentinel-2Level-1C can be processed into Level-2A data using Sen2Cor. Level-2A data is atmospherically corrected, geometrically corrected, and orthorectified, and can be widely used in vegetation monitoring, soil analysis, water resources management and other fields.

[0038] The present invention processes and downloads Level-2A data through the Google Earth Engine (GEE) cloud platform, including: (1) defining a function for masking clouds and cumulus clouds in Sentinel-2 remote sensing images, the function uses the QA60 band containing image quality assessment information to screen out image sets between a period of time, excludes images with cloud cover exceeding 5%, obtains an image set after filtering and masking, and takes an average value; (2) in order to meet research needs, the B1 and B5 to B12 bands in the Level-2A data are resampled to a resolution of 10m using cubic convolution interpolation in the Google Earth Engine (GEE) cloud platform, and the resolution of the remaining bands can be maintained at the original 10m. These processing processes help to ensure the consistency and accuracy of the data, so that the characteristics of soil heavy metals can be explored more accurately in subsequent analysis.

[0039] At the same time, the soil in the corresponding study area is sampled and analyzed on the spot, the sampling location and sampling time of these samples are recorded, and the performance indicators of these samples are measured, including Cu content, pH value, humidity and organic matter content, etc., to form a calibration data set and a verification data set. This embodiment takes the study of a certain area in South China as an example. 242 soil surface samples are collected in this area in South China. The content of Cu in the soil samples is determined by reference to the standard method "Determination of Copper, Zinc, Lead, Nickel and Chromium in Soil and Sediment Flame Atomic Absorption Spectrophotometry" (HJ491-2019), wherein the minimum value of soil Cu content is 3mg / kg, the maximum value is 271mg / kg, the standard deviation is 29.7mg / kg, and the coefficient of variation is 87.6%. Among them, the coefficient of variation is high, indicating that there is a large degree of discreteness in the soil Cu content in the sample. At the same time, the average value of soil Cu content in the study area is 33.9mg / kg, which is significantly higher than the background value of surface soil in Guangdong Province and China, and is 3 times the background value of surface soil in Guangdong Province.

[0040] S2. Construct a decision tree: Construct a decision tree based on soil pH, humidity, and organic matter content, and divide the samples into range groups with different pH intervals, humidity intervals, and organic matter content intervals through the decision tree. Classify the samples according to the pH, humidity, and organic matter content values, and construct a heavy metal Cu content inversion model based on the classification, so that the samples can be divided into the corresponding range group inversion model according to their properties, reducing the impact of different soil pH values, humidity, and organic matter content on the final prediction results and improving the inversion accuracy.

[0041] Specifically, the soil samples are divided into different intervals according to pH value, humidity and organic matter content. The interval data of pH value, humidity and organic matter content can be determined according to actual conditions, such as being divided into two intervals, three intervals, four intervals, etc., mainly based on the number of samples, the sensitive relationship between Cu content and indicators, and the indicator range of the study area. If different indicator values ​​have a greater impact on Cu content, the number of intervals can be increased. If the indicator range of the study area is wide (such as a large humidity range distribution), the number of intervals of the corresponding indicator can be appropriately increased. As one implementation method, each indicator in this embodiment is divided into three intervals. Specifically, for pH value, all soil samples are divided into three data sets, such as Figure 2 , the samples are divided into three pH intervals: [pHmin, pH1), [pH1, pH2], and [pH2, pHmax]. For moisture, all soil samples are divided into three data sets, such as Figure 3 , including three humidity intervals: [WETNESSmin, WETNESS1), [WETNESS1, WETNESS2), and [WETNESS2, WETNESSmax]. For organic matter content, all soil samples were divided into three data sets, including three organic matter intervals: [SOMmin, SOM1), [SOM1, SOM2), and [SOM2, SOMmax].

[0042] In order to achieve a reasonable division of each index interval, it is necessary to find the appropriate corresponding pH1 and pH2, WETNESS1 and WETNESS2, SOM1 and SOM2, that is, to find the optimal value of the interval node, and then obtain the optimal division interval based on these optimal node values, so that the decision tree can allocate samples to the corresponding intervals according to these nodes and facilitate the establishment of subsequent models. The calculation of these optimal node values ​​and the establishment method of the corresponding decision tree are described in detail below.

[0043] Taking soil pH as an example, firstly, soil pH is taken as a factor and the pH decision tree correlation coefficient R is used as pHAs the evaluation index, two pH optimal nodes pH1 and pH2 are obtained, thereby obtaining three pH intervals [pHmin, pH1), [pH1, pH2) and [pH2, pHmax]. Then, these three pH intervals are used as classification targets to construct a pH decision tree. The pH decision tree is used to divide samples into corresponding pH intervals according to pH values.

[0044] Wherein, pHmin represents the minimum pH value observed in the soil sample, pHmax represents the maximum pH value observed in the soil sample, and pHmin<pH1<pH2<pHmax ， Bi represents the spectral reflectance of the i-th band of the Sentinel-2 satellite, where i = 1, 2, 3, 4, 5, 6, 7, 8, 8A, 9, 11, 12, corresponding to the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 bands of the Sentinel-2 satellite; S pH represents the soil Cu content at a specified soil pH value, mg / kg; Rn represents the correlation coefficient R between the spectral reflectance of the i-th band in the n-th branch of the Sentinel-2 satellite and the soil Cu content at a specified pH value, where the n-th branch refers to one of the three pH value intervals, n = 1, 2, 3, and the correlation coefficient X i and Y i are the i-th observation values ​​of X and Y respectively, and are the average values ​​of X and Y, respectively, k is the number of samples, where X and Y represent the spectral reflectance and the soil Cu content at a specified pH value, respectively.

[0045] Specifically, the method for obtaining the optimal pH nodes pH1 and pH2 is as follows: select two pH values, pH1 and pH2, respectively, pHmin<pH1<pH2<pHmax, as shown in the figure, to form three pH intervals [pHmin, pH1), [pH1, pH2) and [pH2, pHmax]. The step sizes of pH1 and pH2 are determined by the number of samples and the distribution of pH values. The samples are assigned to different intervals according to their pH values. The process is repeated iteratively, and it is ensured that the number of samples in each interval is not less than the set value 70. The evaluation index R ph When the maximum value is reached, the corresponding decision tree is the best, and the optimal nodes pH1 and pH2 are obtained.

[0046] In this embodiment, according to the minimum pH value pHmin of the soil sample being 3.63, the maximum pH value pHmax being 9.25, the step size of pH1 being 0.1, and the step size of pH2 being 0.05, the samples are actually divided, and the threshold combination of the pH value decision tree is shown in Table 1.

[0047] Table 1 Threshold combinations of pH value decision tree

[0048]

[0049]

[0050] Figure 4 The R values ​​of different bands and different pH threshold combinations in Table 1 are shown. pH value, Figure 3 It shows that when pH1 and pH2 are 5.85 and 7.75 respectively, the R pH All of them can reach the maximum value, that is, the correlation coefficient between the spectral reflectance of all bands and the soil Cu content at the specified soil pH value can reach the maximum value. Therefore, the soil samples are divided into three pH intervals, namely [3.63, 5.85), [5.85, 7.75) and [7.75, 9.25].

[0051] Taking soil moisture WETNESS as an example, soil moisture WETNESS is used as a factor and the humidity decision tree correlation coefficient R WETNESS As the evaluation index, two humidity optimal nodes WETNESS1 and WETNESS2 are obtained, thereby obtaining three humidity intervals [WETNESSmin, WETNESS1), [WETNESS1, WETNESS2), and [WETNESS2, WETNESSmax]. Then, these three humidity intervals are used as classification targets to construct a humidity decision tree. The humidity decision tree is used to divide the samples into corresponding humidity intervals according to the humidity.

[0052] Among them, WETNESSmin represents the minimum moisture value observed in the soil sample, WETNESSmax represents the maximum moisture value observed in the soil sample, WETNESSmin<WETNESS1<WETNESS2<WETNESSmax, Bi represents the spectral reflectance of the i-th band of the Sentinel-2 satellite, where i = 1, 2, 3, 4, 5, 6, 7, 8, 8A, 9, 11, 12, corresponding to the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 bands of the Sentinel-2 satellite; S WETNESSrepresents the soil Cu content at a specified soil moisture level, mg / kg; Rn represents the correlation coefficient R between the spectral reflectance of the i-th band in the n-th branch of the Sentinel-2 satellite and the soil Cu content at a specified humidity, where the n-th branch refers to one of the three intervals of humidity, n = 1, 2, 3, and the correlation coefficient X i and Y i are the i-th observation values ​​of X and Y respectively, and are the average values ​​of X and Y, respectively, k is the number of samples, where X and Y represent the spectral reflectance and the soil Cu content at a specified humidity, respectively.

[0053] Considering the inconsistency between the time when the remote sensing images were acquired and the moisture content of the samples collected on site, this embodiment performs a tasseled cap transformation on the surface reflectance spectral data of the Sentinel-2 satellite. This transformation reduces the dimensionality of the multispectral band data into three main components: brightness, greenness, and humidity, in order to extract the humidity index T CW The decision tree is constructed by using the Tasseled cap Wetness (Tasseled cap Wetness) formula to reduce the potential impact of soil moisture on the inversion results.

[0054] Humidity index T CW The formula is:

[0055] T CW =0.0649B1+0.1363B2+0.2802B3+0.3072B4+0.5288B5+0.1379B6-0.0001B7-0.0807B8-0.1389B8A

[0056] -0.0302B9-0.4064B11-0.5602B12

[0057] Among them, B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, and B12 are the reflectivities of B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, and B12 bands of Sentinel-2 satellite respectively.

[0058] Specifically, the method for obtaining the optimal humidity nodes WETNESS1 and WETNESS2 is as follows: select two humidity values, WETNESS1 and WETNESS2, WETNESSmin<WETNESS1<WETNESS2<WETNESSmax, and form three pH intervals [WETNESSmin, WETNESS1), [WETNESS1, WETNESS2), and [WETNESS2, WETNESSmax]. The step sizes of WETNESS1 and WETNESS2 are determined by the number of samples and humidity distribution. The samples are assigned to different intervals according to their humidity. The process is repeated iteratively, and it is ensured that the number of samples in each interval is not less than the set value 70. The evaluation index R WETNESS When the maximum value is reached, the corresponding decision tree is the best, and the nodes corresponding to this decision tree are the humidity optimal nodes WETNESS1 and WETNESS2.

[0059] In this embodiment, according to the minimum moisture value WETNESSmin of the soil sample being -0.1115, the maximum moisture value WETNESSmax being 9.25, the step size of WETNESSa being 0.0015, and the step size of WETNESSb being 0.001, the samples are actually divided, and the threshold combination of the moisture decision tree is shown in Table 2.

[0060] Table 2 Threshold combinations of humidity decision tree

[0061]

[0062] Figure 5 The R values ​​of different humidity threshold combinations in different bands are shown in Table 2. WETNESS value, Figure 4 Display, such as Figure 4 As shown in the figure, when WETNESS1 and WETNESS2 are -0.0785 and -0.062 respectively, except for the B5 and B11 bands, R WETNESS The maximum value cannot be reached, and the correlation coefficients between the spectral reflectance of other bands and the soil Cu content under specified soil moisture can all reach the maximum value. Therefore, the soil samples are divided into three moisture intervals, namely [-0.1115, -0.0785), [-0.0785, -0.062) and [-0.062, 0.0296].

[0063] Taking organic matter content SOM as a factor and the organic matter decision tree correlation coefficient R SOMAs the evaluation index, two optimal nodes of organic matter content, SOM1 and SOM2, were obtained, thereby obtaining three organic matter content intervals [SOMmin, SOM1), [SOM1, SOM2), and [SOM2, SOMmax]. Then, these three organic matter content intervals were used as classification targets to construct an organic matter decision tree, which was used to divide samples according to their organic matter content.

[0064] Among them, SOMmin represents the minimum moisture value observed in the soil sample, SOMmax represents the maximum moisture value observed in the soil sample, SOMmin<SOM1<SOM2<SOMmax, Bi represents the spectral reflectance of the i-th band of the Sentinel-2 satellite, where i = 1, 2, 3, 4, 5, 6, 7, 8, 8A, 9, 11, 12, corresponding to the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 bands of the Sentinel-2 satellite; S SOM represents the soil Cu content at a specified soil organic matter content, mg / kg; Rn represents the correlation coefficient R between the spectral reflectance of the i-th band in the n-th branch of the Sentinel-2 satellite and the soil Cu content at a specified organic matter content, where the n-th branch refers to one of the three intervals of pH value, n = 1, 2, 3, and the correlation coefficient X i and Y i are the i-th observation values ​​of X and Y respectively, and are the average values ​​of X and Y, respectively, k is the number of samples, where X and Y represent the spectral reflectance and the soil Cu content at a specified organic matter content, respectively.

[0065] The method for obtaining the optimal nodes SOM1 and SOM2 for organic matter content refers to the method for the aforementioned pH value and humidity value.

[0066] S3. Screening spectral index: By constructing two-dimensional spectral index form and three-dimensional spectral index form for the original bands of Sentinel-2, M feature variables are obtained using the original bands of Sentinel-2, two-dimensional spectral index and three-dimensional spectral index; then the continuous projection algorithm SPA is used to screen out the top N most important bands, thereby reducing the number of features and reducing the collinearity between them; then the Boruta algorithm is used to further refine these N bands. The Boruta algorithm is a feature selection method based on random forests. The Boruta algorithm introduces shadow features, which are generated by disrupting the order of the original features (N bands), adding noise and randomization (these shadow features are used to simulate randomly selected features and compared with the original features). By calculating and comparing the importance scores of each original feature and its corresponding shadow feature, the truly important key features are identified; the most suitable spectral index form is found by iterating the feature screening algorithm of the above SPA+Boruta on the original bands of Sentinel-2, the two-dimensional spectral index and the three-dimensional spectral index. This method not only comprehensively considers all the features that affect the prediction target, but also effectively identifies the important features in the variable selection process instead of just finding the optimal subset.

[0067] The Boruta algorithm is a feature screening method based on random forests. It generates shadow features to evaluate the importance of each feature, thereby helping to identify which features are truly important. The core of the Boruta algorithm is to determine which features are related to the dependent variable by comparing the importance of the original features and the randomly generated features (shadow features). Shadow features are generated by randomly shuffling the values ​​of the original features. The shadow feature generation process includes randomly shuffling each original feature and appending the shadow feature to the original feature data frame to form a new data frame. By training a model, the importance of the original features and shadow features are compared. If the importance of a feature is significantly higher than its corresponding shadow feature, then the feature is considered important. Conversely, if the importance of a feature is significantly lower than its corresponding shadow feature, then the feature is marked as irrelevant.

[0068] Studies have shown that spectral indexes respond to soil heavy metal Cu content. The present invention selects commonly used two-dimensional spectral index forms, including difference index (DI), ratio index (RI) and normalized difference index (NDIs), as well as common three-dimensional spectral index forms, and finds the most suitable spectral index form by iterating all possible spectral indexes (including original band data). Table 3 details all spectral index forms and calculation formulas.

[0069] Table 3 Two-dimensional and three-dimensional spectral indices

[0070]

[0071] Among them, the three-dimensional spectral indices TBI1, TBI2, TBI3, TBI4, and TBI5 represent the spectral indices of the first, second, third, fourth, and fifth transformation forms, respectively, and Ri, Rj, and Rk are the reflectances of any three bands selected from all the bands of the Sentinel-2 satellite.

[0072] This embodiment uses Sentinel-2 original bands and the above-mentioned two-dimensional and three-dimensional spectral indices to obtain 7008 characteristic variables. In the variable selection process, the continuous projection algorithm (SPA) method is used to select the top 100 most important bands. Subsequently, the Boruta algorithm is applied to further refine these 100 bands. Taking pH value as an example, the most suitable spectral index for each pH interval is finally obtained in Table 4.

[0073] Table 4 Comparison of feature band selection results under different branch conditions

[0074]

[0075] Table 4 shows the results of the characteristic band selection with and without considering the decision tree. It is worth noting that in all cases, no changes in the two-dimensional spectral index were selected, indicating that the prediction ability of the two-dimensional spectral index is not as good as that of the three-dimensional spectral index under the current conditions. In addition, the selected characteristic bands are mainly concentrated on the three-dimensional spectral indexes TBI2 and TBI4, indicating that the changes in these indices are highly correlated with the inversion of soil heavy metal Cu and provide more information for the inversion of soil heavy metal Cu, and the forms of TBI2 and TBI4 corresponding to each branch interval are also different.

[0076] S4. Construction and training of inversion model: Taking the key features as independent variables, the gradient boosting tree XGBoost algorithm is used to construct an inversion model for Cu content in soil, and three pH intervals, three humidity intervals, and three organic matter content inversion models are established through decision trees.

[0077] The Gradient Boosting Tree (XGBoost) algorithm is a machine learning algorithm based on the Gradient Boosting Decision Tree (GBDT) with properties such as efficiency, flexibility, and portability. Based on the upgrade ensemble strategy, XGBoost introduces a series of optimization measures to improve the performance and efficiency of the model. The core idea of ​​XGBoost is to build a series of weak learners and gradually correct previous prediction errors to approach the true target function. In each iteration, XGBoost calculates the residuals of the current model and trains a new decision tree based on these residuals. The new decision tree focuses on optimizing the samples that the previous model failed to predict correctly, thereby gradually reducing the prediction error.

[0078] The present invention combines the Optuna framework with the XGBoost algorithm, aiming to invert the Cu content concentration in the soil through the following steps: (1) defining an objective function to evaluate the performance of the model and return the corresponding evaluation index, in which the search space of the hyperparameters is clearly defined; (2) when performing a separate experiment, the data in the data set are classified using a pH decision tree, a moisture decision tree, and an organic matter decision tree, a calibration data set is used to train a regression model, and then a validation data set is used for prediction and the root mean square error (RMSE) is calculated to evaluate the prediction accuracy; (3) through multiple experiments, the optimal parameter configuration of the model is determined according to the value of the loss function. In this embodiment, the XGBoost algorithm needs to adjust 9 key parameters, as shown in Table 5, and the other parameters are set to the default values.

[0079] Table 5 XGBoost hyperparameter value range

[0080]

[0081]

[0082] Using R 2 The accuracy of the soil copper content regression model was evaluated by three accuracy indicators: R (coefficient of determination), RMSE (root mean square error) and RPD (residual prediction deviation). 2 The value is between 0 and 1. The closer it is to 1, the better the model fits the data and the higher the model stability. RMSE represents the difference between the model prediction value and the actual observation value, which is used to measure the prediction accuracy of the model. The smaller the RMSE, the smaller the difference between the prediction value and the actual observation value, and the higher the prediction accuracy of the model. RPD is the ratio of prediction deviation to standard deviation, which can evaluate the prediction ability of the model. Between 0 and 1.40, it means that the model has poor prediction ability, between 1.40 and 2.00, it means that the prediction effect is average, and above 2.00, it means that the model has good prediction ability. The calculation method of each indicator is as follows:

[0083]

[0084] Where n represents the number of soil samples, y i and represent the measured value and predicted value of the i-th soil sample in the validation dataset, respectively, and Represents the average Cu content measured.

[0085] S5, prediction of soil Cu content: the soil to be tested is classified according to its pH value, humidity, and organic matter content through pH decision tree, humidity decision tree, and organic matter decision tree. The inversion model with the best model performance is selected from the inversion model corresponding to the pH value interval, the inversion model corresponding to the humidity interval, and the inversion model corresponding to the organic matter content interval, and the soil spectral index screened in the corresponding branch is input into the inversion model, and the inversion model outputs the Cu content value in the soil, that is, the influencing factor corresponding to the model with the best inversion result is the influencing factor with the least influence on the spectral reflectance. Selecting this model for prediction can improve the accuracy of the prediction, and the factor interval corresponding to the model is also the factor interval to which the sample belongs. This method can comprehensively detect the different influencing factors of the sample, so that each test sample can find the best inversion model and obtain the most accurate prediction result. The spectral index of the soil to be tested in the input model is the most suitable spectral index form of the inversion model found in step S3. In this embodiment, the most suitable three-dimensional spectral indexes of the inversion model for each pH value interval are TBI2 and TBI4, but the TBI2 and TBI4 corresponding to each branch interval are not the same. During inversion, TBI2 and TBI4 of the corresponding pH value interval need to be input into the corresponding model according to the pH value corresponding to the soil.

[0086] This embodiment takes the inversion model of different pH value intervals as an example, and compares it with the inversion model without building a decision tree and performing pH value partitioning (full pH range), demonstrating the performance of the model for each pH interval. The evaluation results are shown in Table 6.

[0087] Table 6 Evaluation of soil copper content inversion results at different soil pH values

[0088]

[0089] Table 6 shows the evaluation results of soil copper content inversion without building a decision tree and under different soil pH conditions. When the decision tree was not built, the inversion effect of soil copper content was poor, and the R 2 The coefficient of determination, RMSE, and RPD did not perform well, and basically lacked predictive power. When considering the decision tree, the XGBoost model had a lower R value in the validation set when the soil pH value was between 5.85 and 7.75. 2 The XGBoost model has an excellent performance in this range and can accurately predict the copper content of soil in this pH range. This result also shows that the decision tree can effectively improve the sensitivity of the spectrum to the copper content of soil, thereby significantly improving the inversion accuracy. In addition, when the pH value is lower than 5.85, the R 2 When the pH value is higher than 7.75, R 2The value is 0.5, which indicates that within these pH ranges, the relationship between spectral reflectance and soil copper content is weak or even misleading, which in turn affects the accuracy of soil copper content prediction within these pH ranges. For soils within these pH ranges, the inversion models corresponding to the soil moisture range and the inversion models corresponding to the organic matter content range can be comprehensively compared, and the inversion model with better prediction performance can be selected for copper content prediction.

[0090] The above detailed description is a specific description of a feasible embodiment of the present invention. The embodiment is not intended to limit the patent scope of the present invention. Any equivalent implementation or modification that does not deviate from the present invention should be included in the patent scope of this case.

Claims

1. A method for inverting the Cu content of agricultural soil using multispectral data, characterized in that: The steps include: S1. Data collection and processing: Acquire Level-2A remote sensing image data of the study area processed by Sen2Cor from the Sentinel-2 satellite, and process and download it through the Google Earth Engine cloud platform to obtain surface reflectance spectral data after filtering and masking, atmospheric correction, geometric correction and orthorectification correction; conduct field sampling and analysis of the soil in the study area, test the performance indicators of the sampled samples including Cu content, pH value, humidity and organic matter content, and form calibration data sets and validation data sets; S2. Construct a decision tree: take soil pH as a factor and pH value as a decision tree correlation coefficient R ph As the evaluation index, two pH value optimal nodes pH1 and pH2 are obtained, thereby obtaining three pH intervals [pHmin, pH1), [pH1, pH2) and [pH2, pHmax]; Taking these three pH intervals as classification targets, a pH decision tree was constructed. The pH decision tree was used to divide samples into corresponding pH intervals according to pH values; pHmin represents the minimum pH value observed in soil samples, pHmax represents the maximum pH value observed in soil samples, pHmin<pH1<pH2<pHmax, Taking soil moisture WETNESS as the factor and humidity decision tree correlation coefficient R WETNESS As the evaluation index, two humidity optimal nodes WETNESS1 and WETNESS2 are obtained, thereby obtaining three humidity intervals [WETNESSmin, WETNESS1), [WETNESS1, WETNESS2), and [WETNESS2, WETNESSmax]. Taking these three humidity intervals as classification targets, a humidity decision tree is constructed, and the humidity decision tree is used to divide the samples into corresponding humidity intervals according to humidity; among them, WETNESSmin represents the minimum humidity value observed in the soil sample, WETNESSmax represents the maximum humidity value observed in the soil sample, and WETNESSmin<WETNESS1<WETNESS2<WETNESSmax ， Taking organic matter content SOM as a factor and the organic matter decision tree correlation coefficient R SOM As the evaluation index, two optimal nodes of organic matter content, SOM1 and SOM2, were obtained, thereby obtaining three organic matter content intervals [SOMmin, SOM1), [SOM1, SOM2), and [SOM2, SOMmax]. Taking these three organic matter content intervals as classification targets, an organic matter decision tree was constructed, which was used to divide the samples into corresponding organic matter content intervals according to the organic matter content. Among them, SOMmin represents the minimum moisture value observed in the soil sample, SOMmax represents the maximum moisture value observed in the soil sample, SOMmin<SOM1<SOM2<SOMmax, Where Bi represents the spectral reflectance of the i-th band of the Sentinel-2 satellite, where i = 1, 2, 3, 4, 5, 6, 7, 8, 8A, 9, 11, 12, corresponding to the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, B12 bands of the Sentinel-2 satellite; S pH、 S WET and S SOM represents the soil Cu content at the specified soil pH value, humidity and organic matter content, mg / kg; Rn represents the correlation coefficient R between the spectral reflectance of the i-th band in the n-th branch of the Sentinel-2 satellite and the soil Cu content at the specified pH value, humidity or organic matter content, where the n-th branch refers to one of the three intervals of pH value, humidity or organic matter content, n = 1, 2, 3, and the correlation coefficient X i and Y i are the i-th observation values ​​of X and Y respectively, and are the average values ​​of X and Y, respectively, where X represents the spectral reflectance, Y represents the soil Cu content at a specified soil pH, moisture or organic matter content, and k is the number of samples; S3. Screening spectral index: By constructing a two-dimensional spectral index form and a three-dimensional spectral index form for the original bands of Sentinel-2, using the original bands of Sentinel-2, the two-dimensional spectral index and the three-dimensional spectral index, M feature variables are obtained; then the continuous projection algorithm SPA is used to screen out the first N most important bands; then the Boruta algorithm is used to further refine these N bands. The Boruta algorithm is a feature selection method based on random forests. It introduces shadow features, calculates and compares the importance scores of each original feature and its corresponding shadow feature, and identifies the truly important key features; through the SPA and Boruta feature screening algorithms, the most suitable spectral index form is found; S4, construction and training of inversion model: using the key features as independent variables, using the gradient boosting tree XGBoost algorithm to construct an inversion model for Cu content in soil, and performing classification inversion using decision trees constructed with pH, ​​moisture, and organic matter content as branching criteria; using pH value decision trees, moisture decision trees, and organic matter decision trees to classify the data in the data set; using the calibration data set to train the inversion model, using the validation data set to predict and test the performance of the model; determining the optimal parameter configuration for each model based on the value of the loss function; S5. Prediction of soil Cu content: The soil to be tested is classified according to its pH value, humidity, and organic matter content through the pH value decision tree, humidity decision tree, and organic matter decision tree. The inversion model with the best model performance is selected from the inversion models corresponding to the pH value interval, the inversion model corresponding to the humidity interval, and the inversion model corresponding to the organic matter content interval, and the soil spectral index screened in the corresponding branch is input into the inversion model, and the inversion model outputs the Cu content value in the soil.

2. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S4, the Optuna framework is used in model training. The performance of each model is evaluated by defining an objective function and returning the corresponding evaluation indicators. The search space of hyperparameters is defined in the objective function.

3. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: The two-dimensional spectral index includes a difference index DI, a ratio index RI and a normalized difference index NDIs, DI (Ri, Rj)) = Ri-Rj, The three-dimensional spectral index includes TBI1, TBI2, TBI3, TBI4, and TBI5, which represent the spectral index of the first, second, third, fourth, and fifth transformation forms, respectively. TBI3=(Ri-Rj)-(Rj-Rk), Where Ri, Rj and Rk are the reflectances of any three bands selected from all the bands of the Sentinel-2 satellite.

4. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S1, the method for processing and downloading through the Google Earth Engine cloud platform includes: defining a function for masking clouds and cumulus clouds in Sentinel-2 remote sensing images in Google Earth Engine, the function uses the QA60 band containing image quality assessment information to screen out an image set over a period of time, excludes images with cloud cover exceeding 5%, obtains a filtered and masked image set and takes its average value.

5. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S1, the method for processing and downloading through the Google Earth Engine cloud platform also includes: using cubic convolution interpolation in the Google Earth Engine cloud platform to resample the B1 band and B5-B12 bands in the Level-2A data to a resolution of 10m.

6. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S2, when constructing the humidity decision tree, the surface reflectance spectral data of the Sentinel-2 satellite is transformed by tassel cap, and the multi-spectral band data is reduced to three main components: brightness, greenness and humidity, so as to extract the humidity index T CW Formula to build decision tree; Humidity index T CW The formula is: Among them, B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, and B12 are the reflectances of the B1, B2, B3, B4, B5, B6, B7, B8, B8A, B9, B11, and B12 bands respectively.

7. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S2, the method for obtaining the optimal pH nodes pH1 and pH2 is as follows: select two pH values, pHa and pHb, respectively, pHmin<pHa<pHb<pHmax, to form three pH intervals [pHmin, pHa), [pHa, pHb) and [pHb, pHmax], assign samples to corresponding intervals according to their pH values, repeat the process iteratively, and ensure that the number of samples in each interval is not less than the set value. The evaluation index R pH When the maximum value is reached, the corresponding decision tree is the best, and the corresponding pHa is pH1, and pHb is pH2.

8. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S2, the method for obtaining the optimal humidity nodes WETNESS1 and WETNESS2 is as follows: select two humidity values, WETNESSa and WETNESSb, WETNESSmin<WETNESSa<WETNESSb<WETNESSmax, to form three humidity intervals [WETNESSmin, WETNESSa), [WETNESSa, WETNESSb), [WETNESSb, WETNESSmax], assign samples to corresponding intervals according to their humidity, repeat the process iteratively, and ensure that the number of samples in each interval is not less than the set value. The evaluation index R WETNESS When the maximum value is reached, the corresponding decision tree is the best. At this time, the corresponding WETNESSa is WETNESS1, and WETNESSb is WETNESS2.

9. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: In step S4, the hyperparameters that need to be adjusted in the XGBoost algorithm include: num_boost_round, eta, max_depth, subsample, colsample_bytree, gamma, min_child_weight, rel_lambda, and alpha. The value ranges of each hyperparameter are set as follows: num_boost_round is [500,2000], eta is [0.01,0.3], max_depth is [3,10], subsample is [0.5,1], colsample_bytree is [0.5,1], gamma is [0,1], min_child_weight is [1,10], rel_lambda is [0,1], and alpha is [0,1].

10. The method for inverting the Cu content of agricultural soil using multispectral data according to claim 1, characterized in that: The coefficient of determination R 2 , root mean square error RMSE and residual prediction deviation RPD were used to evaluate the accuracy of each soil Cu content inversion model. Where n represents the number of soil samples, y i and represent the measured value and model predicted value of the i-th soil sample in the validation dataset, respectively. Represents the measured average value of Cu content.

Images