A method for estimating the sampling amount of a table based on a monitoring binocular camera and a close-up monocular camera

Abstract

The application obtains a method capable of estimating the sampling amount of soil or other similar texture sampling objects in real time based on monitoring the image information of binocular cameras and close-up monocular cameras, wherein the method comprises: 1) three-dimensional reconstruction point cloud data of a sampling point and classification processing of corresponding point cloud data of soil in different loose states; 2) evaluating the loose volume ratio of soil in different loose states at the sampling point, and then formulating the corresponding equivalent volume proportionality coefficient; 3) comparing the point cloud data before and after sampling and combining the equivalent volume proportionality coefficient of each part of soil to calculate the sampling amount estimation result. The core point of the application is to divide soil in different loose states (turned-out soil, sunken soil and natural soil) to solve the key problem faced by sampling amount estimation. As a completely visual sampling amount estimation method, the application can be widely applied in the fields of agricultural automation, mechanization and intelligentization and celestial body detection.

Description

Technical Field

[0001] This invention relates to the field of soil characteristic measurement and point cloud processing, and in particular to a method for estimating the amount of surface sampling based on a monitoring binocular camera and a close-up monocular camera. Background Technology

[0002] Existing soil characteristic measurement schemes documented in the literature can be mainly divided into two categories: wide-area measurement and sample measurement. The former generally involves remote observation of the topographic and geomorphological features of large-scale soil areas, while the latter often involves measuring samples in a laboratory setting with specific instruments and conditions. Wide-area measurements, for example, utilize ground and aerial photogrammetry to assess the extent of soil erosion over large areas, and also employ remote observation by drones combined with orthophoto mapping and digital elevation models to measure the geomorphological features of large-area soil cavities and channels formed by natural erosion. Furthermore, wide-area soil characteristic measurements play an important role in extraterrestrial exploration (such as lunar and Mars exploration). Sample measurements, on the other hand, can measure the volume, composition, and particle size of soil samples using more detailed detection methods. Examples include directly measuring the physical dimensions of soil clods, indirectly measuring soil volume through drainage methods, and using multifractal analysis and laser diffraction to detect the distribution of soil particle diameters in dry soil. The soil characteristic measurement needs addressed by this invention differ from the two existing types mentioned above: this invention assesses local soil volume changes within a sampling area rather than large-scale, wide-area topography; and it achieves real-time soil characteristic measurement (sampling volume estimation) rather than uniform measurement after sampling. In summary, this invention does not contain any prior art with similar background and objectives, but rather is specifically designed for a novel soil characteristic measurement need.

[0003] This invention calculates the sampled volume solely based on visual information. It first uses visual information to reconstruct a 3D point cloud of the sampled area, and then performs volume calculations based on the point cloud. Therefore, the original data addressed by this invention consists of two sets of point cloud data—before and after sampling—that do not constitute a closed volume. By comparing the differences between the two sets of point cloud data, the volume values ​​of each part of the point cloud change caused by sampling are obtained. In this functional direction, the open-source software CloudCompare can already perform various calculations and analyses based on point clouds. Its point cloud mesh subtraction and Poisson surface fitting functions are also suitable for some of the calculation requirements of this invention. Therefore, this invention uses the underlying code of CloudCompare as a reference to build the most basic point cloud volume calculation function, while all other innovative core functions are independently designed and completed by this invention. The programming development of this invention is based on the Windows 10 system and Visual Studio 2019 software, using C / C++ as the programming language. The main special libraries used are the Point Cloud Library (PCL) and the Computational Geometry Algorithms Library (CGAL). Summary of the Invention

[0004] The present invention was made in view of the above circumstances, and its purpose is to design and implement a method for estimating the table sampling quantity based on image information from a monitoring binocular camera and a close-up monocular camera under unattended conditions.

[0005] To achieve the above objectives, this invention specifically addresses a key issue in image-based estimation of sampling volume: generally, the sampler disrupts the natural soil condition during sampling, resulting in varying soil compaction in different parts of the sampling area. Ignoring this compaction issue leads to significant errors in the calculation results. This invention divides the soil damaged by sampling into three parts: 1) natural soil, i.e., the part unaffected by sampling and retaining its original soil characteristics; 2) overturned soil, i.e., the soil accumulated outside the sampling pit due to the extra volume insertion of the sampling shovel and changes in soil porosity; 3) indented soil, i.e., the loose soil that slips into the sampling pit due to the loss of static equilibrium caused by sampling disturbance.

[0006] Based on the above classification of soils in different compaction states, the present invention designs a sampling volume estimation scheme as follows. First, a binocular camera reconstructs the three-dimensional point cloud of the sampling area before and after sampling based on the principle of binocular vision. The difference between the heights of the two is used to judge each point cloud. If the value is greater than 0, the height of that part of the soil increases, corresponding to the overturned soil outside the pit; if the value is less than 0, the height of that part of the soil decreases, corresponding to the part inside the pit. The part inside the pit includes the natural soil that was sampled and the subsidence soil part. By comparing the actual point cloud with the point cloud of the ideal sampling trajectory, the subsidence soil part above the ideal depression is obtained. A close-up monocular camera provides soil images of the sampling area. After image standardization, the loose soil and solid soil parts in the image are divided by a one-time manual operation. Then, the porosity of the soil in each sub-image is calculated based on features such as texture and color. Furthermore, the equivalent volume ratio coefficient of each part of the soil in different compaction states is obtained by combining the theoretical relationship between porosity and volume. In addition, the soil quality of the sampling area can be detected by using soil image features captured by a close-up camera or spectral information from a geological spectrometer. Then, the loose volume ratio of the soil in the sampling area can be obtained by looking up a table of known soil loose volume ratios. Finally, the sampling volume estimation result can be obtained by multiplying each part of the point cloud data—natural soil, subsidence soil, and overturned soil—with their corresponding equivalent volume ratio coefficients and then superimposing the results.

[0007] Based on the above sampling estimation scheme, the framework of the calculation procedure designed in this invention is as follows:

[0008] Step 1: Import the global point cloud before and after sampling;

[0009] Step 2: Automatically extract the neighborhood of the sampling pits in the global point cloud;

[0010] Step 3: Remove outliers from the neighborhood point cloud of the sampling pit before and after sampling to prevent outliers from interfering with subsequent volume calculations.

[0011] Step 4: Perform plane fitting on the point cloud before sampling. Using the standard of rotating the plane to be parallel to the xOy plane and translating the center to the centroid position of the plane, correct the direction of the point cloud coordinates and the origin position before and after sampling.

[0012] Step 5: Compute the sampling volume of the points before and after sampling based on the aforementioned standardized processing.

[0013] Step 6: Display and save the volume calculation results and related information.

[0014] Furthermore, in order to achieve the above objectives, the present invention is designed with the following special program functions:

[0015] In step 1, the grid parameters required for subsequent calculations can be automatically calculated based on the characteristics of point cloud data, such as point density. This mainly includes two functions: automatically obtaining the coordinates of the extreme points in the point cloud and automatically determining the grid step size and the number of horizontal and vertical grids.

[0016] Step 2 first obtains a coarse screening of the neighborhood of the sampling pit based on the three-dimensional coordinates of the ideal sampling point (xyz), and then further obtains a fine screening of the neighborhood using extreme value search, 8-connected region search, etc.

[0017] Step 3 can be selected from two outlier removal strategies proposed in this invention: a two-level screening and removal method based on z-value and a screening and removal method based on point cloud distance.

[0018] Step 4 uses least squares fitting, rotation matrix, and translation matrix to perform coordinate correction on the point cloud before and after sampling.

[0019] Step 5 uses two methods to calculate the volume difference caused by the changes in point cloud before and after sampling: direct calculation using the grid method and calculation after fitting the Poisson surface. In addition, the outward soil part is screened out by subtracting the height of the point cloud before sampling from the point cloud after sampling and judging the positive or negative of the result. The volume of the inward soil is evaluated by subtracting the ideal sampling volume from the volume in the pit of the reconstructed point cloud.

[0020] Step 6 includes colorizing and saving the grid filling results before and after sampling, as well as the grid method calculation results, according to the standard JET color bar.

[0021] Invention Effects

[0022] According to the present invention, estimation results of sample volume can be obtained based on image information from a monitoring binocular camera and a close-up monocular camera under unattended conditions. Through 17 sets of actual soil sampling experiments, the test results show that the relative error of the estimation results for approximately 76% of the sampling pits is below 15%, and the relative errors for approximately 65% ​​and 35% of the sampling pits are below 10% and 5%, respectively. The minimum absolute error can reach 1.6g, and the minimum relative error can reach 1.1%. Excluding abnormal pits, the average absolute error is 7.9g, and the root mean square error is 9.4g. This demonstrates that the present invention can provide valuable sample volume estimation results based solely on image information, representing a novel method for soil characteristic measurement (sample volume estimation). Furthermore, the innovative two-level screening and removal method based on z-value and the screening and removal method based on point cloud distance designed in this invention can be applied as general-purpose point cloud preprocessing methods to other point cloud processing problems. Attached Figure Description

[0023] Figure 1 This is a flowchart of the overall scheme for the method of estimating the sample size.

[0024] Figure 2 This provides the overall framework for the sample size estimation procedure.

[0025] Figure 3 This is a flowchart illustrating the overall process of automatically capturing the neighborhood of the sampling pit in step 2 of the specific implementation.

[0026] Figure 4 The image shows the fine sieving effect of the sampling pit neighborhood obtained in step 2.6 of the specific implementation method, where the left image is the sampling pit neighborhood point cloud before sampling and the right image is the sampling pit neighborhood point cloud after sampling.

[0027] Figure 5 This is a flowchart of the first implementation method under step 3 of the specific implementation plan—a two-level filtering and removal method based on z-values.

[0028] Figure 6 This is a flowchart of the second implementation method under step 3 of the specific implementation plan—the filtering and removal method based on point cloud distance.

[0029] Figure 7 This is a schematic diagram illustrating the principle of using the vector inner product to distinguish between some outliers and sparse valid points in step 3(2).4 of the specific implementation method.

[0030] Figure 8 The flowchart for step 3(2).5.1 of the specific implementation is shown.

[0031] Figure 9 This is a schematic diagram illustrating the principle of distinguishing outliers and sparse valid points using projection length in step 3(2).6 of the specific implementation method.

[0032] Figure 10 This is a flowchart of the first implementation method—direct calculation using the grid method—in step 5 of the specific implementation.

[0033] Figure 11 The diagram shows the triangular mesh required for deriving the calculation method of the centroid coordinate interpolation in step 5(1).3 of the specific implementation method.

[0034] Figure 12 The example of an ideal sampling trajectory point cloud in step 5(1).4.2 of the specific implementation method is shown in green, where the left figure is the observation view above the side of the soil point cloud and the right figure is the observation view below the side of the soil point cloud.

[0035] Figure 13 This is a schematic diagram of the soil division method of each part of the sampling pit neighborhood used in the derivation of the sampling quantity estimation formula in step 5(1).6 of the specific implementation method.

[0036] Figure 14This is a flowchart of the processing of the second implementation method—Poisson surface fitting method—in step 5 of the specific implementation method.

[0037] Figure 15 The following is an example of the point cloud results of grid filling before and after sampling in step 6 of the specific implementation method. The left image is the point cloud before sampling, and the right image is the point cloud after sampling. The color of the point cloud is more reddish the farther away from the soil plane, and more blue the closer to the soil plane.

[0038] Figure 16 This is an example of the point cloud of the mesh result after subtraction in step 6 of the specific implementation method, where the larger the absolute value of the subtraction, the more red the color and the smaller the absolute value, the more blue the color.

[0039] Figure 17 This is an example of the Poisson surface reconstruction result in step 6 of the specific implementation method. Detailed Implementation

[0040] Hereinafter, examples of methods for implementing the technology of the present invention will be described in detail with reference to the accompanying drawings.

[0041] Figure 1 This is a flowchart illustrating the overall sampling estimation scheme described in the above-described invention. Based on this overall scheme flowchart, the overall program framework designed in this invention is as follows: Figure 2 As shown below, the program steps and their related technologies, principles, processing flow and other details are described in detail below.

[0042] Step 1: Import point clouds before and after sampling

[0043] Import the global point cloud before and after sampling. Alternatively, you can choose to directly import the neighborhood of the sampled pits after manual extraction (in which case the program will automatically skip step 2 below).

[0044] Step 2: Automatic extraction of the neighborhood of the sampling pit

[0045] The target sampling pit's neighborhood point cloud is automatically extracted from the global point cloud imported in step 1 above, significantly reducing the computational load and time for subsequent sampling estimation and minimizing interference from areas outside the sampling pit's neighborhood. The processing flow for automatically extracting the sampling pit's neighborhood is as follows: Figure 3 As shown below, the specific details of each part of the process are described in the following text.

[0046] Step 2.1: Coarse screening of the sampling pit neighborhood

[0047] The xyz three-dimensional coordinates of the ideal sampling point (given by the sampling requirements), combined with the width of the sampling shovel (determined by the hardware conditions of the sampling device) and with a large margin (which can be customized), directly intercept the xyz coordinate range to obtain a coarse screening of the neighborhood of the sampling pit.

[0048] Step 2.2: SOR (Statistical Outlier Removal) Outlier Filtering and Removal

[0049] The coarse-screened point cloud results are filtered and removed using SOR (Striking Around the Origin) outlier removal to prevent outliers from interfering with subsequent fine-screening of the sampling pit neighborhood. SOR is a classic outlier removal method (included in the PCL point cloud library), which judges outliers based on the average and standard deviation of the distances to neighboring points. The specific calculation basis is as follows:

[0050] d max =d average +nσ d (1)

[0051] In the formula, d max The maximum allowed distance value (if the calculated maximum distance to a point exceeds this value, that point will be identified as an outlier and removed), d average σ is the average distance from a given point to k points (k is user-defined), where n is the user-defined standard deviation multiplier, and σ is the average distance from that point to the given point. d Let $\frac{k}{k}$ be the standard deviation of the distances from a given point to that point.

[0052] Step 2.3: Calculate the volume difference before and after sampling and save the results.

[0053] The volume difference before and after sampling is calculated using the grid method. This method is detailed in step 5 below, implementation method 1. In short, the point cloud is subdivided according to the grid, the grid value is set according to the height of the data points in the grid, the values ​​of the corresponding grids before and after sampling are subtracted, and finally the subtraction results of all grids are superimposed and multiplied by the grid area to obtain the volume difference (integral idea).

[0054] Step 2.4: Coarse screening of neighborhood maximum / minimum points

[0055] Find the maximum and minimum values ​​near the ideal sampling point in the volume difference result point cloud. These two values ​​correspond to the outward peak and the bottom of the pit, respectively. (Note that if the outward outward is severe, some outward outward areas will not be included in the selection area if the maximum value is not considered.)

[0056] Step 2.5: Connected Component Search

[0057] Using the maximum and minimum points as seed points respectively, search for and mark the 8-connected regions based on the given threshold judgment conditions (which can be customized).

[0058] Step 2.6: Fine sieving of the sampling pit neighborhood

[0059] Calculate the boundary coordinates of the connected region, and then truncate the xyz coordinate range of the point cloud based on this boundary, retaining a certain margin (which can be customized). This yields the neighborhood of the sampled pits after fine screening, as shown in the image. Figure 4 As shown.

[0060] Step 3: Outlier Filtering and Removal

[0061] Because outliers near sampling pits in the reconstructed point cloud can cause significant errors in volume calculation, accurate calculation of the sampling volume requires the program to have a reliable and stable outlier filtering and removal function. The program designed in this invention provides two outlier filtering and removal methods: a two-level filtering method based on z-values ​​(first embodiment) and a filtering method based on point cloud distance (second embodiment). It is recommended to use the first embodiment when the sampling and monitoring camera shooting environment is good, and to use the second embodiment when the sampling and monitoring camera shooting environment is moderate or poor.

[0062] [First Implementation Method – Two-Tier Filtering Removal Method Based on z-Value]

[0063] The processing flow of this implementation method is as follows: Figure 5 As shown, its characteristics are: fast program execution speed, able to handle most sampling and outlier cases (a special case is when the point cloud has a significant slope, causing outliers to extend in a direction almost parallel to the xOy plane, in which case the z value of the outlier does not change much, so the method fails), and in extreme cases, it can determine whether it is worth executing based on the magnitude of the removal.

[0064] Step 3(1).1: First-level outlier filtering and removal based on z-value distance

[0065] Step 3(1).1.1: Search for the data point with the largest z-value in the point cloud.

[0066] Step 3(1).1.2: Count the number of data points within a certain z-coordinate range below this point.

[0067] Step 3(1).1.3: Further judgment and search based on the number of data points.

[0068] If the number of data points is less than the set threshold (which can be customized), the highest point is considered an outlier (because outliers are generally caused by factors above the soil, such as sampler occlusion, so only outliers above the point cloud are considered), the point is marked, and the search continues for the next highest point until it exceeds the set threshold (this filtering process can be referenced). Figure 5 (Illustrative image on the left side of the middle section).

[0069] Step 3(1).2: SOR method trial processing

[0070] Select the SOR parameter with a stronger removal force (to avoid the problem of incomplete removal of some point clouds, which can be customized) to perform trial processing on the point cloud.

[0071] Step 3(1).3: Determine the effect of the trial treatment and decide whether to retain the secondary screening removal results of SOR.

[0072] When the SOR trial processing causes a significant drop in the highest point of the point cloud, and the drop exceeds the threshold (which can be customized), it is considered that the first-level filtering in step 3(1).1 is not clean, and the point cloud result after the second-level filtering by SOR is retained.

[0073] [Second Implementation Method - Point Cloud Distance-Based Filtering and Removal Method]

[0074] The processing flow of this implementation method is as follows: Figure 6 As shown, its characteristics are: it can avoid sparse valid points being regarded as outliers for removal to the greatest extent, the outlier removal effect is good, and the program execution speed is slightly slower than the first implementation method, but the overall execution time is still at a relatively fast level (generally about 2 seconds in the experimental test of this invention).

[0075] Step 3(2).1: Extract candidate outlier points using the SOR method

[0076] The SOR method is used to filter out outlier candidate points in the point cloud to be processed, and these candidate points are listed. The candidate points include "true" outliers and "false" outliers, which are actually sparse and valid points. The latter refers to data points that have few nearby data points and are misidentified as outliers by the SOR method. They generally appear at outward peaks or steep pit walls in the reconstructed point cloud (the essential reason is that the number of pixels at such locations in the stereo image is small, resulting in low reconstruction accuracy).

[0077] Step 3(2).2: KNN (K-Nearest Neighbor) search

[0078] The candidate points are traversed, and KNN (Knowledge Neighbors) is used to search for data points near each candidate point (which can greatly improve the search speed and is included in the PCL point cloud library).

[0079] Step 3(2).3: Find the vector pointing from the candidate point to its neighboring points.

[0080] For each neighboring point of a candidate point, calculate the vector formed by the candidate point and the candidate point itself, and then take the dot product of each of these vectors with the vectors formed by the candidate point's other neighboring points. The dot product is the vector dot product, and its formula is as follows:

[0081]

[0082] In the formula, Let (x, y, z) be a vector, and (x, y, z) be its corresponding three-dimensional coordinates, with subscripts 1 and 2 corresponding to each. To reduce computation, the program in this invention first calculates the first neighboring vector, and then calculates the inner product of each new neighboring vector with the first neighboring vector.

[0083] Step 3(2).4: Inner product judgment and filtering

[0084] The algorithm checks if there are neighboring points with an inner product less than 0. If so, the candidate point is considered a sparse valid point and added back to the point cloud; otherwise, it proceeds to the next layer of filtering. The principle behind this decision is explained in detail below. The relevant formulas for the existence of inner products are as follows:

[0085]

[0086] In the formula, For vector magnitude, The angle between two vectors. For example... Figure 7 As shown, when there are outliers far from the effective surface (i.e., a quasi-continuous surface composed of closely spaced effective points), the angle between vectors pointing from candidate points to neighboring points is less than 90°, and the inner product of these vectors is greater than 0. However, for sparse effective points and outliers close to the effective surface, there is generally at least one set of angles greater than 90° between their many three-dimensional vectors, and the inner product is less than 0. Therefore, outliers far from the effective surface can be screened based on whether there is an inner product less than 0, which greatly reduces the computational load for the second screening in step 3(2).5.

[0087] Step 3(2).5: Calculate the projection length index of the remaining candidate points

[0088] Calculate the projection length of any neighbor vector of the remaining candidate point after step 3(2).4 onto the nearest valid surface normal vector of the candidate point.

[0089] Step 3(2).5.1: Calculate the distance from the candidate point to the nearest valid surface.

[0090] The calculation process is as follows: Figure 8 As shown, the calculation method can be visualized as follows: using candidate points as the center of a circle, the radius of the circle is continuously increased by a fixed step size until a certain step size increase causes a large number of valid data points to flood in, at which point the circle touches the valid surface (this state corresponds to...). Figure 8 The inner circle), and then continue to increase the step size by one unit, within this increment range (corresponding to Figure 8 If a large number of valid data points flow into the area between the inner and outer circles, it proves that the effective surface has indeed been reached.

[0091] Step 3(2).5.1.1: Distance information statistics of candidate points

[0092] The distance information from the candidate point to all other data points is calculated at intervals with an accuracy of 1 (corresponding to 1 mm in reality) (this statistical accuracy can be customized and adjusted). For example, there are 5 data points in the range of 0 to 1, and 8 data points in the range of 1 to 2.

[0093] Step 3(2).5.1.2: Preliminary calculation of the possible distance values ​​between the candidate points and the effective surface.

[0094] The data points are judged sequentially in ascending order of distance. The position where the number of data points in the first range of d0 to (d0+1) is greater than the set threshold (which can be customized) is detected, and d0 is recorded as the possible distance value of the candidate point from the effective surface.

[0095] Step 3(2).5.1.3: Determine whether the calculated distance value is accurate and reliable.

[0096] Perform incremental judgment, that is, judge whether the number of data points in the range (d0+1) to (d0+2) (referred to as the incremental range or incremental area) is greater than the set threshold (which can be customized). If so, d0 is accepted as the distance of the candidate point from the effective surface and the loop is exited. Otherwise, it means that the range d0 to (d0+1) has not reached the effective surface (because the density of the points on the 1mm scale is not continuous) and the large number in the range is due to interference from other relatively concentrated outliers. The judgment result of d0 is rejected and the sequential judgment process in step 3(2).5.1.2 above continues.

[0097] Step 3(2).5.2: Calculate the incremental region neighbor vector of the candidate point

[0098] Take the vector pointing from the candidate point to any point within the increment range. This vector is named the incremental region neighbor vector for ease of explanation in subsequent calculation operations.

[0099] Step 3(2).5.3: Calculate the normal vector of the effective surface within the increment range.

[0100] Two vectors are formed by any four points within the increment range. and The cross product yields the normal vector representation of the effective surface within the increment range. The calculation formula is as follows:

[0101]

[0102] Step 3(2).5.4: Calculate the projection length index

[0103] The projection length L of the candidate point is calculated using the following formula:

[0104]

[0105] Step 3(2).6: Determine and filter based on projection length

[0106] The algorithm determines whether the projected length L of a candidate point is less than a set threshold (which can be customized). If it is, the candidate point is considered a sparse and valid point and added back to the point cloud. Otherwise, it is identified as an outlier and removed. The principle behind this determination is as follows: Figure 9 As shown, for sparse valid points, since they are not outside the valid surface, the calculated projection length is generally small; while for outliers outside the valid surface, their projection length is large. Therefore, by setting a projection length threshold, sparse valid points and outliers close to the valid surface can be distinguished. It should be noted that although projection length judgment can directly filter out sparse valid points, its calculation process is relatively complex. Therefore, using inner product judgment for primary filtering is essential and can significantly reduce the overall computation time.

[0107] Step 4: Point cloud coordinate correction

[0108] The point cloud volume calculation process of this invention uses the z-value as the height information by default. Therefore, a coordinate correction method is designed to make the height direction of the point cloud as parallel as possible to the z-axis direction before and after sampling.

[0109] Step 4.1: Soil plane fitting

[0110] Calculate the least squares fitting plane for the point cloud before sampling (the point cloud before sampling is generally relatively flat, so its plane fitting effect is relatively reliable) and make its origin coincide with the centroid of the point cloud before sampling.

[0111] Step 4.2: Calculate the rotation and translation matrices.

[0112] Based on the normal vector of the fitted plane and the position of the origin, we can obtain the rotation matrix that rotates the z-axis of the point cloud before sampling to be parallel to the normal vector, and the translation matrix that translates the origin to the centroid position (these can be solved using the built-in calculation functions in the point cloud library PCL).

[0113] Step 4.3: Perform the same coordinate correction operation on the point cloud before and after sampling.

[0114] By using the rotation and translation matrices obtained from the point cloud before sampling, the same coordinate rotation and translation are performed on the point cloud before and after sampling to achieve coordinate correction of the point cloud.

[0115] Step 5: Point cloud volume calculation

[0116] The program of this invention uses two methods to calculate the volume difference caused by the change of point cloud before and after sampling: direct calculation using the grid method (first embodiment) and calculation after fitting the Poisson surface (second embodiment). The sampling quantity estimation results obtained by the two methods can be used as mutual reference to increase the reliability of the results and provide more information to users.

[0117] When the changes in soil compaction caused by the sampling described above are not considered, the results of the two implementation methods should be similar. If there is a large deviation between the two, it is generally due to the interference of stones in the sampling area or the abnormal estimation of the sampling volume caused by factors such as the accuracy of the point cloud reconstruction. In this case, it is necessary to observe the status of each key point cloud shown and saved in step 6 and conduct further problem analysis and investigation.

[0118] When considering the different densities and looseness states of natural soil, overturned soil, and subsidence soil, the first implementation method (i.e., direct calculation using the grid method) should be used to estimate the sampling volume. The second implementation method does not have the function of sieving the soil point cloud of each part. Therefore, its result is only a measure of the degree to which the influence of looseness variation is considered. If the two results are still similar, it means that the looseness variation of soil in each part of the sampling pit is not obvious or the volume of soil with large looseness variation is negligible compared to the total sampling volume.

[0119] [Implementation Method 1 - Direct Calculation Using the Mesh Method]

[0120] The overall computational process of the mesh method is as follows: Figure 10 As shown, the idea is to project the point cloud in a certain direction and simplify it into infinitesimal elements, then process and calculate each element separately and finally superimpose them.

[0121] Step 5(1).1: Create a mesh based on the set parameters and point cloud features

[0122] A mesh is created using user-defined parameters such as the plane, minimum coordinate points, number of horizontal and vertical grid points, and step size. To simplify the setting process for some parameters, the program includes a function to automatically calculate some mesh parameters based on features such as point density in the point cloud data. Specifically, this includes:

[0123] 1) Automatically obtain the coordinates of the extreme points in the point cloud: The program iterates through all data points in the point cloud before and after sampling and compares them to obtain the maximum and minimum values ​​of the xyz coordinates.

[0124] 2) Automatically determine grid step size and number of horizontal and vertical grid points: i) The program counts the number of data points in the point cloud, obtains the range of point cloud data based on the coordinates of the extreme points, and calculates the average point cloud data density in the x and y directions; ii) The grid step size is automatically set to a multiple of the average point cloud data density by the user-specified scaling factor; iii) The number of horizontal and vertical grid points (corresponding to the x and y coordinates respectively) under the given grid step size is calculated based on the coordinates of the extreme points.

[0125] Step 5(1).2: Data points are projected onto the grid and filled.

[0126] This section references the underlying code of the aforementioned CloudCompare open-source software, which performs statistical analysis on the data points within each grid and assigns corresponding values ​​to each grid according to the user-defined projection strategy. Optional projection strategies include projecting based on the minimum, average, maximum, and median values ​​of all data points within the grid.

[0127] Step 5(1).3: Fill blank grid

[0128] This section references the underlying code of the aforementioned CloudCompare open-source software, which implements the filling of blank grids according to user-defined strategies. Optional blank grid filling strategies include filling with the minimum, average, or maximum values ​​of all grids, or not filling, filling with user-defined values, or interpolation. Interpolation is primarily based on Delaunay triangulation and barycentric coordinate interpolation. The specific implementation process of interpolation is as follows: 1) Filter out all valid points in the point cloud; 2) Construct a Delaunay triangulation based on all valid points; 3) Traverse all triangulations and filter out all invalid points in the point cloud within the grid range; 4) Calculate the relevant parameters for calculating the barycentric coordinates at the invalid points based on the projected x and y coordinates of the invalid points and the corner points of the Delaunay triangulation (assuming x and y are position information, and z is height information), and then perform a weighted summation of the z values ​​of the Delaunay triangulation corner points based on these parameters to finally obtain the interpolation at the invalid points (i.e., the z value of the barycentric coordinates at that point). The specific calculation method for barycentric coordinate interpolation is described below. Figure 11 In the model, let the vertices of the triangle be A, B, and C, and use x... * y * z * Let α, β, and γ represent their coordinate values, respectively. For any point P inside the triangle, there uniquely exist α, β, and γ that simultaneously satisfy:

[0129] α+β+γ=1 (6)

[0130] (x P ,y P )=α(x A ,y A )+β(xB ,y B )+γ(x C ,y C (7)

[0131] By transforming the above equation, we can see that:

[0132]

[0133] The above equation shows that a stable state can be achieved when objects of mass α, β, and γ are suspended at points A, B, and C respectively, with P as the support center. This is the so-called "center of gravity coordinates". From equation (7), we can obtain the following system of equations:

[0134]

[0135] Simplifying the above equation using equation (6), we get:

[0136]

[0137] Rewriting the above system of equations in matrix form, we get:

[0138]

[0139] From the above formula, we can see that if the x and y coordinates of points A, B, C, and P are known, the values ​​of β and γ can be obtained, and then the value of α can be obtained according to formula (6). Finally, the height z at point P is calculated according to the following formula. P The centroid coordinate interpolation can then be obtained:

[0140] z P =αz A +βz B +γz C (12)

[0141] Step 5(1).4: Screening and classification of soils in different compaction states and calculation of their volume

[0142] As mentioned earlier, this invention classifies the soil within the sampling pit's vicinity into overturned soil, natural soil, and indented soil based on their compaction state. The specific sieving methods for each type of soil are detailed below.

[0143] Step 5(1).4.1: Screening of the point cloud corresponding to the overturned soil

[0144] Traverse all grids, subtract the pre-sampling point cloud from the sampled point cloud, and the grids with positive results correspond to the outward soil portion of the point cloud. Superimpose the subtraction results of all grids of this type and multiply by the grid area to obtain the outward soil volume V. outside .

[0145] Step 5(1).4.2: Screening of point cloud corresponding to subsidence soil

[0146] Screening of subsidence soil requires the use of an ideal sampling shovel trajectory. This information can be obtained by continuously photographing the sampling shovel during the sampling process and then fitting the image, or by recording the sampler's pose data during the digging process. Generally, the lowest point of the ideal sampling shovel trajectory point cloud should be close to the lowest point of the actual sampling pit, while the ideal sampling point cloud in other areas of the sampling pit should be lower than the reconstructed soil point cloud (because the upper end of the subsidence soil is reconstructed as the bottom of the actual sampling pit). An example of an ideal sampling shovel trajectory point cloud is shown below. Figure 12 The point cloud is shown in the green area.

[0147] Step 5(1).4.2.1: Calculate the ideal sampling volume

[0148] The point cloud before sampling and the point cloud of the ideal sampling trajectory are respectively meshed (the steps are the same as the meshing method described above), and the values ​​of each mesh are subtracted. Then, the subtraction results of all meshes are superimposed and finally multiplied by the mesh area to obtain the ideal sampling volume V. ideal .

[0149] Step 5(1).4.2.2: Calculate the volume within the pit in the reconstructed point cloud.

[0150] Traverse all grids, subtract the unsampled point cloud from the sampled point cloud, and the grids with negative results correspond to the point cloud inside the pit. Superimpose the subtraction results of all grids of this type and multiply by the grid area to obtain the reconstructed point cloud pit volume V. inside .

[0151] Step 5(1).4.2.3: Estimate the volume of the subsidence soil

[0152] The difference between the ideal sampling volume and the volume within the reconstructed point cloud pit can be used to estimate the volume of the subsided soil. The specific calculation formula is as follows:

[0153] V slide ≈V ideal -V inside (13)

[0154] In the formula, V slide This represents the volume of the subsidence soil.

[0155] Step 5(1).5: Determination of the volume ratio of loose nuts

[0156] like Figure 1As shown, a close-up camera can be used to capture detailed images of the soil in the sampling area. The loose soil and compacted soil in the images are then manually divided, and the image features of each image region are calculated (based on previous literature, indicators such as mean, variance, mean of contrast, and variance of contrast can be used). The ratio of key distinguishing features is then used to estimate the loose-compact volume ratio. Alternatively, a close-up camera or geospectroscopy can be used to detect the soil type in the sampling area, and the loose-compact volume ratio of that soil type can be obtained by referring to previous literature. The specific implementation and verification of the above approach still require future research; only its possibility is proposed here. The user obtains the loose-compact volume ratio α of the soil in the sampling area. set Then, custom settings can be made for it in the program of this invention.

[0157] Step 5(1).6: Formulation and application of equivalent volume scaling factor

[0158] Using the natural soil density as the standard, let α be the equivalent volume ratio coefficient of the outward-turned soil portion. outside The equivalent volume ratio of the subsidence soil portion is α. slide According to Figure 13 Based on the classification of excavated soil, natural soil, and indented soil shown, the following formula can be used to estimate the sampling volume:

[0159] V sample =V inside +α slide V slide -α outside V outside (14)

[0160] In the formula, V sample The sampling volume.

[0161] Whether the subsidence soil is caused by soil slippage or by falling during the sampling and collection process, its state is always loose. Therefore, its equivalent volume ratio can always be estimated using the loose volume ratio determined by the user in step 5(1).5, that is:

[0162] α slide =α set (15)

[0163] Unlike subsidence soil, outward-turned soil exhibits different states of compaction. If the outward turning is significant, it is generally caused by the final accumulation of soil pushed forward during the sampling process. In this case, because this portion of the pushed soil has already undergone "loosening" due to the sampling operation, its state is relatively loose, similar to subsidence soil. Conversely, when the outward turning is less pronounced, it is generally caused by the insertion of extra volume from the sampling shovel, resulting in the bulging of nearby soil. In this case, because this portion of soil has not undergone "loosening" and its bulge originates from the compression of deeper soil layers, its state is not as loose as subsidence soil, but rather close to the compacted state of natural soil.

[0164] To quantify the different eversion conditions described above, users can conduct shallow and deep sampling experiments. Analysis of the experimental data yields specific indicators that can classify the degree of eversion. One discriminant indicator obtained from experiments conducted in the target sampling environment is the external-to-internal volume ratio β:

[0165]

[0166] The larger the outward volume relative to the volume within the reconstructed point cloud pit, the larger the β value, and vice versa. Therefore, the outward-to-inward volume ratio β can be used as an indicator to measure the degree of outward movement. Based on experience from multiple sampling experiments, this invention sets the equivalent volume ratio coefficient of the outward movement as follows:

[0167]

[0168] In the formula, α is in the range of 0.53 < β < 0.57. outside The value of is a linear interpolation based on the magnitude of β. Users can customize this part in the program of this invention, adjusting the specific division index and interpolation method according to the sampling environment. If the user's sampling is relatively stable and uniform, the equivalent volume ratio coefficient α of the overturned soil portion can be simplified. outside It is a user-defined constant.

[0169] Step 5(1).7: Sample Size Estimation

[0170] Combining equations (13) and (14), we can obtain:

[0171] V sample =(1-α) slide V inside +α slide V ideal -α outside V outside (18)

[0172] The above formula is the final calculation formula for the sampling volume estimation. From the first two terms of this formula, it can also be seen that the intervention method of the ideal sampling volume can be understood from another perspective as a weighted average of the volume within the pits of the reconstructed point cloud, where the weight of the ideal sampling volume is α. slide The weight of the reconstructed point cloud pit volume is (1-α) slide ).

[0173] The outward soil volume V obtained in step 5(1).4 outside and the volume of the indented soil V slide and the equivalent volume ratio α of the outward-turned soil portion obtained in step 5(1).6 outside The equivalent volume ratio α of the subsidence soil portion slide Substituting into equation (18), the final calculation realizes the estimation of sampling volume considering different soil compaction conditions in different parts.

[0174] [Second Implementation Method – Poisson Surface Fitting]

[0175] The overall calculation process for point cloud volume calculation based on Poisson surface fitting is as follows: Figure 14 As shown, the specific steps are described below.

[0176] Step 5(2).1: Global smoothing

[0177] Global smoothing is performed on the point clouds before and after sampling. The program of this invention uses the moving least squares method built into the point cloud library PCL, and the relevant parameters can be customized and adjusted by the user.

[0178] Step 5(2).2: Poisson surface fitting

[0179] The Poisson surface fitting is performed on the point cloud before and after smoothing. The program of this invention uses the Poisson surface fitting function built into the point cloud library PCL, and the relevant parameters can be customized and adjusted by the user.

[0180] Step 5(2).3: Calculate the volume of the fitted surface before and after sampling using the mesh method.

[0181] This part is the same as described in the first embodiment above, except that the input is replaced with the fitted surface before and after sampling. Then, the infinitesimal element is projected onto the xOy plane, the z value is taken as the height information and the difference is calculated. The difference result is multiplied by the area of ​​the infinitesimal element and the infinitesimal elements are superimposed to obtain the point cloud volume calculation result based on the fitted surface before and after sampling.

[0182] Step 6: Display and save the calculation results

[0183] In addition to the real-time output of key parameters at each stage of the above steps, the program of this invention will uniformly output the sampling estimation result and relevant indicators of the reliability of the result after the calculation is completed (such as the average number of nearest neighbors of the grid and the proportion of grids containing valid data points in the same grid before and after sampling when calculating directly using the grid method). Furthermore, if the grid method is used for direct calculation, the user can choose to display or save the point cloud of the grid-filled results before and after sampling (e.g., ...). Figure 15 (as shown) and the point cloud of the mesh result after subtraction (as shown) Figure 16 As shown), if the Poisson surface fitting method is used, the Poisson surface reconstruction results before and after sampling can be displayed or saved (e.g., Figure 17 (As shown).

Claims

1. A method for estimating sample size based on a surveillance binocular camera and a close-up monocular camera, the technical features of which include three parts, as follows: The first part involves the 3D reconstruction of point cloud data at sampling points and its classification processing. Specifically: First, based on the principle of binocular stereo vision, 3D point cloud data before and after sampling are obtained from image information from a monitoring binocular camera. Second, the soil condition near the sampling points is reassessed after sampling, and the soil is subdivided into natural soil, overturned soil, and subsidence soil, considering the impact of changes in soil compaction caused by the sampling operation on the estimation of sample volume. Finally, the point cloud data corresponding to different compaction states of soil are classified using point cloud data features and ideal sampling trajectories. The second part involves determining the equivalent volume ratio coefficients of soils in different compaction states at sampling points. Specifically, this involves: first, manually dividing and extracting regions of different compaction states from soil images captured by a close-up monocular camera; second, performing digital image processing on each sub-image and using image features to reasonably evaluate the porosity of soils in different compaction states; and finally, based on the ideal theoretical relationship between porosity and volume, obtaining the equivalent volume ratio coefficients between soils in different compaction states from the porosity. The third part is the estimation of the sampling volume and the presentation of related results. Specifically, the volume of soil in different compaction states is obtained by subtracting the point cloud data before and after sampling. The volume of each soil in different compaction states is multiplied by the equivalent volume ratio coefficient of its corresponding compaction state and then calculated by superposition to obtain the sampling volume estimation result. Finally, the parameters, images and key information of the point cloud related to the results are displayed and presented.

2. The sampling quantity estimation method based on a surveillance binocular camera and a close-up monocular camera according to claim 1, wherein the first part of its technical features—the 3D reconstructed point cloud data of the sampling points and its classification processing—is specifically implemented as follows: First, the binocular camera reconstructs the global 3D point cloud before and after sampling based on the principle of binocular vision. The sampling pit neighborhood in the global point cloud is automatically cropped. Then, outliers in the sampling pit neighborhood point cloud are filtered and removed to prevent outliers from interfering with subsequent volume calculations. Secondly, a plane fitting is performed on the point cloud before sampling. The standard is to make the origin of the fitting plane coincide with the centroid of the point cloud before sampling and to rotate the fitting plane to be parallel to the xOy plane. The direction and origin position of the point cloud coordinates before and after sampling are corrected. Then, the corresponding point cloud data at the same position before and after sampling are subtracted in the height direction. Then, the difference in height at each location is used to determine the soil height. If the height of the point cloud after sampling at a certain location is greater than the height of the point cloud before sampling at that location, the soil height at that location has increased due to the sampling operation. Therefore, the point cloud data at that location is assigned to the outward turning part of the pit caused by sampling. Conversely, if the height of the point cloud after sampling at a certain location is less than the height of the point cloud before sampling at that location, the soil height at that location has decreased due to the sampling operation. Therefore, the point cloud data at that location is assigned to the soil pit after sampling. Finally, since the pit includes the natural soil that was sampled and the subsidence soil formed by the sliding of the soil around the sampling pit, the point cloud inside the actual soil pit obtained in the previous step is compared with the sampling trajectory point cloud under the ideal sampling condition, which does not consider the sliding of the soil around the pit into the pit after sampling. The difference between the two is used to obtain the subsidence soil part above the ideal sampling pit, thus realizing the division of the natural soil part and the subsidence soil part inside the actual pit.

3. The method for estimating the sampling volume based on a binocular monitoring camera and a close-up monocular camera according to claim 1, specifically the determination of the equivalent volume ratio coefficient of soil at different loose states at sampling points in the second part of its technical features, is implemented as follows: First, a close-up monocular camera captured magnified images of the soil in the sampling area. Secondly, after the image is standardized, the loose soil and solid soil parts in the image are divided by human screening and judgment to obtain sub-images of each part; Then, the texture and color features of each part of the loose soil and solid soil sub-images are calculated from the image information, and the porosity of each part of the loose soil and solid soil is obtained according to empirical relationships. Finally, based on the theoretical relationship between void ratio and volume, the equivalent volume ratio coefficients of soils in different loose states were obtained.

4. The method for estimating the sampling volume based on a binocular surveillance camera and a close-up monocular camera according to claim 1, wherein the second part of its technical features is the determination of the equivalent volume ratio coefficient of soil in different loose states at the sampling points, another technical route is to test the soil quality of the sampling area and then obtain the corresponding loose volume ratio of the soil by looking up a table: the soil quality test is achieved by performing digital image processing on the soil image taken by the close-up monocular camera and extracting texture and color features, or by directly using a geological spectrometer to obtain the spectral information of the sampling area by scanning and further analyzing to obtain the soil composition of the sampling area.

Images