Cultivated soil pH estimation method based on multi-source remote sensing and spatial weighting

Abstract

The invention discloses a cultivated land soil pH estimation method based on multi-source remote sensing and spatial weighting, which comprehensively utilizes multi-source remote sensing environment information such as optics, radar, terrain, weather and the like, constructs multi-source remote sensing environment factor data, and improves response sensitivity to cultivated land soil pH change. A space weighting mechanism is introduced based on latitude and longitude coordinates of each sampling point, so that the adaptability of the model to cultivated land space heterogeneity is enhanced; the weighted random forest model constructed by taking the enhanced feature matrix as the input of the random forest model improves the spatial continuity and prediction precision of the cultivated land soil pH estimation result; the method disclosed by the invention has good universality and expandability, and is suitable for rapidly estimating the pH of the cultivated land soil in multiple regions and multiple scales.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of digital soil mapping and cultivated soil environment monitoring, and particularly relates to a cultivated soil pH estimation method based on multi-source remote sensing and spatial weighting. BACKGROUND

[0002] Soil pH is an important physicochemical parameter for measuring the acid-base properties of soil, and has important guiding significance for crop growth, soil nutrient availability and agricultural management decision-making. In cultivated land management and precision agriculture, it is particularly crucial to obtain high-precision and high-resolution spatial distribution information of soil pH. With the rapid development of remote sensing technology and geographic spatial information technology, the integration of multi-source remote sensing environmental factors and the use of machine learning models to estimate soil properties have become an important research direction in digital soil mapping. Existing researches mostly use optical remote sensing, terrain and meteorological factors to construct soil pH estimation models, and rarely systematically introduce radar remote sensing data. Radar factors can indirectly reflect soil humidity, texture and surface roughness, and have potential response relationship with pH. In addition, the commonly used modeling methods include support vector machine, random forest and gradient boosting tree, which mainly focus on the statistical correlation between environmental factors and soil properties, and do not fully consider the spatial heterogeneity characteristics of soil. Especially in areas with significant differences in cultivated soil or uneven spatial distribution of samples, the estimation accuracy of the model is easily affected. SUMMARY

[0003] The purpose of the present application is to provide a cultivated soil pH estimation method based on multi-source remote sensing and spatial weighting, which can improve the perception ability of the model for the spatial distribution of soil pH and the adaptability to spatially heterogeneous areas, and is suitable for high-precision estimation and dynamic monitoring of cultivated soil pH at regional scale.

[0004] The above-mentioned purpose of the present application is realized by the following technical means:

[0005] A cultivated soil pH estimation method based on multi-source remote sensing and spatial weighting, comprising the following steps:

[0006] Step 1: arranging a plurality of sampling points in the target area, collecting soil samples of a set depth at each sampling point in the target area based on a global land cover dataset, measuring the actual pH value of the soil samples collected at each sampling point, and recording the latitude and longitude coordinates of each sampling point;

[0007] Step 2: extracting multi-source remote sensing environmental factor data corresponding to the position according to the latitude and longitude coordinates of each sampling point, wherein the multi-source remote sensing environmental factor data includes the mean value of multiple secondary feature data;

[0008] Step 3: normalizing the mean value of each type of secondary feature data of each sampling point extracted in step 2, and constructing an original feature matrix , the actual pH value of the soil sample of each sampling point is constructed as a label matrix ;

[0009] Step 4, based on the latitude and longitude coordinates between the sampling points, a spatially weighted feature matrix is constructed , and the original feature matrix is spatially weighted and fused to obtain an enhanced feature matrix ;

[0010] Step 5, set the loss function to enhance the feature matrix , the random forest model is constructed as a weighted random forest model, the weighted random forest model is trained, and the trained weighted random forest model is obtained;

[0011] Step 6, arrange sampling points in the target area to be measured, and sequentially execute steps 1-5 to obtain the trained weighted random forest model of the target area to be measured, input the enhanced feature matrix of the target area to be measured, and obtain the prediction matrix of the soil pH value of the target area to be measured.

[0012] As described above, the spatially weighted feature matrix is constructed by a Gaussian kernel function in step 4; the candidate interval of the Gaussian kernel function is set in step 5, and is traversed with a set step size, each time the reconstructed enhanced feature matrix is calculated, the weighted random forest model is trained, the root mean square error loss function is used as the loss function, the minimum root mean square error is used as the evaluation index, the optimal value is determined, the model parameters are saved after training, and the trained weighted random forest model is obtained;

[0013] is a parameter set to adjust the spatial influence range.

[0014] As described above, step 2 specifically includes the following steps:

[0015] The multi-source remote sensing environmental factor data includes multiple types of first-level feature data, which are optical vegetation index, radar backscattering coefficient, meteorological factor, and terrain factor;

[0016] The optical vegetation index is calculated by the following method: selecting the image of the target area within a set time range, and respectively removing cloud and cloud shadow interference, respectively calculating atmospheric impedance vegetation index, soil adjustment vegetation index, normalized vegetation index, normalized difference water body index, normalized difference red edge vegetation index, conversion vegetation index, normalized difference cloud index, difference vegetation index, enhanced vegetation index, and ratio vegetation index These two types of feature data are calculated.

[0017] The time mean synthesis of each type of vegetation index at the pixel level is performed to obtain the mean value of each type of vegetation index of each pixel in the set time period, and the mean value of each type of vegetation index of the pixel where each sampling point is located is extracted to constitute the optical vegetation index feature of the corresponding sampling point;

[0018] The radar backscattering coefficient is calculated by selecting radar images, extracting VV polarization, VH polarization, and the ratio of VV polarization and VH polarization in the set time period, and performing time mean synthesis of VV polarization, VH polarization, and the ratio of VV polarization and VH polarization in the set time period at the pixel level after correction processing of all radar images to obtain the mean value of VV polarization, VH polarization, and the ratio of VV polarization and VH polarization of each pixel in the target time period, and then extracting the mean value of VV polarization, VH polarization, and the ratio of VV polarization and VH polarization of the pixel where the sampling point is located to constitute the radar backscattering coefficient of the corresponding sampling point;

[0019] The weather factor is calculated by using a public climate data set to obtain the original weather data of the target region in the set time period, resampling to the set resolution, and aligning with the latitude and longitude coordinates of each sampling point to obtain the resampled weather data of the target region;

[0020] The annual temperature data, annual precipitation data, and annual soil moisture data in the resampled weather data of the target region are extracted respectively, and the mean values are calculated to obtain the annual average temperature, annual average precipitation, and annual average soil moisture of the target region, and the annual average temperature, annual average precipitation, and annual average soil moisture at the location of each sampling point are extracted to constitute the weather factor of the corresponding sampling point;

[0021] The terrain factor is calculated by performing spatial interpolation of the digital elevation model based on the digital elevation model to the set resolution, and extracting the elevation, slope, aspect, and terrain shadow index of the location of each sampling point in the spatially interpolated digital elevation model by a terrain analysis module in a cloud computing platform to constitute the terrain factor of the corresponding sampling point.

[0022] The original feature matrix of step 3 as described above , is the number of sampling points, is the number of types of secondary feature data included in the multi-source remote sensing environmental factor data, and the original feature matrix Each row of data of the original feature matrix is a feature data vector composed of the mean values of each type of secondary feature data of a sampling point.

[0023] The label matrix , is a real number set.

[0024] Step 4 specifically comprises the following steps:

[0025] Step 4.1, calculate the Euclidean distance between each sampling point and other sampling points based on the latitude and longitude coordinates of all sampling points, and construct a spatial distance matrix ;

[0026] Step 4.2, construct a spatial weight matrix using a Gaussian kernel function , and normalize the spatial weight matrix to obtain a normalized spatial weight matrix ;

[0027] Step 4.3, multiply the normalized spatial weight matrix and the original feature matrix by matrix multiplication to obtain a spatial weighted feature matrix , , the first row in the spatial weighted feature matrix is the spatial weighted feature vector of the first sampling point.

[0028] The spatial distance matrix is specifically constructed as follows:

[0029] Based on the latitude and longitude coordinates of all sampling points, calculate the Euclidean distance between each sampling point and other sampling points, and construct the Euclidean distance between each sampling point and other sampling points as the spatial distance vector of the corresponding sampling point according to the sampling point sequence number, and all spatial distance vectors of the sampling points are sequentially constructed into a spatial distance matrix , is the Euclidean distance between sampling point and sampling point , and are the sequence numbers of the sampling points, , , each row in the spatial distance matrix is the Euclidean distance between the first sampling point and other sampling points, and when , the Euclidean distance is 0.

[0030] The spatial weight matrix is calculated based on the following formula:

[0031]

[0032] The spatial weight matrix The normalization processing is performed on each row of the spatial weight matrix , so that the sum of each spatial weight of each sampling point is 1, based on the following formula:

[0033]

[0034] In the formula, is the spatial weight between the sampling point and the sampling point .

[0035] The enhanced feature matrix of step 4 as described above is obtained by splicing the original feature matrix and the spatial weighted feature matrix by column:

[0036]

[0037] The enhanced feature matrix .

[0038] A computer device comprising a memory and a processor, the memory storing a computer program, and the processor implementing the steps of the method as described above when executing the computer program.

[0039] A computer readable storage medium having a computer program stored thereon, the computer program being executed by a processor to implement the steps of the method as described above.

[0040] The present application has the following beneficial effects relative to the prior art:

[0041] The method of the present application comprehensively utilizes optical, radar, terrain, meteorological and other multi-source remote sensing environmental information to improve the response sensitivity to the change of the pH value of the cultivated land soil; the spatial weighting mechanism is introduced to enhance the adaptability of the model to the spatial heterogeneity of the cultivated land; the enhanced feature matrix is used as the input of the random forest model to construct a weighted random forest model, which improves the spatial continuity and prediction accuracy of the cultivated land soil pH estimation result; the method of the present application has good universality and scalability, and is suitable for rapid estimation of the pH value of the cultivated land soil in multiple regions and multiple scales. BRIEF DESCRIPTION OF DRAWINGS

[0042] Figure 1 is a flowchart of the method of the present application;

[0043] Figure 2 is a verification diagram of the cultivated land soil pH value inversion model of a certain region in Hubei Province obtained based on the RF model of the present application;

[0044] Figure 3A verification diagram of a cultivated land soil pH value inversion model of a certain region in Hubei Province based on the SWRF model of the application. DETAILED DESCRIPTION

[0045] In order to facilitate those skilled in the art to understand and implement the present application, the present application will be further described in detail below in conjunction with examples, and the examples described herein are only used to illustrate and explain the present application, and are not a limitation on the present application.

[0046] Example 1

[0047] As shown in the figure, a cultivated land soil pH estimation method based on multi-source remote sensing and spatial weighting includes the following steps: Figure 1

[0048] In this embodiment, the cultivated land soil pH of a certain region in Hubei Province is estimated;

[0049] Step 1, multiple sampling points are arranged in the target region, based on the 10-meter resolution WorldCover global land cover dataset, soil samples of a set depth at each sampling point in the target region are collected, the pH value of the soil sample collected at each sampling point is measured, and the latitude and longitude coordinates of each sampling point are recorded;

[0050] In this embodiment, based on the 10-meter resolution WorldCover global land cover dataset released by the European Space Agency (ESA) in 2022, 526 typical cultivated land sampling points are arranged within the cultivated land range of a certain region, soil samples of the surface layer (0.1-20 cm depth) at each sampling point are collected, the soil pH value data is obtained through laboratory measurement, and the latitude and longitude coordinates of each sampling point are recorded synchronously.

[0051] Step 2, obtain multi-source remote sensing environmental factor data: based on the Google Earth Engine (GEE) cloud computing platform, according to the latitude and longitude coordinates of each sampling point, automatically extract the multi-source remote sensing environmental factor data at its corresponding position, and combine the mean value calculation of different time sequences to ensure the representativeness and timeliness of the factors, which includes the following steps:

[0052] The multi-source remote sensing environmental factor data includes multiple types of first-level feature data, which are optical vegetation index, radar backscatter coefficient, meteorological factor, and terrain factor, respectively;

[0053] ​Step 2.1, extract the optical vegetation index, the optical vegetation index includes multiple secondary feature data, respectively, atmospheric resistance vegetation index (ARVI), soil adjustment vegetation index (SAVI), normalized vegetation index (NDVI), normalized difference water index (NDWI), normalized difference red edge vegetation index (NDRE), transformed vegetation index (TVI), normalized difference cloud index (NDCI), difference vegetation index (DVI), enhanced vegetation index (EVI) and ratio vegetation index (RVI), specifically including the following steps:

[0054] Step 2.1.1, based on the multi-spectral remote sensing data (Sentinel-1 is a earth observation satellite launched by European Space Agency (ESA), mainly used to obtain various data of earth's surface; Level-2A data is corrected by atmosphere, which can be directly used for analysis of surface reflectance) of Sentinel-2 (Level-2A), select all available images of the target area from March 1, 2024 to May 31, 2024, and respectively remove cloud and cloud shadow interference, and calculate the various vegetation indexes of the optical vegetation index;

[0055] Step 2.1.2, respectively, the various vegetation indexes calculated in step 2.1.1 are subjected to pixel-level time mean synthesis (i.e. time mean synthesis of each pixel of various vegetation indexes in the set time range (this embodiment is from March 1, 2024 to May 31, 2024) of all images), the mean value of each pixel of various vegetation indexes in the set time range is obtained, and the mean value of each sampling point is extracted to form the optical vegetation index feature of the corresponding sampling point;

[0056] Step 2.2, extract the radar backscattering coefficient, the radar backscattering coefficient includes multiple secondary feature data, respectively, VV polarization (the polarization direction of the electromagnetic wave emitted and received by the radar is vertical), VH polarization (the polarization direction of the electromagnetic wave emitted by the radar is vertical (Vertical), while the polarization direction of the electromagnetic wave received is horizontal (Horizontal)), and the ratio of VV polarization and VH polarization (by comparing VV and VH polarization signals, more effective ground object classification and vegetation monitoring can be carried out), specifically including the following steps:

[0057] Step 2.2.1, select Sentinel-1 SAR radar images (Sentinel-1 is a European Space Agency (ESA) launched Earth observation satellite; SAR is Synthetic Aperture Radar), extract VV polarization, VH polarization, and the ratio of VV polarization and VH polarization from March 1, 2024 to May 31, 2024, reflecting the surface structure characteristics and soil water content.

[0058] Step 2.2.2, after all the radar images are corrected (the correction in this embodiment includes radiation correction and geometric correction), the VV polarization, VH polarization, and the ratio of VV polarization and VH polarization in the target period (March 1, 2024 to May 31, 2024) are synthesized at the pixel level, and the average value of the VV polarization, VH polarization, and the ratio of VV polarization and VH polarization in the target period is obtained. The average value of the VV polarization, VH polarization, and the ratio of VV polarization and VH polarization of the pixel where the sampling point is located is extracted to form the radar backscattering coefficient of the corresponding sampling point.

[0059] Step 2.3, extract meteorological factors, which include multiple secondary feature data, including annual average temperature, annual average precipitation, and annual average soil moisture, specifically:

[0060] Step 2.3.1, use the TERRACLIMATE public climate dataset (global land surface monthly climate and climate water balance dataset) to obtain the original meteorological data of the target area for the whole year (January 1 to December 31, 2024). Since the original meteorological data is usually 1 km or coarser resolution, it is resampled to 10-meter resolution by cubic interpolation method and aligned with the latitude and longitude coordinates of each sampling point to obtain the resampled meteorological data of the target area.

[0061] Step 2.3.2, extract the annual temperature data, annual precipitation data, and annual soil moisture data from the resampled meteorological data of the target area, respectively, and calculate the average value to obtain the annual average temperature, annual average precipitation, and annual average soil moisture of the target area. Extract the annual average temperature, annual average precipitation, and annual average soil moisture at the location of each sampling point to form the meteorological factors of the corresponding sampling point.

[0062] Step 2.4, extract terrain factors, including multiple types of secondary feature data, including elevation (Elevation), slope (Slope), aspect (Aspect), and hillshade index (Hillshade), specifically: based on SRTM (Space Shuttle Radar Topography Mission) 30m resolution digital elevation model (DEM), spatial interpolation of digital elevation model to 10m resolution, in the cloud computing platform, through the terrain analysis module, the elevation, slope, aspect, and hillshade index of each sampling point in the spatially interpolated digital elevation model are extracted to form the terrain factors of the corresponding sampling point.

[0063] This step processes the multi-source remote sensing environmental factor data through spatial alignment and scale unification (i.e., resampling to 10 meters), providing consistent and complete original feature data for subsequent model construction.

[0064] Step 3, normalize the mean of each type of secondary feature data extracted in step 2, and process the normalized mean of each type of secondary feature data of each sampling point as a feature data vector of the corresponding sampling point, randomly divide the feature data vectors of all sampling points into training data set and validation data set, and finally construct the original feature matrix of the feature data vectors of each sampling point in the training data set , the actual pH value of the soil sample collected by each sampling point constitutes the label vector , is a real set;

[0065] wherein the is the number of sampling points in the training data set, is the number of types of secondary feature data included in the multi-source remote sensing environmental factor data, and each row of data in the original feature matrix is a feature data vector of a sampling point;

[0066] The multi-source remote sensing environmental factor data extracted in step 2 of the present embodiment are all continuous variables, therefore, the Min-Max normalization method is used to normalize the mean of 10 types of vegetation index, the mean of VV polarization, the mean of VH polarization, the mean of the ratio VV / VH, the annual average temperature, the annual average precipitation, the annual average soil moisture, the elevation, the slope, the aspect, and the hillshade index, scaled to the interval [0, 1], which can reduce the interference of variable dimension difference on modeling and improve the stability and efficiency of model training; if it is a categorical variable, then use one-hot encoding for normalization.

[0067] The application obtains the measured data of cultivated land soil pH in the target area and the corresponding spatial position coordinates; integrates the multi-source remote sensing environmental factor data related to soil pH, wherein the introduction of the radar backscattering coefficient enhances the response capability to soil structure and humidity state.

[0068] Step 4, in order to enhance the adaptability of the model to spatial heterogeneity, the application introduces a spatial weighting mechanism, based on the geographical spatial distance (longitude and latitude coordinates) between the sampling points, a spatial weighting feature matrix is constructed, and the original feature matrix is fused by spatial weighting according to the spatial weighting feature matrix ,

[0069] Step 4.1, based on the longitude and latitude coordinates of all sampling points, a spatial distance matrix is constructed, specifically: the spatial distance between each sampling point and other sampling points is calculated by using the Euclidean distance formula, the Euclidean distance between each sampling point and other sampling points is constructed as the spatial distance vector of the corresponding sampling point according to the sampling point serial number, and all spatial distance vectors of the sampling points are constructed as the spatial distance matrix in order according to the sampling point serial number, , , , , , , , , , , , , , ,

[0070] Step 4.2, based on the spatial distance matrix , a spatial weight matrix is constructed by using a Gaussian kernel function; in order to improve the numerical stability of the weight, the spatial weight of each row of the spatial weight matrix is normalized, so that the sum of each spatial weight of each sampling point is 1, and a normalized spatial weight matrix is obtained;

[0071] The spatial weight matrix is calculated based on the following formula:

[0072] (1)

[0073] In the formula, is the Euclidean distance between the sampling point With sampling points Spatial weights between them The parameter representing the range of influence of the set adjustment space;

[0074] Normalization is based on the following formula:

[0075] (2)

[0076] Step 4.3: Construct the spatially weighted feature matrix The normalized spatial weight matrix is ​​obtained through matrix multiplication. With the original feature matrix Multiplying them together yields the spatially weighted eigenvalue matrix. ,Right now:

[0077] (3)

[0078] In the formula, the spatial weighted characteristic matrix The first in The line is the first The spatially weighted feature vector of each sampling point represents the feature vector of each sampling point. The comprehensive feature obtained by taking the center as the reference and performing a weighted average after considering the spatial proximity of all other sampling points includes information from the surrounding sampling points.

[0079] Step 4.4: Construct an enhanced feature matrix that integrates spatial perception information. : The original feature matrix With spatial weighted characteristic matrix Concatenate the columns to obtain the enhanced feature matrix. Based on the following formula:

[0080] (4)

[0081] Enhanced feature matrix It also includes the original feature information of each sampling point and the weighted information of its spatial neighborhood, thereby enhancing the model's ability to express geospatial variations.

[0082] Step 5: Enhance the feature matrix As input to the random forest model, the random forest model is constructed into a weighted random forest model (SWRF model), and then a loss function is set (the label matrix of the measured soil pH value is calculated). The weighted random forest model is trained by analyzing the loss between the predicted soil pH value and the predicted value matrix output by the model, and then determining the optimal model. The values ​​are then saved, and the model parameters are obtained to obtain the trained weighted random forest model;

[0083] Training a weighted random forest model specifically involves determining the optimal Gaussian kernel function. The value is determined using a grid search combined with cross-validation, within the set range. Within the candidate interval, traverse the interval with a set step size, repeating steps 4.2 to 4.4 (i.e., reconstructing the enhanced feature matrix) for each iteration. Train a weighted random forest model and determine the optimal model by minimizing the root mean square error (RMSE). value;

[0084] This embodiment is set in Within the candidate interval [0.01, 0.5], the optimal value is obtained by traversing the interval with a step size of 0.01. The value is 0.08.

[0085] Model evaluation: Construct the original feature matrix from the feature data vectors of each sampling point in the validation dataset. Then, perform step 4 to obtain the enhanced feature matrix corresponding to the validation dataset. A random forest model (RF model) was used as a control group to verify the original feature matrix corresponding to the dataset. This serves as input to the random forest model; it is used to validate the augmented feature matrix corresponding to the dataset. This is the input for the random forest model; the model parameters are set consistently, including: number of decision trees: 500; maximum tree depth: 42; other parameters use default values.

[0086] Step 6: Set up sampling points in the target area to be tested, and execute steps 1 to 4 in sequence to obtain the enhanced feature matrix of the target area to be tested. Then, perform step 5 to obtain the trained weighted random forest model of the target region to be tested, and then use the enhanced feature matrix of the target region to obtain the model. Input the data to obtain the prediction matrix of soil pH values ​​for the target area.

[0087] To verify the effectiveness of the spatial adaptive weighting mechanism, the pH of the topsoil (0.1-20cm) in a certain area of ​​Hubei Province, predicted based on multi-source remote sensing environmental factor data, was used as the reference standard. The pH values ​​were calculated from R... 2 The effectiveness of the method of this invention in predicting the pH of topsoil (0.1-20cm) in a certain area of ​​Hubei Province was evaluated using three evaluation indicators: coefficient of determination (RCD), RMSE, and MAE (mean absolute error).

[0088]

[0089] in, R represents 2the proportion of improved ratio, or the proportion of RMSE error reduction, or the proportion of MAE error reduction, the precision of the cultivated land surface layer (0.1-20cm) soil pH in a certain area of Hubei Province predicted by the SWRF model of the application, the precision of the cultivated land surface layer (0.1-20cm) soil pH in a certain area of Hubei Province predicted by the RF model.

[0090] Table 1 is a comparison result table of the prediction of the RF model and the SWRF model

[0091]

[0092] The results show that the precision of the cultivated land surface layer (0.1-20cm) soil pH in a certain area of Hubei Province predicted by the method of the application is significantly improved compared with the precision predicted by the RF model, and the modeling precision R 2 is improved by 30.9% (from 0.463 to 0.606), the RMSE error is reduced by 14.3% (from 0.791 to 0.678), and the MAE error is reduced by 17.5% (from 0.618 to 0.510), as shown in Figure 2 and Figure 3

[0093] In terms of the three evaluation indexes, the SWRF model proposed in the application is superior to the RF model, which shows that the introduction of the spatial weight matrix significantly improves the prediction precision and spatial adaptability of the model.

[0094] In one embodiment, a computer device is also provided, including a memory and a processor, the memory storing a computer program, and the processor implementing the steps in the above method embodiments when executing the computer program.

[0095] In one embodiment, a computer readable storage medium is provided, which stores a computer program, and the computer program implements the steps in the above method embodiments when executed by a processor.

[0096] In one embodiment, a computer program product is provided, including a computer program, and the computer program implements the steps in the above method embodiments when executed by a processor.

[0097] It should be noted that the embodiments described in the application are only examples illustrating the spirit of the application. Those skilled in the art to which the application belongs can make various modifications or supplements to the described embodiments or replace them with similar ways without departing from the spirit of the application or exceeding the scope defined by the appended claims.​

Claims

1. A method for estimating cultivated soil pH based on multi-source remote sensing and spatial weighting, characterized in that, The method comprises the following steps: Step 1, arranging a plurality of sampling points in a target area, collecting soil samples of a set depth of each sampling point in the target area based on a global land cover dataset, measuring the actual pH value of the soil samples collected at each sampling point, and recording the longitude and latitude coordinates of each sampling point; Step 2, extracting multi-source remote sensing environmental factor data of the corresponding position according to the longitude and latitude coordinates of each sampling point, wherein the multi-source remote sensing environmental factor data comprises mean values of a plurality of secondary feature data; Step 3, normalize the mean value of each type of secondary feature data of each sampling point extracted in step 2, and construct into an original feature matrix The actual pH value of the soil sample of each sampling point is constructed into a label matrix ; Step 4, based on the latitude and longitude coordinates between the sampling points, a spatially weighted feature matrix is constructed , and the original feature matrix is spatially weighted and fused to obtain an enhanced feature matrix ; Step 5, setting a loss function to enhance the feature matrix The random forest model is constructed as a weighted random forest model for the input of the random forest model, the weighted random forest model is trained, and the trained weighted random forest model is obtained. Step 6, layout sampling points in the target area to be measured, sequentially execute steps 1~5 to obtain the trained weighted random forest model of the target area to be measured, input the enhanced feature matrix of the target area to be measured obtain the prediction matrix of the soil pH value of the target area to be measured.

2. The method according to claim 1, wherein, The step 4 is to construct a spatially weighted feature matrix by a Gaussian kernel function The step 5 is to set a Gaussian kernel function The step 6 is to traverse the candidate interval with a set step length, and calculate a reconstructed enhanced feature matrix each time The step 7 is to train a weighted random forest model, adopt a root mean square error loss function as a loss function, and determine an optimal value by taking the root mean square error as a judgment index, save model parameters after training, and obtain the trained weighted random forest model Parameter for the set adjustment space impact range. 3.The method of claim 2, wherein, The step 2 specifically comprises the following steps: The multi-source remote sensing environmental factor data comprises a plurality of primary feature data, which are respectively optical vegetation index, radar backscatter coefficient, meteorological factor, and terrain factor; The optical vegetation index is calculated by the following method: selecting images of a set time range in the target area, and respectively removing cloud and cloud shadow interference, and respectively calculating atmospheric impedance vegetation index, soil-adjusted vegetation index, normalized difference water index, normalized difference red edge vegetation index, conversion vegetation index, normalized difference cloud index, difference vegetation index, enhanced vegetation index, and ratio vegetation index, which are secondary feature data; Pixel-level time mean value synthesis is performed on each type of vegetation index to obtain the mean value of each type of vegetation index of each pixel in the set time period, and the mean value of each type of vegetation index of the pixel where each sampling point is located is extracted to form the optical vegetation index feature of the corresponding sampling point; The radar backscatter coefficient is calculated by the following method: selecting radar images, extracting VV polarization, VH polarization, and the ratio of VV polarization and VH polarization in the set time period, and performing pixel-level time mean value synthesis on the VV polarization, VH polarization, and the ratio of VV polarization and VH polarization in the set time period after correction processing of all radar images to obtain the mean value of VV polarization, VH polarization, and the ratio of VV polarization and VH polarization of each pixel in the target time period, and then extracting the mean value of VV polarization, the mean value of VH polarization, and the mean value of the ratio of VV polarization and VH polarization of the pixel where the sampling point is located to form the radar backscatter coefficient of the corresponding sampling point; The meteorological factor is calculated by the following method: obtaining original meteorological data of the target area in the set time period by using a public climate dataset, resampling to a set resolution, and aligning with the longitude and latitude coordinates of each sampling point to obtain resampled meteorological data of the target area; Extracting annual temperature data, annual precipitation data, and annual soil moisture data in the resampled meteorological data of the target area, respectively, and calculating the mean values to obtain the annual average temperature, the annual average precipitation, and the annual average soil moisture of the target area, and extracting the annual average temperature, the annual average precipitation, and the annual average soil moisture of the position where each sampling point is located to form the meteorological factor of the corresponding sampling point; The terrain factor is calculated by: based on a digital elevation model, spatially interpolating the digital elevation model to a set resolution, extracting the elevation, slope, aspect, and terrain shadow index of each sampling point in the spatially interpolated digital elevation model by a terrain analysis module in a cloud computing platform to form the terrain factor of the corresponding sampling point.

4. The method according to claim 2, wherein, The original feature matrix of step 3 , is the number of sampling points, is the number of secondary feature data included in the multi-source remote sensing environmental factor data, the original feature matrix Each row of data of the original feature matrix is a feature data vector composed of the mean values of various types of secondary feature data of a sampling point. The label matrix , is the set of real numbers.

5. The method according to claim 2, wherein, The step 4 specifically includes the following steps: Step 4.1, calculate the Euclidean distance between each sampling point and other sampling points based on the latitude and longitude coordinates of all sampling points, and construct a spatial distance matrix ; Step 4.2, constructing a spatial weight matrix using a Gaussian kernel function , and normalizing the spatial weight matrix to obtain a normalized spatial weight matrix ; Step 4.3: The normalized spatial weight matrix is ​​obtained through matrix multiplication. With the original feature matrix Multiplying them together yields the spatially weighted eigenvalue matrix. , Spatial weighted feature matrix The first in The line is the first Spatial weighted feature vector of each sampling point.

6. The method according to claim 5, wherein, The spatial distance matrix is constructed in particular by Calculate the Euclidean distance between each sampling point and other sampling points based on the latitude and longitude coordinates of all sampling points. Construct a spatial distance vector for each sampling point according to its sampling point index. Then, construct a spatial distance matrix from all sampling point spatial distance vectors in sequence according to their sampling point indices. , Sampling points With sampling points The Euclidean distance between them and These are all the serial numbers of the sampling points. , Spatial distance matrix each The line is the first The Euclidean distance between each sampling point and other sampling points, when When the Euclidean distance is 0, the distance between the two sides is zero.

7. The method according to claim 6, wherein, The spatial weight matrix is calculated based on the following equation: ； normalizing the spatial weight matrix normalizing each row of the spatial weight matrix each sample point to 1 based on the following equation: ； wherein is the spatial weight between sample point and sample point . 8.The method of claim 2, wherein, the enhanced feature matrix of step 4 to obtain the original feature matrix the spatially weighted feature matrix concatenating by column: ； enhanced feature matrix . 9.A computer device, comprising a memory and a processor, wherein the memory stores a computer program, and the computer device is configured to perform the method according to any one of claims 1-8 when the computer program is executed by the processor. The processor implements the steps of the method of any one of claims 1 to 8 when executing the computer program.

10. A computer-readable storage medium having stored thereon a computer program, characterized in that, The computer program, when executed by the processor, implements the steps of the method of any one of claims 1 to 8.

Images