Method and system for three-dimensional reconstruction of grain pile surface based on millimeter wave radar

Abstract

This invention discloses a method and system for three-dimensional reconstruction of grain pile surfaces based on millimeter-wave radar. The method comprises: S1) three-dimensional point cloud acquisition and preprocessing; S2) point cloud density analysis and missing region identification; S3) local slope field and dynamic angle of repose estimation; S4) surface completion and optimization under physical constraints; and S5) surface meshing and volume calculation. The system includes a millimeter-wave radar three-dimensional scanning module, a point cloud preprocessing module, a dynamic angle of repose estimation module, a physically constrained surface reconstruction module, and a volume calculation module. It faithfully reproduces measured data while physically and reasonably completing obscured areas, avoiding excessive smoothing or artificial assumptions, and significantly improving the accuracy of volume estimation.

Description

Technical Field

[0001] This invention belongs to the field of grain storage measurement and three-dimensional reconstruction technology, specifically involving a method and system for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar. Background Technology

[0002] Accurate measurement of grain reserves is a crucial foundation for national food security, storage scheduling, and asset management. The accurate estimation of grain pile volume and reserves directly impacts inventory management and financial auditing, and is also critical for operational aspects such as grain condition monitoring, ventilation control, and grain storage safety assessment. In actual storage environments, to ensure grain quality and storage safety, silos are typically equipped with structures such as temperature-measuring cables, central pressure-reducing pipes, vertical ventilation ducts, and operating platforms. While these facilities function, they also physically obstruct the grain pile surface, leading to data blind spots in traditional measurement methods and consequently affecting the accuracy of 3D reconstruction and volume calculation. Therefore, achieving high-precision and robust 3D reconstruction of the grain pile surface in complex silo environments with pervasive dust and severe obstruction is a pressing technical challenge that needs to be addressed in this field.

[0003] Currently, grain pile volume measurement technologies mainly fall into two categories: contact and non-contact. Contact methods, such as manual depth sounding, capacitive sensing, and pressure / weight sensing, are simple to implement but easily affected by material characteristics, installation location, and environmental factors, making it difficult to obtain complete three-dimensional surface morphology. With the development of sensing technology, non-contact optical methods, represented by laser ranging, structured light, and depth cameras, have become a research hotspot and are gradually being applied to grain pile surface reconstruction. However, inside silos with high dust concentrations, optical signals are subject to severe scattering and absorption, resulting in sparse point clouds, increased noise, and even complete data loss. Especially when there are obstructions such as temperature measuring cables and vertical pipes, strip-shaped or regional "holes" often appear in the collected point clouds. Although some studies have attempted to complete the data through multi-viewpoint registration or interpolation based on simple geometric models (such as cones), the repetitive structure and weak texture features inside the silo easily lead to registration failure or reconstruction results deviating from the actual physical morphology.

[0004] Compared to optical sensors, millimeter-wave radar exhibits unique advantages in warehouse environment sensing due to its strong penetration capabilities through dust and fog, its insensitivity to lighting conditions, and its ability to operate in all weather conditions. Existing millimeter-wave radar-based grain warehouse measurement methods mainly include single-point / single-line ranging, compressed sensing and signal enhancement, machine learning, and deep learning. Among these, the use of frequency-modulated continuous wave (FMCW) radar for liquid level measurement has verified the high-precision potential of millimeter-wave radar; and some progress has been made in classifying storage tank environments by analyzing radar curve characteristics. However, the above methods can only acquire a limited number of measurement points and cannot reconstruct a complete three-dimensional geometric surface. Compressed sensing-based grain warehouse measurement methods enhance radar echo quality through sparse reconstruction, but this method relies on the assumption of signal sparsity, limiting its performance under complex grain surface morphology. While the grain level estimation method utilizing radar spectrogram features and combining machine / deep learning has verified the measurement potential of millimeter-wave radar, it still has significant limitations: single-point measurement cannot reconstruct three-dimensional geometry; compressed sensing methods have insufficient performance under complex surface morphology; and data-driven methods rely heavily on a large amount of labeled data, have weak generalization ability when conditions such as grain pile slope and humidity change, and generally lack modeling of the physical properties of particulate materials.

[0005] In recent years, millimeter-wave radar 3D reconstruction and enhancement technologies have developed rapidly. Examples include using radar SLAM to stitch together multi-frame point clouds, proposing cross-modal supervised radar point cloud enhancement models, and applying diffusion models to millimeter-wave point cloud super-resolution. These works demonstrate the potential of millimeter-wave radar to approach the 3D perception accuracy of optical sensors. However, directly applying these technologies to grain storage scenarios still faces fundamental challenges: purely data-driven methods ignore the inherent physical laws governing the accumulation behavior of particulate materials. The surface morphology of grain piles is strictly constrained by the "angle of repose," which is the maximum slope angle that bulk materials can form while remaining stable under their own weight. This physical law provides strong prior knowledge for the reconstruction of occluded areas. Existing methods fail to effectively utilize this domain knowledge, potentially resulting in physically unreasonable surface morphologies (e.g., excessively steep slopes) in areas with missing data.

[0006] In summary, existing methods have significant shortcomings at the sensor level, reconstruction algorithm level, and application level for the specific scenario of grain silos. At the sensor level, traditional millimeter-wave radar has low angular resolution, making it difficult to acquire dense 3D point clouds. At the reconstruction algorithm level, existing point cloud completion methods mostly employ pure geometric interpolation (ignoring physical constraints) or pure data-driven methods (relying on training data distribution), lacking effective fusion of physical models and observational data. At the application level, the unique internal occlusion structure of grain silos and the angle of repose characteristics of granular materials have not yet been systematically modeled and utilized. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a method and system for three-dimensional reconstruction of grain pile surfaces based on millimeter-wave radar. This method integrates prior physical knowledge, overcomes limitations caused by occlusion and low angular resolution, and achieves high-precision three-dimensional reconstruction of grain pile surfaces.

[0008] The specific plan is as follows: The method for three-dimensional reconstruction of grain pile surfaces based on millimeter-wave radar includes the following steps. S1): 3D point cloud acquisition and preprocessing; a narrow-beam millimeter-wave radar combined with a 3D mechanical scanning turntable is used to acquire the original 3D point cloud of the grain pile surface, and spatial point cloud filtering is performed based on the geometric prior of the grain silo to remove outliers, resulting in a preprocessed set of observation points. P raw ; S2): Point cloud density analysis and missing region identification; calculation of the observed point cloud set. P raw The local spatial density of each point in the cloud is used to divide the point cloud space into reliable observation regions Ω based on the density distribution. obs and missing region Ω miss ; S3): Local slope field and dynamic angle of repose estimation; in the reliable observation area Ω obs Within the grain pile, the surface normal vector and slope angle are calculated through local surface fitting to construct a local slope field. Based on the statistical characteristics of the slope field, the dynamic angle of repose θ is adaptively estimated. r ; S4): Surface completion and optimization under physical constraints; constructing an energy functional model that integrates data fidelity terms, surface smoothing terms, and angle of repose constraint terms. E(S) With the reliable observation area Ω obs Using the data as the anchor point, the dynamic angle of repose θ r As a physical constraint, the energy functional model is solved through numerical optimization in the missing region Ω. miss Generate a three-dimensional surface height field that conforms to physical laws. S(x,y) ; S5): Surface meshing and volume calculation; optimizing the continuous height field. S(x,y) The grain pile is discretized into a triangular mesh model, and the total volume of the grain pile is calculated based on the triangular mesh model. V .

[0009] In step S1), the narrow-beam millimeter-wave radar employs a frequency-modulated continuous wave system, compressing the beamwidth through a cascaded waveguide slot antenna and Fresnel lens; the three-dimensional mechanical scanning turntable operates according to a preset azimuth step size. Dth and pitch step Perform a scan and measure the radial distance at each scan angle. R ijAnd convert it to Cartesian coordinates.

[0010] radial distance R ij The following method is used to obtain the discrete sequence by sampling the beat signal, obtaining the preliminary spectrum by using Fast Fourier Transform, and then refining the spectrum within the target frequency band by using Chirp-Z Transform. The radial distance is calculated based on the linear relationship between the beat frequency and the distance. The Chirp-Z Transform samples along a spiral path in the z-plane. By setting the initial sampling point A and the sampling interval control parameter W, the distance is calculated within the frequency band of interest. f L ，f H Internal implementation M Uniform sampling improves frequency resolution to d ( f CZT )=( f H - f L ) / M .

[0011] Local spatial density in step S2) r k The calculation formula is r k = N k / (4 / 3 πr 3 ),in, N k Indicates that in p k Centered on, with radius r The number of points within the spherical neighborhood; p k It represents any point in the observation point cloud, and its neighborhood radius. r Average point spacing d avg Associated, set as r=α · d avg ,in α The value is [2,3]; density threshold r th The lower quartile of the distribution of all point cloud density values ​​satisfies r k ≥r th The points are included in the reliable observation area Ω obs The rest are assigned to the missing region Ω. miss .

[0012] The method for estimating the local slope angle in step S3) is as follows: for the reliable observation area Ω obs Each point inside p k The surface normal vector is estimated using principal component analysis with its neighborhood points. n k Surface normal vector n k The local slope angle is obtained by constructing the local covariance matrix and finding the eigenvector corresponding to the minimum eigenvalue. α k Defined as normal vector n k and vertical direction e z =(0,0,1) T The included angle between them: α k =arccos (| n k · e z | / ‖ n k || ); The dynamic angle of repose θ r The estimation method is as follows: the reliable observation area Ω obs The horizontal plane is divided into grids, and the slope angle within each grid cell is calculated. α k After removing abnormally large slope values ​​corresponding to the warehouse wall structure, the statistical measure is used to calculate the weighted average of the slope angles of all grids, with the number of points in each grid as the weight, as the dynamic angle of repose θ. r The estimated value; the dynamic angle of repose θ r It is used to adaptively reflect the influence of changes in grain variety, moisture content, and accumulation history on the angle of repose.

[0013] The energy functional model constructed in step S4) E(S) Represented as: E(S)=E data (S)+λE smooth (S)+μ E repose (S) ,in, E data (S) For data fidelity, constrain the reconstructed surface within the reliable observation region Ω. obs The internal structure is consistent with the measured point cloud. E smooth (S) For smoothing terms, constrain the geometric smoothness of the reconstructed surface; E repose(S) As the angle of repose constraint term, the constrained reconstructed surface is in the missing region Ω. miss The local slope within the slope does not exceed the dynamic angle of repose θ r ; l and m These are the weighting coefficients; The data fidelity item E data (S) Represented as: E data (S) =∑ pk∈Ωobs oh k [ S(x k ,y k )-z k ] 2 ;in, oh k =p k / r max As weight, r k For point p k Local density, r max This represents the maximum density value. The smoothing item E smooth (S) Represented as: ,in β The weights are second-order smoothing weights; The angle of repose constraint E repose (S) Represented as: E repose (S)= ∑ (i,j)∈Ωmiss Ψ(‖▽ S ij ||-tan i r (l) ); Where Ψ(·) is a differentiable soft threshold penalty function; Defined as: when u ≤0, Ψ( u ) = 0; when 0 < u ≤γ,Ψ( u ) =u 2 / 2γ; when u >γ,Ψ(u ) =u -γ / 2; γ is the preset threshold.

[0014] The method for numerically optimizing the energy functional model in step S4) is as follows: the alternating direction multiplier method is used, and auxiliary variables are introduced. u=▽S The original problem is transformed into a constrained optimization problem. An augmented Lagrangian function is constructed, and then the height field is updated iteratively. S Auxiliary variables u and Lagrange multipliers or The convergence process continues until the convergence condition is met; the convergence condition considers the original residual, the dual residual, and the rate of energy change simultaneously, and occurs when all three are below a set threshold. σ= 10 -4 Stop iterating when the time comes.

[0015] The method for volume calculation in step S5) is as follows: The optimized continuous height field... S(x,y) Discretize the data into a regular grid on the horizontal plane and connect them to form a triangular grid model. Total volume of grain pile V By accumulating the values ​​of each triangular facet and the reference plane at the bottom of the silo, z=z ref The volume of the prism formed is obtained by the following formula: The matrix is ​​composed of triangular facets Δ k The coordinates of the three vertices constitute the equation.

[0016] A three-dimensional reconstruction system for grain pile surface based on millimeter-wave radar, including The millimeter-wave radar 3D scanning module includes a narrow-beam millimeter-wave radar and a 2D mechanical scanning gimbal, used to acquire the original 3D point cloud of the grain pile surface. The point cloud preprocessing module is used to perform spatial filtering on the original point cloud based on the geometric prior of the grain warehouse, and remove outliers; the density analysis and missing region identification module is used to calculate the local spatial density of the point cloud and divide the reliable observation area and the missing region. The dynamic angle of repose estimation module is used to estimate the slope field through local surface fitting within a reliable observation area and adaptively calculate the dynamic angle of repose. The physical constraint surface reconstruction module is used to construct an energy functional model that integrates data fidelity, smoothing, and repose angle constraints, and generates a complete three-dimensional surface height field through numerical optimization. The volume calculation module is used to discretize the height field into a triangular mesh and calculate the total volume of the grain pile.

[0017] The narrow-beam millimeter-wave radar in the three-dimensional scanning module of the millimeter-wave radar adopts a waveguide slot antenna and Fresnel lens cascade design to compress the beam width, and adopts a frequency-modulated continuous wave system to achieve range dimension super-resolution measurement through Chirp-Z transform.

[0018] Compared with existing technologies for three-dimensional reconstruction and volume measurement of grain pile surfaces, the millimeter-wave radar-based three-dimensional reconstruction method for grain pile surfaces proposed in this invention has the following significant advantages and beneficial effects: (1) Robust reconstruction driven by physical information overcomes the limitations of data-driven methods. Existing methods rely on a large amount of labeled data or pure geometric interpolation, which has weak generalization ability and easily ignores the physical characteristics of particulate materials, resulting in unreasonable slope of the shading area. This invention is the first to explicitly embed the physical constraint of the angle of repose into the energy functional model and dynamically estimate the actual angle of repose to guide surface completion, so that the reconstruction result is both geometrically continuous and meets the static stability. Even when training data is scarce or the grain type, moisture content and stacking method change, it still maintains high-precision physical consistency reconstruction, and the generalization ability and engineering practicality are significantly improved.

[0019] (2) Dynamic angle of repose adaptive estimation to adapt to changing storage conditions. Existing methods use fixed empirical values ​​as the upper limit of slope, which cannot reflect the actual differences in the angle of repose under different grain types, moisture contents, and stacking methods. This invention obtains the slope field through local surface fitting within a reliable observation area and dynamically infers the true angle of repose based on a statistical weighting strategy. No manual calibration or pre-stored material parameters are required; the system automatically adapts to changes in storage conditions, significantly improving scene adaptability and measurement accuracy.

[0020] (3) Energy functional optimization that integrates data fidelity, smoothing, and physical constraints achieves physically reasonable surface completion. Existing methods often employ simple interpolation or global shape assumptions to address point cloud voids caused by temperature measuring cables, pressure reducing pipes, etc., which easily introduce systemic biases. This invention constructs an energy functional model containing data fidelity terms, smoothing terms, and repose angle constraint terms, and uses ADMM for efficient solution. The data fidelity terms are weighted according to the local density of the point cloud, the smoothing terms ensure natural transitions, and the repose angle constraint terms limit the upper limit of the slope of the missing region using a soft threshold penalty function. While faithfully reflecting the measured data, it physically and reasonably completes the occluded areas, avoiding excessive smoothing or artificial assumptions, and significantly improving the accuracy of volume estimation. Attached Figure Description

[0021] Figure 1 This is a comparison of point cloud filtering before and after different operating conditions: warehouse entry, warehouse closing, and warehouse exit.

[0022] Figure 2 This is a spatial distribution map of point cloud density under three operating conditions: warehouse entry, warehouse closing, and warehouse exit.

[0023] Figure 3The spatial distribution of the local angle of repose for an initial grain load of 1216 kg under the conditions of grain entering and leaving the warehouse.

[0024] Figure 4 The iterative convergence curves of the energy functional optimization process under different grain loading conditions are shown.

[0025] Figure 5 A quantitative error table is provided for nine sets of experiments under three working conditions and three grain loading amounts.

[0026] Figure 6 This is a comparison table of relative errors in volume estimation.

[0027] Figure 7 The error distribution heatmap and profile curves for each method under the outbound working condition are shown.

[0028] Figure 8 This is a graph showing the sensitivity analysis of parameters. Detailed Implementation

[0029] The technical solutions in the embodiments of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the implementation of the present invention, and not all of it. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0030] This invention belongs to the field of grain storage measurement and three-dimensional reconstruction technology. Specifically, it relates to a millimeter-wave radar three-dimensional reconstruction method for grain pile surface based on physical information and dynamic repose angle constraints. It is applicable to high-precision three-dimensional reconstruction and volume measurement of grain pile surface under conditions of dust interference and internal structural obstruction in closed storage environments such as silos.

[0031] This embodiment uses a scaled-down cylindrical steel silo model as the experimental platform. The silo consists of six sections connected by prefabricated slots, with an inner diameter of 1.207 m, a depth of 2.672 m, and a flat bottom. To simulate the working environment above the silo, a horizontal grid-like steel frame is erected above the silo to mount and fix sensors. The experimental material is dry bulk wheat with a standard bulk density of 750 kg / m³, and the average static angle of repose measured using the standard tilting method is 24.6° ± 1.2°.

[0032] The radar measurement system employs a narrow-beam millimeter-wave radar with the following core parameters: center frequency 80 GHz, modulation bandwidth 3 GHz (theoretical range resolution 5 cm), and modulation period 5 ms. It utilizes a two-stage beam compression system consisting of a waveguide slot array and a Fresnel lens, achieving a half-power beamwidth of 2° and a maximum range of 50 m. The radar is mounted on a high-precision, three-dimensional electronically controlled pan-tilt unit on the steel frame of the warehouse roof, with a pointing accuracy of 0.01°. Omnidirectional scanning is achieved through program control, with an azimuth range of 0° to 360°, an elevation range of -90° to 90°, and an angle step of 0.1°.

[0033] To obtain a high-precision evaluation benchmark, a multi-site data fusion strategy was adopted: radar sensors were placed sequentially on four precisely calibrated fixed stations on the top of the silo to scan and obtain four sets of independent point clouds. Then, the iterative nearest point algorithm was used to register and fuse the multi-site cloud data to generate a complete and detailed three-dimensional model of the grain pile surface as the ground truth, with an average registration error of less than 2 mm.

[0034] This embodiment achieves three-dimensional reconstruction and volume measurement of the grain pile surface according to the following steps.

[0035] S1): 3D point cloud acquisition and preprocessing; a narrow-beam millimeter-wave radar combined with a 3D mechanical scanning turntable is used to acquire the original 3D point cloud of the grain pile surface, and spatial point cloud filtering is performed based on the geometric prior of the grain silo to remove outliers, resulting in a preprocessed set of observation points. P raw ; In step S1), the narrow-beam millimeter-wave radar employs a frequency-modulated continuous wave system, compressing the beamwidth through a cascaded waveguide slot antenna and Fresnel lens; the three-dimensional mechanical scanning turntable operates according to a preset azimuth step size. Dth and pitch step Perform a scan and measure the radial distance at each scan angle. R ij And convert it to Cartesian coordinates.

[0036] radial distance R ij The following method is used to obtain the discrete sequence by sampling the beat signal, obtaining the preliminary spectrum by using Fast Fourier Transform, and then refining the spectrum within the target frequency band by using Chirp-Z Transform. The radial distance is calculated based on the linear relationship between the beat frequency and the distance. The Chirp-Z Transform samples along a spiral path in the z-plane. By setting the initial sampling point A and the sampling interval control parameter W, the distance is calculated within the frequency band of interest. f L ，f H Internal implementation M Uniform sampling improves frequency resolution to d (f CZT )=( f H - f L ) / M .

[0037] Specifically, it includes: (1) Radar scanning and ranging Start the 3D electronically controlled pan-tilt unit and operate it in preset azimuth step sizes. Dth =0.1 ° and pitch step =0.1 ° Perform point-by-point scanning. At each scanning angle Below, the millimeter-wave radar transmits a linear frequency modulated continuous wave signal, with the starting frequency of the transmitted signal being... f c = 80GHz, modulation bandwidth B =3GHz, modulation period T c =5ms; After receiving the echo, it is mixed with the transmitted signal to obtain the beat signal. The beat signal is then numerically sampled to obtain a discrete sequence. x(n) First, a preliminary spectrum is obtained using a Fast Fourier Transform (FFT). Then, a Chirp-Z Transform is applied within the frequency band corresponding to the target distance to refine the spectrum. The Chirp-Z Transform parameters are set as follows: initial sampling point A=1, sampling interval control parameter... W =1, in the frequency band of interest[[ f L ，f H [Inside] M= 4096 points uniform sampling; based on beat frequency f IF With distance R linear relationship f IF = 2 BR / cT c The radial distance was calculated. R ij .

[0038] (2) Coordinate transformation Based on the measured radial distance R ij and current scanning angle Convert polar coordinates to Cartesian coordinates in the radar coordinate system using the following formula: Traverse all scanning angles, azimuth step count. N θ =3600, pitch steps Obtain the original point cloud set P raw In this embodiment, approximately 80,000 original points were obtained under different working conditions.

[0039] (3) Geometric constraint filtering Based on the prior geometric model of the silo, the cylindrical silo wall radius is 0.6035 m, and the silo bottom plane height z ref =0, for each point p ij Validity assessment: If the point is located inside the cylindrical surface of the warehouse wall and its height is within a reasonable range, it is retained; otherwise, it is removed as an outlier. For example... Figure 1 As shown, the comparison before and after filtering is presented under different operating conditions: warehouse entry, warehouse closing, and warehouse exit. The outlier removal rate reaches 57%~63%, the average number of effective point clouds is approximately 31806 points, and the surface coverage reaches 68.03%. The preprocessed point cloud is denoted as... ,in N ≈31806.

[0040] S2): Point cloud density analysis and missing region identification; calculation of the observed point cloud set. P raw The local spatial density of each point in the cloud is used to divide the point cloud space into reliable observation regions Ω based on the density distribution. obs and missing region Ω miss ; Local spatial density in step S2) r k The calculation formula is r k = N k / (4 / 3 πr 3 ),in, N k Indicates that in p k Centered on, with radius r The number of points within the spherical neighborhood; p k It represents any point in the observation point cloud, and its neighborhood radius. r Average point spacing d avg Associated, set as r=α · d avg ,in α The value is [2,3]; density threshold r th The lower quartile of the distribution of all point cloud density values ​​satisfies r k ≥r thThe points are included in the reliable observation area Ω obs The rest are assigned to the missing region Ω. miss .

[0041] The specific steps for point cloud density analysis and missing region identification are as follows: (1) Calculate local density In this embodiment, d avg g≈0.86cm, set the neighborhood radius r=α · d avg ,Pick α =3, therefore... r =2.85cm. For each point p k The search is centered on it with a radius of [missing information]. r The number of points in the spherical neighborhood N k According to the formula r k = N k / (4 / 3 πr 3 Calculate the local spatial density.

[0042] (2) Divide the area like Figure 2 As shown, the spatial distribution of point cloud density is illustrated under three operating conditions: grain entering the warehouse, leveling the warehouse, and grain exiting the warehouse. Light blue indicates missing areas. The results show that under the grain entering condition, low-density areas are mainly concentrated at the top of the conical grain pile and the steep edges; under the leveling condition, the density distribution is generally uniform, but density decreases at local edges near the warehouse wall; under the grain exit condition, the point cloud density decreases significantly in the concave slope area formed by unloading, presenting clear strip-shaped missing areas. The density distribution map highly matches the morphological characteristics of the actual grain pile.

[0043] S3): Local slope field and dynamic angle of repose estimation; in the reliable observation area Ω obs Within the grain pile, the surface normal vector and slope angle are calculated through local surface fitting to construct a local slope field. Based on the statistical characteristics of the slope field, the dynamic angle of repose θ is adaptively estimated. r ; The method for estimating the local slope angle in step S3) is as follows: for the reliable observation area Ω obs Each point inside p k The surface normal vector is estimated using principal component analysis with its neighborhood points. n k Surface normal vector n k It is obtained by constructing the local covariance matrix and finding the eigenvector corresponding to the minimum eigenvalue; Local covariance matrix ,in Let the neighborhood centroid be denoted by . Find . C k The eigenvector corresponding to the smallest eigenvalue is used as the normal vector. n k ; Local slope angle α k Defined as normal vector n k and vertical direction e z =(0,0,1) T The included angle between them: α k =arccos (| n k · e z | / ‖ n k || ); The dynamic angle of repose θ r The estimation method is as follows: the reliable observation area Ω obs The horizontal plane is divided into grids, and the slope angle within each grid cell is calculated. α k After removing abnormally large slope values ​​corresponding to the warehouse wall structure, the statistical measure is used to calculate the weighted average of the slope angles of all grids, with the number of points in each grid as the weight, as the dynamic angle of repose θ. r The estimated value; the dynamic angle of repose θ r It is used to adaptively reflect the influence of changes in grain variety, moisture content, and accumulation history on the angle of repose.

[0044] In this embodiment, Ω obs The horizontal plane is divided into 10cm × 10cm grids. The slope angle of all points within each grid is calculated, and obvious outliers, such as points on the silo wall with an angle greater than 45°, are removed. The weighted average of the slope angles of all grids is calculated using the number of points within each grid as the weight, and this average is taken as the dynamic angle of repose θ of the current grain pile. r .

[0045] like Figure 3 As shown, the spatial distribution of the local angle of repose of the initial 1216 kg grain load under loading and unloading conditions is illustrated. The value fluctuates between 22.5° and 27.3°, reflecting the local accumulation differences in the grain pile caused by the loading and unloading process. Under both conditions, the weighted average dynamic angle of repose is 24.19° and 24.53°, respectively, with deviations from the static angle of repose measurement of 24.9° for the experimental wheat being less than 3%, thus verifying the rationality and accuracy of the strategy of dynamically inferring the angle of repose from the observation data.

[0046] S4): Surface completion and optimization under physical constraints; constructing an energy functional model that integrates data fidelity terms, surface smoothing terms, and angle of repose constraint terms. E(S) With the reliable observation area Ω obs Using the data as the anchor point, the dynamic angle of repose θ r As a physical constraint, the energy functional model is solved through numerical optimization in the missing region Ω. miss Generate a three-dimensional surface height field that conforms to physical laws. S(x,y) ; An energy functional model is constructed to formulate the surface reconstruction problem as a height field defined on the horizontal domain Ω. S(x,y) The energy minimization problem. The energy functional model constructed in step S4). E(S) Represented as: E(S)=E data (S)+λE smooth (S)+μE repose (S) ,in, E data (S) For data fidelity, constrain the reconstructed surface within the reliable observation region Ω. obs The internal structure is consistent with the measured point cloud. E smooth (S) For smoothing terms, constrain the geometric smoothness of the reconstructed surface; E repose (S) As the angle of repose constraint term, the constrained reconstructed surface is in the missing region Ω. miss The local slope within the slope does not exceed the dynamic angle of repose θ r ; l and m These are the weighting coefficients; The data fidelity item E data (S) Represented as: E data (S) =∑ pk∈Ωobs oh k [ S(x k ,y k )-z k ] 2 ;in, oh k =p k / r max As weight, rk For point p k Local density, r max The maximum density value is used; different observation points are assigned confidence levels based on the point cloud density. The smoothing item E smooth (S) Represented as: ,in β The weights are second-order smoothing weights; in this embodiment, we take... β =1.

[0047] The angle of repose constraint E repose (S) Represented as: E repose (S)= ∑ (i,j)∈Ωmiss Ψ(‖▽ S ij ||-tan i r (l) ); Where Ψ(·) is a differentiable soft threshold penalty function; Defined as: when u ≤0, Ψ( u ) = 0; when 0 < u ≤γ,Ψ( u ) =u 2 / 2γ; when u >γ,Ψ( u ) =u -γ / 2; γ is a preset threshold. In this embodiment, γ is 0.2, and the smoothing weight is... l =2.0, physical constraint weight m =20.0.

[0048] The method for numerically optimizing the energy functional model in step S4) is as follows: the alternating direction multiplier method is used, and auxiliary variables are introduced. u=▽S The original problem is transformed into a constrained optimization problem. An augmented Lagrangian function is constructed, and then the height field is updated iteratively. S Auxiliary variables u and Lagrange multipliers or The convergence process continues until the convergence condition is met; the convergence condition considers the original residual, the dual residual, and the rate of energy change simultaneously, and occurs when all three are below a set threshold. σ= 10 -4 Stop iterating when the time comes.

[0049] Construct the augmented Lagrangian function: minS,u E data (S) +λ E smooth ( u )+μ E repose ( u Iterative update steps: Update S Solving the linear least squares problem using the conjugate gradient method; updating... u Solving the proximal operator problem involves soft thresholding operations; updating multipliers. or Penalty parameters r Press min( r max 1.05 r Gradually increase, initially r =1.0; r max =100. The convergence condition is that the original residual <10. -4 Dual residual <10 -4 And the relative rate of change of energy is <10 -4 .

[0050] like Figure 4 As shown, the iterative convergence curves of the energy functional optimization process under different grain loading conditions are presented. The results show that, in typical scenarios, such as the grain loading condition of 1216 kg, the convergence condition can be met after 130-160 iterations using the ADMM solver.

[0051] Optimization yields continuous height field S(x,y) Its domain covers the entire horizontal cross-section of the silo, and the silo diameter is 1.207m.

[0052] S5): Surface meshing and volume calculation; optimizing the continuous height field. S(x,y) The grain pile is discretized into a triangular mesh model, and the total volume of the grain pile is calculated based on the triangular mesh model. V .

[0053] The method for volume calculation in step S5) is as follows: The optimized continuous height field... S(x,y) Discretize the height field using a regular grid on the horizontal plane. S(x,y) On the horizontal plane d x =d y = 1 cm Discretize the regular grid and connect adjacent grid points to form a triangular grid model. ; N t There are approximately 28,000 triangular facets.

[0054] Using the bottom plane z=z ref Using the reference plane as a reference, the total volume of the grain pile is obtained by summing the volumes of the prism formed by each triangular facet and the reference plane. V , The matrix is ​​composed of triangular facets Δ k The coordinates of the three vertices constitute the equation, i.e. ( x i ,y i ,z i ) is a triangular facet Δ k The coordinates of the three vertices are given. In this embodiment, for the grain loading condition of 1216kg, the calculated volume is 1.628 m³, and the relative error between the actual value of 1.621 m³ and the actual value is 0.42%.

[0055] A three-dimensional reconstruction system for grain pile surface based on millimeter-wave radar, including The millimeter-wave radar 3D scanning module includes a narrow-beam millimeter-wave radar and a 2D mechanical scanning gimbal, used to acquire the original 3D point cloud of the grain pile surface. The point cloud preprocessing module is used to perform spatial filtering on the original point cloud based on the geometric prior of the grain warehouse, and remove outliers; the density analysis and missing region identification module is used to calculate the local spatial density of the point cloud and divide the reliable observation area and the missing region. The dynamic angle of repose estimation module is used to estimate the slope field through local surface fitting within a reliable observation area and adaptively calculate the dynamic angle of repose. The physical constraint surface reconstruction module is used to construct an energy functional model that integrates data fidelity, smoothing, and repose angle constraints, and generates a complete three-dimensional surface height field through numerical optimization. The volume calculation module is used to discretize the height field into a triangular mesh and calculate the total volume of the grain pile.

[0056] The narrow-beam millimeter-wave radar in the three-dimensional scanning module of the millimeter-wave radar adopts a waveguide slot antenna and Fresnel lens cascade design to compress the beam width, and adopts a frequency-modulated continuous wave system to achieve range dimension super-resolution measurement through Chirp-Z transform.

[0057] To verify the superiority of the method of this invention, the classical inverse distance weighted interpolation method (IDW) and the smooth reconstruction method based on moving least squares (MLS) were compared on the same dataset. Figure 5As shown, the quantitative errors of nine sets of experiments under three working conditions and three grain loading amounts are summarized, including root mean square error (RMSE), mean error (MAE), and maximum error (MaxError). The results show that the method of this invention achieves the lowest error in all tests, with an average RMSE reduction of 57.4% compared to IDW and 32.8% compared to MLS. Figure 6 As shown, the comparison of relative errors in volume estimation is illustrated. The relative error of the method described in this invention is less than 1% (maximum 0.63%) under all operating conditions, while the error ranges of the IDW and MLS methods are 1.79%~3.52% and 1.17%~2.73%, respectively. Especially in the outbound operating condition (768 kg) where data is severely lacking, the relative error of the method described in this invention is only 0.63%, significantly outperforming the comparative methods. Figure 7 The error distribution heatmaps and profile curves of each method under the outbound operating condition are presented. IDW and MLS produce large errors and physically unreasonable oscillations in the missing regions, while the surface generated by the method of this invention has a smooth transition and the profile is highly consistent with the true value.

[0058] like Figure 8 Parameter robustness analysis shows that when smoothing weights l ∈[0.5,5], physical constraint weight m When the resolution is in the range [10,50], the reconstruction error RMSE remains around 12.0 mm with a fluctuation of <20%, indicating that the algorithm is not sensitive to parameters and has good robustness. In terms of computational efficiency, the complete process takes an average of about 180 seconds, using an Intel Core i7-12700K, 16GB RAM, and MATLAB R2022b. The surface optimization stage accounts for about 80% of the computation, which meets the requirements of engineering applications.

[0059] Those skilled in the art should understand that the above parameters can be appropriately adjusted according to the actual grain warehouse size, radar system parameters, and grain varieties, all of which fall within the protection scope of this invention.

[0060] This embodiment describes the invention in detail using a scaled-down silo as an example, but the invention is equally applicable to actual large grain silos. In practical applications, the radar is installed at the center or off-center of the silo roof. The scanning range of the three-dimensional gimbal is adjusted according to the silo size, with an azimuth angle of 0° to 360° and an elevation angle determined by the silo height. Since actual grain silos are much taller, typically 10 to 30 meters, the radar ranging range needs to be increased accordingly. However, the Chirp-Z transform super-resolution ranging method used in this invention can still guarantee centimeter-level accuracy. The geometric priors of obstructions inside the grain silo, such as temperature measuring cables and pressure reducing pipes, such as their location and size, can be pre-calibrated and used for point cloud filtering and missing area identification. The dynamic estimation method for the angle of repose is directly applicable to actual grain piles without modification. Therefore, this invention has good scalability and practical application value.

[0061] The technical means disclosed in this invention are not limited to those disclosed in the above embodiments, but also include technical solutions composed of any combination of the above technical features. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of this invention, and these improvements and modifications are also considered within the scope of protection of this invention.

Claims

1. A method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar, characterized in that: Includes the following steps, S1): 3D point cloud acquisition and preprocessing; a narrow-beam millimeter-wave radar combined with a 3D mechanical scanning turntable is used to acquire the original 3D point cloud of the grain pile surface, and spatial point cloud filtering is performed based on the geometric prior of the grain silo to remove outliers, resulting in a preprocessed set of observation points. P raw ; S2): Point cloud density analysis and missing region identification; Calculate the set of observation points P raw The local spatial density of each point in the cloud is used to divide the point cloud space into reliable observation regions Ω based on the density distribution. obs and missing region Ω miss ; S3): Local slope field and dynamic angle of repose estimation; in the reliable observation area Ω obs Within the grain pile, the surface normal vector and slope angle are calculated through local surface fitting to construct a local slope field. Based on the statistical characteristics of the slope field, the dynamic angle of repose θ is adaptively estimated. r ; S4): Surface completion and optimization under physical constraints; constructing an energy functional model that integrates data fidelity terms, surface smoothing terms, and angle of repose constraint terms. E(S) With the reliable observation area Ω obs Using the data as the anchor point, the dynamic angle of repose θ r As a physical constraint, the energy functional model is solved through numerical optimization in the missing region Ω. miss Generate a three-dimensional surface height field that conforms to physical laws. S(x,y) ; S5): Surface meshing and volume calculation; optimizing the continuous height field. S(x,y) The grain pile is discretized into a triangular mesh model, and the total volume of the grain pile is calculated based on the triangular mesh model. V .

2. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 1, characterized in that: In step S1), the narrow-beam millimeter-wave radar employs a frequency-modulated continuous wave system, compressing the beamwidth through a cascaded waveguide slot antenna and Fresnel lens; the three-dimensional mechanical scanning turntable operates according to a preset azimuth step size. Δθ and pitch step Perform a scan and measure the radial distance at each scan angle. R ij And convert it to Cartesian coordinates.

3. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 2, characterized in that: radial distance R ij The following method is used to obtain the discrete sequence by sampling the beat signal, obtaining the preliminary spectrum by using Fast Fourier Transform, and then refining the spectrum within the target frequency band by using Chirp-Z Transform. The radial distance is calculated based on the linear relationship between the beat frequency and the distance. The Chirp-Z Transform samples along a spiral path in the z-plane. By setting the initial sampling point A and the sampling interval control parameter W, the distance is calculated within the frequency band of interest. f L ，f H Internal implementation M Uniform sampling improves frequency resolution to δ ( f CZT )=( f H - f L ) / M .

4. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 1, characterized in that: Local spatial density in step S2) ρ k The calculation formula is: ρ k = N k / (4 / 3 πr 3 ),in, N k Indicates in p k Centered on a cluster of observation points P raw The radius of each point in the middle is r The number of points within the spherical neighborhood; p k It represents any point in the observation point cloud, and its neighborhood radius. r Average point spacing d avg Associated, set as r=α·d avg ,in α The value is [2,3]; density threshold ρ th The lower quartile of the distribution of all point cloud density values ​​satisfies ρ k ≥ρ th The points are included in the reliable observation area Ω obs The rest are assigned to the missing region Ω. miss .

5. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 1, characterized in that: The method for estimating the local slope angle in step S3) is as follows: for the reliable observation area Ω obs Each point within p k The surface normal vector is estimated using principal component analysis with its neighborhood points. n k Surface normal vector n k The local slope angle is obtained by constructing the local covariance matrix and finding the eigenvector corresponding to the minimum eigenvalue. α k Defined as normal vector n k and vertical direction e z =(0,0,1) T The included angle between them: α k =arccos (| n k · e z | / ‖ n k ||); The dynamic angle of repose θ r The estimation method is as follows: the reliable observation area Ω obs The horizontal plane is divided into grids, and the slope angle within each grid cell is calculated. α k After removing abnormally large slope values ​​corresponding to the warehouse wall structure, the statistical measure is used to calculate the weighted average of the slope angles of all grids, with the number of points in each grid as the weight, as the dynamic angle of repose θ. r The estimated value; the dynamic angle of repose θ r It is used to adaptively reflect the influence of changes in grain variety, moisture content, and accumulation history on the angle of repose.

6. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 1, characterized in that: The energy functional model constructed in step S4) E(S) Represented as: E(S) = E data (S)+λE smooth (S)+μE repose (S) ,in, E data (S) For data fidelity, constrain the reconstructed surface within the reliable observation region Ω. obs The internal structure is consistent with the measured point cloud. E smooth (S) For smoothing terms, constrain the geometric smoothness of the reconstructed surface; E repose (S) As the angle of repose constraint term, the constrained reconstructed surface is in the missing region Ω. miss The local slope within the slope does not exceed the dynamic angle of repose θ r ; λ and μ These are the weighting coefficients; The data fidelity item E data (S) Represented as: E data (S) =∑ pk∈Ωobs ω k [ S(x k ,y k )-z k ] 2 ;in, ω k =ρ k / ρ max As weight, ρ k For point p k Local density, ρ max This represents the maximum density value. The smoothing item E smooth (S) Represented as: in β The weights are second-order smoothing weights; The angle of repose constraint E repose (S) Represented as: E repose (S)= ∑ (i,j)∈Ωmiss Ψ(‖▽ S ij ||-tan θ r (l) ); where Ψ(·) is a differentiable soft threshold penalty function, Defined as: when u ≤0, Ψ( u ) = 0; when 0 < u ≤γ,Ψ( u ) =u 2 / 2γ; when u >γ,Ψ( u ) =u -γ / 2; γ is the preset threshold.

7. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 1, characterized in that: The method for numerically optimizing the energy functional model in step S4) is as follows: the alternating direction multiplier method is used, and auxiliary variables are introduced. u=▽S The original problem is transformed into a constrained optimization problem. An augmented Lagrangian function is constructed, and then the height field is updated iteratively. S Auxiliary variables u and Lagrange multipliers η The convergence process continues until the convergence condition is met; the convergence condition considers the original residual, the dual residual, and the rate of energy change simultaneously, and occurs when all three are below a set threshold. σ= 10 -4 Stop iterating when the time comes.

8. The method for three-dimensional reconstruction of grain pile surface based on millimeter-wave radar according to claim 1, characterized in that: The method for volume calculation in step S5) is as follows: The optimized continuous height field... S(x,y) Discretize the data into a regular grid on the horizontal plane and connect them to form a triangular grid model. Total volume of grain pile V By accumulating the values ​​of each triangular facet and the reference plane at the bottom of the silo, z=z ref The volume of the prism formed is obtained by the following formula: The matrix is ​​composed of triangular facets Δ k The coordinates of the three vertices constitute the equation.

9. A three-dimensional reconstruction system for grain pile surface based on millimeter-wave radar, comprising the three-dimensional reconstruction method for grain pile surface based on millimeter-wave radar as described in any one of claims 1 to 8, characterized in that: It also includes The millimeter-wave radar 3D scanning module includes a narrow-beam millimeter-wave radar and a 2D mechanical scanning gimbal, used to acquire the original 3D point cloud of the grain pile surface. The point cloud preprocessing module is used to perform spatial filtering on the original point cloud based on the geometric prior of the grain warehouse, and to remove outliers. The density analysis and missing region identification module is used to calculate the local spatial density of point clouds and delineate reliable observation areas and missing regions; The dynamic angle of repose estimation module is used to estimate the slope field through local surface fitting within a reliable observation area and adaptively calculate the dynamic angle of repose. The physical constraint surface reconstruction module is used to construct an energy functional model that integrates data fidelity, smoothing, and repose angle constraints, and generates a complete three-dimensional surface height field through numerical optimization. The volume calculation module is used to discretize the height field into a triangular mesh and calculate the total volume of the grain pile.

10. The three-dimensional reconstruction system for grain pile surface based on millimeter-wave radar according to claim 9, characterized in that: The narrow-beam millimeter-wave radar in the three-dimensional scanning module of the millimeter-wave radar adopts a waveguide slot antenna and Fresnel lens cascade design to compress the beam width, and adopts a frequency-modulated continuous wave system to achieve range dimension super-resolution measurement through Chirp-Z transform.

Images