A method and system for analyzing and evaluating the structure of farmland soil

Abstract

This application relates to a method and system for analyzing and evaluating farmland soil structure. The method includes: reconstructing a first point cloud using a multi-view stereo matching algorithm, and simultaneously acquiring a second point cloud containing color information; preprocessing and extracting targets from both point clouds; introducing adaptive voxel downsampling into the second point cloud; dynamically iteratively estimating the optimal voxel side length based on the downsampling rate, bounding box volume, and total number of points; compressing data using the voxel centroid while preserving geometric features; using the first optimized point cloud as a benchmark, performing coarse registration with side length constraints using an improved RANSAC algorithm based on fast point feature histograms; completing fine registration using a KD-Tree accelerated iterative nearest-point algorithm to generate a fused point cloud; extracting elevation values ​​and removing overlapping points; calculating the root mean square height and autocorrelation length; quantifying vertical undulation and horizontal spatial correlation; and finally, classifying the soil structure into smooth, transitional, and rough levels according to preset thresholds and outputting an evaluation report.

Description

Technical Field

[0001] This application relates to the field of soil analysis technology, and in particular to a method and system for analyzing and evaluating the structure of farmland soil. Background Technology

[0002] Soil surface micro-topography is a core physical parameter characterizing the structural state of farmland soil. Its quantitative index, soil roughness, directly determines surface hydrological processes, soil erosion intensity, seed contact and germination environment, and the quality of agricultural machinery operations. Accurate, rapid, and non-destructive acquisition of farmland soil geomorphological information is a key prerequisite for realizing smart agricultural management such as precision seeding, variable irrigation, and conservation tillage. Research shows that soil roughness has become an important auxiliary parameter for soil moisture retrieval, crop growth simulation, and the evaluation of high-standard farmland construction.

[0003] Current soil micro-topography reconstruction technology is mainly based on photogrammetry methods that combine motion-restoring structures and multi-view stereo matching. It uses consumer-grade digital cameras to collect multi-view two-dimensional image sequences and utilizes feature point matching and bundle adjustment to solve three-dimensional geometric information. It has the advantages of low equipment cost, realistic texture restoration, and high detail resolution.

[0004] Regarding the aforementioned techniques, photogrammetry is extremely sensitive to lighting conditions. In the shadow areas of large soil clods after farmland plowing, under low light conditions in the early morning or late evening, or under diffuse light conditions on cloudy days, the number of feature points extracted is sharply reduced, and large areas of voids are easily generated in the reconstructed point cloud. In particular, the reconstruction integrity of the depression areas on the undulating and rough soil surface is seriously insufficient, which affects the accuracy of soil structure analysis.

[0005] Based on this, this application provides a method and system for analyzing and evaluating the structure of farmland soil. Summary of the Invention

[0006] To address the issue that photogrammetry methods are extremely sensitive to lighting conditions, the number of feature points extracted is drastically reduced in the shadow areas of large soil clods after farmland plowing, in low light conditions in the early morning or late evening, or under diffuse light conditions on cloudy days. This results in large areas of voids in the reconstructed point cloud, especially in the reconstruction of depressions on undulating and rough soil surfaces, which seriously compromises the accuracy of soil structure analysis. This application provides a method and system for analyzing and evaluating farmland soil structure.

[0007] Firstly, this application provides a method for analyzing and evaluating the structure of farmland soil, which adopts the following technical solution: including: The first point cloud data and the second point cloud data of the soil surface of the farmland to be tested are acquired by the imaging device. The first point cloud data has absolute coordinates and actual physical scale, and the second point cloud data contains color information. The first point cloud data and the second point cloud data are denoised, filtered and target region extracted respectively to obtain the first optimized point cloud and the second optimized point cloud. The voxel side length is adaptively estimated according to the preset downsampling rate. The second optimized point cloud is dimensionality reduced and compressed by voxel core downsampling to obtain the second optimized point cloud with reduced number of points but preserved geometric features. For the first optimized point cloud and the second optimized point cloud, coarse registration based on geometric feature matching and fine registration based on distance minimization are used to solve the rotation matrix and translation vector between the two point clouds. The second optimized point cloud is then transformed into the coordinate system of the first optimized point cloud to generate a fused point cloud that simultaneously possesses the detail accuracy of the first point cloud data and the integrity of the second point cloud data. Elevation values ​​of all points are extracted from the fused point cloud, and overlapping points are removed to form an elevation dataset. The root mean square height of the vertical undulation is calculated based on the elevation dataset, and the autocorrelation length of the horizontal spatial correlation length is calculated based on the autocorrelation function of the horizontal elevation sequence. Based on the root mean square height and autocorrelation length, and in accordance with the preset grading threshold, the farmland soil structure is divided into different roughness levels, generating structural analysis results.

[0008] Preferably, acquiring the first point cloud data and the second point cloud data of the farmland soil surface using the imaging device includes: A multi-view, high-overlapping sequence of two-dimensional images of the soil surface is acquired using a digital camera. A three-dimensional point cloud is reconstructed from the two-dimensional image sequence using a multi-view stereo matching algorithm based on motion recovery structure. The reconstructed three-dimensional point cloud is then scaled and georeferenced using the world coordinates of ground control points to obtain the first point cloud data with absolute coordinates and actual physical scale. A ToF camera is used to simultaneously acquire RGB and depth images of the soil surface. The depth image is mapped to the viewpoint of the RGB camera to generate an RGB-D image. Based on the camera intrinsic parameter model, the pixel coordinates of the RGB-D image are converted into three-dimensional world coordinates to obtain the second point cloud data containing color information.

[0009] Preferably, the step of denoising, filtering, and extracting target regions from the first point cloud data and the second point cloud data respectively to obtain a first optimized point cloud and a second optimized point cloud, and adaptively estimating the voxel side length according to a preset downsampling rate, and performing dimensionality reduction and compression on the second optimized point cloud through voxel core downsampling to obtain a second optimized point cloud with reduced point count but preserved geometric features, includes: KD-Tree spatial indexes are constructed for the first point cloud data and the second point cloud data respectively. For each point in the point cloud, a specified number of neighboring points are searched, the Euclidean distance from the target point to each neighboring point is calculated, and the average value is obtained. After traversing all points, the global mean and standard deviation of all average distances are calculated. Based on the global mean and the preset standard deviation multiple threshold, outliers that do not meet the interval requirements are identified and removed, resulting in the first intermediate point cloud and the second intermediate point cloud after outlier removal processing. Using the coordinate range of four preset ground control points in three-dimensional space as the boundary, a three-dimensional bounding box is constructed. Dimension threshold constraints are applied to the first intermediate point cloud and the second intermediate point cloud respectively. Point clouds that fall inside the three-dimensional bounding box are retained, and background and edge noise points are removed to obtain the first and second extracted point clouds that only contain the surface of the farmland soil to be tested. For each point in the first extracted point cloud and the second extracted point cloud, the Euclidean distance from each point in its neighborhood to the target point is used as the independent variable. The Gaussian weight of each neighboring point is calculated using a Gaussian function. The three-dimensional coordinates of each neighboring point are weighted and averaged. The position coordinates of the current point are updated with the weighted average value to obtain a first smooth point cloud and a second smooth point cloud with unchanged geometric features but smooth surface. The first smooth point cloud is set as the first optimized point cloud. For the second smoothed point cloud, input the preset target downsampling rate, calculate its bounding box volume and the total number of points in the point cloud, and estimate the initial voxel side length based on the target downsampling rate, bounding box volume and the total number of points in the point cloud; The three-dimensional space of the second smooth point cloud is divided into a voxel grid according to the initial voxel side length. The actual number of points contained in each non-empty voxel is counted, and the average density of the non-empty voxels is calculated. The average density is compared with the target number of points per voxel, where the target number of points per voxel is the reciprocal of the target downsampling rate. If the average density is less than the target number of voxels, the voxel side length is increased by a preset growth factor; if the average density is not less than the target number of voxels, the voxel side length is decreased by a growth factor; the process is iterated until the average density approaches the target number of voxels to obtain the optimal voxel side length. The second smooth point cloud is re-divided into voxel grids based on the optimal voxel side length. For each non-empty voxel, the centroid coordinates of all points within the target non-empty voxel are calculated. The centroid coordinates are used to replace all the original points within the target non-empty voxel to complete the dimensionality reduction and compression, resulting in the second optimized point cloud with a reduced number of points but preserved geometric features.

[0010] Preferably, the step of calculating the rotation matrix and translation vector between the first optimized point cloud and the second optimized point cloud based on coarse registration using geometric feature matching and fine registration based on distance minimization, and transforming the second optimized point cloud to the coordinate system of the first optimized point cloud to generate a fused point cloud that simultaneously possesses the detail accuracy of the first point cloud data and the integrity of the second point cloud data, includes: For the first optimized point cloud and the second optimized point cloud, the one with more points is used as the original point cloud and the other is used as the target point cloud; the elevation range of the two point clouds is obtained, and two elevation thresholds are calculated according to the first segmentation coefficient and the second segmentation coefficient. The regions with elevation values ​​less than the first threshold or greater than the second threshold are retained, and the flat regions in the middle are removed to obtain the original optimized point cloud and the target optimized point cloud after elevation screening. The normal vectors of each point in the filtered point cloud are estimated, a local coordinate system is established for each point, and the three angular features of the point and its neighboring points are calculated. The three angular features are divided into 11 statistical intervals and normalized, and then spliced ​​to generate a 33-dimensional fast point feature histogram descriptor. At least three points are randomly sampled from the original optimized point cloud. The corresponding points are searched in the target optimized point cloud based on the nearest neighbor matching of the FPFH descriptor. Polygons are constructed with the sampled points and the corresponding points respectively. The dissimilarity vectors of the corresponding side lengths are calculated. Each element is the absolute value of the side length difference divided by the larger value of the two corresponding side lengths. Based on the preset side length difference threshold, if any dissimilarity vector element is greater than the side length difference threshold, it is judged as a mismatch and resampling is performed. The corresponding point pairs that meet the threshold condition are retained. The initial rotation matrix and translation vector are solved by singular value decomposition to obtain the coarse registration transformation matrix and transform the target optimized point cloud. Based on the coarsely registered target optimized point cloud, a KD-Tree spatial index is constructed for the original optimized point cloud. For each point in the target optimized point cloud, the corresponding point with the closest Euclidean distance in the original optimized point cloud is searched using the KD-Tree. An optimization model is constructed with the weighted sum of the squared Euclidean distances between corresponding points as the objective function. The optimal rotation matrix and translation vector that minimize the objective function are solved iteratively. After each iteration, the point cloud coordinates are updated and the correspondence is re-established until the number of iterations reaches the maximum value or the error change is less than the convergence threshold, and the final transformation matrix is ​​obtained. The second optimized point cloud is transformed to the coordinate system of the first optimized point cloud according to the final transformation matrix to generate the fused point cloud; the bidirectional Hausdorff distance between the two point clouds is calculated, and if it is less than a preset threshold, the registration is deemed qualified and the fused point cloud is output; otherwise, the parameters are adjusted and the registration is repeated.

[0011] Preferably, the step of extracting the elevation values ​​of all points from the fused point cloud, removing overlapping points to form an elevation dataset; calculating the root mean square height of the vertical undulation based on the elevation dataset, and calculating the autocorrelation length of the horizontal spatial correlation length based on the autocorrelation function of the horizontal elevation sequence, includes: Traverse all points in the fused point cloud, extract the Z-axis coordinate value of each point as the elevation value, and construct the original elevation dataset. Using the spatial hash grid method, based on a preset distance threshold, points whose Euclidean distance is less than the distance threshold are identified as overlapping points and removed, retaining a single representative point to obtain the deduplicated elevation dataset; The total number of elevation points in the elevation dataset is counted, the arithmetic mean of all elevation values ​​is calculated, and the root mean square height is calculated. The root mean square height represents the degree of vertical undulation of the soil surface. The larger the value, the more uneven the surface, and the smaller the value, the flatter the surface. Elevation sequences are extracted from the fused point cloud along the horizontal X-axis at a preset sampling interval, or along the horizontal Y-axis at a preset sampling interval; a normalized autocorrelation function is calculated based on the elevation value and its displacement, where the displacement is equal to the sampling interval multiplied by a natural number index; the autocorrelation function values ​​corresponding to different displacements are calculated iteratively to generate an autocorrelation function curve; The autocorrelation function curve is smoothly fitted using the cubic spline interpolation method. The displacement value corresponding to the autocorrelation function value being equal to the reciprocal of the natural constant is calculated. The displacement value is determined as the autocorrelation length. The autocorrelation length characterizes the degree of spatial correlation of farmland soil surface in the horizontal direction. The larger the value, the gentler the surface undulation and the stronger the spatial continuity. The smaller the value, the more intense the surface undulation and the more significant the spatial heterogeneity. High-order sequences are extracted along the X-axis, Y-axis, and at a 45-degree angle to the X-axis, respectively. The autocorrelation length in each direction is calculated, and the arithmetic mean of the three is taken as the final autocorrelation length.

[0012] Preferably, the step of extracting high-order intervals from the fused point cloud along the horizontal X-axis at a preset sampling interval, or along the horizontal Y-axis at a preset sampling interval; calculating a normalized autocorrelation function based on the elevation value and its displacement, where the displacement is equal to the sampling interval multiplied by a natural number index; iterating through and calculating the autocorrelation function values ​​corresponding to different displacements; and generating an autocorrelation function curve includes: From the fused point cloud, discrete sampling is performed along the target horizontal direction at preset sampling intervals to obtain an elevation value sequence arranged in spatial order, denoted as the elevation sequence; the length of the elevation sequence is the total number of sampling points N. Define the sequence index displacement j, where j is a natural number and its value ranges from 1 to N-1; define the actual physical displacement x, where the actual physical displacement x is equal to the sequence index displacement j multiplied by the preset sampling interval; iterate through each integer value of the sequence index displacement j from 1 to N-1 and calculate the normalized autocorrelation function value corresponding to each displacement. For each sequence index displacement j, starting from the first position of the high-order sequence, traverse from index i from 1 to Nj, taking the i-th elevation value and the (i+j)-th elevation value in sequence, where i is a natural number from 1 to N, multiply the two, sum all the products, and use this as the numerator; sum the squares of all N elevation values ​​in the high-order sequence, and use this as the denominator; divide the numerator by the denominator to obtain the normalized autocorrelation function value corresponding to the sequence index displacement j. A discrete point set is constructed with the actual physical displacement x as the abscissa and the corresponding normalized autocorrelation function value as the ordinate. The discrete point set is then smoothly fitted using cubic spline interpolation to generate a continuous autocorrelation function curve. On the autocorrelation function curve, solve for the actual physical displacement value corresponding to the normalized autocorrelation function value being equal to the reciprocal of the natural constant; determine the corresponding actual physical displacement value as the autocorrelation length in the target horizontal direction.

[0013] Preferably, the method further includes: Based on the values ​​of root mean square height and autocorrelation length, and in accordance with the preset soil roughness classification threshold, the soil structure is classified into one of smooth, transitional, or rough types; a structural analysis report containing root mean square height, autocorrelation length, and soil structure evaluation level is generated and output.

[0014] Secondly, this application discloses a farmland soil structure analysis and evaluation device, which adopts the following technical solution, including: The point cloud imaging module is used to acquire first point cloud data and second point cloud data of the soil surface of the farmland to be measured through imaging equipment. The first point cloud data has absolute coordinates and actual physical scale, and the second point cloud data contains color information. The point cloud optimization module is used to denoise, filter and extract target regions from the first point cloud data and the second point cloud data respectively to obtain the first optimized point cloud and the second optimized point cloud. It also adaptively estimates the voxel side length according to the preset downsampling rate and performs dimensionality reduction and compression on the second optimized point cloud through voxel core downsampling to obtain the second optimized point cloud with reduced number of points but preserved geometric features. The point cloud fusion module is used to solve the rotation matrix and translation vector between the first optimized point cloud and the second optimized point cloud based on coarse registration of geometric feature matching and fine registration based on distance minimization, and transform the second optimized point cloud into the coordinate system of the first optimized point cloud to generate a fused point cloud that has both the detail accuracy of the first point cloud data and the integrity of the second point cloud data. The elevation calculation module is used to extract the elevation values ​​of all points from the fused point cloud, remove overlapping points to form an elevation dataset; calculate the root mean square height of the vertical undulation based on the elevation dataset, and calculate the autocorrelation length of the horizontal spatial correlation length based on the autocorrelation function of the horizontal elevation sequence. The structural analysis module is used to classify farmland soil structure into different roughness levels based on the root mean square height and autocorrelation length, and compare them with preset grading thresholds, and generate structural analysis results.

[0015] Thirdly, this application also provides a control device, the device comprising: It includes a memory and a processor, wherein the memory stores a computer program that can be loaded by the processor and executed, such as the above-described method for analyzing and evaluating farmland soil structure.

[0016] Fourthly, this application also provides a computer-readable storage medium storing a computer program that can be loaded by a processor and executed as described above for the analysis and evaluation of farmland soil structure.

[0017] In summary, this application employs a dual-source acquisition strategy to obtain a first point cloud and a second point cloud on the surface of the farmland soil to be tested. The first point cloud is reconstructed from a multi-view stereo matching algorithm based on structure-of-motion reconstructing, and its absolute geographic scale is assigned by ground control points. The second point cloud is based on RGB and depth images simultaneously acquired by a time-of-flight depth camera, which are then mapped to RGB-D images and converted into a dense point cloud containing color information. Subsequently, statistical filtering is applied to both point clouds to remove outliers, pass-through filtering is applied to extract the target region within the bounding box of the control points, and Gaussian filtering is applied to smooth surface noise. For the second point cloud, adaptive voxel downsampling is further introduced: the optimal voxel side length is dynamically estimated based on the preset downsampling rate, bounding box volume, and total number of points, and the original point set is replaced by the centroid within the voxel, achieving data compression while preserving geometric features. In the registration stage, firstly, smooth areas are segmented and removed based on elevation thresholds, while undulating feature points are retained to reduce computational redundancy. Fast point feature histogram descriptors are extracted, and a modified RANSAC algorithm is used to introduce side length difference constraints for coarse registration, obtaining an initial transformation matrix. Then, KDTree is used to accelerate nearest neighbor search, and fine registration is completed through an iterative nearest point algorithm, transforming the second point cloud to the coordinate system of the first point cloud, generating a fused point cloud that combines photogrammetric detail accuracy with depth camera integrity. Next, elevation values ​​are extracted from the fused point cloud, and after removing overlapping points, the root mean square height is calculated to represent vertical undulation. Elevation sequences are extracted at equal intervals along the horizontal direction, and the autocorrelation length is solved by fitting the sequence autocorrelation function to represent horizontal spatial correlation. Multiple directional means can be selected to characterize anisotropy. Finally, based on preset thresholds, the soil structure is divided into smooth, transitional, and rough levels, and an evaluation report is output. By leveraging the complementary advantages of heterogeneous point clouds, adaptive compression, and geometric constraint registration, the problem of a single sensor being unable to simultaneously achieve detailed accuracy and surface integrity was solved. This enabled high-precision three-dimensional quantitative analysis of farmland soil micro-topography, providing a data foundation for soil erosion assessment and the construction of high-standard farmland. Attached Figure Description

[0018] Figure 1 This is a flowchart illustrating a method for analyzing and evaluating the structure of farmland soil.

[0019] Figure 2 This is a structural block diagram of a farmland soil structure analysis and evaluation device. Detailed Implementation

[0020] The following combination Figures 1-2 This application will be described in further detail.

[0021] Reference Figure 1 The embodiments of this application include at least steps S10 to S50.

[0022] S10 acquires first and second point cloud data of the soil surface of the farmland to be tested through a camera. The first point cloud data has absolute coordinates and actual physical scale, while the second point cloud data contains color information.

[0023] S20, the first point cloud data and the second point cloud data are denoised, filtered and target region extracted respectively to obtain the first optimized point cloud and the second optimized point cloud. The voxel side length is adaptively estimated according to the preset downsampling rate. The second optimized point cloud is dimensionality reduced and compressed by voxel core downsampling to obtain the second optimized point cloud with reduced number of points but preserved geometric features.

[0024] S30: For the first optimized point cloud and the second optimized point cloud, based on coarse registration using geometric feature matching and fine registration based on distance minimization, solve the rotation matrix and translation vector between the two point clouds, and transform the second optimized point cloud to the coordinate system of the first optimized point cloud to generate a fused point cloud that simultaneously possesses the detail accuracy of the first point cloud data and the integrity of the second point cloud data.

[0025] S40: Extract the elevation values ​​of all points from the fused point cloud, remove overlapping points to form an elevation dataset; calculate the root mean square height of the vertical undulation based on the elevation dataset, and calculate the autocorrelation length of the horizontal spatial correlation length based on the autocorrelation function of the horizontal elevation sequence.

[0026] S50, based on the root mean square height and autocorrelation length, and in accordance with the preset grading threshold, classifies farmland soil structure into different roughness levels and generates structural analysis results.

[0027] Specifically, the system reconstructs a first point cloud with absolute geographic scale using a structure-on-motion motion reconstruction-multi-view stereo matching algorithm, while simultaneously acquiring a second point cloud containing color information using a time-of-flight depth camera, forming heterogeneous data input. Then, statistical filtering, pass-through filtering, and Gaussian filtering are applied to the two point clouds respectively to remove noise and extract target regions. Adaptive voxel downsampling is introduced into the second point cloud—the optimal voxel side length is dynamically and iteratively estimated based on a preset downsampling rate, point cloud bounding box volume, and total number of points, replacing the original point set with the voxel centroid, achieving data compression while preserving key geometric features. Based on this, using the first optimized point cloud as a benchmark, a fast point feature histogram is used... An improved RANSAC algorithm is used for coarse registration with side length constraints to obtain an initial transformation matrix. Then, the KDTree algorithm is used to accelerate the iterative nearest point algorithm to complete fine registration, transforming the second optimized point cloud to the coordinate system of the first optimized point cloud, generating a fused point cloud that combines photogrammetric detail accuracy with the complete coverage capability of a depth camera. Elevation values ​​are extracted from the fused point cloud and overlapping points are removed. The root mean square height is calculated to quantify the degree of vertical undulation. Equally spaced elevation sequences are extracted along the horizontal direction and fitted with an autocorrelation function. The autocorrelation length is solved to characterize the horizontal spatial correlation characteristics. Finally, the soil structure is classified into evaluation levels such as smooth, transitional, and rough according to the preset grading threshold, and an analysis report is output.

[0028] In some embodiments, step S10 specifically includes the following steps: using a digital camera to acquire a multi-view, high-overlapping sequence of two-dimensional images of the soil surface, reconstructing a three-dimensional point cloud from the two-dimensional image sequence using a multi-view stereo matching algorithm based on motion recovery structure, and scaling and georeferencing the reconstructed three-dimensional point cloud using the world coordinates of ground control points to obtain first point cloud data with absolute coordinates and actual physical scale. A ToF camera is used to simultaneously acquire RGB and depth images of the soil surface. The depth image is then mapped to the viewpoint of the RGB camera to generate an RGB-D image. Based on the camera intrinsic parameter model, the pixel coordinates of the RGB-D image are converted into three-dimensional world coordinates to obtain second point cloud data containing color information.

[0029] Specifically, a multi-view, high-overlap image sequence was acquired using a digital camera. A 3D point cloud was reconstructed based on a structure-of-motion (COMO) multi-view stereo matching algorithm. Scale scaling and georeferenced registration were performed using world coordinates of ground control points to obtain a first point cloud with absolute coordinates and actual physical scale, ensuring geometric accuracy and traceability of the true scale. Secondly, RGB and depth images were simultaneously acquired using a time-of-flight depth camera. The depth image was mapped to the RGB camera viewpoint to generate an RGB-D image. Based on a camera intrinsic model, pixel coordinates were converted to 3D world coordinates to obtain a second point cloud containing color information, achieving rapid and complete surface coverage in weakly textured environments. The two data streams output high-precision geometric point clouds and high-completeness color point clouds, forming a complementary heterogeneous point cloud pair. This provides a data foundation with significantly different features for subsequent adaptive preprocessing and rigid registration fusion.

[0030] In some embodiments, step S20 specifically includes the following steps: constructing KD-Tree spatial indexes for the first point cloud data and the second point cloud data respectively; for each point in the point cloud, searching a specified number of neighboring points, calculating the Euclidean distance from the target point to each neighboring point and obtaining the average value; after traversing all points, calculating the global mean and standard deviation of all average distances; based on the global mean and a preset standard deviation multiple threshold, identifying and removing outliers that do not meet the interval requirements, obtaining the first intermediate point cloud and the second intermediate point cloud after outlier removal processing; constructing a three-dimensional bounding box with the coordinate range of four preset ground control points in three-dimensional space as the boundary, applying dimensional threshold constraints to the first intermediate point cloud and the second intermediate point cloud respectively, retaining the point cloud falling inside the three-dimensional bounding box, removing background and edge noise points, obtaining the first extracted point cloud and the second extracted point cloud that only contain the surface of the farmland soil to be tested. For each point in the first and second extracted point clouds, the Euclidean distance from each point in its neighborhood to the target point is used as the independent variable. The Gaussian weight of each neighboring point is calculated using a Gaussian function. The 3D coordinates of each neighboring point are then weighted and averaged. The position coordinates of the current point are updated using the weighted average value. This results in a first smooth point cloud and a second smooth point cloud with unchanged geometric features but smooth surface. The first smooth point cloud is set as the first optimized point cloud. For the second smooth point cloud, a preset target downsampling rate is input, and its bounding box volume and the total number of points in the point cloud are calculated. The initial voxel side length is estimated based on the target downsampling rate, bounding box volume, and the total number of points in the point cloud. The 3D space of the second smooth point cloud is divided into a voxel grid according to the initial voxel edge length. The actual number of points contained in each non-empty voxel is counted, and the average density of the non-empty voxels is calculated. The average density is compared with the target number of points per voxel, which is the reciprocal of the target downsampling rate. If the average density is less than the target number of points per voxel, the voxel edge length is increased by a preset growth factor. If the average density is not less than the target number of points per voxel, the voxel edge length is decreased by the growth factor. The adjustment is iterated until the average density approaches the target number of points per voxel, and the optimal voxel edge length is obtained. The second smooth point cloud is re-divided into a voxel grid based on the optimal voxel edge length. For each non-empty voxel, the centroid coordinates of all points in the target non-empty voxel are calculated. The centroid coordinates are used to replace all the original points in the target non-empty voxel to complete the dimensionality reduction and compression, resulting in a second optimized point cloud with a reduced number of points but preserved geometric features.

[0031] For the second smoothed point cloud, input the preset target downsampling rate, calculate its bounding box volume and the total number of points in the point cloud, and estimate the initial voxel side length based on the target downsampling rate, bounding box volume, and total number of points in the point cloud, as shown in the following formula: ; in, It is the initial voxel side length estimate, representing the side length of each cube cell when the voxel mesh is first generated; is the preset target downsampling rate, a dimensionless pure number ranging from 0.01 to 0.5, representing the ratio of the number of points in the downsampled point cloud to the number of points in the original point cloud. v is the bounding box volume of the second smoothed point cloud, in cubic millimeters, the volume of the cuboid bounded by the maximum and minimum values ​​of the point cloud in the X, Y, and Z directions. N is the total number of points in the second smoothed point cloud.

[0032] Specifically, the system constructs a KD-Tree index through statistical filtering, eliminating outliers based on the neighborhood average distance and global standard deviation to effectively suppress scattered noise. Then, a 3D bounding box is constructed using ground control points as boundaries for direct filtering, accurately extracting the target soil region and removing background and edge interference. Subsequently, Gaussian filtering is used to perform weighted average smoothing of the point cloud coordinates, eliminating surface burrs while preserving geometric features. Addressing the large data volume of the second point cloud, adaptive voxel downsampling is further introduced: the initial voxel side length is dynamically estimated based on the target downsampling rate, bounding box volume, and total number of points. Through iterative adjustment, the average density of non-empty voxels approximates the target number per voxel, and the original point set is replaced by the voxel centroid, achieving point cloud dimensionality reduction and compression while maximizing the preservation of local geometric features. Thus, the smoothed first point cloud is directly set as the first optimized point cloud, and the second point cloud, after multi-level filtering and adaptive compression, is output as the second optimized point cloud, providing lightweight data input with good feature preservation for subsequent registration of heterogeneous point clouds.

[0033] In some embodiments, step S30 specifically includes the following steps: For the first optimized point cloud and the second optimized point cloud, the one with more points is used as the original point cloud, and the other is used as the target point cloud; the elevation range of the two point clouds is obtained, and two elevation thresholds are calculated according to the first segmentation coefficient and the second segmentation coefficient. Regions with elevation values ​​less than the first threshold or greater than the second threshold are retained, and the flat areas in the middle are removed to obtain the original optimized point cloud and the target optimized point cloud after elevation screening; the normal vector of each point in the screened point cloud is estimated, a local coordinate system is established for each point, and three angular features between the point and its neighboring points are calculated. The three angular features are divided into 11 statistical intervals and normalized. A 33-dimensional fast point feature histogram descriptor is generated by splicing the data. At least three points are randomly sampled from the original optimized point cloud, and corresponding points are searched in the target optimized point cloud based on nearest neighbor matching of the FPFH descriptor. Polygons are constructed using the sampled points and corresponding points respectively, and the dissimilarity vectors of corresponding side lengths are calculated. Each element of the vector is the absolute value of the side length difference divided by the larger of the two corresponding side lengths. Based on a preset side length difference threshold, if any dissimilarity vector element is greater than the side length difference threshold, it is determined to be a mismatch and resampling is performed. The corresponding point pairs that meet the threshold condition are retained, and the initial rotation matrix and translation vector are solved by singular value decomposition to obtain the coarse registration transformation matrix and transform the target optimized point cloud. Based on the coarsely registered target optimized point cloud, a KD-Tree spatial index is constructed for the original optimized point cloud. For each point in the target optimized point cloud, the KD-Tree searches for the closest Euclidean distance corresponding point in the original optimized point cloud. An optimization model is constructed with the weighted sum of the squared Euclidean distances between corresponding points as the objective function. The optimal rotation matrix and translation vector that minimize the objective function are iteratively solved. After each iteration, the point cloud coordinates are updated and the correspondence is re-established until the number of iterations reaches the maximum value or the error change is less than the convergence threshold, thus obtaining the final transformation matrix. The second optimized point cloud is transformed to the coordinate system of the first optimized point cloud according to the final transformation matrix to generate a fused point cloud. The bidirectional Hausdorff distance between the two point clouds is calculated. If it is less than a preset threshold, the registration is deemed qualified and the fused point cloud is output; otherwise, the parameters are adjusted and re-registered.

[0034] The formula for calculating the centroid coordinates of all points within the non-empty voxel of the target is as follows: ; in, It is the centroid coordinate of the k-th non-empty voxel, a three-dimensional vector used to replace all the original points within that voxel. This represents the number of original point clouds contained within the k-th non-empty voxel, expressed in units of points. This represents the three-dimensional coordinates of the j-th point.

[0035] Specifically, based on the elevation undulation characteristics of the soil surface, the system uses dual-threshold segmentation to eliminate smooth areas and retain significant feature points, effectively reducing registration computational redundancy and improving feature recognition. Subsequently, it estimates the normal vectors of the filtered point cloud and extracts 33-dimensional fast point feature histogram descriptors, establishing robust local feature representations for geometric matching. In the coarse registration stage, an improved RANSAC algorithm is used, introducing side length difference constraints—a polygon is constructed using sampled points and corresponding points, the vectors of dissimilarity of corresponding side lengths are calculated and compared with a preset threshold, retaining only matching point pairs that satisfy rigid geometric consistency. The initial transformation matrix is ​​solved through singular value decomposition, achieving robust initial alignment under large pose differences. In the fine registration stage, using the coarse registration result as initial values, a KD-Tree spatial index is constructed on the original optimized point cloud to accelerate nearest neighbor search. The iterative nearest point algorithm optimizes the weighted sum of the squared Euclidean distances between corresponding point pairs, iteratively solving for the optimal rotation matrix and translation vector until convergence. Finally, the second optimized point cloud is transformed to the coordinate system of the first optimized point cloud to generate a fused point cloud. The registration quality is quantitatively evaluated using bidirectional Hausdorff distance, forming a closed loop of adaptive parameter adjustment. Thus, through the joint collaboration of geometric feature selection, rigid constraint coarse registration, and fine registration, high-precision fusion of heterogeneous point clouds is achieved.

[0036] In some embodiments, step S40 specifically includes the following steps: traversing all points in the fused point cloud, extracting the Z-axis coordinate value of each point as the elevation value, and constructing an original elevation dataset; using the spatial hash grid method, based on a preset distance threshold, determining points with Euclidean distance less than the distance threshold as overlapping points and removing them, retaining a single representative point, and obtaining a deduplicated elevation dataset; counting the total number of elevation points in the elevation dataset, calculating the arithmetic mean of all elevation values, and calculating the root mean square height, which characterizes the degree of undulation of the soil surface in the vertical direction. The larger the value, the more uneven the surface, and the smaller the value, the flatter the surface; extracting elevation sequence from the fused point cloud along the horizontal X-axis direction at a preset sampling interval, or extracting elevation sequence along the horizontal Y-axis direction at a preset sampling interval; calculating the normalized autocorrelation function based on the elevation value and its displacement, where the displacement is equal to the sampling interval multiplied by the natural number index, traversing and calculating the autocorrelation function values ​​corresponding to different displacements, and generating an autocorrelation function curve; The autocorrelation function curve was smoothly fitted using cubic spline interpolation. The displacement value corresponding to the autocorrelation function value being equal to the reciprocal of the natural constant was calculated, and the displacement value was determined as the autocorrelation length. The autocorrelation length characterizes the spatial correlation of farmland soil surface in the horizontal direction. The larger the value, the gentler the surface undulation and the stronger the spatial continuity. The smaller the value, the more severe the surface undulation and the more significant the spatial heterogeneity. High-order sequences were extracted along the X-axis, Y-axis and the direction at a 45-degree angle to the X-axis, and the autocorrelation length in each direction was calculated. The arithmetic mean of the three was taken as the final autocorrelation length.

[0037] The root mean square height is calculated using the following formula: ; Wherein, RMSH represents the root mean square height, which characterizes the degree of vertical undulation of the soil surface; a larger value indicates a more uneven surface. N is the total number of elevation points in the deduplicated elevation dataset. This represents the elevation value (Z-axis coordinate) of the i-th point, in millimeters. This represents the arithmetic mean of all elevation values.

[0038] Specifically, the system traverses and merges point clouds to extract Z-axis elevation values. It then uses a spatial hash grid method to remove overlapping points with Euclidean distances less than a threshold, forming a deduplicated elevation dataset to avoid bias in statistical results caused by repeated sampling. Next, the total number of statistical points and their arithmetic mean are used to calculate the root mean square height, quantifying the intensity of vertical undulations in the soil surface. For horizontal characterization, elevation sequences are extracted along the X-axis, Y-axis, and 45-degree direction at preset sampling intervals. A normalized autocorrelation function is calculated based on the sequence index displacement, and autocorrelation curves are generated by traversing each displacement. Cubic spline interpolation is used to smoothly fit discrete points, and the actual displacement corresponding to the decay of the autocorrelation function value to the reciprocal of the natural constant is determined as the autocorrelation length, used to characterize the spatial correlation and surface smoothness in the horizontal direction. Finally, the arithmetic mean of the multi-directional autocorrelation lengths is fused to effectively characterize the anisotropic features of soil roughness. This achieves a lossless transformation from merged point clouds to vertical and horizontal dual-dimensional statistical parameters, providing a quantitative basis for soil structure classification and evaluation.

[0039] Furthermore, considering the specific calculation process of the autocorrelation length, the corresponding processing steps are as follows: Discretely sample the fused point cloud along the target horizontal direction at preset sampling intervals to obtain a sequence of elevation values ​​arranged in spatial order, denoted as the elevation sequence; the length of the elevation sequence is the total number of sampling points N; define the sequence index displacement j, where j is a natural number ranging from 1 to N-1; define the actual physical displacement x, which is equal to the sequence index displacement j multiplied by the preset sampling interval; traverse the sequence index displacement j from 1 to... For each integer value of N-1, calculate the normalized autocorrelation function value corresponding to each displacement. For each sequence index displacement j, starting from the first position of the high-order sequence, traverse from index i from 1 to Nj, taking the i-th elevation value and the (i+j)-th elevation value in sequence, where i is a natural number from 1 to N, multiply the two, and sum all the products as the numerator. Sum the squares of all N elevation values ​​in the high-order sequence as the denominator. Divide the numerator by the denominator to obtain the normalized autocorrelation function value corresponding to the sequence index displacement j. A discrete point set is constructed with the actual physical displacement x as the abscissa and the corresponding normalized autocorrelation function value as the ordinate. The discrete point set is then smoothly fitted using cubic spline interpolation to generate a continuous autocorrelation function curve. On the autocorrelation function curve, the actual physical displacement value corresponding to the normalized autocorrelation function value being equal to the reciprocal of the natural constant is calculated. The corresponding actual physical displacement value is then determined as the autocorrelation length in the target horizontal direction.

[0040] For each sequence index shift j, the normalized autocorrelation function value is calculated as follows: ; in, represents the normalized autocorrelation function value when the sequence index shift is j. It is a dimensionless pure number with a value range of [-1, 1], representing the degree of autocorrelation at the spatial interval j×Δx. N represents the length of the high-order sequence (total number of sampling points). This represents the elevation value of the i-th position in the high-order sequence.

[0041] Specifically, the system performs equally spaced discrete sampling of the fused point cloud along the target horizontal direction at a preset sampling interval, generating a high-order sequence arranged in spatial order with a sequence length of N. Then, it defines the sequence index displacement j (from 1 to N-1) and the actual physical displacement x (j × sampling interval), establishing a mapping relationship between the displacement and the autocorrelation function value. For each displacement j, it iterates from index i to Nj, multiplies the i-th and (i+j)-th elevation values ​​and sums them as the numerator, uses the sum of squares of all elevation values ​​as the denominator, and divides to obtain the normalized autocorrelation function value corresponding to that displacement. After iterating through all j values, a discrete point set is constructed with the actual physical displacement x as the abscissa and the autocorrelation function value as the ordinate. Cubic spline interpolation is used to smoothly fit the discrete point set, generating a continuous autocorrelation function curve. The actual physical displacement corresponding to the autocorrelation function value decaying to the reciprocal of the natural constant on the curve is calculated, and this is determined as the autocorrelation length in that direction. This step transforms the calculation of the autocorrelation function from an abstract formula into an executable sequence operation. Through equal-interval sampling, sliding window accumulation, interpolation fitting, and threshold solving, it achieves a high-precision quantitative characterization of the horizontal spatial correlation of the soil surface.

[0042] In some embodiments, step S50 specifically includes the following steps: based on the values ​​of root mean square height and autocorrelation length, and in accordance with a preset soil roughness grading threshold, classify the soil structure into one of smooth, transitional, or rough types; generate and output a structural analysis report containing root mean square height, autocorrelation length, and soil structure evaluation level.

[0043] The implementation principle of the farmland soil structure analysis and evaluation method in this application is as follows: A first point cloud and a second point cloud of the farmland soil surface to be tested are acquired through a dual-source acquisition strategy. The first point cloud is reconstructed from a multi-view image sequence of a digital camera based on a motion recovery structure-multi-view stereo matching algorithm and assigned an absolute geographic scale by ground control points. The second point cloud is based on RGB and depth images simultaneously acquired by a time-of-flight depth camera, which are then mapped to RGB-D images and converted into a dense point cloud containing color information. Subsequently, statistical filtering is performed on both point clouds to remove outliers, pass-through filtering is used to extract the target area within the bounding box of the control points, and Gaussian filtering is used to smooth surface noise. For the second point cloud, adaptive voxel downsampling is further introduced: the optimal voxel side length is dynamically estimated based on a preset downsampling rate, bounding box volume, and total number of points, and the original point set is replaced by the centroid within the voxel, achieving data compression while preserving geometric features. In the registration stage, firstly, smooth areas are segmented and removed based on elevation thresholds, while undulating feature points are retained to reduce computational redundancy. Fast point feature histogram descriptors are extracted, and a modified RANSAC algorithm is used to introduce side length difference constraints for coarse registration, obtaining an initial transformation matrix. Then, KDTree is used to accelerate nearest neighbor search, and fine registration is completed through an iterative nearest point algorithm, transforming the second point cloud to the coordinate system of the first point cloud, generating a fused point cloud that combines photogrammetric detail accuracy with depth camera integrity. Next, elevation values ​​are extracted from the fused point cloud, and after removing overlapping points, the root mean square height is calculated to represent vertical undulation. Elevation sequences are extracted at equal intervals along the horizontal direction, and the autocorrelation length is solved by fitting the sequence autocorrelation function to represent horizontal spatial correlation. Multiple directional means can be selected to characterize anisotropy. Finally, based on preset thresholds, the soil structure is divided into smooth, transitional, and rough levels, and an evaluation report is output. By leveraging the complementary advantages of heterogeneous point clouds, adaptive compression, and geometric constraint registration, the problem of a single sensor being unable to simultaneously achieve detailed accuracy and surface integrity was solved. This enabled high-precision three-dimensional quantitative analysis of farmland soil micro-topography, providing a data foundation for soil erosion assessment and the construction of high-standard farmland.

[0044] Figure 1 This is a flowchart illustrating a method for analyzing and evaluating farmland soil structure in one embodiment. It should be understood that, although... Figure 1 The steps in the flowchart are shown sequentially as indicated by the arrows, but these steps are not necessarily executed in the order indicated by the arrows; unless explicitly stated herein, there is no strict order requirement for the execution of these steps, and they can be executed in other orders; and Figure 1At least some of the steps in the process may include multiple sub-steps or multiple stages. These sub-steps or stages are not necessarily completed at the same time, but can be executed at different times. The execution order of these sub-steps or stages is not necessarily sequential, but can be executed in turn or alternately with other steps or at least some of the sub-steps or stages of other steps.

[0045] Based on the same technical concept, referring to Figure 2 This application also provides a farmland soil structure analysis and evaluation device, which adopts the following technical solution: the device includes: The point cloud imaging module is used to acquire first and second point cloud data of the soil surface of the farmland to be tested through imaging equipment. The first point cloud data has absolute coordinates and actual physical scale, and the second point cloud data contains color information. The point cloud optimization module is used to denoise, filter and extract target regions from the first point cloud data and the second point cloud data respectively to obtain the first optimized point cloud and the second optimized point cloud. It also adaptively estimates the voxel side length according to the preset downsampling rate, and performs dimensionality reduction and compression on the second optimized point cloud through voxel core downsampling to obtain a second optimized point cloud with reduced number of points but preserved geometric features. The point cloud fusion module is used to solve the rotation matrix and translation vector between the first optimized point cloud and the second optimized point cloud based on coarse registration of geometric feature matching and fine registration based on distance minimization. It also transforms the second optimized point cloud into the coordinate system of the first optimized point cloud to generate a fused point cloud that has both the detail accuracy of the first point cloud data and the integrity of the second point cloud data. The elevation calculation module is used to extract the elevation values ​​of all points from the fused point cloud, remove overlapping points to form an elevation dataset; it calculates the root mean square height of the vertical undulation based on the elevation dataset, and calculates the autocorrelation length of the horizontal spatial correlation length based on the autocorrelation function of the horizontal elevation sequence. The structural analysis module is used to classify farmland soil structure into different roughness levels based on the root mean square height and autocorrelation length, and compare them with preset grading thresholds, and generate structural analysis results.

[0046] In some embodiments, the point cloud shooting module is specifically used to acquire a multi-view, high-overlapping sequence of two-dimensional images of the soil surface using a digital camera, reconstruct a three-dimensional point cloud from the two-dimensional image sequence using a multi-view stereo matching algorithm based on motion recovery structure, and scale and georeference register the reconstructed three-dimensional point cloud using the world coordinates of ground control points to obtain first point cloud data with absolute coordinates and actual physical scale. A ToF camera is used to simultaneously acquire RGB and depth images of the soil surface. The depth image is then mapped to the viewpoint of the RGB camera to generate an RGB-D image. Based on the camera intrinsic parameter model, the pixel coordinates of the RGB-D image are converted into three-dimensional world coordinates to obtain second point cloud data containing color information.

[0047] In some embodiments, the point cloud optimization module is specifically used to construct KD-Tree spatial indexes for the first point cloud data and the second point cloud data respectively. For each point in the point cloud, a specified number of neighboring points are searched, the Euclidean distance from the target point to each neighboring point is calculated and the average value is obtained. After traversing all points, the global mean and standard deviation of all average distances are calculated. Based on the global mean and the preset standard deviation multiple threshold, outliers that do not meet the interval requirements are identified and removed, resulting in the first intermediate point cloud and the second intermediate point cloud after outlier removal processing. Using the coordinate range of four preset ground control points in three-dimensional space as the boundary, a three-dimensional bounding box is constructed. Dimension threshold constraints are applied to the first intermediate point cloud and the second intermediate point cloud respectively. The point cloud that falls inside the three-dimensional bounding box is retained, and background and edge noise points are removed to obtain the first and second extracted point clouds that only contain the surface of the farmland soil to be tested. For each point in the first extracted point cloud and the second extracted point cloud, the Euclidean distance from each point in its neighborhood to the target point is used as the independent variable. The Gaussian weight of each neighboring point is calculated using a Gaussian function. The three-dimensional coordinates of each neighboring point are weighted and averaged. The position coordinates of the current point are updated with the weighted average value to obtain the first smooth point cloud and the second smooth point cloud with unchanged geometric features but smooth surface. The first smooth point cloud is set as the first optimized point cloud. For the second smooth point cloud, input the preset target downsampling rate, calculate its bounding box volume and the total number of points in the point cloud, and estimate the initial voxel side length based on the target downsampling rate, bounding box volume and the total number of points in the point cloud; The three-dimensional space of the second smooth point cloud is divided into a voxel grid according to the initial voxel side length. The actual number of points contained in each non-empty voxel is counted, and the average density of the non-empty voxels is calculated. The average density is compared with the target number of points per voxel, where the target number of points per voxel is the reciprocal of the target downsampling rate. If the average density is less than the target number of voxels, the voxel side length is increased by a preset growth factor; if the average density is not less than the target number of voxels, the voxel side length is decreased by a growth factor; the process is iterated until the average density approaches the target number of voxels to obtain the optimal voxel side length. The second smooth point cloud is re-divided into voxel grids based on the optimal voxel side length. For each non-empty voxel, the centroid coordinates of all points in the target non-empty voxel are calculated. The centroid coordinates are used to replace all the original points in the target non-empty voxel to complete the dimensionality reduction and compression, resulting in a second optimized point cloud with a reduced number of points but preserved geometric features.

[0048] In some embodiments, the point cloud fusion module is specifically used to select the first optimized point cloud and the second optimized point cloud as the original point cloud and the other as the target point cloud; obtain the elevation range of the two point clouds; calculate two elevation thresholds according to the first segmentation coefficient and the second segmentation coefficient; retain the region where the elevation value is less than the first threshold or greater than the second threshold; remove the middle flat region; and obtain the original optimized point cloud and the target optimized point cloud after elevation screening. The normal vectors of each point in the filtered point cloud are estimated, a local coordinate system is established for each point, and the three angular features of the point and its neighboring points are calculated. The three angular features are divided into 11 statistical intervals and normalized, and then spliced ​​to generate a 33-dimensional fast point feature histogram descriptor. At least three points are randomly sampled from the original optimized point cloud. The corresponding points are searched in the target optimized point cloud based on the nearest neighbor matching of the FPFH descriptor. Polygons are constructed with the sampled points and the corresponding points respectively. The dissimilarity vectors of the corresponding side lengths are calculated. Each element is the absolute value of the side length difference divided by the larger value of the two corresponding side lengths. Based on the preset side length difference threshold, if any dissimilarity vector element is greater than the side length difference threshold, it is judged as a mismatch and resampling is performed. The corresponding point pairs that meet the threshold condition are retained. The initial rotation matrix and translation vector are solved by singular value decomposition to obtain the coarse registration transformation matrix and transform the target optimized point cloud. Based on the coarsely registered target optimized point cloud, a KD-Tree spatial index is constructed for the original optimized point cloud. For each point in the target optimized point cloud, the corresponding point with the closest Euclidean distance in the original optimized point cloud is searched using the KD-Tree. An optimization model is constructed with the weighted sum of the squared Euclidean distances between corresponding points as the objective function. The optimal rotation matrix and translation vector that minimize the objective function are solved iteratively. After each iteration, the point cloud coordinates are updated and the correspondence is re-established until the number of iterations reaches the maximum value or the error change is less than the convergence threshold, and the final transformation matrix is ​​obtained. The second optimized point cloud is transformed to the coordinate system of the first optimized point cloud based on the final transformation matrix to generate a fused point cloud. The bidirectional Hausdorff distance between the two point clouds is calculated. If it is less than a preset threshold, the registration is deemed qualified and the fused point cloud is output. Otherwise, the parameters are adjusted and the registration is repeated.

[0049] In some embodiments, the elevation calculation module is specifically used to traverse all points in the fused point cloud, extract the Z-axis coordinate value of each point as the elevation value, and construct the original elevation dataset. Using the spatial hash grid method, based on a preset distance threshold, points whose Euclidean distance is less than the distance threshold are identified as overlapping points and removed, while retaining a single representative point to obtain a deduplicated elevation dataset. The total number of elevation points in the elevation dataset is counted, the arithmetic mean of all elevation values ​​is calculated, and the root mean square height is calculated. The root mean square height characterizes the degree of vertical undulation of the soil surface. The larger the value, the more uneven the surface, and the smaller the value, the flatter the surface. Elevation sequences are extracted from the fused point cloud along the horizontal X-axis at a preset sampling interval, or along the horizontal Y-axis at a preset sampling interval; a normalized autocorrelation function is calculated based on the elevation value and its displacement, where the displacement is equal to the sampling interval multiplied by a natural number index; the autocorrelation function values ​​corresponding to different displacements are calculated iteratively to generate an autocorrelation function curve; The autocorrelation function curve was smoothly fitted using the cubic spline interpolation method. The displacement value corresponding to the autocorrelation function value being equal to the reciprocal of the natural constant was calculated. The displacement value was determined as the autocorrelation length. The autocorrelation length characterizes the degree of spatial correlation of farmland soil surface in the horizontal direction. The larger the value, the gentler the surface undulation and the stronger the spatial continuity. The smaller the value, the more intense the surface undulation and the more significant the spatial heterogeneity. High-order sequences are extracted along the X-axis, Y-axis, and at a 45-degree angle to the X-axis, respectively. The autocorrelation length in each direction is calculated, and the arithmetic mean of the three is taken as the final autocorrelation length.

[0050] In some embodiments, the elevation calculation module is specifically used to perform equally spaced discrete sampling from the fused point cloud along the target horizontal direction at a preset sampling interval to obtain an elevation value sequence arranged in spatial order, denoted as the elevation sequence; the length of the elevation sequence is the total number of sampling points N. Define the sequence index displacement j, where j is a natural number and its value ranges from 1 to N-1; define the actual physical displacement x, where the actual physical displacement x is equal to the sequence index displacement j multiplied by the preset sampling interval; iterate through each integer value of the sequence index displacement j from 1 to N-1 and calculate the normalized autocorrelation function value corresponding to each displacement. For each sequence index displacement j, starting from the first position of the high-order sequence, traverse from index i from 1 to Nj, taking the i-th elevation value and the (i+j)-th elevation value in sequence, where i is a natural number from 1 to N, multiply the two, sum all the products, and use this as the numerator; sum the squares of all N elevation values ​​in the high-order sequence, and use this as the denominator; divide the numerator by the denominator to obtain the normalized autocorrelation function value corresponding to the sequence index displacement j. A discrete point set is constructed with the actual physical displacement x as the abscissa and the corresponding normalized autocorrelation function value as the ordinate. The discrete point set is then smoothly fitted using cubic spline interpolation to generate a continuous autocorrelation function curve. On the autocorrelation function curve, solve for the actual physical displacement value corresponding to the normalized autocorrelation function value being equal to the reciprocal of the natural constant; determine the corresponding actual physical displacement value as the autocorrelation length in the target horizontal direction.

[0051] In some embodiments, the structural analysis module is specifically used to classify the soil structure into one of smooth, transitional, or rough types based on the values ​​of the root mean square height and autocorrelation length, and in accordance with a preset soil roughness grading threshold; and to generate and output a structural analysis report containing the root mean square height, autocorrelation length, and soil structure evaluation level.

[0052] This application also discloses a control device.

[0053] Specifically, the control device includes a memory and a processor, with the memory storing a computer program that can be loaded by the processor and executed using the aforementioned farmland soil structure analysis and evaluation method.

[0054] This application also discloses a computer-readable storage medium.

[0055] Specifically, the computer-readable storage medium stores a computer program that can be loaded by a processor and executed, such as the above-described method for analyzing and evaluating farmland soil structure. The computer-readable storage medium includes, for example, various media capable of storing program code, such as a USB flash drive, a portable hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk.

[0056] The above are all preferred embodiments of this application, and are not intended to limit the scope of protection of this application. Therefore, all equivalent changes made in accordance with the structure, shape and principle of this application should be covered within the scope of protection of this application.

Claims

1. A method for analyzing and evaluating the structure of farmland soil, characterized in that, include: The first point cloud data and the second point cloud data of the soil surface of the farmland to be tested are acquired by the imaging device. The first point cloud data has absolute coordinates and actual physical scale, and the second point cloud data contains color information. The first point cloud data and the second point cloud data are denoised, filtered and target region extracted respectively to obtain the first optimized point cloud and the second optimized point cloud. The voxel side length is adaptively estimated according to the preset downsampling rate. The second optimized point cloud is dimensionality reduced and compressed by voxel core downsampling to obtain the second optimized point cloud with reduced number of points but preserved geometric features. For the first optimized point cloud and the second optimized point cloud, coarse registration based on geometric feature matching and fine registration based on distance minimization are used to solve the rotation matrix and translation vector between the two point clouds. The second optimized point cloud is then transformed into the coordinate system of the first optimized point cloud to generate a fused point cloud that simultaneously possesses the detail accuracy of the first point cloud data and the integrity of the second point cloud data. Elevation values ​​of all points are extracted from the fused point cloud, and overlapping points are removed to form an elevation dataset. The root mean square height of the vertical undulation is calculated based on the elevation dataset, and the autocorrelation length of the horizontal spatial correlation length is calculated based on the autocorrelation function of the horizontal elevation sequence. Based on the root mean square height and autocorrelation length, and in accordance with the preset grading threshold, the farmland soil structure is divided into different roughness levels, generating structural analysis results.

2. The method for analyzing and evaluating farmland soil structure according to claim 1, characterized in that, The acquisition of the first and second point cloud data of the soil surface of the farmland to be tested through the imaging device includes: A multi-view, high-overlapping sequence of two-dimensional images of the soil surface is acquired using a digital camera. A three-dimensional point cloud is reconstructed from the two-dimensional image sequence using a multi-view stereo matching algorithm based on motion recovery structure. The reconstructed three-dimensional point cloud is then scaled and georeferenced using the world coordinates of ground control points to obtain the first point cloud data with absolute coordinates and actual physical scale. A ToF camera is used to simultaneously acquire RGB and depth images of the soil surface. The depth image is mapped to the viewpoint of the RGB camera to generate an RGB-D image. Based on the camera intrinsic parameter model, the pixel coordinates of the RGB-D image are converted into three-dimensional world coordinates to obtain the second point cloud data containing color information.

3. The method for analyzing and evaluating farmland soil structure according to claim 1, characterized in that, The process of denoising, filtering, and extracting target regions from the first and second point cloud data respectively to obtain a first optimized point cloud and a second optimized point cloud, adaptively estimating the voxel side length according to a preset downsampling rate, and performing dimensionality reduction and compression on the second optimized point cloud through voxel core downsampling to obtain a second optimized point cloud with reduced point count but preserved geometric features, includes: KD-Tree spatial indexes are constructed for the first point cloud data and the second point cloud data respectively. For each point in the point cloud, a specified number of neighboring points are searched, the Euclidean distance from the target point to each neighboring point is calculated, and the average value is obtained. After traversing all points, the global mean and standard deviation of all average distances are calculated. Based on the global mean and the preset standard deviation multiple threshold, outliers that do not meet the interval requirements are identified and removed, resulting in the first intermediate point cloud and the second intermediate point cloud after outlier removal processing. Using the coordinate range of four preset ground control points in three-dimensional space as the boundary, a three-dimensional bounding box is constructed. Dimension threshold constraints are applied to the first intermediate point cloud and the second intermediate point cloud respectively. Point clouds that fall inside the three-dimensional bounding box are retained, and background and edge noise points are removed to obtain the first and second extracted point clouds that only contain the surface of the farmland soil to be tested. For each point in the first extracted point cloud and the second extracted point cloud, the Euclidean distance from each point in its neighborhood to the target point is used as the independent variable. The Gaussian weight of each neighboring point is calculated using a Gaussian function. The three-dimensional coordinates of each neighboring point are weighted and averaged. The position coordinates of the current point are updated with the weighted average value to obtain a first smooth point cloud and a second smooth point cloud with unchanged geometric features but smooth surface. The first smooth point cloud is set as the first optimized point cloud. For the second smoothed point cloud, input the preset target downsampling rate, calculate its bounding box volume and the total number of points in the point cloud, and estimate the initial voxel side length based on the target downsampling rate, bounding box volume and the total number of points in the point cloud; The three-dimensional space of the second smooth point cloud is divided into a voxel grid according to the initial voxel side length. The actual number of points contained in each non-empty voxel is counted, and the average density of the non-empty voxels is calculated. The average density is compared with the target number of points per voxel, where the target number of points per voxel is the reciprocal of the target downsampling rate. If the average density is less than the target number of voxels, the voxel side length is increased by a preset growth factor; if the average density is not less than the target number of voxels, the voxel side length is decreased by a growth factor; the process is iterated until the average density approaches the target number of voxels to obtain the optimal voxel side length. Based on the optimal voxel side length, the second smooth point cloud is re-divided into voxel grids. For each non-empty voxel, the centroid coordinates of all points within the target non-empty voxel are calculated. The centroid coordinates are used to replace all the original points within the target non-empty voxel to complete dimensionality reduction and compression, resulting in the second optimized point cloud with a reduced number of points but unchanged geometric features.

4. The method for analyzing and evaluating farmland soil structure according to claim 3, characterized in that, The process involves calculating the rotation matrix and translation vector between the first and second optimized point clouds using coarse registration based on geometric feature matching and fine registration based on distance minimization. The second optimized point cloud is then transformed into the coordinate system of the first optimized point cloud to generate a fused point cloud that simultaneously possesses the detailed precision of the first point cloud data and the completeness of the second point cloud data. This includes: For the first optimized point cloud and the second optimized point cloud, the one with more points is used as the original point cloud and the other is used as the target point cloud; the elevation range of the two point clouds is obtained, and two elevation thresholds are calculated according to the first segmentation coefficient and the second segmentation coefficient. The regions with elevation values ​​less than the first threshold or greater than the second threshold are retained, and the flat regions in the middle are removed to obtain the original optimized point cloud and the target optimized point cloud after elevation screening. The normal vectors of each point in the filtered point cloud are estimated, a local coordinate system is established for each point, and the three angular features of the point and its neighboring points are calculated. The three angular features are divided into 11 statistical intervals and normalized, and then spliced ​​to generate a 33-dimensional fast point feature histogram descriptor. At least three points are randomly sampled from the original optimized point cloud. The corresponding points are searched in the target optimized point cloud based on the nearest neighbor matching of the FPFH descriptor. Polygons are constructed with the sampled points and the corresponding points respectively. The dissimilarity vectors of the corresponding side lengths are calculated. Each element is the absolute value of the side length difference divided by the larger value of the two corresponding side lengths. Based on the preset side length difference threshold, if any dissimilarity vector element is greater than the side length difference threshold, it is judged as a mismatch and resampling is performed. The corresponding point pairs that meet the threshold condition are retained. The initial rotation matrix and translation vector are solved by singular value decomposition to obtain the coarse registration transformation matrix and transform the target optimized point cloud. Based on the coarsely registered target optimized point cloud, a KD-Tree spatial index is constructed for the original optimized point cloud. For each point in the target optimized point cloud, the corresponding point with the closest Euclidean distance in the original optimized point cloud is searched using the KD-Tree. An optimization model is constructed with the weighted sum of the squared Euclidean distances between corresponding points as the objective function. The optimal rotation matrix and translation vector that minimize the objective function are solved iteratively. After each iteration, the point cloud coordinates are updated and the correspondence is re-established until the number of iterations reaches the maximum value or the error change is less than the convergence threshold, and the final transformation matrix is ​​obtained. The second optimized point cloud is transformed to the coordinate system of the first optimized point cloud according to the final transformation matrix to generate the fused point cloud; the bidirectional Hausdorff distance between the two point clouds is calculated, and if it is less than a preset threshold, the registration is deemed qualified and the fused point cloud is output; otherwise, the parameters are adjusted and the registration is repeated.

5. The method for analyzing and evaluating farmland soil structure according to claim 1, characterized in that, The process of extracting elevation values ​​from all points in the fused point cloud, removing overlapping points to form an elevation dataset; calculating the root mean square height of the vertical undulation based on the elevation dataset, and calculating the autocorrelation length of the horizontal spatial correlation length based on the autocorrelation function of the horizontal elevation sequence, includes: Traverse all points in the fused point cloud, extract the Z-axis coordinate value of each point as the elevation value, and construct the original elevation dataset. Using the spatial hash grid method, based on a preset distance threshold, points whose Euclidean distance is less than the distance threshold are identified as overlapping points and removed, retaining a single representative point to obtain the deduplicated elevation dataset; The total number of elevation points in the elevation dataset is counted, the arithmetic mean of all elevation values ​​is calculated, and the root mean square height is calculated. The root mean square height represents the degree of vertical undulation of the soil surface. The larger the value, the more uneven the surface, and the smaller the value, the flatter the surface. Elevation sequences are extracted from the fused point cloud along the horizontal X-axis at a preset sampling interval, or along the horizontal Y-axis at a preset sampling interval; a normalized autocorrelation function is calculated based on the elevation value and its displacement, where the displacement is equal to the sampling interval multiplied by a natural number index; the autocorrelation function values ​​corresponding to different displacements are calculated iteratively to generate an autocorrelation function curve; The autocorrelation function curve is smoothly fitted using the cubic spline interpolation method. The displacement value corresponding to the autocorrelation function value being equal to the reciprocal of the natural constant is calculated. The displacement value is determined as the autocorrelation length. The autocorrelation length characterizes the degree of spatial correlation of farmland soil surface in the horizontal direction. The larger the value, the gentler the surface undulation and the stronger the spatial continuity. The smaller the value, the more intense the surface undulation and the more significant the spatial heterogeneity. High-order sequences are extracted along the X-axis, Y-axis, and at a 45-degree angle to the X-axis, respectively. The autocorrelation length in each direction is calculated, and the arithmetic mean of the three is taken as the final autocorrelation length.

6. The method for analyzing and evaluating farmland soil structure according to claim 5, characterized in that, The process involves extracting elevation data from the fused point cloud along the horizontal X-axis at a preset sampling interval, or along the horizontal Y-axis at a preset sampling interval; calculating a normalized autocorrelation function based on the elevation value and its displacement, where the displacement is equal to the sampling interval multiplied by a natural number index; iterating through and calculating the autocorrelation function values ​​corresponding to different displacements; and generating an autocorrelation function curve, including: From the fused point cloud, discrete sampling is performed along the target horizontal direction at preset sampling intervals to obtain an elevation value sequence arranged in spatial order, denoted as the elevation sequence; the length of the elevation sequence is the total number of sampling points N. Define the sequence index displacement j, where j is a natural number and its value ranges from 1 to N-1; define the actual physical displacement x, where the actual physical displacement x is equal to the sequence index displacement j multiplied by the preset sampling interval; iterate through each integer value of the sequence index displacement j from 1 to N-1 and calculate the normalized autocorrelation function value corresponding to each displacement. For each sequence index displacement j, starting from the first position of the high-order sequence, traverse from index i from 1 to Nj, taking the i-th elevation value and the (i+j)-th elevation value in sequence, where i is a natural number from 1 to N, multiply the two, sum all the products, and use this as the numerator; sum the squares of all N elevation values ​​in the high-order sequence, and use this as the denominator; divide the numerator by the denominator to obtain the normalized autocorrelation function value corresponding to the sequence index displacement j. A discrete point set is constructed with the actual physical displacement x as the abscissa and the corresponding normalized autocorrelation function value as the ordinate. The discrete point set is then smoothly fitted using cubic spline interpolation to generate a continuous autocorrelation function curve. On the autocorrelation function curve, solve for the actual physical displacement value corresponding to the normalized autocorrelation function value being equal to the reciprocal of the natural constant; determine the corresponding actual physical displacement value as the autocorrelation length in the target horizontal direction.

7. The method for analyzing and evaluating farmland soil structure according to claim 1, characterized in that, The method further includes: Based on the values ​​of root mean square height and autocorrelation length, and in accordance with the preset soil roughness classification threshold, the soil structure is classified into one of smooth, transitional, or rough types; a structural analysis report containing root mean square height, autocorrelation length, and soil structure evaluation level is generated and output.

8. A device for analyzing and evaluating the structure of farmland soil, characterized in that, The device includes: The point cloud imaging module is used to acquire first point cloud data and second point cloud data of the soil surface of the farmland to be measured through imaging equipment. The first point cloud data has absolute coordinates and actual physical scale, and the second point cloud data contains color information. The point cloud optimization module is used to denoise, filter and extract target regions from the first point cloud data and the second point cloud data respectively to obtain the first optimized point cloud and the second optimized point cloud. It also adaptively estimates the voxel side length according to the preset downsampling rate and performs dimensionality reduction and compression on the second optimized point cloud through voxel core downsampling to obtain the second optimized point cloud with reduced number of points but preserved geometric features. The point cloud fusion module is used to solve the rotation matrix and translation vector between the first optimized point cloud and the second optimized point cloud based on coarse registration of geometric feature matching and fine registration based on distance minimization, and transform the second optimized point cloud into the coordinate system of the first optimized point cloud to generate a fused point cloud that has both the detail accuracy of the first point cloud data and the integrity of the second point cloud data. The elevation calculation module is used to extract the elevation values ​​of all points from the fused point cloud, remove overlapping points to form an elevation dataset; calculate the root mean square height of the vertical undulation based on the elevation dataset, and calculate the autocorrelation length of the horizontal spatial correlation length based on the autocorrelation function of the horizontal elevation sequence. The structural analysis module is used to classify farmland soil structure into different roughness levels based on the root mean square height and autocorrelation length, and compare them with preset grading thresholds, and generate structural analysis results.

9. A control device, characterized in that, The device includes: A memory and a processor, wherein the memory stores a computer program that can be loaded by the processor and executed as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer program is stored that can be loaded by a processor and executed as described in any one of claims 1 to 7.

Images