Farmland soil quality intelligent evaluation method and system based on big data

Abstract

The invention discloses a big data-based farmland soil quality intelligent evaluation method and system, and the method comprises the steps: collecting farmland data, carrying out the preprocessing, obtaining soil data and remote sensing data, carrying out the cooperative Kriging inversion through the remote sensing data and the soil data, obtaining surface topographic data, and carrying out the calculation of the surface topographic data. The method comprises the following steps: constructing a three-dimensional partitioning model through a hierarchical Bayesian interpolation method based on a partitioning model according to surface topographic data, combining the first and second three-dimensional partitioning models by using B-spline modeling to obtain a three-dimensional model, carrying out particle size distribution fractal to obtain fractal features, and calculating a soil quality index by using a random forest based on the fractal features and soil data. And inputting the fractal features, the soil data and the soil quality indexes into a fractal neural network for training to obtain a soil quality intelligent evaluation model. According to the method, fractal feature calculation and model training are utilized, and the accuracy of farmland soil quality evaluation and the model generalization ability are improved.

Description

Technical Field

[0001] The present invention relates to the field of intelligent assessment of soil quality, and particularly to an intelligent assessment method and system for farmland soil quality based on big data. Background Art

[0002] With the advancement of agricultural modernization, significant progress has been made in the technology for assessing farmland soil quality, but there are still some deficiencies. Although rapid laboratory testing technologies such as electrochemical analysis and spectral analysis have emerged in recent years, they rely on manual sampling and laboratory analysis, which are inefficient and costly. On the other hand, although remote sensing and digital soil technologies can quickly obtain large-area data, there is still a need to further improve the details and accuracy of soil properties. In addition, there are also limitations in data integration and intelligent applications in existing technologies. For example, it is difficult to integrate soil data in space and time, and it is difficult to achieve efficient sharing. These problems limit the accuracy and coverage area of soil quality assessment, and there is an urgent need to develop a more efficient and intelligent assessment system. Summary of the Invention

[0003] The purpose of the present invention is to provide an intelligent assessment method and system for farmland soil quality based on big data.

[0004] To achieve the above purpose, the present invention is implemented according to the following technical solutions:

[0005] The first aspect of the present invention provides an intelligent assessment method for farmland soil quality based on big data, including:

[0006] S1 Collect farmland data for preprocessing to obtain soil data and remote sensing data;

[0007] S2 Use the remote sensing data and soil data for co-kriging inversion to obtain surface terrain data, calculate the box dimension and slope value based on the surface terrain data, and divide the surface terrain data using the box dimension and slope value to obtain the first terrain area and the second terrain area;

[0008] S3 Use regular grid partitioning for the first terrain area and irregular grid partitioning for the second terrain area, and construct a three-dimensional partition model through hierarchical Bayesian interpolation based on the regular grid partitioning and irregular grid partitioning;

[0009] S4 Divide the terrain transition area based on the boundaries of the first terrain area and the second terrain area, interpolate the terrain transition area using the inverse distance weighted method, and perform nonlinear correction on the interpolation using the radial basis function method to obtain the first and second three-dimensional partition models. Merge the first and second three-dimensional partition models through B-spline modeling to obtain a three-dimensional model;

[0010] S5 Perform fractal of the particle size distribution on the three-dimensional model to obtain fractal features, and calculate soil quality indicators using random forest based on the fractal features and soil data;

[0011] S6 inputs the fractal features, soil data, and soil quality indicators into the fractal neural network for training to obtain an intelligent soil quality assessment model.

[0012] As a further method, the method for collecting farmland data for preprocessing to obtain soil data and remote sensing data includes:

[0013] Collect soil samples at multiple depths in the farmland and remote sensing image data, and perform cleaning and duplicate removal to obtain soil data and remote sensing data. The soil data includes soil bulk density, soil porosity, water-stable aggregate content, soil texture, soil permeability, soil compaction, sampling depth, soil particle size distribution, water-stable aggregate size distribution, and particle size quality.

[0014] As a further method, the method for using remote sensing data and soil data for co-kriging inversion to obtain surface terrain data includes:

[0015] Align the longitude and latitude of soil data points based on the geographical coordinates of remote sensing data points to obtain an initial surface topographic map. Calculate the semi-variogram parameters of soil data points and remote sensing data points respectively according to the initial surface topographic map. The semi-variogram parameters include nugget value, sill value, range, and partial sill value. Calculate the main semi-variogram function of soil data points and the auxiliary semi-variogram function of remote sensing data based on the semi-variogram parameters, calculate the cross semi-variogram function according to the main semi-variogram function and the auxiliary semi-variogram function, obtain a coefficient matrix based on the three semi-variogram functions, use soil data as the main variable and remote sensing data as the auxiliary variable to construct a co-kriging model. The co-kriging model formula is:

[0016]

[0017] Where Z(s) is the surface terrain data or soil attribute spatial data to be predicted at spatial location s, which is the output prediction value of the model. n is the number of soil data points, m is the global mean, that is, the average value of soil data Z s (s i ) or remote sensing data Z r (s j ). The average value of Z s (s i ) is the value of the i-th soil data point, Z r (s j ) is the value of the j-th remote sensing data point. K is the coefficient matrix constructed by the semi-variogram functions and cross semi-variogram functions of soil data points and remote sensing data points, is the inverse matrix element, which is used to determine the weight of soil data point i on the predicted value Z(s), is the inverse matrix element, which is used to determine the weight of remote sensing data point j on the predicted value Z(s), fq It is the q-th element of vector f, used for multiplying with the elements of the inverse matrix, and ∈(s) is the random error term;

[0018] Interpolate the initial surface topographic map area outside the soil data points according to the co-kriging model, and inversely obtain the surface topographic data.

[0019] As a further method, the method of dividing the surface topographic data using the box dimension and slope value to obtain the first topographic area and the second topographic area includes:

[0020] Convert the surface topographic data into a regular grid through a geographic information system, calculate the box dimension for the grid cells of the regular grid, calculate the slope value based on the elevation values of adjacent grid cells, statistically analyze the joint distribution of the box dimension and slope value of all grid cells, determine the box dimension threshold and slope value threshold based on the 95% quantile of the joint distribution points. If the grid cell satisfies that the box dimension is less than the box dimension threshold and the slope value is less than the slope value threshold, it is divided into the first topographic area with simple terrain; if the grid cell satisfies that the box dimension is greater than the box dimension threshold or the slope value is greater than the slope value threshold, it is divided into the second topographic area with complex terrain.

[0021] As a further method, the method of using regular grid blocks for the first topographic area, using irregular grid blocks for the second topographic area, and constructing a three-dimensional block model based on the regular grid blocks and irregular grid blocks through the hierarchical Bayesian interpolation method includes:

[0022] Use a grid size matching the resolution of the surface topographic data for regular grid blocks in the first topographic area; divide the second topographic area into irregular grid blocks through a triangulation algorithm, determine the prior probability distributions of the soil data and surface topographic data by block based on the hierarchical Bayesian interpolation model, construct a likelihood function for the soil data and surface topographic data of the regular grid blocks, irregular grid blocks, and adjacent blocks, multiply the likelihood function with the prior probability distribution according to Bayes' theorem and normalize to obtain the posterior distribution, and iteratively update the parameters of the hierarchical Bayesian interpolation model based on the posterior distribution until the maximum number of iterations or the optimal parameters are reached and then end, to obtain a three-dimensional block model containing soil property and surface topographic data characteristics.

[0023] As a further method, the method of using the inverse distance weighting method to interpolate the terrain transition area and using the radial basis function method for nonlinear correction of the interpolation to obtain the first and second three-dimensional partition models includes:

[0024] Based on the first topographic area, extend 2 grid blocks on both sides centered on the boundary of the three-dimensional block model to form a strip area; check whether there are physical gaps between the three-dimensional block models in the second topographic area. For the three-dimensional block models with physical gaps, extend 3 grid blocks on both sides centered on the model boundary to form a strip area. All the strip areas of the first and second topographic areas are used as the topographic transition area;

[0025] Determine the points to be interpolated in the topographic transition area, search for the known data points within the width range of the strip area of the points to be interpolated, select the correction base points according to the boundary of the topographic transition area, and perform interpolation and correction on the points to be interpolated through the inverse distance weighted-radial basis function combined interpolation formula. Based on the corrected topographic transition area and the three-dimensional block model, obtain the first and second three-dimensional partition models;

[0026] The inverse distance weighted-radial basis function combined interpolation formula is:

[0027]

[0028] where (x,y) are the horizontal and vertical coordinates of the point to be interpolated, n is the total number of known data points around the point to be interpolated, i is the index for traversing the known data points, and the value range is from 1 to n, which is used to process each known data point in turn, c pi is the attribute value of the i-th known data point, a pi is the horizontal coordinate of the i-th known data point, b pi is the vertical coordinate of the i-th known data point, α is the weight parameter, m is the total number of correction base points, j is the index for traversing the correction base points, γ j is the coefficient corresponding to the j-th correction base point, a qj the horizontal coordinate of the j-th correction base point, b qj the vertical coordinate of the j-th correction base point, φ is the radial basis function, k is the mean value of the attribute values c of the known data points pi which is used to adjust the offset of the entire function.

[0029] As a further method, the method for merging the first and second three-dimensional partition models through B-spline modeling to obtain a three-dimensional model includes:

[0030] Determine the first control vertices based on the grid block nodes passing through the boundary of the first three-dimensional partition model using cubic uniform B-splines, calculate the first basis functions according to the recurrence formula of cubic uniform B-spline basis functions using the first control vertices, obtain the first surface point coordinates through the cubic uniform B-spline surface equation using the first control vertices and the first basis functions, connect and interpolate the first surface point coordinates to obtain the first boundary terrain surface, determine the second control vertices based on the grid block nodes passing through the boundary of the second three-dimensional partition model using non-uniform rational B-splines, calculate the second basis functions according to the recurrence formula of non-uniform rational B-spline basis functions using the second control vertices, obtain the second surface point coordinates through the non-uniform rational B-spline surface equation using the second control vertices and the second basis functions, connect and interpolate the second surface point coordinates to obtain the second boundary terrain surface.

[0031] Extend 2.5 grid blocks on both sides of the surface points along the intersection line of the first boundary terrain surface and the second boundary terrain surface to obtain a buffer zone with a width of 5 grid blocks, and use the buffer zone to connect the first boundary terrain surface and the second boundary terrain surface through the non-uniform rational B-spline surface splicing algorithm to obtain a three-dimensional model.

[0032] As a further method, the method for performing grain size distribution fractal on the three-dimensional model to obtain fractal characteristics and calculating soil quality indicators using random forest based on the fractal characteristics and soil data includes:

[0033] Extract the soil particle size distribution, water-stable aggregate size distribution, and particle size quality of the soil data based on the three-dimensional model, extract the number of particles at different particle sizes according to the soil particle size distribution and water-stable aggregate size distribution, calculate the Minkowski dimension of the number of particles at different particle sizes as the overall fractal dimension, calculate the generalized fractal dimension using different particle sizes and the corresponding particle size quality, calculate the aggregate sensitivity index according to the water-stable aggregate size distribution and particle size quality, and perform fractal using the grain size distribution fractal formula based on the overall fractal dimension, generalized fractal dimension, and aggregate sensitivity index to obtain fractal characteristics including the fractal dimension and characteristic particle size. The grain size distribution fractal formula is:

[0034]

[0035] where SFC is the fractal characteristic, FD q is the generalized fractal dimension, q is the particle weight, ASI is the aggregate sensitivity index, ASI max is the maximum value of the aggregate sensitivity index, ε is a very small constant to prevent the denominator from being zero, PDV is the particle size variance, PDA is the average particle size, GFD is the overall fractal dimension, GFD max is the maximum value of the overall fractal dimension, and δ is a very small constant to prevent the denominator from being zero;

[0036] The fractal features and soil data are used as the input data set to train a random forest model. The bootstrap sampling method is used to draw multiple samples from the input data set with replacement to train decision trees, and the decision trees are used to predict soil quality indicators. The soil quality indicators are obtained based on the multi-decision tree voting mechanism.

[0037] As a further method, the method of inputting the fractal features, soil data, and soil quality indicators into a fractal neural network for training to obtain a soil quality intelligent evaluation model includes:

[0038] The fractal features, soil data, and soil quality indicators are subjected to Z-score standardization and organized into a fractal neural network data set. Based on the fractal neural network data set, a training set, a validation set, and a test set are obtained. The number of input layer nodes is set according to the feature dimension of the fractal neural network data set. The ReLU function is used as the activation function for the hidden layer, and the output layer is set to 1 node. The mean squared error is used as the loss function, and the Xavier method is used to initialize the weights of the fractal neural network to obtain a fractal neural network model. The training set is input into the fractal neural network model for training, the gradients of the model parameters are calculated based on the backpropagation algorithm, and the model parameters are iteratively updated according to the gradients. The validation set is used to evaluate the performance of the model based on the mean squared error until the mean squared error is less than 0.01% or the maximum number of training times is reached to obtain a soil quality intelligent evaluation model. The test set is input into the soil quality intelligent evaluation model to obtain an evaluation result, and the mean absolute error, root mean square error, and determination coefficient between the evaluation result and the soil quality indicators are used to measure the evaluation result.

[0039] The second aspect of the present invention provides a farmland soil quality intelligent evaluation system based on big data, including:

[0040] A data acquisition module: used to collect farmland data for preprocessing to obtain soil data and remote sensing data;

[0041] A block modeling module: used to perform co-kriging inversion using remote sensing data and soil data to obtain surface terrain data, calculate the box dimension and slope value based on the surface terrain data, divide the surface terrain data using the box dimension and slope value to obtain a first terrain area and a second terrain area, perform regular grid partitioning on the first terrain area, perform irregular grid partitioning on the second terrain area, and construct a three-dimensional block model based on the regular grid partitioning and irregular grid partitioning through hierarchical Bayesian interpolation;

[0042] A three-dimensional modeling module: based on the boundaries of the first terrain area and the second terrain area, a terrain transition area is divided, the inverse distance weighting method is used to interpolate the terrain transition area, and the radial basis function method is used to perform nonlinear correction on the interpolation to obtain first and second three-dimensional partition models. Based on the first and second three-dimensional partition models, merging is performed through B-spline modeling to obtain a three-dimensional model;

[0043] Evaluation and modeling module: Conduct fractal analysis of the grain size distribution of the 3D model to obtain fractal features. Based on the fractal features and soil data, use random forest to calculate soil quality indicators. Input the fractal features, soil data, and soil quality indicators into the fractal neural network for training to obtain an intelligent soil quality assessment model.

[0044] Compared with the prior art, the embodiments of the present invention have at least the following advantages or beneficial effects:

[0045] (1) The present invention integrates remote sensing data and soil data through co-kriging inversion technology, and constructs a 3D block model by combining hierarchical Bayesian interpolation method, realizing refined modeling of complex terrain. Adopt a differential modeling strategy for different terrain areas, and combine the transitional zone processing method corrected by inverse distance weighting and radial basis function to ensure smooth transition of the 3D model in the terrain continuous area and capture details in the mutation area, improving the spatial interpolation accuracy of surface terrain and soil attribute data, and providing a reliable data basis for subsequent fractal feature calculation and quality assessment.

[0046] (2) The present invention conducts fractal analysis of the grain size distribution of the soil structure, combines indicators such as Minkowski dimension and aggregate sensitivity index, comprehensively reflects the characteristics of soil particle size distribution and aggregate stability, enables the fractal features to more accurately characterize the soil structure characteristics, provides input features with high correlation for the random forest model, and improves the reliability of soil quality indicator calculation.

[0047] (3) The present invention constructs an evaluation model combining fractal neural network and random forest, improves the adaptability of the model to different terrains and soil types, uses hierarchical Bayesian interpolation and B-spline surface splicing technology to integrate multi-source data features, effectively balances the weights of fractal features and soil data, and significantly optimizes indicators such as mean absolute error and root mean square error of the model in soil quality assessment, providing efficient and reliable decision support for precision agriculture management. Description of the Drawings

[0048] Figure 1 It is a flowchart of the steps of an intelligent method for evaluating farmland soil quality based on big data in an embodiment of the present invention. Detailed Embodiments

[0049] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0050] Refer to Figure 1As shown in the figure, the present invention provides an intelligent evaluation method for farmland soil quality based on big data, including:

[0051] S1 Collect farmland data for preprocessing to obtain soil data and remote sensing data;

[0052] In actual evaluation, 500 hectares of farmland is used as the research area, which includes 350 hectares of plain cultivated area, 100 hectares of sloping cultivated land transition area, and 50 hectares of erosion gully area. 200 soil sampling points in the 0-20 cm and 20-40 cm soil layers are arranged according to the terrain zoning, including 140 in the plain area, 50 in the sloping cultivated area, and 10 in the erosion gully area. The soil data is measured by the ring knife method, laser particle size analyzer, etc. The satellite images with a resolution of 10 m of Sentinel-2 satellite and the remote sensing images with a resolution of 5 cm of unmanned aerial vehicle are obtained synchronously, and NDVI and EVI are calculated. After cleaning and deduplication of the above data, a soil data set containing 200 groups of data is finally obtained. For example, the data of a certain point in the plain cultivated area: bulk density 1.15 g / cm 3 、porosity 52%, water-stable aggregates > 0.25 mm 75%, sand grains 13% / silt grains 65% / clay grains 22%, permeability 15 mm / h, compactness 180 kPa, NDVI = 0.81, EVI = 0.79, and corresponding remote sensing data.

[0053] S2 Use remote sensing data and soil data for co-kriging inversion to obtain surface terrain data, calculate the box dimension and slope value based on the surface terrain data, and divide the surface terrain data using the box dimension and slope value to obtain the first terrain area and the second terrain area;

[0054] It should be noted that traditional soil quality evaluation usually uses the form of sparse sampling of soil or feature extraction of remote sensing data. Using co-kriging inversion of remote sensing data and soil data can fuse soil sampling data and remote sensing image data, combine the attribute information of limited sampling points with the spatial details of large-area remote sensing data, generate continuous and high-precision terrain data through co-inversion, and can distinguish different geomorphic units more clearly, avoiding the limitations of simply relying on sampling points or a single data source.

[0055] In actual evaluation, the soil data and remote sensing data are geographically coordinated and unified into WGS84 and UTM50N projections through ArcGIS to generate an initial surface topographic map, and a preliminary surface organic matter raster map with a resolution of 10 m is generated by the inverse distance weighted method with soil organic matter content. 45 soil samples in the plain cultivated area are selected, and the soil organic matter is used as the main variable to calculate the main semi-variogram parameters, and the semi-variogram parameters are obtained, including nugget value 0.8 g 2 / kg 2 、sill value 55 g 2 / kg 2, range 200m; NDVI is used as an auxiliary variable to calculate the auxiliary variogram parameters, including nugget value 0.01, sill value 0.15, and range 300m. At the same time, the cross-semivariogram is calculated to construct the coefficient matrix, and the weights are determined by matrix inversion. For example, at a certain point in the plain area, the soil weight is 0.6 and the remote sensing weight is 0.4. After obtaining the weights, the co-kriging model is applied to interpolate the areas without soil samples such as the transition area of sloping farmland and the erosion gully area to generate surface terrain data with a resolution of 5m. For example, for a certain pixel in the plain area, the organic matter is 24g / kg, NDVI = 0.78, and the elevation is 52m. Specifically, 30 samples are reserved from 200 samples and not involved in the inversion for accuracy verification, including 20 in the plain, 8 in the sloping cultivated area, and 2 in the erosion gully area. The verification results show that the mean absolute error between the interpolated data and the actual soil data is 1.5g / kg, the root mean square error is 2.0g / kg, and the determination coefficient R 2 = 0.92, which is 15% higher than the single kriging method.

[0056] In actual evaluation, based on the 5m-resolution surface terrain data, the box dimension is calculated using the grid covering method. Different grid sizes such as 10m, 20m... 100m are used for covering, and the number of grids N(δ) containing effective terrain data is counted. The absolute value of the slope is obtained by linearly regressing and fitting the relationship between lnN(δ) and lnδ. For example, the box dimension of a certain area in the plain cultivation area is 2.6; using the slope tool in ArcGIS, the slope value is calculated by the differential method of a 3×3 pixel window, and the slope of a certain point in the erosion gully area is 25°; the joint distribution of the box dimension and slope value of all grid cells is counted, and the box dimension threshold and slope value threshold are determined based on the 95% quantile of the joint binning. The area with box dimension ≥ 2.5 and slope value ≤ 15° is divided into the first terrain area. For example, in most areas of the plain cultivation, the box dimension is 2.5 - 2.7 and the slope is 5° - 12°. The area with box dimension < 2.5 or slope value > 15° is divided into the second terrain area. For example, in some areas of the transition area of sloping farmland, the box dimension is 2.3 - 2.4 and the slope is 18° - 22°.

[0057] S3 divides the first terrain area into regular grid blocks and the second terrain area into irregular grid blocks, and constructs a three-dimensional block model through hierarchical Bayesian interpolation based on the regular grid blocks and irregular grid blocks;

[0058] It should be noted that in S2, the slope and box dimension of the farmland terrain have divided the surface terrain into the first terrain area with relatively gentle terrain and little terrain change and the second terrain area with a larger slope and more serious terrain erosion. Therefore, in subsequent modeling, the first terrain area is divided into regular grid blocks and the second terrain area is divided into irregular grid blocks.

[0059] In the actual evaluation, for the first topographic area, the plain area, a grid size of 100m×100m matching the resolution of the surface topographic data is used, and regular grid partitioning is carried out through the ArcGIS Pro fishnet tool; for the second topographic area, the sloping cultivated area, it is divided into irregular grid blocks with side lengths of 50 - 150m using the QGIS triangulation algorithm to adapt to the terrain. Based on the hierarchical Bayesian interpolation model, the prior probability distribution of the soil data in the plain area with its organic matter being N(25, 3 2 ) and the surface topographic data with its elevation being U(200m, 250m) is determined. For the sloping cultivated area, the prior distribution of its compactness being N(1.5MPa, 0.2 2 ) is determined. A likelihood function is constructed for the soil data and surface topographic data of the regular grid blocks, irregular grid blocks, and adjacent blocks. According to Bayes' theorem, the likelihood function is multiplied by the prior probability distribution and normalized to obtain the posterior distribution. Then, based on the posterior distribution, the parameters of the hierarchical Bayesian interpolation model are iteratively updated through WinBUGS software and JAGS software until the maximum number of iterations of 500 times is reached or the optimal parameters are obtained, and then it ends. Finally, a three-dimensional block model containing soil properties and surface topographic data characteristics is obtained.

[0060] For the first topographic area, a grid size of 100m×100m matching the resolution of the surface topographic data is adopted, and regular grid partitioning is carried out through the ArcGIS Pro fishnet tool to ensure that the grid evenly covers the simple terrain area; for the second topographic area, it is divided into irregular grid blocks with side lengths of 50 - 150m using the QGIS triangulation algorithm to adapt to the complex terrain area, such as areas with a slope of 18 - 22° in the sloping cultivated land transition area or areas with a box dimension > 2.5 in the erosion gully area. Based on the hierarchical Bayesian interpolation model, it is determined that the soil data in the first topographic area, such as the organic matter content, follows a normal distribution N(25, 3 2 ) with a mean of 25g / kg and a variance of 3 2 ) and the surface topographic data. For example, the elevation follows a uniform distribution U(200m, 250m) with a prior probability distribution. The soil compactness in the second topographic area follows a normal distribution N(1.5MPa, 0.2 2 ) with a mean of 1.5MPa and a variance of 0.2 2 ) as the prior distribution. A likelihood function is constructed for the soil data (bulk density, porosity, etc.) and surface topographic data (elevation, slope, etc.) of the regular grid blocks, irregular grid blocks, and adjacent blocks. According to Bayes' theorem, the likelihood function is multiplied by the prior probability distribution and normalized to obtain the posterior distribution. Then, based on the posterior distribution, the parameters of the hierarchical Bayesian interpolation model are iteratively updated through WinBUGS software and JAGS software. The maximum number of iterations is set to 500 times, and it is terminated in advance when the parameter change rate < 0.1% for 50 consecutive iterations. Finally, a soil property is obtained (for example, the bulk density in the plain area is 1.15 ± 0.08g / cm3 , porosity of 48.5±4.2% in sloping farmland areas) and surface topography data characteristics (elevation 50-60m in plain areas, slope 15°-25° in sloping farmland areas).

[0061] S4 divides a terrain transition zone based on the boundary between the first terrain zone and the second terrain zone, interpolates the terrain transition zone using an inverse distance weighted method, performs nonlinear correction on the interpolation using a radial basis function method to obtain first and second three-dimensional partition models, and merges the first and second three-dimensional partition models through B-spline modeling to obtain a three-dimensional model;

[0062] It's important to explain that the terrain transition zone, defined as the region between the first and second terrain zones, can have complex topographic features and require special processing. Delineating the terrain transition zone better captures the details of terrain changes, providing a foundation for subsequent interpolation and correction. The inverse distance weighted method (IDW) estimates the values of unknown points based on the values of known points and distance weights, generating smoother and more reasonable terrain data and improving model accuracy. The radial basis function (RBF), a nonlinear correction method, further optimizes interpolation results, better adapting to nonlinear terrain changes and improving model accuracy.

[0063] In the actual assessment, for the first terrain area, two regular grid blocks were extended to both sides of the three-dimensional block model boundary to form a 200-m wide strip area. For the second terrain area with physical gaps (≥5m elevation mutations), three irregular grid blocks were extended to both sides to form a 300-450m wide strip area. A total of 85 hectares of transition area (accounting for 17%) were delineated. Taking the interpolation point P (116.360°E, 39.915°N) in the transition area as an example, eight known data points within 300m of it were searched. For example, for point A in the plain area (116.340°E, 39.910°N, organic matter 24g / kg) and point B in the sloping area (116.370°E, 39.920°N, compaction 1.6MPa), five terrain feature points at the boundary were selected as correction base points, and the inverse distance weighted-radial basis function combination formula was used, where the weight parameter α = 0.6 and the radial basis function was a Gaussian kernel φ(r) = exp(-(r / 50m) 2Interpolate and correct by combining the known point attribute values with the correction base point coefficients. Finally, the organic matter content of point P is obtained as 21 g / kg, and the first and second three-dimensional partition models are obtained. For the boundary grid nodes of the first three-dimensional partition model, cubic uniform B-splines are used, and the control vertices are determined at 20 m intervals, with 5 on each side. The coordinates of the surface points are calculated through the recurrence formula of the basis function to generate the first boundary surface with an elevation fluctuation ≤ 1 m. For the boundary nodes of the second terrain area, non-uniform rational B-splines (NURBS) are used. In the area with a slope > 15°, the control vertices are encrypted at 5 m intervals and the weight factor W = 2 is set to generate the second boundary surface that accurately captures the 30° steep slope turning. Extend 2.5 grid blocks to both sides along the intersection line of the two boundary surfaces, and its total width of 500 m forms a buffer zone. In the buffer zone, the NURBS surface splicing algorithm is used. By adjusting the weights of the control vertices, for example, the vertex weight in the plain area is 1 and the vertex weight of the steep slope in the sloping cultivated land area is 3, the continuous transition of the surface curvature is realized, and the merged three-dimensional model is obtained, with an elevation error < 0.5 m in the transition area.

[0064] S5 Perform fractal analysis on the particle size distribution of the three-dimensional model to obtain fractal characteristics, and calculate the soil quality index using random forest based on the fractal characteristics and soil data.

[0065] In the actual evaluation, based on the soil data of the three-dimensional model, the particle size distribution of soil particles, the particle size distribution of water-stable aggregates, and the particle size quality of 200 samples are extracted. For example, the mass proportion of the particle size distribution of soil particles in sample A01 in the plain cultivated area is 22% for clay (< 0.002 mm), 65% for silt (0.002 - 0.05 mm), and 13% for sand (0.05 - 2 mm). The mass proportion of the particle size distribution of water-stable aggregates in sample B05 in the transition area of the sloping cultivated land is 12% for > 5 mm, 28% for 2 - 5 mm, 40% for 0.25 - 2 mm, and 20% for < 0.25 mm. Then calculate the fractal characteristics. Perform double logarithmic regression on the particle size distribution of soil particles with particle size intervals of [0.001, 0.01, 0.1, 1, 10] mm and the corresponding particle numbers. Through the fitting equation logN(r) = (3 - GFD)logr + C, the GFD of the plain area sample is obtained as 2.75 (determination coefficient R 2 = 0.92), and calculate the generalized fractal dimension FD q Take the particle weight q = 0, 1, 2, 3, and calculate the fractal dimension under different weights. For example, when q = 0, FD0 is the capacity dimension. Through the particle size - mass distribution fitting logm(r) = (3 - FD0)logr, FD0 = 2.68 is obtained, and calculate the aggregate sensitivity index ASI, defined as where the proportion of 0.25 mm aggregates is 80, the particle size variance PDV is 0.08, and ASI = 1000. Set ASI max to 1500, GFD max to 3.0, and δ to 10-6 , the average particle size PDA is 1.2 mm. Substituting it into the fractal formula of particle size distribution, the calculated fractal characteristic SFC = 0.82. Subsequently, an input data set is constructed based on the fractal characteristics, integrating 6-dimensional fractal characteristics GFD = 2.75, FD0 = 2.68, FD1 = 2.55, FD2 = 2.42, FD3 = 2.30, ASI = 1000, a total of 6 dimensions. The soil data includes soil bulk density 1.3 g / cm 3 , porosity 45%, organic matter content 20 g / kg and other 14-dimensional attributes to form 20-dimensional features. The input data set is divided into a training set (140 samples), a validation set (40 samples), and a test set (20 samples) according to 7:2:1. The random forest model is used for training, setting the number of decision trees to 100 and the maximum depth to 15. The Bootstrap self-sampling method is adopted, and 80% of the samples are drawn with replacement for training each time. The hyperparameters are optimized through grid search, obtaining the best learning rate of 0.1 and the minimum number of samples in leaf nodes of 5. The training outputs the soil fertility index (SFI) and the soil erosion resistance index (SEI). For example, sample A01 gets SFI = 85 and SEI = 78, and sample B05 gets SFI = 65 and SEI = 60. Finally, the model performance is verified. The test set results show that the mean absolute error of SFI is 3.2, the root mean square error is 4.8, R 2 = 0.91, the mean absolute error of SEI is 4.5, the root mean square error is 5.2, R 2 = 0.89. The fractal characteristics are positively correlated with the soil organic matter content r = 0.78 and negatively correlated with the soil erosion modulus r = -0.65.

[0066] S6 Inputs the fractal characteristics, soil data, and soil quality indicators into the fractal neural network for training to obtain a soil quality intelligent evaluation model.

[0067] In actual evaluation, the fractal characteristics, soil data, and soil quality indicators are subjected to Z-score standardization. For example, the mean value of soil bulk density is 1.2, and after standardization, it is -0.5. After sorting, the fractal neural network data set is obtained and divided into a training set, a validation set, and a test set. When constructing the fractal neural network model, the input layer is set with 2 nodes (consistent with the fractal feature dimension), and the number of hidden layer nodes is dynamically optimized according to the overall fractal dimension (GFD). For example, when GFD > 2.8, the first hidden layer is set to 48 nodes (complex soil particle distribution), when GFD < 2.5, it is set to 24 nodes, and the second hidden layer is adjusted proportionally. For example, 48 nodes correspond to 24 nodes, and 24 nodes correspond to 12 nodes. The activation function of the hidden layer adopts the ReLU function, and the loss function is based on ) and adds a fractal feature regularization term (FD j is the true fractal feature value, is the predicted value, and j represents different fractal features. The total loss function L = MSE + L fractal , the model weights are initialized by Xavier. The training set is input into the model, and the Ada optimizer is used. The learning rate is set to 0.001. Based on the backpropagation algorithm, the parameter gradients are calculated. The performance is evaluated by calculating the total loss on the validation set every 5 epochs. The training stops when the total loss on the validation set is less than 0.01% or the maximum number of training times, 500 times, is reached. The test set is input into the model, and the mean absolute error, root mean square error, and coefficient of determination are calculated to evaluate the model performance.

[0068] In this embodiment, the method for collecting farmland data for preprocessing to obtain soil data and remote sensing data includes:

[0069] Soil samples at multiple depths of the farmland and remote sensing image data are collected, cleaned, and duplicate data is removed to obtain soil data and remote sensing data. The soil data includes soil bulk density, soil porosity, water-stable aggregate content, soil texture, soil permeability, soil compaction, sampling depth, soil particle size distribution, water-stable aggregate size distribution, and particle size quality.

[0070] In this embodiment, the method for using remote sensing data and soil data for co-kriging inversion to obtain surface terrain data includes:

[0071] Based on the geographical coordinates of the remote sensing data points, the soil data points are aligned in longitude and latitude to obtain an initial surface topographic map. The semi-variogram parameters of the soil data points and remote sensing data points are calculated respectively according to the initial surface topographic map. The semi-variogram parameters include nugget value, sill value, range, and partial sill value. The main semi-variogram function of the soil data points and the auxiliary semi-variogram function of the remote sensing data are calculated based on the semi-variogram parameters. The cross-semi-variogram function is calculated according to the main semi-variogram function and the auxiliary semi-variogram function. A coefficient matrix is obtained based on the three semi-variogram functions. Taking the soil data as the main variable and using the remote sensing data as the auxiliary variable, a co-kriging model is constructed. The co-kriging model formula is:

[0072]

[0073] where Z(s) is the surface terrain data or soil property spatial data to be predicted at the spatial position s, which is the output predicted value of the model. n is the number of soil data points, m is the global mean, that is, the average value of the soil data Z s (s i ) or the remote sensing data Z r (s j ), Z s (s i ) is the value of the i-th soil data point, and Z r (s j) is the value of the j-th remote sensing data point, K is the coefficient matrix constructed by the semivariogram and cross-semivariogram of soil data points and remote sensing data points, is the inverse matrix element, used to determine the weight of soil data point i on the predicted value Z(s), is the inverse matrix element, used to determine the weight of remote sensing data point j on the predicted value Z(s), f q is the q-th element of vector f, used to multiply with the inverse matrix element, ∈(s) is the random error term;

[0074] Interpolate the initial surface topographic map area outside the soil data points according to the co-kriging model, and inversely obtain the surface topographic data.

[0075] In this embodiment, the method of using the box dimension and slope value to divide the surface topographic data to obtain the first topographic area and the second topographic area includes:

[0076] Convert the surface topographic data into a regular grid through a geographic information system, calculate the box dimension for the grid cells of the regular grid, calculate the slope value based on the elevation values of adjacent grid cells, statistically analyze the joint distribution of the box dimension and slope value of all grid cells, determine the box dimension threshold and slope value threshold based on the 95% quantile of the joint distribution points. If the grid cell satisfies that the box dimension is less than the box dimension threshold and the slope value is less than the slope value threshold, it is divided into the first topographic area with simple terrain. If the grid cell satisfies that the box dimension is greater than the box dimension threshold or the slope value is greater than the slope value threshold, it is divided into the second topographic area with complex terrain.

[0077] In this embodiment, the method of using regular grid partitioning for the first topographic area, using irregular grid partitioning for the second topographic area, and constructing a three-dimensional partitioning model based on regular grid partitioning and irregular grid partitioning through hierarchical Bayesian interpolation method includes:

[0078] Use a grid size matching the resolution of the surface topographic data for regular grid partitioning of the first topographic area; divide the second topographic area into irregular grid partitions through a triangulation algorithm. Determine the prior probability distributions of soil data and surface topographic data by partitioning based on the hierarchical Bayesian interpolation model. Construct likelihood functions for the soil data and surface topographic data of regular grid partitions, irregular grid partitions, and adjacent partitions. Multiply the likelihood function with the prior probability distribution according to Bayes' theorem and normalize to obtain the posterior distribution. Iteratively update the parameters of the hierarchical Bayesian interpolation model based on the posterior distribution until the maximum number of iterations or the optimal parameters are reached and then end to obtain a three-dimensional partitioning model containing soil property and surface topographic data characteristics.

[0079] In this embodiment, the method of using the inverse distance weighting method to interpolate the terrain transition area and using the radial basis function method for nonlinear correction of the interpolation to obtain the first and second three-dimensional partitioning models includes:

[0080] Based on the first topographic area, extend 2 grid blocks to both sides centered on the boundary of the three-dimensional block model to form a strip area; check whether there are physical gaps between the three-dimensional block models in the second topographic area, and for the three-dimensional block models with physical gaps, extend 3 grid blocks to both sides centered on the model boundary to form a strip area, and take all the strip areas of the first and second topographic areas as the topographic transition area;

[0081] Determine the points to be interpolated in the topographic transition area, search for the known data points within the width range of the strip area of the points to be interpolated, select the correction base points according to the boundary of the topographic transition area, and perform interpolation and correction on the points to be interpolated through the inverse distance weighted-radial basis function combined interpolation formula, and obtain the first and second three-dimensional partition models based on the corrected topographic transition area and the three-dimensional block model.

[0082] The inverse distance weighted-radial basis function combined interpolation formula is:

[0083]

[0084] where (x, y) are the horizontal and vertical coordinates of the point to be interpolated, n is the total number of known data points around the point to be interpolated, i is the index for traversing the known data points, and the value range is from 1 to n, which is used to process each known data point in turn, c pi is the attribute value of the i-th known data point, a pi is the horizontal coordinate of the i-th known data point, b pi is the vertical coordinate of the i-th known data point, α is the weight parameter, m is the total number of correction base points, j is the index for traversing the correction base points, γ j is the coefficient corresponding to the j-th correction base point, a qj is the horizontal coordinate of the j-th correction base point, b qj is the vertical coordinate of the j-th correction base point, φ is the radial basis function, k is the mean value of the known data point attribute value c pi which is used to adjust the offset of the entire function.

[0085] In this embodiment, the method for merging the first and second three-dimensional partition models through B-spline modeling to obtain a three-dimensional model includes:

[0086] Determine the first control vertices based on the grid block nodes passing through the boundary of the first three-dimensional partition model using cubic uniform B-splines, calculate the first basis functions according to the recurrence formula of cubic uniform B-spline basis functions using the first control vertices, obtain the coordinates of the first surface points using the first control vertices and the first basis functions through the cubic uniform B-spline surface equation, connect and interpolate the coordinates of the first surface points to obtain the first boundary terrain surface, determine the second control vertices based on the grid block nodes passing through the boundary of the second three-dimensional partition model using non-uniform rational B-splines, calculate the second basis functions according to the recurrence formula of non-uniform rational B-spline basis functions using the second control vertices, obtain the coordinates of the second surface points using the second control vertices and the second basis functions through the non-uniform rational B-spline surface equation, connect and interpolate the coordinates of the second surface points to obtain the second boundary terrain surface.

[0087] Extend 2.5 grid blocks on both sides of the surface points along the intersection line of the first boundary terrain surface and the second boundary terrain surface to obtain a buffer zone with a width of 5 grid blocks, and use the buffer zone to connect the first boundary terrain surface and the second boundary terrain surface through the non-uniform rational B-spline surface splicing algorithm to obtain a three-dimensional model.

[0088] In this embodiment, the method for performing grain size distribution fractal on the three-dimensional model to obtain fractal features and calculating soil quality indicators using random forest based on the fractal features and soil data includes:

[0089] Extract the soil particle size distribution, water-stable aggregate size distribution and particle size quality from the soil data of the three-dimensional model, extract the number of particles at different particle sizes according to the soil particle size distribution and water-stable aggregate size distribution, calculate the Minkowski dimension as the overall fractal dimension for the number of particles at different particle sizes, calculate the generalized fractal dimension using different particle sizes and the corresponding particle size quality, calculate the aggregate sensitivity index according to the water-stable aggregate size distribution and particle size quality, and perform fractal using the grain size distribution fractal formula based on the overall fractal dimension, generalized fractal dimension and aggregate sensitivity index to obtain fractal features including the fractal dimension and characteristic particle size. The grain size distribution fractal formula is:

[0090]

[0091] Where SFC is the fractal feature, FD q is the generalized fractal dimension, q is the particle weight, ASI is the aggregate sensitivity index, ASI max is the maximum value of the aggregate sensitivity index, ε is a very small constant to prevent the denominator from being zero, PDV is the particle size variance, PDA is the average particle size, GFD is the overall fractal dimension, GFD max is the maximum value of the overall fractal dimension, and δ is a very small constant to prevent the denominator from being zero;

[0092] The fractal features and soil data are used as input datasets to train a random forest model. The bootstrap sampling method is used to draw multiple samples from the input datasets with replacement to train decision trees, and the decision trees are used to predict soil quality indicators. Based on the voting mechanism of multiple decision trees, soil quality indicators are obtained.

[0093] In this embodiment, the method for inputting the fractal features, soil data, and soil quality indicators into a fractal neural network for training to obtain a soil quality intelligent evaluation model includes:

[0094] The fractal features, soil data, and soil quality indicators are subjected to Z-score standardization and organized into a fractal neural network dataset. Based on the fractal neural network dataset, a training set, a validation set, and a test set are obtained. The number of input layer nodes is set according to the feature dimension of the fractal neural network dataset. The ReLU function is used as the activation function for the hidden layer, and the output layer is set to 1 node. The mean squared error is used as the loss function, and the Xavier method is used to initialize the weights of the fractal neural network to obtain a fractal neural network model. The training set is input into the fractal neural network model for training, and the gradients of the model parameters are calculated based on the backpropagation algorithm. The model parameters are iteratively updated according to the gradients. The validation set is used to evaluate the performance of the model according to the mean squared error until the mean squared error is less than 0.01% or the maximum number of training times is reached to obtain a soil quality intelligent evaluation model. The test set is input into the soil quality intelligent evaluation model to obtain an evaluation result, and the mean absolute error, root mean square error, and coefficient of determination of the evaluation result and the soil quality indicators are used to measure the evaluation result.

[0095] The second aspect of the present invention also provides a farmland soil quality intelligent evaluation system based on big data, including:

[0096] A data acquisition module: used to collect farmland data for preprocessing to obtain soil data and remote sensing data;

[0097] A block modeling module: used to perform co-kriging inversion using remote sensing data and soil data to obtain surface terrain data, calculate the box dimension and slope value based on the surface terrain data, divide the surface terrain data using the box dimension and slope value to obtain a first terrain area and a second terrain area, perform regular grid partitioning on the first terrain area, perform irregular grid partitioning on the second terrain area, and construct a three-dimensional block model based on the regular grid partitioning and irregular grid partitioning through hierarchical Bayesian interpolation;

[0098] A three-dimensional modeling module: based on the boundaries of the first terrain area and the second terrain area, divide a terrain transition area, perform interpolation on the terrain transition area using the inverse distance weighting method, and perform nonlinear correction on the interpolation using the radial basis function method to obtain first and second three-dimensional partition models, and merge them through B-spline modeling based on the first and second three-dimensional partition models to obtain a three-dimensional model;

[0099] Evaluation and modeling module: Perform grain size distribution fractal on the three-dimensional model to obtain fractal features. Use random forest to calculate soil quality indicators based on the fractal features and soil data. Input the fractal features, soil data, and soil quality indicators into the fractal neural network for training to obtain an intelligent soil quality evaluation model.

[0100] The above content is only an example and explanation of the structure of the present invention. Those skilled in the art of this technology can make various modifications or supplements to the described specific embodiments or use similar methods for substitution. As long as they do not deviate from the structure of the invention or exceed the scope defined by this claims, they should fall within the protection scope of the present invention.

Claims

1. An intelligent evaluation method for farmland soil quality based on big data, characterized in that, It includes the following steps: S1 Collect farmland data for preprocessing to obtain soil data and remote sensing data; S2 Use the remote sensing data and soil data for co-Kriging inversion to obtain surface terrain data, calculate the box dimension and slope value based on the surface terrain data, and divide the surface terrain data using the box dimension and slope value to obtain the first terrain area and the second terrain area; S3 Use regular grid partitioning for the first terrain area and irregular grid partitioning for the second terrain area, and construct a three-dimensional partition model based on the regular grid partitioning and irregular grid partitioning through hierarchical Bayesian interpolation; S4 Divide the terrain transition area based on the boundaries of the first terrain area and the second terrain area, interpolate the terrain transition area using the inverse distance weighting method, and perform nonlinear correction on the interpolation using the radial basis function method to obtain the first and second three-dimensional partition models. Merge them through B-spline modeling based on the first and second three-dimensional partition models to obtain a three-dimensional model; S5 Perform fractal of the grain size distribution on the three-dimensional model to obtain fractal characteristics, and calculate the soil quality index using random forest based on the fractal characteristics and soil data; S6 Input the fractal characteristics, soil data, and soil quality index into the fractal neural network for training to obtain a soil quality intelligent evaluation model.

2. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, wherein, The method for collecting farmland data for preprocessing to obtain soil data and remote sensing data includes: Collect soil samples at multiple depths of the farmland and remote sensing image data, and perform cleaning and duplicate removal to obtain soil data and remote sensing data. The soil data includes soil bulk density, soil porosity, water-stable aggregate content, soil texture, soil permeability, soil compaction, sampling depth, soil particle size distribution, water-stable aggregate particle size distribution, and particle size quality.

3. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, characterized in that The method for using the remote sensing data and soil data for co-Kriging inversion to obtain surface terrain data includes: Align the longitude and latitude of the soil data points based on the geographical coordinates of the remote sensing data points to obtain an initial surface topographic map. Calculate the semi-variogram parameters of the soil data points and remote sensing data points respectively according to the initial surface topographic map. The semi-variogram parameters include nugget value, sill value, range, and partial sill value. Calculate the main semi-variogram function of the soil data points and the auxiliary semi-variogram function of the remote sensing data based on the semi-variogram parameters, calculate the cross semi-variogram function according to the main semi-variogram function and the auxiliary semi-variogram function, obtain the coefficient matrix based on the three semi-variogram functions, use the soil data as the main variable, and use the remote sensing data as the auxiliary variable to construct a co-Kriging model. The co-Kriging model formula is: Among them, Z(s) is the predicted surface terrain data or soil property spatial data at spatial position s, which is the output prediction value of the model. n is the number of soil data points, and m is the global mean, that is, the mean value of soil data Z s (s i ) or remote sensing data Z r (s j ). The value of Z s (s i ) is the value of the i-th soil data point, and the value of Z r (s j ) is the value of the j-th remote sensing data point. K is the coefficient matrix constructed by the semivariogram and cross-semivariogram of soil data points and remote sensing data points. is the inverse matrix element, which is used to determine the weight of soil data point i on the predicted value Z(s). is the inverse matrix element, which is used to determine the weight of remote sensing data point j on the predicted value Z(s). f q is the q-th element of vector f, which is used to multiply with the inverse matrix element. ∈(s) is the random error term; Interpolate the area of the initial surface topographic map outside the soil data points according to the co-Kriging model to inversely obtain the surface terrain data.

4. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, wherein, The method for dividing the surface terrain data using the box dimension and slope value to obtain the first terrain area and the second terrain area includes: Convert the surface terrain data into regular grids through a geographic information system, calculate the box dimension of the grid cells of the regular grids, calculate the slope value based on the elevation values of adjacent grid cells, statistically analyze the joint distribution of the box dimension and slope value of all grid cells, determine the box dimension threshold and slope value threshold based on the 95% quantile of the joint distribution points. If the grid cell satisfies that the box dimension is less than the box dimension threshold and the slope value is less than the slope value threshold, it is divided into the first terrain area with simple terrain. If the grid cell satisfies that the box dimension is greater than the box dimension threshold or the slope value is greater than the slope value threshold, it is divided into the second terrain area with complex terrain.

5. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, wherein, The method of using regular grid blocks for the first terrain area, using irregular grid blocks for the second terrain area, and constructing a three-dimensional block model through hierarchical Bayesian interpolation based on regular grid blocks and irregular grid blocks includes: Use a grid size matching the resolution of the surface terrain data to perform regular grid blocks on the first terrain area; divide the second terrain area into irregular grid blocks through a triangulation algorithm, determine the prior probability distributions of soil data and surface terrain data through block-based hierarchical Bayesian interpolation models, construct likelihood functions for the soil data and surface terrain data of regular grid blocks, irregular grid blocks, and adjacent blocks, multiply the likelihood function by the prior probability distribution according to Bayes' theorem and normalize to obtain the posterior distribution, and iteratively update the parameters of the hierarchical Bayesian interpolation model based on the posterior distribution until the maximum number of iterations or optimal parameters are reached and then end to obtain a three-dimensional block model containing soil attribute and surface terrain data characteristics.

6. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, characterized in that, The method of using the inverse distance weighted method to interpolate the terrain transition area and using the radial basis function method to perform nonlinear correction on the interpolation to obtain the first and second three-dimensional partition models includes: Based on the first terrain area, extend 2 grid blocks to both sides centered on the boundary of the three-dimensional block model to form a strip area; check whether there are physical gaps between the three-dimensional block models of the second terrain area, and for the three-dimensional block models with physical gaps, extend 3 grid blocks to both sides centered on the model boundary to form a strip area, and use all the strip areas of the first terrain area and the second terrain area as the terrain transition area; Determine the interpolation points for the terrain transition area, search for known data points within the width range of the strip area of the interpolation points, select correction base points according to the boundary of the terrain transition area, perform interpolation and correction on the interpolation points through the inverse distance weighted-radial basis function combined interpolation formula, and obtain the first and second three-dimensional partition models based on the corrected terrain transition area and the three-dimensional block model; The inverse distance weighted-radial basis function combined interpolation formula is: where (x, y) are the abscissa and ordinate of the point to be interpolated, n is the total number of known data points around the point to be interpolated, i is the index for traversing the known data points, with a value range from 1 to n, used to process each known data point in sequence, c pi is the attribute value of the i-th known data point, a pi is the abscissa of the i-th known data point, b pi is the ordinate of the i-th known data point, α is the weight parameter, m is the total number of correction base points, j is the index for traversing the correction base points, γ j is the coefficient corresponding to the j-th correction base point, a qj the abscissa of the j-th correction base point, b qj the ordinate of the j-th correction base point, φ is the radial basis function, k is the mean value of the attribute values c pi of the known data points, used to adjust the offset of the entire function.

7. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, characterized in that The method of merging through B-spline modeling based on the first and second three-dimensional partition models to obtain a three-dimensional model includes: Determine the first control vertices based on the grid block nodes passing through the boundary of the first three-dimensional partition model by cubic uniform B-spline, calculate the first basis functions according to the first control vertices through the recurrence formula of cubic uniform B-spline basis functions, obtain the first surface point coordinates through the cubic uniform B-spline surface equation using the first control vertices and the first basis functions, connect and interpolate the first surface point coordinates to obtain the first boundary terrain surface, determine the second control vertices based on the grid block nodes passing through the boundary of the second three-dimensional partition model by non-uniform rational B-spline, calculate the second basis functions according to the second control vertices through the recurrence formula of non-uniform rational B-spline basis functions, obtain the second surface point coordinates through the non-uniform rational B-spline surface equation using the second control vertices and the second basis functions, connect and interpolate the second surface point coordinates to obtain the second boundary terrain surface. Extend 2.5 grid blocks on both sides of the surface points along the intersection line of the first boundary terrain surface and the second boundary terrain surface to obtain a buffer zone with a width of 5 grid blocks, and use the buffer zone to connect the first boundary terrain surface and the second boundary terrain surface through the non-uniform rational B-spline surface splicing algorithm to obtain a three-dimensional model.

8. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, characterized in that The method of performing grain size distribution fractal on the three-dimensional model to obtain fractal features and calculating soil quality indicators using random forest based on the fractal features and soil data includes: Extract the soil particle size distribution, water-stable aggregate size distribution and particle size quality from the soil data of the three-dimensional model, extract the number of particles at different particle sizes according to the soil particle size distribution and water-stable aggregate size distribution, calculate the Minkowski dimension as the overall fractal dimension for the number of particles at different particle sizes, calculate the generalized fractal dimension using different particle sizes and the corresponding particle size quality, calculate the aggregate sensitivity index according to the water-stable aggregate size distribution and particle size quality, perform fractal using the grain size distribution fractal formula based on the overall fractal dimension, generalized fractal dimension and aggregate sensitivity index to obtain fractal features including fractal dimension and characteristic particle size, and the grain size distribution fractal formula is: Among them, SFC is the fractal feature, FD q is the generalized fractal dimension, q is the particle weight, ASI is the aggregate sensitivity index, and ASI max is the maximum value of the aggregate sensitivity index, ε is the minimum constant to prevent the denominator from being zero, PDV is the particle size variance, PDA is the average particle size, GFD is the overall fractal dimension, GFD max is the maximum value of the overall fractal dimension, and δ is a minimum constant to prevent the denominator from being zero; Input the fractal features and soil data as the input data set into the random forest model for training, use the bootstrap sampling method to draw multiple samples with replacement from the input data set to train decision trees, use the decision trees to predict the soil quality indicators, and obtain the soil quality indicators based on the multi-decision tree voting mechanism.

9. The intelligent evaluation method for farmland soil quality based on big data according to claim 1, wherein The method of inputting the fractal features, soil data, and soil quality indicators into the fractal neural network for training to obtain the soil quality intelligent evaluation model includes: The fractal features, soil data, and soil quality indicators are Z-score standardized and organized into a fractal neural network dataset. Based on the fractal neural network dataset, training sets, validation sets, and test sets are obtained. The number of input layer nodes is set according to the feature dimension of the fractal neural network dataset. The ReLU function is used as the activation function for the hidden layer, and the output layer is set to 1 node. The mean square error is used as the loss function, and the Xavier method is used to initialize the weights of the fractal neural network to obtain the fractal neural network model. The training set is input into the fractal neural network model for training, and the gradients of the model parameters are calculated based on the backpropagation algorithm. The model parameters are iteratively updated according to the gradients. The validation set is used to evaluate the model performance based on the mean square error until the mean square error is less than 0.01% or the maximum number of training times is reached, obtaining the intelligent soil quality assessment model. The test set is input into the intelligent soil quality assessment model to obtain the assessment results, and the mean absolute error, root mean square error, and determination coefficient between the assessment results and the soil quality indicators are used to measure the assessment results.

10. An intelligent evaluation system for farmland soil quality based on big data, which is used to execute the method described in any one of claims 1-9, characterized in that, Including: Data acquisition module: used to collect farmland data for preprocessing to obtain soil data and remote sensing data; Block modeling module: used to perform co-kriging inversion using remote sensing data and soil data to obtain surface terrain data, calculate the box dimension and slope value based on the surface terrain data, divide the surface terrain data using the box dimension and slope value to obtain the first terrain area and the second terrain area, perform regular grid block division on the first terrain area, perform irregular grid block division on the second terrain area, and construct a three-dimensional block model based on the regular grid block division and irregular grid block division through hierarchical Bayesian interpolation; Three-dimensional modeling module: divide the terrain transition area based on the boundaries of the first terrain area and the second terrain area, perform interpolation on the terrain transition area using the inverse distance weighting method, and perform nonlinear correction on the interpolation using the radial basis function method to obtain the first and second three-dimensional partition models. The first and second three-dimensional partition models are merged through B-spline modeling to obtain the three-dimensional model; Evaluation modeling module: perform fractal of the grain size distribution on the three-dimensional model to obtain fractal features, calculate soil quality indicators using random forest based on the fractal features and soil data, and input the fractal features, soil data, and soil quality indicators into the fractal neural network for training to obtain the intelligent soil quality assessment model.

Images