A Precise Monitoring Method for Warehouse Materials Based on the Fusion of LiDAR and Millimeter-Wave Radar

Abstract

This invention relates to the field of multi-sensor fusion technology, and more particularly to a method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar. The method includes the following steps: simultaneously scanning a material pile using lidar and millimeter-wave radar in both modes to obtain penetration echo data and surface measurement data; dividing the material pile into depth layers using the penetration echo data and calculating the interlayer dielectric parameters of the depth layers; detecting density gradient anomalies based on the interlayer dielectric parameters to obtain a layered density feature map; determining the reliability of the point cloud from the surface measurement data, and using the layered density feature map to fill in missing areas and reconstruct the surface contour to obtain a fused surface contour; calculating the support stress distribution on the layered density feature map and the fused surface contour to obtain the vertical pressure distribution field and support strength matrix. This invention significantly improves the accuracy of calculating the storage capacity of stored materials through deep layered analysis and coupling analysis of internal and external structures.

Description

Technical Field

[0001] This invention relates to the field of multi-sensor fusion technology, and in particular to a method for precise monitoring of warehouse materials based on the fusion of lidar and millimeter-wave radar. Background Technology

[0002] Existing material storage monitoring technologies, relying on single sensors, cannot adapt to the unique environment of material warehouses. High-concentration dust generated during material processing and storage severely interferes with lidar signals, causing data loss or distortion. While alternative technologies such as infrared and ultrasonic sensors also have measurement blind spots in dusty environments, resulting in low reliability of the monitoring system and an inability to continuously provide complete information on the status of the material pile. Traditional material level measurement technologies can only acquire surface contour information and cannot penetrate the material surface to detect the internal structure. When uneven material accumulation or long-term storage leads to internal voids and loose areas, there is a lack of effective identification and early warning mechanisms, and collapse accidents often occur suddenly without warning, causing casualties and equipment damage. Existing technologies typically use a single bulk density coefficient to convert volume to weight, ignoring the density differences between different types of materials (refined materials, coarse materials, etc.). Furthermore, they do not consider the vertical density gradient changes caused by compaction effects. In mixed storage conditions, measurement errors can reach as high as 15-20%, seriously affecting inventory management and cost control.

[0003] In summary, existing technologies suffer from problems such as poor sensor environmental adaptability, lack of internal structure detection capabilities, and insufficient measurement accuracy, which urgently need to be addressed. Summary of the Invention

[0004] Therefore, it is necessary to provide a method for precise monitoring of warehouse materials based on the fusion of lidar and millimeter-wave radar to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, a method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar includes the following steps:

[0006] Step S1: Perform dual-mode synchronous scanning of the material pile using laser-millimeter-wave radar to obtain a spatiotemporally aligned point cloud, which includes penetration echo data and surface measurement data; use the penetration echo data to divide the material pile into depth layers and calculate the interlayer dielectric parameters of the depth layers; perform density gradient anomaly detection based on the interlayer dielectric parameters to obtain a layered density feature map;

[0007] Step S2: Determine the reliability of the point cloud of the surface measurement data, use the layered density feature map to fill in the missing areas and reconstruct the surface contour to obtain the fused surface contour;

[0008] Step S3: Spatial registration of the layered density feature map with the fused surface profile, and calculation of the support stress distribution of each spatial unit to obtain the vertical pressure distribution field and support strength matrix; identification of voxel weak regions of the support strength matrix based on the vertical pressure distribution field, and weak region connectivity analysis to obtain weak region connectivity data; quantitative assessment of collapse risk based on weak region connectivity data to obtain stress risk distribution map.

[0009] Step S4: Identify the material type marking data of the material pile based on the stress risk distribution map; determine the bulk density coefficient mapping table based on the material type marking data; dynamically correct the zoned bulk density of the material pile using the bulk density coefficient mapping table to obtain the corrected storage quantity data.

[0010] This invention achieves non-destructive penetration detection of the internal structure of a material silo in a dusty environment through dual-modal synchronous scanning, overcoming the limitations of a single sensor. LiDAR provides millimeter-level surface accuracy, while millimeter-wave radar provides centimeter-level penetration capability; their complementary advantages form a comprehensive perception. Spatiotemporal aligned point clouds ensure data consistency across time and space, laying the foundation for subsequent analysis. Depth layer segmentation and interlayer dielectric parameter calculation enable quantitative characterization of the material's internal physical properties, while density gradient anomaly detection accurately identifies internal voids and loose areas—structural hazards difficult to detect using traditional methods. This step, through electromagnetic wave physical property analysis, visualizes the internal state of the sealed warehouse, providing early risk identification capabilities and preventing safety accidents caused by changes in internal structure. Fine surface contour reconstruction effectively solves the problem of missing laser data in high-concentration dust environments, improving the reliability of the monitoring system in harsh environments. A point cloud quality assessment mechanism automatically identifies low-reliability areas, and millimeter-wave data supplementation technology effectively fills measurement blind spots. Dual-source data fusion and stitching employs a distance-weighted averaging method, ensuring data continuity in transition areas and avoiding abrupt changes in height at the stitching points. The surface smoothing reconstruction algorithm generates a continuous surface model with millimeter-level accuracy, adaptively adjusting the window size for areas with large curvature changes while preserving detailed features. The fused surface profile provides a complete and accurate description of the external geometry of the material pile, offering a reliable basis for volume calculation and internal / external structural coupling analysis, enabling continuous and stable operation of the monitoring system in harsh environments. Internal / external structural coupling analysis establishes the correlation between the internal density distribution of the material and its surface morphology, realizing the transformation from static structural features to dynamic mechanical analysis. Spatial registration ensures strict correspondence between internal and external data, and the calculation of support stress distribution simulates the actual pressure transmission process, considering stress dispersion effects. Weak zone connectivity analysis aggregates discrete weak points into physically meaningful structural units, enhancing the analysis's relevance through shape factor and overburden impact assessment. Quantitative collapse risk assessment comprehensively considers four key factors: weak zone volume, overburden pressure, surrounding support capacity, and safety factor, transforming qualitative judgments into quantitative indicators. Risk level classification and 3D visualization enable operators to intuitively identify dangerous locations, providing a scientific basis for safety management decisions. This step achieves accurate detection of internal voids in materials and quantitative early warning of collapse risks, fundamentally improving the level of warehouse safety. The zoned dynamic density correction overcomes the limitations of traditional single density coefficient calculation methods, significantly improving measurement accuracy under mixed storage conditions. Density zoning and material type identification enable accurate classification of different types of materials, and the density coefficient determination fully considers the influence of three key factors: stacking time, degree of compaction, and depth. The time compaction factor reflects the natural compaction effect of materials as stacking time increases, the depth compaction correction considers the change in pressure gradient in the vertical direction, and the standardized density coefficient ensures consistent comparison between different material types.Precise calculation of partition volume combined with bulk density coefficient enables accurate calculation of partition quality. The calibration mechanism ensures the reliability of the results by comparing with historical records, reducing measurement error from 15-20% in traditional methods to less than 3%, providing accurate data support for inventory management, cost control, and production planning, while improving the accuracy of material quality management.

[0011] Therefore, this invention provides a method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar, achieving data complementarity and automatic switching in dusty environments. It innovatively employs a layered analysis of penetration depth and a coupled analysis process of internal and external structures, enabling accurate identification of internal voids and quantitative assessment of collapse risk. Simultaneously, it introduces a dynamic bulk density correction mechanism based on spatial partitioning, applying differentiated bulk density coefficients to materials of different densities, significantly improving the accuracy of storage volume calculation. This multi-dimensional fusion method provides a novel technical path for precise monitoring of stored materials. Attached Figure Description

[0012] Figure 1 This is a flowchart illustrating the steps of a method for precise monitoring of warehouse materials based on the fusion of lidar and millimeter-wave radar.

[0013] The objectives, features, and advantages of this invention will be further explained in conjunction with the embodiments and with reference to the accompanying drawings. Detailed Implementation

[0014] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0015] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0016] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0017] To achieve the above objectives, please refer to Figure 1 This invention provides a method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar, comprising the following steps:

[0018] Step S1: Perform dual-mode synchronous scanning of the material pile using laser-millimeter-wave radar to obtain a spatiotemporally aligned point cloud, which includes penetration echo data and surface measurement data; use the penetration echo data to divide the material pile into depth layers and calculate the interlayer dielectric parameters of the depth layers; perform density gradient anomaly detection based on the interlayer dielectric parameters to obtain a layered density feature map;

[0019] In this embodiment of the invention, dual-modal synchronous scanning employs a hardware trigger to achieve millisecond-level synchronous acquisition of LiDAR and millimeter-wave radar. Coordinate transformation matrices and timestamp verification ensure strict spatiotemporal alignment of the dual-source data. The LiDAR provides 1.2 million high-precision surface data points, while the millimeter-wave radar provides a 140GHz penetration echo sequence, forming a spatiotemporally aligned point cloud containing surface geometry and internal structure information. Depth layer segmentation is based on echo time interval analysis, setting a 0.5 nanosecond time threshold to determine independent reflection layers using a formula. Calculate the actual thickness and divide the material into 5-8 layers of varying thickness. The interlayer dielectric parameter is determined by the reflection coefficient. Calculation, where and The first Layer and first The relative permittivity of the layers, considering the propagation distance correction factor. Density gradient anomaly detection calculates the interlayer variation rate. A normal variation baseline is established based on the data of the central area of ​​the material pile. Abnormal points that exceed the baseline range are marked. The void type is determined according to the dielectric constant value. Aggregated abnormal areas are formed through spatial clustering. A density level value (0-10) is assigned to construct a hierarchical density feature map.

[0020] Step S2: Determine the reliability of the point cloud of the surface measurement data, use the layered density feature map to fill in the missing areas and reconstruct the surface contour to obtain the fused surface contour;

[0021] In this embodiment of the invention, surface contour fine reconstruction extracts surface measurement data from the spatiotemporally aligned point cloud using lidar, and point cloud quality is assessed through reflection intensity analysis. A baseline reflection intensity value of 150 is set; points with an intensity below 60 are marked as low-reliability points. The effective point ratio for each grid cell is calculated. ,when Units with a density below 50% are marked as missing units. A region growing algorithm is used to aggregate adjacent missing units, defining regions with an area exceeding 100 cm² as missing data areas, and boundary point sets are extracted. For each missing region, millimeter-wave surface data is extracted from the hierarchical density feature map, processed into a 5cm × 5cm grid, and subjected to median filtering to form a supplementary height dataset. Within a 20cm wide transition zone, a distance-weighted average method is used to achieve smooth fusion of the two-source data, as shown in the following formula:

[0022] ;

[0023] in This represents the height value after merging. This indicates the altitude measured by the lidar. This indicates the altitude measured by millimeter-wave radar. and These are the weighting coefficients for laser data and millimeter-wave data, respectively. Finally, a local quadratic surface fitting method is used for smooth surface reconstruction, and the fitting model is... The window size is adaptively adjusted for edge regions, ultimately generating a fused surface profile with a resolution of 5mm.

[0024] Step S3: Spatial registration of the layered density feature map with the fused surface profile, and calculation of the support stress distribution of each spatial unit to obtain the vertical pressure distribution field and support strength matrix; identification of voxel weak regions of the support strength matrix based on the vertical pressure distribution field, and weak region connectivity analysis to obtain weak region connectivity data; quantitative assessment of collapse risk based on weak region connectivity data to obtain stress risk distribution map.

[0025] In this embodiment of the invention, the internal and external structural coupling analysis spatially registers the layered density feature map with the fused surface contour, requiring a registration accuracy of 5mm. Three-dimensional rigid body transformation is used to eliminate positional deviations. The material pile space is divided into a regular voxel mesh of 10cm × 10cm × 10cm, and material density values ​​(0-750kg / m³) are assigned according to density levels. Vertical pressure transmission calculation is based on voxel gravity accumulation, considering the pressure dispersion effect of a Gaussian distribution, where the calculation formula is as follows:

[0026] ;

[0027] in Indicates position Vertical pressure at the location, This indicates the mass of the voxel at that location. This indicates that the acceleration due to gravity is 9.8 m / s². This represents the pressure transmission weighting coefficient. The location of adjacent voxels in the upper layer is indicated, and the support strength assessment establishes a mapping relationship between density level and bearing capacity (0-90000Pa). Weak zone connectivity analysis screens weak voxels with a safety margin SM < 0, establishing a 26-neighborhood weak point adjacency graph. Isolated weak points with fewer than 8 voxels are filtered and merged. Characteristic parameters such as geometric center, volume, and shape factor are calculated to assess the impact of the overburden layer. Quantitative assessment of collapse risk extracts the surrounding support strength Ssurr, accumulates the overburden pressure Ptotal, dynamically determines the safety factor SF based on volume and shape, and applies the risk formula. Calculate the standard risk index, where Indicates the risk index. Represents the volume of the weak region (m³). This represents the total overburden pressure (Pa). Indicates the strength of the surrounding support (Pa). The safety factor is represented. Risk levels are divided into four levels: 0-30 (low), 30-60 (medium), 60-80 (high), and 80-100 (extremely high), generating a three-dimensional stress risk distribution map.

[0028] Step S4: Identify the material type marking data of the material pile based on the stress risk distribution map; determine the bulk density coefficient mapping table based on the material type marking data; dynamically correct the zoned bulk density of the material pile using the bulk density coefficient mapping table to obtain the corrected storage quantity data;

[0029] In this embodiment of the invention, the dynamic correction of bulk density by zone is based on the stress risk distribution map to divide density regions, setting clear standards for low-density regions (<300kg / m³), medium-density regions (300-500kg / m³), and high-density regions (>500kg / m³), and gradient analysis is used to determine the region boundaries. Material type identification is performed according to a density-material type comparison table. >520kg / m³ is considered refined feed. <430kg / m³ is considered coarse material, and 430-520kg / m³ is considered mixed material. The bulk density coefficient is determined by extracting the baseline bulk density from the material type marking data (600kg / m³ for fine material, 400kg / m³ for coarse material, and 500kg / m³ for mixed material), and then calculating the time-based compaction factor based on the storage time. ,in For maximum compaction gain (0.3 for fine materials, 0.5 for coarse materials), The compaction rate is 0.1 m / day. Indicates the stockpiling time (days). Deep compaction correction divides the stockpile into a surface layer (…). ), middle layer ( ) and bottom layer ( ), calculate the overall compaction coefficient Standardized bulk density coefficient Volume calculations were performed using the 5cm voxel method to accurately measure the volume of each section. The mass calculation formula is as follows. ,in Indicates the zoned correction mass (kg). This represents the volume of the partition (m³). This indicates that the reference density of concentrate is 600 kg / m³. This indicates the bulk density coefficient of the material. This represents the depth correction factor. Finally, the quality of each partition is summarized, and a calibration factor is introduced when the calculation deviation exceeds ±3%. Make corrections to generate corrected storage data that includes total storage, distribution ratio, and timestamps.

[0030] Preferably, the dual-modal synchronous scanning of the material pile in step S1 includes:

[0031] The dual radar response signals of the material pile are collected simultaneously by lidar and millimeter-wave radar.

[0032] The coordinate system is transformed and the timestamp is checked to obtain the checked response signal;

[0033] Radar point cloud data fusion is performed on the calibrated response signal to obtain a spatiotemporally aligned point cloud.

[0034] In this embodiment, dual-modal synchronous scanning first achieves millisecond-level synchronous acquisition of data from the lidar and millimeter-wave radar via a hardware trigger. This trigger employs a high-precision timing circuit to generate a rectangular pulse signal with a width of 10 microseconds. The rising edge of the pulse simultaneously triggers both radar devices. The lidar is configured in line-scan mode with a scanning frequency set to 10Hz. The single-scan angle range is 0° to 180°, with an angular resolution of 0.25° and a range resolution of 5mm. Each scan acquires high-density surface point cloud data of 1.2 million points. The millimeter-wave radar operates in the 140GHz band, using FMCW (Frequency Modulated Continuous Wave) mode with a bandwidth of 6GHz. Its scanning cycle is strictly synchronized with the lidar, achieving a penetration depth of 3 meters and a depth resolution of 10mm. Each transmission and reception consists of 16 echo signals forming an echo sequence. A unified timestamp is added to the trigger time of both radar response data, achieving microsecond-level accuracy.

[0035] The coordinate system transformation is implemented using a rigid body transformation matrix. Based on the installation position relationship of the two radars (horizontal spacing 50cm, vertical height difference 30cm, azimuth deviation 5°), a unified Cartesian coordinate system is established. The original millimeter-wave radar data is in polar coordinate form (r, θ, φ), which is converted to Cartesian coordinates (x, y, z) in the unified coordinate system using the coordinate transformation matrix T. The transformation matrix T includes a rotation matrix R and a translation vector t, satisfying the equation P' = R × P + t, where P is the original coordinate point and P' is the transformed coordinate point. The rotation matrix R is determined by three Euler angles (roll = 2°, pitch = 3°, yaw = 5°), and the translation vector t is (50, -10, 30)cm. After the transformation, all measurement data are unified into a global coordinate system with the warehouse center as the origin. Next, timestamp verification is performed to detect differences in timestamps between the two radar data packets. When the difference exceeds 1 millisecond, linear interpolation correction is applied to the sampling points. The correction formula is: t_corrected = t_original + Δt, where Δt is the time deviation. After time synchronization is completed, each spatial location has dual radar measurement data with strictly aligned timestamps.

[0036] Radar point cloud data fusion is achieved by constructing a unified data structure. This data structure includes: three-dimensional spatial coordinates (x, y, z), laser reflection intensity values ​​I_laser (range 0-255), millimeter-wave echo intensity I_mmw (range 0-100dB), millimeter-wave penetration depth array D_penetration (containing depth values ​​of multiple reflection points), and measurement timestamp t. For each spatial location point, the lidar provides precise surface coordinates and reflection intensity, while the millimeter-wave radar provides penetration depth information and internal multi-layer echo intensity for the corresponding location. When laser data is missing at a point (reflection intensity below the threshold of 40), the surface reflection point of the millimeter-wave data is automatically used as the coordinate value for that location. The global point cloud density is 1000 points / square meter, and the total number of points is maintained at around 1.2 million, forming a complete spatiotemporally aligned point cloud containing fine surface features and internal penetration information. This point cloud data achieves strict alignment of the two radars in the temporal and spatial dimensions, laying the foundation for subsequent deep layering analysis and internal and external structural coupling analysis.

[0037] Preferably, step S1, which involves dividing the material pile into depth layers and calculating the interlayer dielectric parameters of the depth layers, includes:

[0038] Penetration echo data from millimeter-wave radar is extracted from spatiotemporally aligned point clouds, and echo sequence separation is performed to form multi-layer echo pulse sequences;

[0039] When the time interval between two adjacent echoes in a multi-layer echo pulse sequence is greater than the preset echo time threshold, it is determined to be an independent reflection layer. Based on the reflection layer, the material pile is divided into multiple depth layers of varying thicknesses from the surface to the bottom, and the boundary coordinates of the depth layers are recorded.

[0040] The interlayer dielectric parameters of the material stack are calculated based on the depth layer boundary coordinates.

[0041] In this embodiment of the invention, millimeter-wave radar penetration echo data is extracted from spatiotemporally aligned point clouds. For each spatial scan point, multi-layer echo information stored in its D_penetration array is extracted. Time-domain analysis is used to process the echo data, setting a signal strength threshold of -60dBm. When the echo signal strength is higher than this threshold, it is recorded as a valid echo point. The time-domain echo signal is converted to the frequency domain using a Fast Fourier Transform (FFT) to identify the dominant frequency component and filter out noise interference. Envelope detection is performed on the processed echo signal to extract peak point positions; each peak corresponds to a reflection interface within the material. These peak points are organized into an echo pulse sequence in chronological order, where each element contains the echo arrival time t_i and the corresponding echo intensity I_i. For each scan point position, an independent echo pulse sequence is generated, forming a multi-layer echo distribution matrix in three-dimensional space.

[0042] Independent reflection layers are determined based on multi-layer echo pulse sequences. An echo time threshold of 0.5 nanoseconds is set. When the time interval between two adjacent echoes, Δt = t_i + 1 - t_i, is greater than this threshold, it is considered a reflection interface between layers of different depths. The propagation speed v of millimeter waves in the material is set to 0.6c (c is the speed of light in vacuum, 2.998 × 10^8 m / s). The actual thickness of each layer is calculated using the formula d = v × Δt / 2, where v is the propagation speed of electromagnetic waves in the medium, Δt is the echo time difference, and the division by 2 is due to the round-trip propagation of the signal. Based on the calculation results, the material pile is divided into multiple layers of varying thicknesses from the surface to the bottom. Generally, the material pile is divided into 5 to 8 depth layers, with the surface layer typically 15-20 cm thick, the middle layer 20-40 cm thick, and the bottom layer 30-50 cm thick. Record the upper and lower boundary coordinates of each depth layer, with the upper boundary coordinates as (x, y, z_upper) and the lower boundary coordinates as (x, y, z_lower), and construct the depth layer boundary coordinate set.

[0043] The interlayer dielectric parameters of the material stack are calculated based on the boundary coordinates of the depth layers. The relative permittivity of each layer is calculated using the ratio of the intensity of the incident wave to the reflected wave, based on the principle of electromagnetic wave interface reflection. For the [missing information], [missing information] Layer and First The interface of the layer, the reflectivity ,in and The first Layer and first The relative permittivity of the layer. In practical applications, it is used to measure the intensity of incident waves. and reflected wave intensity Through the reflection coefficient formula Inversely solve for the permittivity. Consider the attenuation of electromagnetic waves during propagation and introduce a propagation distance correction factor. ,in The attenuation coefficient is... The relative permittivity ranges from 1.5 to 2.5 for the coarse material region, from 3.0 to 4.5 for the fine material region, and from 5.0 to 7.0 for regions with high moisture content. By calculating the difference in permittivity between adjacent layers, a complete table of interlayer permittivity parameters is established, laying the foundation for subsequent density gradient anomaly detection.

[0044] Preferably, step S1, which involves detecting density gradient anomalies based on interlayer dielectric parameters, includes:

[0045] The dielectric change rate of the interlayer dielectric parameters is calculated to obtain the interlayer change rate matrix;

[0046] Data from the central region of the material pile is selected from the interlayer change rate matrix to establish a normal change benchmark and obtain the stratification benchmark range.

[0047] Identify outlier locations in the inter-layer rate of change matrix based on the stratified baseline range;

[0048] Based on the interlayer dielectric parameters, the void characteristics of the abnormal point location are determined, and a void type classification table is obtained;

[0049] Aggregate abnormal regions based on the cavity type classification table to obtain aggregated abnormal regions;

[0050] Density anomaly marker data is generated based on aggregated anomaly regions;

[0051] Based on the density anomaly marker data, density level values ​​are assigned to the depth layers to obtain the layered density feature map.

[0052] In this embodiment of the invention, density gradient anomaly detection first calculates the rate of change of the interlayer dielectric parameter table. For each spatial measurement point (x, y), the dielectric constant sequence {ε_1, ε_2, ..., ε_n} of each depth layer is extracted, and the rate of change of dielectric constant between adjacent layers is calculated as CR_i = (|ε_i - ε_i+1| / ε_i+1) × 100%, where CR_i represents the rate of change of dielectric constant between the i-th layer and the (i+1)-th layer, and ε_i and ε_i+1 represent the relative dielectric constants of the i-th layer and the (i+1)-th layer, respectively. An interlayer rate of change matrix with dimension m × n is constructed, where m is the number of scan points and n is the number of depth layers minus 1. Then, data from the central region of the material pile is extracted from the interlayer rate of change matrix to establish a normal change benchmark. The central region is defined as all measurement points within a radius less than 40% of the maximum radius of the material pile from the horizontal center point of the material pile. The mean μ_i and standard deviation σ_i of the rate of change at each depth layer within the region are statistically analyzed, and the normal range of change is defined as [μ_i-2σ_i, μ_i+2σ_i]. Reference ranges are established for different depth layers: typically 15%±10% from the surface to the second layer, 10%±8% between intermediate layers, and 5%±5% in deep regions. Based on the established layered reference ranges, the values ​​in the interlayer rate of change matrix are compared point-by-point. When the rate of change CR_i at a measurement point in layer i exceeds the upper limit of the reference range for the corresponding depth, it is marked as an anomaly. Specifically, when the rate of change exceeds 30%, it is directly identified as a strong anomaly. The three-dimensional coordinates (x, y, z_i) of each anomaly point and its corresponding depth layer i are recorded, constructing a set of anomaly point locations.

[0053] For the marked anomaly locations, the interlayer dielectric parameter values ​​are checked to determine the void characteristics. Void types are categorized into three types: when the dielectric constant ε_i < 1.5, it is classified as a complete void, labeled as type 0; when 1.5 ≤ ε_i < 2.5, it is classified as a loose region, labeled as type -1 to -3, with specific values ​​linearly mapped from (2.5 - ε_i); when ε_i ≥ 2.5 and significantly different from adjacent layers, it is classified as a density abrupt change region, labeled as type 1. The spatial coordinates, depth layer, and corresponding void type of each anomaly are recorded to form a void type classification table. Anomaly regions are then aggregated using a spatial proximity clustering algorithm. An aggregation threshold of 20 cm is set; when the horizontal distance between two anomalies is less than this threshold and the depth layer difference is no more than one layer, they are classified into the same aggregation region. For each aggregation region, its center coordinates are calculated as the average of the coordinates of the contained anomalies, and its size is the minimum circumscribed ellipsoidal volume containing all points. The primary anomaly type is the void type with the highest frequency of occurrence. The aggregation results are stored in the aggregation anomaly region table, which includes fields such as region ID, center coordinates, volume size, and main anomaly type.

[0054] Density anomaly marker data is generated based on aggregated anomaly regions. The monitoring space is divided into a regular grid of 10cm × 10cm × 10cm. For each grid cell, it is checked whether it is located within a certain aggregated anomaly region. If it is located within the anomaly region, a corresponding anomaly type marker (0, -1 to -3, or 1) is assigned; if it is not located within the anomaly region, its density level is determined based on the dielectric constant value ε_i at that location. The mapping relationship between dielectric constant and density is set as follows: ε_i = 1.5~2.0 corresponds to density levels 1~2 (very low density), ε_i = 2.0~3.0 corresponds to density levels 3~5 (low density), ε_i = 3.0~4.0 corresponds to density levels 6~8 (medium density), and ε_i = 4.0~7.0 corresponds to density levels 9~10 (high density). The specific density level value is determined by linear interpolation, for example, ε_i = 2.5 corresponds to density level 4. For each depth layer, the corresponding density level value is assigned based on the density anomaly marker data, and a complete layered density feature map is constructed. This feature map is represented in the form of a three-dimensional color image, with different colors representing different density levels: red indicates void areas (level 0), yellow indicates loose areas (levels 1-3), green indicates normal areas (levels 4-7), and blue indicates compacted areas (levels 8-10). This density gradient anomaly detection method achieves non-destructive detection of the internal structure of materials through dielectric parameter analysis, especially the accurate identification of potential voids.

[0055] Preferably, step S2 includes:

[0056] Surface measurement data from lidar are extracted from spatiotemporally aligned point clouds to perform lidar point cloud quality assessment and obtain point cloud reliability distribution.

[0057] The missing regions of the point cloud reliability distribution are identified to obtain the set of missing region boundaries;

[0058] Based on the missing region boundary set, millimeter-wave radar surface data is extracted from the hierarchical density feature map to fill in the missing regions, resulting in a supplementary height dataset.

[0059] The point cloud reliability distribution and the supplementary height dataset are fused and stitched together to obtain a complete height point matrix.

[0060] A smooth surface reconstruction is performed on the complete height lattice to obtain the fused surface profile.

[0061] In this embodiment of the invention, the fine reconstruction of the surface contour first extracts surface measurement data from the spatiotemporally aligned point cloud using lidar, and then performs point cloud quality assessment. The three-dimensional coordinates (x, y, z) and reflection intensity value I_laser of each point are extracted, and data quality is judged by analyzing the reflection intensity. A baseline reflection intensity value I_baseline is set to 150 (range 0-255). When the reflection intensity I_laser of a point is lower than 40% of the baseline value (i.e., 60), it is marked as a low-reliability point, indicating that the point is affected by dust interference. The proportion of valid points within each 10cm×10cm grid is calculated: R_valid = N_valid / N_total, where N_valid is the number of points with reflection intensity higher than the threshold, and N_total is the total number of points in the grid. The R_valid value is mapped to the range of 0-100% to construct a point cloud reliability distribution matrix. Next, connected component analysis is performed on the point cloud reliability distribution to identify missing regions. When the proportion of valid points R_valid of a certain grid cell is lower than 50%, it is marked as a potentially missing cell. A region growing algorithm was used to aggregate adjacent missing units. When the area of ​​a continuous region exceeded 100 cm², it was identified as a missing data region. The α-shape algorithm was used to extract the boundary point set of each missing region, and the three-dimensional coordinates (x_b, y_b, z_b) of the boundary points were recorded to form the boundary set of the missing region. Each missing region was recorded with its ID number, center position, area size, and list of boundary points.

[0062] Millimeter-wave radar surface data is extracted from the hierarchical density feature map based on the boundary set of missing regions. For each missing region, its horizontal range (x_min, x_max, y_min, y_max) is determined, and the first layer (surface layer) data of the hierarchical density feature map is extracted within this range. The millimeter-wave radar surface data contains spatial location (x_m, y_m) and the corresponding surface height value z_m. Since the horizontal resolution of millimeter-wave radar (approximately 5cm) is lower than that of lidar (approximately 0.5cm), gridding is required. The missing regions are divided into regular grids of 5cm × 5cm, and the average height value of the millimeter-wave data is calculated for each grid. To eliminate noise in the millimeter-wave data, a medium-range filter with a window size of 3×3 is applied to filter out outliers. The processed height data is used as a supplementary height dataset to fill in the missing regions of the lidar data.

[0063] Dual-source data fusion and stitching achieves seamless connection between laser point cloud and millimeter-wave supplementary data. To ensure a smooth data transition, a 20cm wide transition zone is set at the boundary of the missing region. Within the transition zone, the fusion height value is calculated using a distance-weighted average method: z_fused=(w_l×z_laser+w_m×z_mmw) / (w_l+w_m), where w_l and w_m are the weighting coefficients for the laser data and millimeter-wave data, respectively. The weighting coefficients are calculated based on the distance d from the point to the boundary of the missing region: w_l=(20-d) / 20 (when d<20cm), w_m=d / 20 (when d<20cm); when d≥20cm, w_l=0, w_m=1 inside the missing region, and w_l=1, w_m=0 outside the missing region. This method achieves a smooth transition between the two data sources, avoiding abrupt height changes at the stitching point. The fusion result is reorganized into a regular height matrix with a resolution of 2cm×2cm, covering the entire surface area of ​​the material pile.

[0064] A smooth surface reconstruction is performed on the complete height lattice to generate a continuous and smooth surface profile. A local quadratic surface fitting method is used to fit each point p(x,y,z) and points within a 5×5 grid around it. During the fitting process, to reduce the influence of noise, different weights are assigned to each point: the center point has a weight of 1.0, and the weight decreases linearly with increasing distance. For areas with large curvature changes (such as the edge of the material pile), the fitting window size is adaptively adjusted, and a 3×3 window is used in the edge region to preserve detailed features. After fitting, spline interpolation is performed on the global surface to generate a uniform grid height field with a resolution of 5 mm. The surface reconstruction accuracy reaches the millimeter level, meeting the requirements for accurate calculation of the material pile volume. The fused surface profile includes three-dimensional coordinates and surface normal vector information, fully describing the external geometry of the material pile and providing basic data for subsequent coupled analysis with the internal structure. The innovation here lies in using millimeter-wave radar data to supplement the measurement blind zone of lidar in dusty environments, realizing the complementary advantages of dual-modal sensors.

[0065] Preferably, step S3 involves spatially registering the layered density feature map with the fused surface profile and calculating the support stress distribution for each spatial unit, including:

[0066] Spatial data registration is performed between the layered density feature map and the fused surface contour to obtain the registered spatial data;

[0067] Based on the fused surface profile and hierarchical density feature map, the material storage space is divided into a regular voxel grid, forming a density-marked voxel grid.

[0068] Vertical pressure transmission is calculated on density-marked voxel meshes to obtain the vertical pressure distribution field;

[0069] The support strength matrix is ​​obtained by evaluating the support strength of the density-marked voxel mesh based on the vertical pressure distribution field.

[0070] In this embodiment of the invention, the internal and external structure coupling analysis first involves spatial data registration between the layered density feature map and the fused surface profile. Using a fixed reference point in the warehouse as a baseline, a three-dimensional rigid body transformation method is employed to eliminate positional discrepancies between the two sets of data. The transformation includes a translation vector t and a rotation matrix R, ensuring that the two sets of data are strictly aligned to a unified coordinate system. The alignment accuracy must be within 5mm; exceeding this error will lead to deviations in subsequent analysis. After implementing the rigid body transformation, the consistency between the surface profile data and the layered depth values ​​of the density feature map is verified, and the root mean square error is calculated. ,when When the thickness is greater than 5mm, fine correction is performed on the local area using bilinear interpolation. After registration, each point of the surface profile ( ) and the corresponding surface points in the density feature map ( Strictly correspondence, and The difference should not exceed 5 mm. The registration result is stored as registration spatial data in a unified three-dimensional coordinate system, including spatial location, surface height, and internal density distribution information.

[0071] The material storage space is divided into a regular voxel grid based on the registration space data, forming a density-marked voxel grid. The voxel size is set to 10cm × 10cm × 10cm, which is sufficiently fine to capture local structural changes without causing excessive computation. Each voxel is assigned a unique three-dimensional index (…). ), corresponding to actual spatial coordinates ( ),in( The origin is set to ( ). For each voxel, it is determined whether it is located inside the material pile: if the height of the voxel's center point is lower than the surface profile height at that location, it is marked as an internal voxel; otherwise, it is marked as an external voxel and does not participate in subsequent calculations. For all internal voxels, the density level value (0-10) at the corresponding location is extracted from the layered density feature map and recorded in the density-marked voxel grid data structure. The density-marked voxel grid is stored using a three-dimensional sparse matrix, and each voxel records its spatial index, center coordinates, depth layer, and corresponding density level.

[0072] The vertical pressure transmission calculation is based on the principle of voxel gravity accumulation. First, a corresponding material density value is assigned to each density level: density level 0 (void) corresponds to 0 kg / m³; density levels 1-3 (loose zone) correspond to 200-350 kg / m³; density levels 4-7 (normal zone) correspond to 350-550 kg / m³; and density levels 8-10 (compacted zone) correspond to 550-750 kg / m³. Specific density values ​​are determined through linear interpolation, such as density level 5 corresponding to 425 kg / m³. The mass of each voxel is calculated: m = ρ × V, where ρ is the material density and V is the voxel volume (0.001 m³). For surface voxels, the vertical pressure they experience is only their own weight: P = m × g, g = 9.8 m / s². For internal voxels, the gravity contribution of all voxels above them is accumulated, considering the stress dispersion effect. Vertical pressure transmission equation: P(i,j,k)=m(i,j,k)×g+∑w(i',j',k'-1)×P(i',j',k'-1), where (i',j',k'-1) are the upper-layer adjacent voxels, and w is the pressure transmission weighting coefficient, which decreases with increasing horizontal distance and follows a Gaussian distribution. Weighting calculation formula: ,in Horizontal distance The standard deviation is given. The calculation results form a three-dimensional vertical pressure distribution field P(i,j,k), with units of Pa.

[0073] The support strength assessment calculates the maximum bearing capacity of voxels based on their density grades. A mapping relationship is established between density grades and support strength: density grade 0 (void) has a support strength of 0 Pa; density grades 1-3 (loose zones) have a support strength of 10000-30000 Pa; density grades 4-7 (normal zones) have a support strength of 30000-60000 Pa; and density grades 8-10 (compacted zones) have a support strength of 60000-90000 Pa. Specific support strength values ​​S(i,j,k) are determined through piecewise linear mapping. The bearing safety factor for each voxel is calculated: SF(i,j,k) = S(i,j,k) / P(i,j,k). When SF < 1, it indicates that the pressure borne by the voxel exceeds its support capacity, posing a structural risk. The calculation results are stored as a three-dimensional support strength matrix, containing the location index, support strength value, vertical pressure value, and safety factor for each voxel. The core innovation of this step lies in transforming the density distribution information inside the material into a mechanical model, enabling coupled analysis of internal voids and external pressure, and providing a mechanical basis for collapse risk assessment. By accurately simulating the vertical pressure transmission process, potential areas with insufficient support were effectively identified, laying the foundation for subsequent connectivity analysis of weak areas.

[0074] Preferably, step S3, which involves identifying voxel weak regions of the support strength matrix based on the vertical pressure distribution field and performing connectivity analysis of these weak regions, includes:

[0075] By combining the support strength matrix with the vertical pressure distribution field to determine the support strength and pressure of voxels, the location of weak voxels can be identified.

[0076] Establish a weakness adjacency graph for the locations of weakness voxels;

[0077] The weak adjacency graph is filtered and merged into a set of effective connected regions.

[0078] The geometric features of the effective connected regions are calculated to obtain the region geometric parameters;

[0079] The overlying influence was assessed based on regional geometric parameters and vertical pressure distribution field, yielding weak-area connectivity data.

[0080] In this embodiment of the invention, the weak region connectivity analysis first compares the support strength matrix and the vertical pressure distribution field voxel by voxel to screen weak voxels. The safety margin SM(i,j,k) = S(i,j,k) - P(i,j,k) is calculated for each voxel, where S(i,j,k) is the support strength and P(i,j,k) is the vertical pressure. When SM < 0, the voxel cannot withstand the pressure from above and is marked as a weak voxel. The three-dimensional index (i,j,k), spatial coordinates (x,y,z), negative safety margin value |SM|, and depth layer k of each weak voxel are recorded to construct a weak voxel location table. Next, a weak adjacency graph is established. For each voxel in the weak voxel location table, its 26 adjacent positions (up, down, left, right, front, back, and diagonal directions) are checked to determine adjacent weak points. An adjacency matrix is ​​used to store the connection relationships. A matrix element A(m,n) = 1 indicates that the m-th weakness is adjacent to the n-th weakness; otherwise, A(m,n) = 0. The adjacency criterion is the Euclidean distance between the centers of the two voxels. (Maximum distance between diagonally adjacent nodes). This method forms an adjacency graph G(V,E) describing the spatial distribution of weaknesses, where V is the set of weaknesses and E is the set of adjacency relationships.

[0081] Perform connected component analysis on the weak point adjacency relationship graph using the breadth - first search (BFS) algorithm. Starting from any unprocessed weak point voxel, gradually expand through the adjacency relationship, and classify all connected weak points into the same connected region. When no further expansion is possible, this connected region is determined to be complete and assigned a unique ID. Repeat this process for all weak points until all are classified. For each connected region, count the number of voxels N it contains. Set the voxel number threshold Nmin = 8, and mark the regions with N < Nmin as isolated weak points, which are not regarded as structural threats. For medium - sized regions with 8 ≤ N < 15, calculate the minimum distance dmin between them and other regions. When dmin < 20 cm (the distance between two voxels), merge them into the nearest large region. When merging, update the voxel list of the connected region and the adjacency relationship graph. Repeat this process until no further merging is possible to form an effective set of connected regions.

[0082] Calculate the geometric features of the effective set of connected regions. Calculate the geometric center coordinates of each region: , and the volume V = N×0.001 m³ (the volume of a single voxel is 0.001 m³). Determine the spatial range of the region: Calculate the maximum and minimum coordinate values in each direction to obtain the bounding box size L. Calculate the shape factor SF = Lmax / Lmin, where Lmax and Lmin are the maximum and minimum values in L respectively. The shape factor reflects the geometric characteristics of the region. SF close to 1 indicates an equiaxed shape, and SF > 3 indicates an elongated shape. Store the geometric parameters corresponding to each region ID to form a regional geometric parameter table. Finally, conduct an overlying influence assessment. For each connected region, determine its top position ztop = max{zi|(i,j,k) ∈ region}. Calculate the overlying material layer thickness h = zsurface - ztop, where zsurface is the height value of the surface profile at this position. Statistically calculate the total pressure of the overlying material layer Ptotal = ∑P(i,j,k), where (i,j,k) are all voxels from the top of the region to the surface. Evaluate the potential influence range: Taking the center of the region as the base point, extend upward at a 45° diffusion angle to the surface and calculate the influence area . The overlying thickness , the total pressure Ptotal, and the influence area [[ID=I0]] are integrated with the regional geometric parameters to form complete weak area connection data.

[0083] Preferably, the quantitative assessment of the collapse risk according to the weak area connection data in step S3 includes:

[0084] Extract the peripheral support force of the weak area connection data using the support strength matrix to obtain the weak area support parameters;

[0085] Perform cumulative calculation of the overlying pressure according to the weak area connection data and the vertical pressure distribution field to obtain the weak area pressure load;

[0086] The safety factor is dynamically determined based on the weak zone support parameters and the weak zone pressure load, resulting in a graded safety factor.

[0087] The standard risk index is calculated based on the graded safety factor.

[0088] Risk levels are classified according to standard risk indices to obtain a graded risk assessment table;

[0089] Spatial risk mapping is performed on the graded risk assessment table to obtain a stress risk distribution map.

[0090] In this embodiment of the invention, the quantitative assessment of collapse risk first involves extracting the surrounding support force from the connectivity data of weak zones. For each connected weak zone, the location of its outer boundary voxel is determined. A boundary voxel is defined as a voxel that belongs to the interior of the weak zone and has at least one adjacent voxel that does not belong to the weak zone. A 3D edge detection algorithm is used to identify the boundary voxels, constructing a boundary voxel set B={b1,b2,...,bn}. For each boundary voxel bi, all non-weak zone voxels within a voxel distance (10cm) around it are searched in the support strength matrix; these voxels constitute the support boundary of the weak zone. The support boundary voxel set is SB={sb1,sb2,...,sbm}. The support strength values ​​{S(sb1),S(sb2),...,S(sbm)} of these support boundary voxels are extracted, and the average value is calculated as the surrounding support strength Ssurr=∑S(sbi) / m of the weak zone. The ID, volume V, center coordinate C, and corresponding surrounding support strength Ssurr of each weak zone are recorded in the weak zone support parameter table. For weak areas with uneven surrounding support, the standard deviation of the support strength is calculated additionally. ,when When / Ssurr>0.3, the surrounding support strength Ssurr is reduced by 20% to reflect the additional risk of uneven support.

[0091] The cumulative overburden pressure calculation is based on the top location of the weak zone and the vertical pressure distribution field. For each weak zone, a set of voxels T={t1,t2,...,tk} is determined at its top location, containing all weak zone voxels with the largest z-coordinate. For each top voxel ti, the pressure values ​​{P(ti,1),P(ti,2),...,P(ti,h)} of all voxels directly above it up to the surface are extracted from the vertical pressure distribution field, where h is the number of overburden voxels. Considering the diffusion effect of pressure transmission, the voxel pressures within a 45° angle diagonally upward are weighted and accumulated. Weighting coefficients are used. According to horizontal distance Sure: ,in Calculate the cumulative pressure Pacc(ti) = ∑w(d) × P(ti,j) borne by each top voxel ti. The total overburden pressure of the weak zone is calculated as Ptotal = ∑Pacc(ti). Record the ID of each weak zone and the corresponding total overburden pressure Ptotal in the weak zone pressure load table as key parameters for assessing collapse risk.

[0092] The safety factor is dynamically determined based on the volume and shape characteristics of the weak zone. First, a base safety factor SF_base is determined: SF_base = 1.5 when volume V < 0.1 m³; SF_base = 2.0 when 0.1 m³ ≤ V < 0.5 m³; and SF_base = 2.5 when V ≥ 0.5 m³. Next, adjustments are made based on the shape factor SF: the safety factor increases by 0.5 when SF > 3 (slender shape); and increases by 1.0 when SF > 5 (extremely slender shape). Finally, the depth factor is considered: the safety factor increases by 0.3 for surface weak zones (<1 m from the surface); and decreases by 0.2 for deep weak zones (>3 m from the surface). The final safety factor SF = SF_base + SF adjustment value. The determined safety factors are stored in a graded safety factor table, corresponding to the weak zone IDs, for subsequent risk index calculations.

[0093] The standard risk index is calculated using the comprehensive risk assessment formula: RI = (V × Ptotal) / (Ssurr × SF), where V is the volume of the weak zone, Ptotal is the total overburden pressure, Ssurr is the strength of the surrounding support, and SF is the safety factor. This formula considers four key factors: the size of the weak zone, the overburden load, the surrounding support capacity, and the safety margin, achieving a comprehensive quantification of subsidence risk. For ease of comparison and classification, the calculation results are normalized, mapping the risk index to a standard range of 0-100: RI_std = 100 × RI / RI_max, where RI_max is the maximum risk index in historical data (set to 10). When the calculated value exceeds 100, it is truncated to 100. The standardized risk index values ​​are recorded in the standard risk index table, corresponding one-to-one with the weak zone ID.

[0094] Risk level classification is based on a standard risk index with clearly defined grading criteria. A risk index of 0-30 indicates low risk, marked as Level 1 (green); 30-60 indicates medium risk, marked as Level 2 (yellow); 60-80 indicates high risk, marked as Level 3 (orange); and 80-100 indicates extremely high risk, marked as Level 4 (red). For weak areas marked as Level 3 and 4, the system automatically generates warning signals. Each weak area records its ID, spatial location, volume, risk index value, and corresponding risk level, forming a graded risk assessment table. This table serves as the core basis for safety management decisions; the higher the risk level, the more stringent the preventative measures required. For Level 4 risk areas, immediate support reinforcement or unloading measures are required.

[0095] Spatial risk mapping maps the risk level information from the graded risk assessment table back to the original 3D spatial grid. For all voxels covered by each weak zone, a corresponding risk level label is assigned. A 20cm wide transition zone is set at the boundary of the weak zone, with the risk level decreasing with distance: d < 5cm, the original level is maintained; 5cm ≤ d < 10cm, the level is reduced by 0.25; 10cm ≤ d < 15cm, the level is reduced by 0.5; and 15cm ≤ d < 20cm, the level is reduced by 0.75. Normal areas outside the weak zone are labeled as level 0 (safe). Complete 3D risk distribution data is generated, including the coordinates and risk level information for each spatial location. The stress risk distribution map is visualized using volumetric rendering technology, with different colors used to mark different risk levels: blue for level 0, green for level 1, yellow for level 2, orange for level 3, and red for level 4. High-risk areas are highlighted by adjusting transparency, allowing operators to visually identify hazardous locations. The innovation of this step lies in transforming the complex multi-parameter risk assessment results into an intuitive three-dimensional risk distribution map, providing a quantitative early warning tool for precise monitoring of material storage, and realizing the accurate assessment and visual representation of potential collapse risks.

[0096] Preferably, step S4 includes:

[0097] Based on the stress risk distribution map, the space of the material pile is divided into regions with different density levels, resulting in a density zoning map;

[0098] Material type identification and classification are performed based on density zoning maps to obtain material type labeling data;

[0099] Calculate the partition volume data based on the density partition map and the fused surface profile;

[0100] Mass distribution correction calculations are performed on the partitioned volume data and bulk density coefficient mapping table to obtain the mass distribution table;

[0101] The storage capacity of the mass distribution table is summarized and calibrated to obtain the corrected storage capacity data.

[0102] In this embodiment of the invention, the dynamic correction of bulk density in zoned areas first involves dividing density regions based on a stress risk distribution map. Density level data from the stress risk distribution map is extracted, and the stockpile space is divided into different density regions. A threshold segmentation method is used to set clear density zoning standards: density levels 0-3 are low-density regions (corresponding to densities <300 kg / m³); density levels 4-7 are medium-density regions (corresponding to densities 300-500 kg / m³); and density levels 8-10 are high-density regions (corresponding to densities >500 kg / m³). At the region boundaries, a gradient analysis algorithm is applied to determine the precise boundaries, calculating the density gradient vector g = ▽ρ of adjacent voxels. When ‖g‖ > 50 kg / m³ / dm, it is identified as a boundary point. To reduce the impact of noise, morphological processing is performed on the preliminary segmentation results, executing opening and closing operations with a radius of 2 voxels to eliminate isolated points and small regions. Finally, a spatial structure map containing region IDs, boundary coordinates, and density classifications is generated, i.e., a density zoning map.

[0103] Material identification and classification are based on density zoning maps and material physical properties. A density-material comparison table is established using material entry records: concentrates (mainly crushed grains or concentrated feed) have a density range of 550-650 kg / m³; coarse materials (mainly straw, hay, etc.) have a density range of 350-450 kg / m³; and mixed materials have a density range of 450-550 kg / m³. For each density zoning, its average density value is calculated. The material type is determined by the judgment rules: when When the concentration is >520 kg / m³, it is classified as a refined feed area; when When the density is less than 430 kg / m³, it is classified as a coarse material zone; when the density is less than or equal to 430 kg / m³, it is classified as a coarse material zone. When the concentration is ≤520 kg / m³, it is classified as a mixed aggregate area. For boundary areas that are difficult to determine directly (within ±15 kg / m³), a spatial continuity rule is adopted: areas adjacent to known fine aggregate areas are classified as fine aggregate, and areas adjacent to known coarse aggregate areas are classified as coarse aggregate. After classification, each area is assigned a unique aggregate type ID: fine aggregate is marked as 1, coarse aggregate as 2, and mixed aggregate as 3, forming aggregate type labeling data.

[0104] The volume calculation for each density zone employs a grid voxel method to accurately measure its volume. The density zone map is overlaid onto the fused surface profile to determine the surface boundary of each zone. Independent volume calculations are performed for each zone, dividing the region into small voxel units with sides of 5 cm. Each voxel is determined to be within its region: a voxel is marked as valid if its center point is below the surface profile and within the zone boundary. The number of valid voxels within each zone is then accumulated. Multiply by the volume of a unit voxel (i.e., 5cm×5cm×5cm), to obtain the precise volume value. For some voxels at the boundary, a volume calculation method involving the intersection of hexahedrons and surface triangular meshes is used to improve the accuracy of volume calculation. The ID, spatial range, material type, and volume value of each density zone are recorded to form zone volume data.

[0105] The corrected mass distribution calculation multiplies the volume of each zone by its corresponding bulk density coefficient. For each zone, a differentiated bulk density coefficient is applied based on its material type and compaction degree. A bulk density coefficient mapping table provides standardized bulk density coefficients for different materials under different compaction conditions: 1.0-1.3 for the fine material zone; 0.65-0.85 for the coarse material zone; and 0.85-0.95 for the mixed material zone. For the same material at different depths, a stratified bulk density coefficient is used: the surface layer (0-1m) uses the baseline value; the middle layer (1-3m) increases by 5-10%; and the bottom layer (>3m) increases by 10-15%. The corrected mass calculation formula for each zone is as follows: ,in For partition volume, The reference density for concentrate is 600 kg / m³. The density coefficient of the material type. This is the depth correction factor. Generate a mass distribution table containing partition ID, volume, material type, bulk density coefficient, and corrected mass.

[0106] Storage aggregation and calibration: The correction quality of all partitions is summed to calculate the total storage. Let the total number of partitions be n, then the total storage is... The calculation results were verified by comparing historical inbound and outbound data. When the calculation deviation exceeded the allowable range of ±3%, a calibration factor was introduced. Make corrections: ,in This is the recorded value. The calibrated storage capacity is... Simultaneously, the distribution ratio of various material types is statistically analyzed, ultimately generating corrected storage data, including total storage volume, weight, distribution ratio, and volume data of each material type, along with a timestamp, providing accurate inventory status at the current moment. The innovation of this method lies in combining internal structural information detected by multi-sensor fusion, enabling the identification and differentiated measurement of areas with different densities. This overcomes the limitations of traditional single-density coefficient calculation methods, significantly improving measurement accuracy under mixed storage conditions, reducing measurement error from the original 15-20% to less than 3%.

[0107] Preferably, step S4, determining the bulk density coefficient mapping table based on the material type marking data, includes:

[0108] Extract the reference bulk density of the material type from the material type labeling data;

[0109] Obtain the material's entry time into the warehouse, and calculate the stacking time factor by combining the material type's benchmark bulk density;

[0110] The comprehensive compaction coefficient is obtained by correcting the depth compaction based on the material type marking data and the time compaction factor.

[0111] The bulk density coefficient is standardized by combining the comprehensive compaction coefficient and the reference bulk density of the material type to obtain the standardized bulk density coefficient.

[0112] In this embodiment of the invention, the bulk density coefficient mapping table first extracts the reference bulk density of each material type from the material type marking data. A standard material type bulk density reference database is established, recording the reference bulk density values ​​of various materials under standard conditions (newly received, uncompacted). The standard bulk density of concentrated materials (including grain flour, concentrate, etc.) is set at 600 kg / m³; the standard bulk density of coarse materials (including hay, straw, etc.) is set at 400 kg / m³; and the standard bulk density of mixed materials is set at 500 kg / m³. For each region in the material type marking data, the corresponding reference bulk density value is matched according to its material type ID (concentrated materials = 1, coarse materials = 2, mixed materials = 3). The system simultaneously records the physical property parameters of different material types, including the bulk density coefficient: the bulk density coefficient of fine materials is set to 1.0 (baseline value), the bulk density coefficient of coarse materials is set to 1.8 (loose structure, easy to compress), and the bulk density coefficient of mixtures is set to 1.3 (medium compressibility). These baseline values ​​are used for subsequent compaction effect calculations.

[0113] The time-based compaction factor is calculated based on the material's arrival time and compaction patterns. Arrival time records for materials in each area are extracted from the warehouse management system, and the time difference between the current monitoring time and the arrival time is calculated. (In days). Long-term storage of materials leads to compaction due to their own weight. A formula for calculating the time-based compaction factor is established: ,in For maximum compaction gain (refined material) Coarse materials Mixture ), Compaction rate coefficient ( / day). This formula describes the exponential approach of the compaction effect to its limit over time. When the storage time exceeds 30 days, the compaction factor tends to stabilize, with a maximum compaction factor of 1.3 for concentrates, 1.5 for coarse materials, and 1.4 for mixtures. Specifically, when the storage time record is missing, a medium compaction state is used by default, i.e. The calculated time compaction factor Stored in relation to the region ID.

[0114] The depth compaction correction takes into account the pressure gradient effect caused by vertical stacking of materials. For each region in the material type marking data, a stratified correction is performed based on its vertical position. The material pile is vertically divided into three layers: surface layer (0-1m from the surface), middle layer (1-3m), and bottom layer (>3m). Different depth correction coefficients are assigned to different depth layers. :surface layer =1.0 (baseline value); middle layer ,in Depth value (meters); bottom layer The maximum value is no more than 1.2. The correction factor increases linearly with depth, reflecting the compaction effect of the upper material's gravity on the lower layer. Differentiated processing is applied to the compression characteristics of different material types: the calculated value is directly used in the concentrate area. Value; Due to its loose structure, the coarse material area exhibits significant compaction, resulting in an enhanced depth correction coefficient. The correction factor for the mixture area is: Time compaction factor With depth correction factor Multiply to obtain the overall compaction coefficient. .

[0115] Standardization of bulk density coefficients converts the comprehensive compaction coefficients of various regions into standardized bulk density coefficients under a unified benchmark. Using the standard bulk density of refined aggregate (600 kg / m³) as the benchmark value, the actual bulk density values ​​of all aggregate types are converted into relative coefficients. Actual bulk density values ​​are calculated as follows: ,in Based on the standard bulk density of the material type, The comprehensive compaction coefficient. Standardized density coefficient calculation formula: For the concentrate area, Range is 1.0-1.3; coarse material zone The range is 0.65-0.85; in the mixed area. The range is 0.85-0.95. To eliminate abrupt changes in bulk density coefficient between adjacent regions, a 10cm wide transition zone is set at the region boundary. Within the transition zone, a distance-weighted average is used to calculate the transition bulk density coefficient. Finally, a complete three-dimensional bulk density coefficient mapping table is constructed, with each record containing spatial coordinates (x, y, z) and the corresponding standardized bulk density coefficient. This density coefficient mapping table enables a refined expression of density for different locations, types of materials, and degrees of compaction in the material silo, providing key parameters for subsequent mass distribution correction calculations. The innovation of this method lies in combining the physical characteristics of the material type, stacking time, and spatial location to establish a dynamically adaptive density coefficient calculation model, effectively solving the measurement error problem caused by the traditional single density coefficient.

[0116] Table 1:

[0117] Material type Stacking time Surface layer (0-1m) Middle layer (1-2m) Middle layer (2-3m) Bottom layer (>3m) concentrate New arrivals (<5 days) 1.00 1.03 1.08 1.13-1.20 concentrate Mid-term (10-20 days) 1.10 1.14 1.19 1.24-1.28 concentrate Long-term (>30 days) 1.20 1.24 1.28 1.30 coarse materials New arrivals (<5 days) 0.67 0.70 0.73 0.75-0.78 coarse materials Mid-term (10-20 days) 0.73 0.76 0.79 0.81-0.83 coarse materials Long-term (>30 days) 0.80 0.82 0.84 0.85 Mixture New arrivals (<5 days) 0.85 0.87 0.89 0.91-0.93 Mixture Mid-term (10-20 days) 0.88 0.90 0.92 0.93-0.94 Mixture Long-term (>30 days) 0.92 0.93 0.94 0.95 .

[0118] Please refer to Table 1 for the standardized bulk density coefficient ( ) Refer to Table 1 (Standardized density coefficient ( The explanation of the comparison table is as follows:

[0119] 1. Standardization Benchmark: The standard bulk density of refined materials (600 kg / m³) is used as the unified benchmark. The actual bulk density of all material types is standardized using a standardization coefficient. Conversion.

[0120] 2. Calculation method: ;

[0121] Standard bulk density for each type of material (600 kg / m³ for fine materials, 400 kg / m³ for coarse materials, and 500 kg / m³ for mixed materials).

[0122] The time compaction factor approaches its maximum value exponentially with increasing storage time.

[0123] : Depth correction factor, which increases linearly with increasing depth.

[0124] 3. Application Instructions:

[0125] Actual mass calculation: ;

[0126] A 10cm transition zone was set at the boundary of the area, and the change in bulk density coefficient was smoothed by distance-weighted average.

[0127] The compaction effect of fine materials is relatively small, while the compaction effect of coarse materials is significant (up to 50%).

[0128] When stored for a long period of time (>30 days), the compaction effect tends to stabilize, and the coefficient no longer increases significantly.

[0129] 4. Advantages: This dynamic bulk density coefficient system takes into account three dimensions: material characteristics, time variation, and spatial location, and achieves accurate measurement in material storage monitoring, reducing the measurement error from the traditional 15-20% to less than 3%.

[0130] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of the equivalents of the application are intended to be included within the invention.

[0131] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar, characterized in that, Includes the following steps: Step S1: Perform dual-mode synchronous scanning of the material pile using laser-millimeter-wave radar to obtain a spatiotemporally aligned point cloud, which includes penetration echo data and surface measurement data; use the penetration echo data to divide the material pile into depth layers and calculate the interlayer dielectric parameters of the depth layers; perform density gradient anomaly detection based on the interlayer dielectric parameters to obtain a layered density feature map; Step S2: Determine the reliability of the point cloud of the surface measurement data, use the layered density feature map to fill in the missing areas and reconstruct the surface contour to obtain the fused surface contour; Step S3: Spatial data registration is performed between the layered density feature map and the fused surface profile to obtain the registered spatial data; the material bin space is divided into regular voxel grids according to the fused surface profile and the layered density feature map to form density-marked voxel grids; vertical pressure transmission calculation is performed on the density-marked voxel grids to obtain the vertical pressure distribution field; the support strength of the density-marked voxel grids is evaluated according to the vertical pressure distribution field to obtain the support strength matrix; By combining the support strength matrix with the vertical pressure distribution field to determine the support strength and pressure of voxels, the location of weak voxels can be identified. A weak point adjacency graph is established for the locations of weak point voxels; the weak point adjacency graph is filtered and merged into small regions to obtain a set of effective connected regions; the regional geometric features of the effective connected regions are calculated to obtain regional geometric parameters; and the overlying influence is assessed based on the regional geometric parameters and the vertical pressure distribution field to obtain weak area connectivity data. The surrounding support force is extracted from the weak zone connectivity data using the support strength matrix to obtain the weak zone support parameters; the overlying pressure is accumulated based on the weak zone connectivity data and the vertical pressure distribution field to obtain the weak zone pressure load; the safety factor is dynamically determined based on the weak zone support parameters and the weak zone pressure load to obtain the graded safety factor; and the standard risk index is calculated based on the graded safety factor. Risk levels are classified according to standard risk indices to obtain a graded risk assessment table; Spatial risk mapping is performed on the graded risk assessment table to obtain a stress risk distribution map; Step S4: Identify the material type labeling data of the material pile based on the stress risk distribution map; Determine the bulk density coefficient mapping table based on the material type marking data; The bulk density of the material pile is dynamically corrected using a bulk density coefficient mapping table to obtain corrected storage quantity data.

2. The method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar according to claim 1, characterized in that, Step S1, which involves simultaneous dual-modal scanning of the material pile, includes: The dual radar response signals of the material pile are collected simultaneously by lidar and millimeter-wave radar. The coordinate system is transformed and the timestamp is checked to obtain the checked response signal; Radar point cloud data fusion is performed on the calibrated response signal to obtain a spatiotemporally aligned point cloud.

3. The method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar according to claim 1, characterized in that, Step S1 involves dividing the material pile into depth layers and calculating the interlayer dielectric parameters of each depth layer, including: Penetration echo data from millimeter-wave radar is extracted from spatiotemporally aligned point clouds, and echo sequence separation is performed to form multi-layer echo pulse sequences; When the time interval between two adjacent echoes in a multi-layer echo pulse sequence is greater than the preset echo time threshold, it is determined to be an independent reflection layer. Based on the reflection layer, the material pile is divided into multiple depth layers of varying thicknesses from the surface to the bottom, and the boundary coordinates of the depth layers are recorded. The interlayer dielectric parameters of the material stack are calculated based on the depth layer boundary coordinates.

4. The method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar according to claim 1, characterized in that, Step S1, which involves density gradient anomaly detection based on interlayer dielectric parameters, includes: The dielectric change rate of the interlayer dielectric parameters is calculated to obtain the interlayer change rate matrix; Data from the central region of the material pile is selected from the interlayer change rate matrix to establish a normal change benchmark and obtain the stratification benchmark range. Identify outlier locations in the inter-layer rate of change matrix based on the stratified baseline range; Based on the interlayer dielectric parameters, the void characteristics of the abnormal point location are determined, and a void type classification table is obtained; Aggregate abnormal regions based on the cavity type classification table to obtain aggregated abnormal regions; Density anomaly marker data is generated based on aggregated anomaly regions; Based on the density anomaly marker data, density level values ​​are assigned to the depth layers to obtain the layered density feature map.

5. The method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar according to claim 1, characterized in that, Step S2 includes: Surface measurement data from lidar are extracted from spatiotemporally aligned point clouds to perform lidar point cloud quality assessment and obtain point cloud reliability distribution. The missing regions of the point cloud reliability distribution are identified to obtain the set of missing region boundaries; Based on the missing region boundary set, millimeter-wave radar surface data is extracted from the hierarchical density feature map to fill in the missing regions, resulting in a supplementary height dataset. The point cloud reliability distribution and the supplementary height dataset are fused and stitched together to obtain a complete height point matrix. A smooth surface reconstruction is performed on the complete height lattice to obtain the fused surface profile.

6. The method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar according to claim 1, characterized in that, Step S4 includes: Based on the stress risk distribution map, the space of the material pile is divided into regions with different density levels, resulting in a density zoning map; Material type identification and classification are performed based on density zoning maps to obtain material type labeling data; Calculate the partition volume data based on the density partition map and the fused surface profile; Mass distribution correction calculations are performed on the partitioned volume data and bulk density coefficient mapping table to obtain the mass distribution table; The storage capacity of the mass distribution table is summarized and calibrated to obtain the corrected storage capacity data.

7. The method for precise monitoring of stored materials based on the fusion of lidar and millimeter-wave radar according to claim 1, characterized in that, Step S4, which involves determining the bulk density coefficient mapping table based on the material type marking data, includes: Extract the reference bulk density of the material type from the material type labeling data; Obtain the material's entry time into the warehouse, and calculate the stacking time factor by combining the material type's benchmark bulk density; The comprehensive compaction coefficient is obtained by correcting the depth compaction based on the material type marking data and the time compaction factor. The bulk density coefficient is standardized by combining the comprehensive compaction coefficient and the reference bulk density of the material type to obtain the standardized bulk density coefficient.

Images