Warehousing material accurate monitoring method based on fusion of laser radar and millimeter wave radar

Abstract

The invention relates to the technical field of multi-sensor fusion, in particular to a storage material accurate monitoring method based on fusion of a laser radar and a millimeter wave radar. The method comprises the following steps: performing bimodal synchronous scanning on a material pile through a laser-millimeter wave radar to obtain penetration echo data and surface measurement data; dividing the material pile into depth layers by using the penetration echo data, and calculating interlayer dielectric parameters of the depth layers; density gradient anomaly detection is carried out according to the interlayer dielectric parameters, and a layered density feature map is obtained; judging the point cloud reliability of the surface measurement data, and performing missing region supplementation and surface contour reconstruction by using the hierarchical density feature map to obtain a fused surface contour; and carrying out support stress distribution calculation on the layered density feature map and the fused surface profile to obtain a vertical pressure distribution field and a support strength matrix. According to the method, the storage capacity calculation precision of material storage is remarkably improved through deep hierarchical analysis and internal and external structure coupling analysis.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of multi-sensor fusion, and in particular to a warehouse material precision monitoring method based on laser radar and millimeter wave radar fusion. BACKGROUND

[0002] The existing material warehouse monitoring technology cannot adapt to the special environment of the material warehouse using a single sensor. High-concentration dust generated during material processing and storage can seriously interfere with the laser radar signal, causing data loss or distortion. Alternative technologies such as infrared and ultrasonic waves also have measurement blind spots in a dusty environment, resulting in low reliability of the monitoring system and the inability to continuously provide complete material pile state information. Traditional material level measurement technology can only obtain surface profile information and cannot penetrate the material surface to detect the internal structure. When the material is unevenly stacked or long-term storage causes internal cavities and loose areas, there is a lack of effective identification and early warning mechanism, and collapse accidents often occur without warning, causing personnel casualties and equipment damage. Existing technologies usually use a single bulk density coefficient to convert volume to weight, ignoring the density differences of different types of materials (fine materials, coarse materials, etc.). And the vertical density gradient change caused by the compaction effect is not considered. In the case of mixed storage, the measurement error can be as high as 15-20%, which seriously affects inventory management and cost control.

[0003] In summary, the existing technology has problems such as poor sensor environmental adaptability, lack of internal structure detection capability, and insufficient measurement accuracy, which need to be solved. SUMMARY

[0004] Therefore, it is necessary to provide a warehouse material precision monitoring method based on laser radar and millimeter wave radar fusion to solve at least one of the above technical problems.

[0005] To achieve the above purpose, a warehouse material precision monitoring method based on laser radar and millimeter wave radar fusion includes the following steps:

[0006] Step S1: Dual-mode synchronous scanning of the material pile is performed by laser-millimeter wave radar to obtain a space-time aligned point cloud, wherein the space-time aligned point cloud includes penetration echo data and surface measurement data; the material pile is divided into depth layers using the penetration echo data, and the interlayer dielectric parameters of the depth layers are calculated; density gradient anomaly detection is performed according to the interlayer dielectric parameters to obtain a layered density feature map;

[0007] Step S2: The point cloud reliability of the surface measurement data is determined, the layered density feature map is used for missing area supplementation and surface profile reconstruction, and a fused surface profile is obtained.

[0008] Step S3: spatially register the hierarchical density feature map with the fused surface profile, and calculate the support stress distribution of each spatial unit to obtain a vertical pressure distribution field and a support strength matrix; identify the voxel weak area of the support strength matrix according to the vertical pressure distribution field, and perform weak area connectivity analysis to obtain weak area connectivity data; perform quantitative evaluation of the collapse risk according to the weak area connectivity data to obtain a stress risk distribution map;

[0009] Step S4: identify the material type marker data of the material pile based on the stress risk distribution map; determine the bulk density coefficient mapping table according to the material type marker data; and dynamically correct the partition bulk density of the material pile by using the bulk density coefficient mapping table to obtain corrected storage quantity data.

[0010] The application realizes the non-destructive penetration detection of the internal structure of the material bin in the dust environment through the dual-mode synchronous scanning, and overcomes the limitations of a single sensor. The laser radar provides millimeter-level surface precision, and the millimeter wave radar provides centimeter-level penetration capability, and the two complement each other to form comprehensive perception. The space-time alignment point cloud ensures the consistency of the data in the time and space dimensions, laying a foundation for subsequent analysis. The deep layer division and the interlayer dielectric parameter calculation realize the quantitative representation of the internal physical characteristics of the material, and the density gradient anomaly detection accurately identifies the internal cavities and loose areas, which are difficult to find in traditional detection methods. This step realizes the visualization of the internal state of the sealed warehouse through electromagnetic wave physical property analysis, provides early risk identification capability, and prevents safety accidents caused by internal structure changes. The fine reconstruction of the surface contour effectively solves the problem of missing laser data in high-concentration dust environments, and improves the reliability of the monitoring system in harsh environments. The point cloud quality evaluation mechanism automatically identifies low-reliability areas, and the millimeter wave data supplement technology realizes effective filling of the measurement blind area. The distance weighted average method is adopted for the dual-source data fusion splicing, ensuring the data continuity of the transition area and avoiding the height mutation at the splicing position. The curved surface smoothing reconstruction algorithm generates a continuous surface model with millimeter-level precision, and adaptively adjusts the window size for areas with large curvature changes to retain the detailed features. The fused surface contour provides a complete and accurate description of the external geometric shape of the material pile, providing a reliable basis for volume calculation and internal-external structure coupling analysis, and realizing the continuous and stable operation of the monitoring system in harsh environments. The internal-external structure coupling analysis establishes the correlation between the internal density distribution and the surface morphology of the material, realizes the transformation from static structural characteristics to dynamic mechanical analysis. Spatial registration ensures the strict correspondence of internal and external data, supports stress distribution calculation, simulates the actual pressure transmission process, and considers the stress dispersion effect. Weak area connectivity analysis aggregates discrete weak points into physically meaningful structural units, and enhances the relevance of the analysis through shape factors and overburden influence evaluation. The quantitative risk assessment of collapse considers four key factors: weak area volume, overburden pressure, surrounding support capacity and safety factor, and converts qualitative judgment into quantitative indicators. Risk level classification and three-dimensional visualization expression enable operators to intuitively identify dangerous positions, providing a scientific basis for safety management decisions. This step realizes accurate detection of internal cavities and quantitative early warning of collapse risk, fundamentally improving the safety level of warehousing. The dynamic correction of the partition bulk density overcomes the limitations of traditional single bulk density coefficient calculation methods, significantly improving the measurement accuracy under mixed storage conditions. Density partitioning and material identification realize accurate classification of different types of materials, and the bulk density coefficient determination fully considers the influence of three key factors: storage time, compaction degree and depth position. The time compaction factor reflects the natural compaction effect of the material with the increase of storage time, the depth compaction correction considers the vertical pressure gradient change, and the standardized bulk density coefficient ensures consistent comparison between different materials.The accurate partition volume calculation is combined with the bulk density coefficient to realize accurate calculation of the partition quality. The calibration mechanism guarantees the result reliability by comparison with historical records, reduces the measurement error from 15-20% of the traditional method to within 3%, provides accurate data support for inventory management, cost control and production planning, and improves the accuracy of material quality management.

[0011] Therefore, the application provides a warehouse material precise monitoring method based on fusion of laser radar and millimeter wave radar, realizes data complementation and automatic switching in a dusty environment, innovatively adopts penetration depth hierarchical analysis and internal and external structure coupling analysis process, can accurately identify internal cavities and quantitatively evaluate collapse risk, and simultaneously introduces a dynamic bulk density correction mechanism based on spatial partition, applies differential bulk density coefficients for different density materials, and significantly improves the storage quantity calculation accuracy. This multi-dimensional fusion method provides a new technical path for warehouse material precise monitoring. BRIEF DESCRIPTION OF DRAWINGS

[0012] Figure 1 FIG. 1 is a step flow diagram of a warehouse material precise monitoring method based on fusion of laser radar and millimeter wave radar.

[0013] The object implementation, functional features and advantages of the application will be further described with reference to the embodiments and the accompanying drawings. DETAILED DESCRIPTION

[0014] The technical method of the application will be described in detail below with reference to the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the application, rather than all the embodiments. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor fall within the protection scope of the application.

[0015] In addition, the accompanying drawings are only schematic diagrams of the application, and are not necessarily drawn to scale. The same reference signs in the drawings represent the same or similar parts, and thus repeated description thereof will be omitted. Some block diagrams shown in the drawings are functional entities, which do not necessarily correspond to physically or logically independent entities. The functional entities can be implemented in software form, or 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," and the like may be used herein to describe various elements, these elements should not be limited by these terms. These terms are used solely to distinguish one element from another. For example, a first element may be referred to as a second element, and similarly, a second element may be referred to as a first element, without departing from the scope of the exemplary embodiments. The term "and / or" as used herein includes any and all combinations of one or more of the listed associated items.

[0017] To achieve this, please refer to Figure 1 The present invention provides a method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar, comprising the following steps:

[0018] Step S1: Perform dual-modal synchronous scanning of the material pile using a laser-millimeter wave radar to obtain a spatiotemporally aligned point cloud, wherein the spatiotemporally aligned point cloud includes penetration echo data and surface measurement data; divide the material pile into depth layers using the penetration echo data, 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 the embodiment of the present invention, dual-mode synchronous scanning uses hardware triggers to achieve millisecond-level synchronous acquisition of laser radar and millimeter-wave radar, and ensures strict alignment of dual-source data in time and space dimensions through coordinate transformation matrix and timestamp calibration. The laser radar provides 1.2 million points of high-precision surface data, and the millimeter-wave radar provides a 140GHz penetration echo sequence, forming a time-space aligned point cloud containing surface geometry and internal structure information. The depth layer division is based on the echo time interval analysis, and a time threshold of 0.5 nanoseconds is set to determine the independent reflection layer. The formula Calculate the actual thickness and divide it into 5-8 layers of varying thickness. The dielectric parameters between layers are calculated by the reflection coefficient. Calculate, where and Respectively Layer and The relative dielectric constant of the layer, taking into account the propagation distance correction factor. Density gradient anomaly detection calculates the inter-layer change rate , a normal change benchmark is established based on the data of the central area of ​​the pile, abnormal points beyond the benchmark range are marked, the cavity type is determined according to the dielectric constant value, aggregated abnormal areas are formed through spatial clustering, and density level values ​​(0-10) are 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 the embodiment of the present application, the surface profile fine reconstruction extracts the surface measurement data of the laser radar from the space-time aligned point cloud, and performs point cloud quality evaluation through reflection intensity analysis. A reflection intensity reference value 150 is set, and when the intensity is lower than 60, it is marked as a low reliability point, and the effective point proportion of each grid unit is calculated When is lower than 50%, it is marked as a missing unit. The region growing algorithm is used to aggregate adjacent missing units, and the area exceeding 100cm² is determined as a data missing area, and the boundary point set is extracted. For each missing area, the millimeter wave surface data is extracted from the hierarchical density feature map, 5cm×5cm gridding processing and median filtering are performed, and a supplementary height data set is formed. In the 20cm wide transition zone, distance weighted average method is used to realize smooth fusion of double source data, and the formula is as follows:

[0022] ;

[0023] Wherein represents the height value after fusion, represents the height measured by laser radar, represents the height measured by millimeter wave radar, and are weight coefficients of laser data and millimeter wave data respectively. Finally, local quadratic surface fitting is used for surface smoothing reconstruction, and the fitting model is obtained, the window size is adaptively adjusted for the edge area, and finally the fusion surface profile with 5mm resolution is generated.

[0024] Step S3: Space registration is performed on the hierarchical density feature map and the fusion surface profile, and the support stress distribution of each space unit is calculated to obtain the vertical pressure distribution field and the support intensity matrix; the voxel weak area of the support intensity matrix is identified according to the vertical pressure distribution field, and weak area connectivity analysis is performed to obtain weak area connectivity data; the weak area connectivity data is used for quantitative evaluation of the collapse risk to obtain a stress risk distribution map;

[0025] In the embodiment of the present application, the internal and external structure coupling analysis performs space registration on the hierarchical density feature map and the fusion surface profile, and the registration accuracy is required to reach 5mm, and three-dimensional rigid body transformation is used to eliminate position deviation. The material pile space is divided into regular voxel grids with a size of 10cm×10cm×10cm, and the material density value (0-750kg / m³) is allocated according to the density level. The vertical pressure transmission calculation is based on voxel gravity accumulation, and the pressure dispersion effect of Gaussian distribution is considered, and the calculation formula is as follows:

[0026] ;

[0027] Wherein represents the vertical pressure at position ; Indicates the quality of the voxel at that position, Indicates the acceleration due to gravity is 9.8m / s², represents the pressure transmission weight coefficient, The position 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 area connectivity analysis screens weak voxels with a safety margin SM < 0, establishes a 26-neighborhood weakness adjacency graph, filters and merges isolated weaknesses with less than 8 voxels, calculates characteristic parameters such as geometric center, volume, and shape factor, and assesses the impact of the overburden layer. Quantitative assessment of collapse risk extracts the peripheral support strength Ssurr, accumulates the overburden pressure Ptotal, dynamically determines the safety factor SF based on the volume and shape, and applies the risk formula Calculate the standard risk index, where represents the risk index, represents the volume of the weak zone (m³), Indicates the total overburden pressure (Pa), Indicates the peripheral support strength (Pa), Indicates the safety factor. Risk levels are divided into four levels: 0-30 (low), 30-60 (medium), 60-80 (high), and 80-100 (very high), generating a three-dimensional stress risk distribution map.

[0028] Step S4: identifying material type identification data of the material pile based on the stress risk distribution map; determining a bulk density coefficient mapping table based on the material type identification data; dynamically correcting the partition bulk density of the material pile using the bulk density coefficient mapping table to obtain corrected storage capacity data;

[0029] In the embodiment of the present invention, the dynamic correction of bulk density is based on the stress risk distribution map to divide the density area, set clear standards for low density area (<300kg / m³), medium density area (300-500kg / m³) and high density area (>500kg / m³), and use gradient analysis to determine the area boundaries. Material type identification is classified according to the density-material type comparison table: >520kg / m³ is fine feed, <430kg / m³ is coarse material, 430-520kg / m³ is mixed material. The bulk density coefficient is determined by extracting the benchmark bulk density from the material type marking data (600kg / m³ for fine material, 400kg / m³ for coarse material, 500kg / m³ for mixed material), and calculating the time compaction factor in combination with the storage time. ,in is the maximum compaction gain (0.3 for fine material and 0.5 for coarse material), is the compaction rate (0.1 / day), Indicates the stockpile time (days). The depth compaction correction divides the stockpile into the surface layer ( ), middle level ( ) and the bottom layer ( ), calculate the integrated compaction coefficient . Normalized bulk density coefficient . Volume calculation uses 5cm voxel method to accurately measure the volume of each partition, and the mass calculation formula , wherein represents the partition corrected mass (kg), represents the partition volume (m³), represents the reference density of concentrate 600 kg / m³, represents the bulk density coefficient of the material, represents the depth correction coefficient. Finally, the mass of each partition is summarized, and when the calculation deviation exceeds ± 3%, a calibration coefficient is introduced for correction to generate corrected storage data containing total storage, distribution ratio and time stamp.

[0030] Preferably, the double-mode synchronous scanning of the material pile in step S1 comprises:

[0031] Synchronously collecting double-radar response signals of the material pile by laser radar and millimeter wave radar;

[0032] Coordinate system transformation and time stamp correction are performed on the double-radar response signals to obtain corrected response signals;

[0033] Radar point cloud data fusion is performed on the corrected response signals to obtain a time-space aligned point cloud.

[0034] In this embodiment, the double-mode synchronous scanning first realizes millisecond-level synchronous collection of laser radar and millimeter wave radar through a hardware trigger. The trigger uses a high-precision timing circuit to generate a 10-microsecond-wide rectangular pulse signal, and the pulse rising edge triggers both radar devices at the same time. The laser radar is configured in line scanning mode, with a scanning frequency of 10Hz, a single scanning angle range of 0° to 180°, an angular resolution of 0.25°, and a distance resolution of 5mm. High-density surface point cloud data of 1.2 million points are collected for each scan. The millimeter wave radar works in the 140GHz frequency band and uses the FMCW (Frequency Modulated Continuous Wave) mode, with a bandwidth of 6GHz. The scanning period is strictly synchronized with the laser radar, and the penetration depth reaches 3 meters with a depth resolution of 10mm. Each transmission receives 16 echo signals to form an echo sequence. The double-radar response data are all added with a unified time stamp of the trigger time, with a precision of microseconds.

[0035] The coordinate system transformation is realized by using a rigid body transformation matrix. According to the installation position relationship (horizontal distance of 50 cm, vertical height difference of 30 cm, and azimuth angle deviation of 5°) of the two radars, a unified Cartesian coordinate system is established. The original data of the millimeter wave radar is in polar coordinate form (r, θ, φ), which is converted into Cartesian coordinates (x, y, z) in the unified coordinate system through the coordinate transformation matrix T. The transformation matrix T includes a rotation matrix R and a translation vector t, which satisfies 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 transformation, all measurement data is unified into a global coordinate system with the warehouse center as the origin. Then, the time stamp is corrected, and the time stamp difference of the two radar data packets is detected. When the difference exceeds 1 millisecond, the sampling points are linearly interpolated and corrected. The correction formula is: t_corrected=t_original+Δt, where Δt is the time deviation. After time synchronization, each spatial position has double radar measurement data with strictly aligned time stamps.

[0036] The radar point cloud data fusion is realized by constructing a unified data structure. The data structure includes: three-dimensional space coordinates (x, y, z), laser reflection intensity value I_laser (range 0-255), millimeter wave echo intensity I_mmw (range 0-100 dB), millimeter wave penetration depth array D_penetration (containing multiple reflection point depth values), and measurement timestamp t. For each spatial position point, the laser radar provides accurate surface coordinates and reflection intensity, and the millimeter wave radar provides penetration depth information and internal multi-layer echo intensity at the corresponding position. When the laser data of a point is missing (reflection intensity is lower than the threshold value 40), the surface reflection point of the millimeter wave data is automatically taken as the coordinate value of the position. The global point cloud density is 1000 points per square meter, and the total number of points is maintained at the order of 1.2 million points, forming a complete time and space aligned point cloud containing surface fine features and internal penetration information. The point cloud data realizes the strict alignment of the two radars in time and space dimensions, laying a foundation for subsequent depth layering analysis and internal and external structure coupling analysis.

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

[0038] The penetration echo data of the millimeter wave radar is extracted from the space-time aligned point cloud, the echo sequence is separated, and a multi-layer echo pulse sequence is formed;

[0039] When the time interval of two adjacent echoes in the multi-layer echo pulse sequence is greater than the preset echo time threshold, it is determined that it is an independent reflection layer, and the material pile is divided into multiple depth layers with different thicknesses from the surface to the bottom according to the reflection layer, and the depth layer boundary coordinates are recorded;

[0040] The interlayer dielectric parameters of the material pile are calculated according to the depth layer boundary coordinates.

[0041] In the embodiment of the application, the penetration echo data of the millimeter wave radar is extracted from the space-time aligned point cloud. For each spatial scanning point, the multi-layer echo information stored in the D_penetration array thereof is extracted. The echo data is processed by using a time domain analysis method, and a signal intensity threshold of -60 dBm is set. When the echo signal intensity is higher than the threshold, it is recorded as an effective echo point. The time domain echo signal is converted into a frequency domain by using a fast Fourier transform, a main frequency component is identified, and noise interference is filtered out. The envelope of the processed echo signal is detected, and the peak point position is extracted. Each peak value corresponds to a reflection interface inside the material. In time sequence, the peak points are organized into an echo pulse sequence, and each element in the sequence contains an echo arrival time t_i and a corresponding echo intensity I_i. For each scanning point position, an independent echo pulse sequence is generated, and a multi-layer echo distribution matrix in three-dimensional space is formed.

[0042] An independent reflection layer is determined according to the multi-layer echo pulse sequence. The echo time threshold is set to 0.5 nanoseconds. When the time interval Δt=t_i+1-t_i of two adjacent echoes is greater than the threshold, it is determined that it is a reflection interface between different depth layers. The propagation speed v of the millimeter wave 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 according to the formula d=v×Δt / 2, wherein v is the propagation speed of the electromagnetic wave in the medium, Δt is the echo time difference, and 2 is divided because the signal propagates back and forth. According to the calculation result, the material pile is divided into multiple depth layers with different thicknesses from the surface to the bottom. Generally, the material pile is divided into 5 to 8 depth layers. The thickness of the surface layer is usually 15-20 cm, the thickness of the middle layer is 20-40 cm, and the thickness of the bottom layer is 30-50 cm. The upper and lower boundary coordinates of each depth layer are recorded. The upper boundary coordinates are (x, y, z_upper), and the lower boundary coordinates are (x, y, z_lower). The depth layer boundary coordinate set is constructed.

[0043] The interlayer dielectric parameters of the material pile are calculated according to the depth layer boundary coordinates. The relative permittivity of each layer of material is calculated by using the principle of electromagnetic wave interface reflection through the intensity ratio of the incident wave and the reflected wave. For the interface between the first layer and the second layer, the reflection coefficient is , wherein and are the first layer and the second layer, respectively. The relative permittivity of the first layer is calculated according to the formula Relative dielectric constant of the layer. In practical applications, the intensity of the incident wave and the intensity of the reflected wave are measured and the reflection coefficient formula is used to calculate the dielectric constant of the layer Inverse dielectric constant. Considering the attenuation of electromagnetic waves during propagation, a propagation distance correction factor is introduced wherein is the attenuation coefficient, is the propagation distance. For the coarse material area, the relative dielectric constant ranges from 1.5 to 2.5; for the fine material area, the dielectric constant ranges from 3.0 to 4.5; for the area with high water content, the dielectric constant can reach 5.0-7.0. By calculating the dielectric constant difference between adjacent layers, a complete interlayer dielectric parameter table is established, laying the foundation for subsequent density gradient anomaly detection.

[0044] Preferably, the density gradient anomaly detection according to the interlayer dielectric parameters in step S1 comprises:

[0045] Calculating the dielectric change rate of the interlayer dielectric parameters to obtain an interlayer change rate matrix;

[0046] Selecting the center area data of the material pile from the interlayer change rate matrix to establish a normal change benchmark and obtain a layered benchmark range;

[0047] Identifying the abnormal point position of the interlayer change rate matrix according to the layered benchmark range;

[0048] Judging the cavity characteristics of the abnormal point position according to the interlayer dielectric parameters to obtain a cavity type classification table;

[0049] Aggregating the abnormal areas according to the cavity type classification table to obtain aggregated abnormal areas;

[0050] Generating density anomaly marker data based on the aggregated abnormal areas;

[0051] Assigning a density level value to the depth layer according to the density anomaly marker data to obtain a layered density feature map.

[0052] ​In the embodiment of the present application, the density gradient anomaly detection first calculates the rate of change of the interlayer dielectric parameter table. For each spatial measurement point (x, y) position, extract the dielectric constant sequence {ε_1, ε_2,..., ε_n} of each depth layer, calculate the dielectric constant change rate CR_i=(|ε_i-ε_i+1| / ε_i+1)×100% between adjacent two layers, wherein CR_i represents the dielectric constant change rate of 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, and an interlayer change rate matrix with a dimension of m×n is constructed, wherein m is the number of scanning points, and n is the number of depth layers minus 1. Then, the center area data of the stockpile is extracted from the interlayer change rate matrix to establish a normal change benchmark. The center area is defined as all measurement points within a range of 40% of the maximum radius of the stockpile from the horizontal center point of the stockpile. The mean value μ_i and the standard deviation σ_i of the change rate of each depth layer in the region are calculated, and the normal change range is set as [μ_i-2σ_i, μ_i+2σ_i]. The benchmark range is established for different depth layers, and the benchmark range of the surface layer to the second layer is usually 15%±10%, the benchmark range between the middle layers is 10%±8%, and the benchmark range of the deep layer area is 5%±5%. According to the established layered benchmark range, the values in the interlayer change rate matrix are compared point by point. When the change rate CR_i of a measurement point in the i-th layer exceeds the upper limit of the benchmark range of the corresponding depth, the measurement point is marked as an abnormal point. In particular, when the change rate exceeds 30%, it is directly determined as a strong abnormal point. The three-dimensional coordinates (x, y, z_i) of each abnormal point and the depth layer i where the abnormal point is located are recorded, and an abnormal point position set is constructed.

[0053] For the marked abnormal point position, the interlayer dielectric parameter value is checked back to determine the cavity characteristics. The cavity type is divided into three categories: when the dielectric constant ε_i<1.5, it is determined as a complete cavity, and the type is marked as 0; when 1.5≤ε_i<2.5, it is determined as a loose area, and the type is marked as -1 to -3, and the specific value is linearly mapped according to (2.5-ε_i); when ε_i≥2.5 and the difference with the adjacent layer is significant, it is determined as a density mutation area, and the type is marked as 1. The spatial coordinates, the depth layer and the corresponding cavity type of each abnormal point are recorded to form a cavity type classification table. Subsequently, the abnormal area is aggregated, and a spatial proximity clustering algorithm is used. The aggregation threshold is set to 20 cm, and when the horizontal distance between two abnormal points is less than the threshold and the depth layer difference is not more than 1 layer, they are classified into the same aggregation area. For each aggregation area, the center coordinates are the average value of the abnormal point coordinates contained therein, the range size is the volume of the smallest circumscribed ellipsoid containing all points, and the main abnormal type is the cavity type with the highest occurrence frequency. The aggregation result is stored in the aggregated abnormal area table, including the fields of area ID, center coordinates, volume size, main abnormal type, etc.

[0054] The density anomaly label data is generated based on the aggregated abnormal area. The monitoring space is divided into regular grids of 10 cm x 10 cm x 10 cm, and for each grid unit, it is checked whether it is located in a certain aggregated abnormal area. If it is located in the abnormal area, the corresponding abnormal type label (0, -1 to -3 or 1) is given; if it is not in the abnormal area, the density level of the position is determined according to the dielectric constant value ε_i. The mapping relationship between dielectric constant and density is set: ε_i=1.5~2.0 corresponds to density level 1~2 (extremely low density), ε_i=2.0~3.0 corresponds to density level 3~5 (low density), ε_i=3.0~4.0 corresponds to density level 6~8 (medium density), and ε_i=4.0~7.0 corresponds to density level 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 given according to the density anomaly label data, and a complete layered density feature map is constructed. The feature map is represented in the form of a three-dimensional color map, with different colors representing different density levels, red representing the hollow area (level 0), yellow representing the loose area (level 1-3), green representing the normal area (level 4-7), and blue representing the compact area (level 8-10). This density gradient anomaly detection method realizes non-destructive detection of the internal structure of the material through dielectric parameter analysis, especially accurate identification of potential voids.

[0055] Preferably, step S2 comprises:

[0056] The surface measurement data of the laser radar is extracted from the spatio-temporal aligned point cloud, the quality of the laser point cloud is evaluated, and the point cloud reliability distribution is obtained;

[0057] The point cloud reliability distribution is identified for data missing area, and a missing area boundary set is obtained;

[0058] The millimeter wave radar surface layer data is extracted from the layered density feature map according to the missing area boundary set to supplement the missing area, and a supplemented height data set is obtained;

[0059] The point cloud reliability distribution and the supplemented height data set are fused and spliced to obtain a complete height point array;

[0060] The complete height point array is subjected to surface smoothing reconstruction to obtain a fused surface contour.

[0061] In the embodiment of the present application, the surface layer profile fine reconstruction first extracts the surface measurement data of the laser radar from the space-time aligned point cloud, and performs point cloud quality evaluation. The three-dimensional coordinates (x, y, z) and the reflection intensity value I laser of each point are extracted, and the data quality is judged by analyzing the reflection intensity. The reflection intensity reference value I baseline is set to 150 (range 0-255), and when the reflection intensity I laser of the point is lower than 40% (i.e. 60) of the reference value, it is marked as a low reliability point, indicating that the point is disturbed by dust. The valid point ratio in 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%, and the point cloud reliability distribution matrix is constructed. Then, the point cloud reliability distribution is analyzed to identify the missing area. When the valid point ratio R valid of a certain grid unit is lower than 50%, it is marked as a potential missing unit. The region growing algorithm is used to aggregate adjacent missing units, and when the continuous area exceeds 100cm², it is determined as a data missing area. The alpha-shape algorithm is used to extract the boundary point set of each missing area, and the three-dimensional coordinates (x_b, y_b, z_b) of the boundary points are recorded to form the missing area boundary set. Each missing area records its ID number, center position, area size and boundary point list.

[0062] The millimeter wave radar surface layer data is extracted from the hierarchical density feature map according to the missing area boundary set. For each missing area, determine its horizontal range (x_min, x_max, y_min, y_max), and extract the first layer (surface layer) data in the hierarchical density feature map within the range. The millimeter wave radar surface layer data includes spatial position (x_m, y_m) and corresponding surface height value z_m. Since the horizontal resolution of the millimeter wave radar (about 5cm) is lower than that of the laser radar (about 0.5cm), grid processing is needed. The missing area is 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 median filter is applied with a window size of 3×3 to filter out outliers. The processed height data is used as a supplementary height data set to fill in the missing areas of the laser data.

[0063] The dual-source data fusion stitching realizes the seamless connection of the laser point cloud and the millimeter wave supplementary data. To ensure smooth transition of data, a transition zone with a width of 20 cm is set at the boundary of the missing area. In the transition zone, the distance weighted average method is used to calculate the fused height value: z_fused=(w_l×z_laser+w_m×z_mmw) / (w_l+w_m), wherein w_l and w_m are weight coefficients of the laser data and the millimeter wave data respectively. The weight coefficients are calculated according to the distance d of the point to the boundary of the missing area: w_l=(20-d) / 20 (when d<20cm), w_m=d / 20 (when d<20cm); when d≥20cm, w_l=0, w_m=1 in the missing area, and w_l=1, w_m=0 outside the missing area. The smooth transition of the two data sources is realized by this method, and the height mutation at the stitching position is avoided. The fusion result is reorganized into a regular height point array with a resolution of 2cm×2cm, covering the entire stockpile surface area.

[0064] The complete height point array is subjected to curved surface smoothing reconstruction 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 the points within a 5x5 grid range around it. In the fitting process, in order to reduce the influence of noise, different weights are given to each point: the weight of the center point is 1.0, and the weight decreases linearly with the increase of distance. For areas with large curvature changes (such as the edges of the stockpile), the size of the fitting window is adjusted adaptively, and a 3x3 window is used in the edge area to preserve the detailed features. After the fitting is completed, the global surface is subjected to spline interpolation to generate a uniform grid height field with a resolution of 5mm. The surface reconstruction accuracy reaches the millimeter level, meeting the demand for accurate calculation of the volume of the stockpile. The fused surface profile contains three-dimensional coordinate and surface normal vector information, which completely describes the external geometric shape of the material pile and provides basic data for subsequent coupling analysis with the internal structure. The innovation point here is to use millimeter wave radar data to supplement the measurement blind area of the laser radar in the dust environment, realizing the complementary advantages of dual-mode sensors.

[0065] Preferably, the spatial registration of the layered density feature map and the fused surface profile in step S3, and the calculation of the support stress distribution of each spatial unit include:

[0066] The spatial data registration of the layered density feature map and the fused surface profile is performed to obtain registered spatial data;

[0067] According to the fused surface profile and the layered density feature map, the space of the material bin is divided into a regular voxel grid to form a density-labeled voxel grid;

[0068] The vertical pressure transmission calculation is performed on the density-labeled voxel grid to obtain a vertical pressure distribution field;

[0069] The support strength of the density-marked voxel grid is evaluated according to the vertical pressure distribution field to obtain the support strength matrix.

[0070] In the embodiment of the present invention, the internal and external structure coupling analysis first performs spatial data registration of the layered density feature map and the fused surface contour. Taking the fixed reference point of the warehouse as the benchmark, the three-dimensional rigid body transformation method is used to eliminate the position deviation between the two data. The transformation includes the translation vector t and the rotation matrix R, so that the two sets of data are strictly aligned to a unified coordinate system. The alignment accuracy is required to be within 5mm. Exceeding this error will lead to subsequent analysis deviations. After implementing the rigid body transformation, the consistency of the surface contour data and the depth value of the density feature map layer is checked, and the root mean square error is calculated. ,when When the distance is >5mm, fine correction is performed on the local area by bilinear interpolation. After registration, each point of the surface contour ( ) and the surface points at the corresponding positions in the density feature map ( ) strictly corresponds to, and The difference does not exceed 5mm. The registration results are stored as registration space data in a unified three-dimensional coordinate system, including spatial position, surface height and internal density distribution information.

[0071] The warehouse space is divided into a regular voxel grid based on the registered spatial data to form a density-labeled voxel grid. The voxel size is set to 10cm×10cm×10cm, which is fine enough to capture local structural changes without causing excessive computational overhead. Each voxel is assigned a unique 3D index ( ), corresponding to the actual space coordinates ( ),in( ) is the coordinate origin. For each voxel, determine whether it is located inside the pile: if the height of the voxel center point is lower than the surface contour 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) of the corresponding location is extracted from the layered density feature map and recorded in the density-labeled voxel grid data structure. The density-labeled voxel grid is stored in 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 voxel-based gravity accumulation principle. First, assign a corresponding material density value to each density level: density level 0 (void) corresponds to 0 kg / m³; density level 1-3 (loose area) corresponds to 200-350 kg / m³; density level 4-7 (normal area) corresponds to 350-550 kg / m³; density level 8-10 (compacted area) corresponds to 550-750 kg / m³. The specific density value is determined by linear interpolation, such as density level 5 corresponds to 425 kg / m³. Calculate the mass of each voxel: m = p x V, where p is the material density and V is the voxel volume 0.001 m³. For surface voxels, the vertical pressure they receive is only their own weight: P = m x g, g = 9.8 m / s². For internal voxels, accumulate the gravity contribution of all voxels above them, considering the stress dispersion effect. The vertical pressure transmission equation is: P(i,j,k) = m(i,j,k) x g + ∑w(i',j',k'-1) x P(i',j',k'-1), where (i',j',k'-1) is the adjacent voxel above, w is the pressure transmission weight coefficient, which decreases with increasing horizontal distance and obeys Gaussian distribution. The weight calculation formula is: where is the horizontal distance, is the standard deviation. The calculation result forms a three-dimensional vertical pressure distribution field P(i,j,k) with a unit of Pa.

[0073] The support strength evaluation calculates the maximum bearing capacity of each voxel according to its density level. Establish a mapping relationship between density level and support strength: density level 0 (void) has a support strength of 0 Pa; density level 1-3 (loose area) has a support strength of 10000-30000 Pa; density level 4-7 (normal area) has a support strength of 30000-60000 Pa; density level 8-10 (compacted area) has a support strength of 60000-90000 Pa. Determine the specific support strength value S(i,j,k) through piecewise linear mapping. Calculate the bearing safety factor of each voxel: 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, and there is a structural risk. Store the calculation result as a three-dimensional support strength matrix, containing the position index, support strength value, vertical pressure value and safety factor of each voxel. The core innovation of this step is to convert the density distribution information inside the material into a mechanical model, realizing the coupling 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 of insufficient support are effectively identified, laying a foundation for subsequent weak zone connectivity analysis.

[0074] Preferably, the step S3 of identifying the voxel weak zone of the support strength matrix according to the vertical pressure distribution field and performing weak zone connectivity analysis comprises:

[0075] The support strength matrix is compared with the vertical pressure distribution field by voxel, and the support strength and the bearing pressure are filtered to screen out the weak voxel position;

[0076] The weak point adjacent relationship graph is established for the weak point voxel position;

[0077] The weak point adjacent relationship graph is filtered and merged in a small area to obtain an effective connected region set;

[0078] The region geometric feature calculation is performed on the effective connected region set to obtain a region geometric parameter;

[0079] The overlying influence evaluation is performed according to the region geometric parameter and the vertical pressure distribution field to obtain weak zone connectivity data.

[0080] In the embodiment of the application, the weak zone connectivity analysis firstly compares the support strength matrix and the vertical pressure distribution field by voxel, and screens out the weak voxel. The safety margin SM(i,j,k)=S(i,j,k)-P(i,j,k) of each voxel is calculated, wherein S(i,j,k) is the support strength, and P(i,j,k) is the vertical pressure. When SM<0, the voxel cannot bear the pressure above, and is marked as a weak voxel. The three-dimensional index (i,j,k), the spatial coordinates (x,y,z), the safety margin negative value |SM| and the depth layer k of each weak voxel are recorded to construct a weak voxel position table. Then, the weak point adjacent relationship graph is established, and for each voxel in the weak voxel position table, 26 adjacent positions (up, down, left, right, front, back and diagonal directions) are checked to determine the adjacent weak points. The adjacent matrix is used to store the connection relationship, and the matrix element A(m,n)=1 indicates that the mth weak point is adjacent to the nth weak point, otherwise A(m,n)=0. The adjacent judgment standard is that the Euclidean distance between the center points of two voxels (maximum distance of diagonal adjacent). The adjacent graph G(V,E) describing the spatial distribution relationship of the weak points is formed by this method, wherein V is the weak point set, and E is the adjacent relationship set.

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

[0082] The set of valid connected regions is calculated for geometric features. Calculate the geometric center coordinates of each region: , the volume V=N×0.001m³ (single voxel volume is 0.001m³). 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 of the region, where Lmax and Lmin are the maximum and minimum values of L, respectively. The shape factor reflects the geometric properties of the region, and SF close to 1 indicates an equiaxed shape, and SF>3 indicates an elongated shape. Store each geometric parameter corresponding to the region ID to form a region geometric parameter table. Finally, perform an overburden influence assessment, and for each connected region, determine its top position ztop=max{zi|(i,j,k)∈region}. Calculate the overburden layer thickness h=zsurface-ztop, where zsurface is the height value of the surface profile at that position. Calculate the total pressure Ptotal=∑P(i,j,k) of the overburden layer, where (i,j,k) is all voxels from the top of the region to the surface. Evaluate the potential influence range: extend upwards to the surface at a 45° diffusion angle with the region center as the base point, and calculate the influence area . Integrate the overburden thickness , total pressure Ptotal and influence area with the region geometric parameters to form a complete weak zone connectivity data.

[0083] Preferably, the step S3 of performing a quantitative assessment of the collapse risk according to the weak zone connectivity data comprises:

[0084] Using the support strength matrix, the weak zone connectivity data is extracted to obtain the weak zone support parameters;

[0085] According to the weak zone connectivity data and the vertical pressure distribution field, the overburden pressure accumulation calculation is performed to obtain the weak zone pressure load;

[0086] According to the weak area support parameters and the weak area pressure load, the safety factor is dynamically determined to obtain a hierarchical safety factor;

[0087] According to the hierarchical safety factor, a standard risk index is calculated;

[0088] According to the standard risk index, a risk level is divided to obtain a hierarchical risk assessment table;

[0089] The hierarchical risk assessment table is mapped to a space risk to obtain a stress risk distribution map.

[0090] In the embodiment of the application, the collapse risk quantitative evaluation first extracts the peripheral support force of the weak area connection data. For each connected weak area, the position of its outer boundary voxel is determined. The boundary voxel is defined as: a voxel belonging to the weak area inside and at least one adjacent voxel not belonging to the weak area. A 3D edge detection algorithm is used to identify the boundary voxel, and a boundary voxel set B={b1, b2,...,bn} is constructed. For each boundary voxel bi, search for all non-weak area voxels within a voxel distance (10 cm) around it in the support strength matrix, and these voxels constitute the support boundary of the weak area. 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 peripheral support strength Ssurr of the weak area. The number ID, volume V, center coordinate C and corresponding peripheral support strength Ssurr of each weak area are recorded in the weak area support parameter table. For the weak area with uneven peripheral support, the standard deviation of the support strength is additionally calculated When When Ssurr>0.3, the peripheral support strength Ssurr is reduced by 20% to reflect the additional risk of uneven support.

[0091] The overburden pressure accumulation calculation is based on the weak area top position and the vertical pressure distribution field. For each weak area, the top position voxel set T={t1, t2,...,tk} is determined, which contains all the weak area 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 above it up to the surface are extracted in the vertical pressure distribution field, where h is the number of overburden layer voxels. Considering the diffusion effect of pressure transmission, the voxel pressures within a 45° range above the diagonal are weighted and accumulated. The weight coefficient According to the horizontal distance Determination: wherein The cumulative pressure Pacc(ti) borne by each top voxel ti is calculated: Pacc(ti) =∑w(d)×P(ti,j). The total overburden pressure of each weak zone is calculated as Ptotal=∑Pacc(ti). The ID of each weak zone and the corresponding total overburden pressure value Ptotal are recorded in the weak zone pressure load table as key parameters for assessing the collapse risk.

[0092] The safety factor is dynamically determined according to the volume size and shape characteristics of the weak zone. First, the basic safety factor SF_base is determined: when the volume V < 0.1 m³, SF_base = 1.5; when 0.1 m³≤V < 0.5 m³, SF_base = 2.0; when V≥0.5 m³, SF_base = 2.5. Then, it is adjusted according to the shape factor SF: when SF>3 (elongated shape), the safety factor increases by 0.5; when SF>5 (extremely elongated), the safety factor increases by 1.0. Finally, the depth factor is considered: the safety factor of surface weak zones (distance from surface <1 m) increases by 0.3; the safety factor of deep weak zones (distance from surface >3 m) decreases by 0.2. The final safety factor SF = SF_base + SF adjustment value. The determined safety factor is stored in correspondence with the weak zone ID as a graded safety factor table for subsequent risk index calculation.

[0093] The standard risk index calculation uses a comprehensive risk assessment formula: RI=(V×Ptotal) / (Ssurr×SF), where V is the weak zone volume, Ptotal is the total overburden pressure, Ssurr is the surrounding support strength, and SF is the safety factor. This formula considers four key factors: weak zone size, overburden load, surrounding support capacity, and safety margin, achieving comprehensive quantification of collapse risk. For ease of comparison and grading, the calculation results are normalized to map the risk index to the standard interval 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 value is recorded in the standard risk index table, corresponding one-to-one with the weak zone ID.

[0094] Risk level classification sets clear grading standards based on the standard risk index. Risk index 0-30 is low risk, marked as level 1 (green); 30-60 is medium risk, marked as level 2 (yellow); 60-80 is high risk, marked as level 3 (orange); and 80-100 is extremely high risk, marked as level 4 (red). For weak zones marked as levels 3 and 4, the system automatically generates a warning signal. Each weak zone records its ID, spatial location, volume size, risk index value, and corresponding risk level, forming a graded risk assessment table. This table becomes the core basis for safety management decisions, and the higher the risk level, the more stringent preventive measures need to be taken. For level 4 risk areas, immediate support reinforcement or unloading measures are required.

[0095] The spatial risk mapping maps the risk level information in the hierarchical risk assessment table back to the original three-dimensional spatial grid. For all voxels covered by each weak zone, the corresponding risk level label is assigned. A transition zone with a width of 20 cm is set at the boundary of the weak zone, and the risk level decreases with distance: the original level is maintained when d<5 cm; it is reduced by 0.25 level when 5 cm≤d<10 cm; it is reduced by 0.5 level when 10 cm≤d<15 cm; it is reduced by 0.75 level when 15 cm≤d<20 cm. For normal areas that are not weak zones, they are labeled as level 0 (safe). A complete three-dimensional risk distribution data is generated, containing the coordinates and risk level information of each spatial position. The stress risk distribution map is visualized using volume rendering technology, and different risk levels are marked with different colors: level 0 is blue, level 1 is green, level 2 is yellow, level 3 is orange, and level 4 is red. The high-risk area is highlighted by adjusting the transparency, so that the operator can intuitively identify the dangerous position. The innovation of this step is to convert the complex multi-parameter risk assessment results into an intuitive three-dimensional risk distribution map, providing a quantitative warning tool for material storage precision monitoring, and realizing the precise assessment and visualization of potential collapse risk.

[0096] Preferably, step S4 comprises:

[0097] Based on the stress risk distribution map, the material pile space where the material pile is located is divided into different density level regions to obtain a density zoning map;

[0098] According to the density zoning map, material species identification and classification are performed to obtain material species labeling data;

[0099] According to the density zoning map and the fused surface contour, partition volume data are calculated;

[0100] The partition volume data and the bulk density coefficient mapping table are subjected to quality distribution correction calculation to obtain a quality distribution table;

[0101] The quality distribution table is subjected to storage quantity summary and calibration to obtain calibrated storage quantity data.

[0102] In the embodiment of the present application, the partition bulk density dynamic correction is first based on the stress risk distribution map to divide the density region. The density level data in the stress risk distribution map is extracted, and the space of the material pile is divided into different density regions. Threshold segmentation method is adopted, and the density partition standard is set: density levels 0-3 are low density region (corresponding to density <300 kg / m³); density levels 4-7 are medium density region (corresponding to density 300-500 kg / m³); and density levels 8-10 are high density region (corresponding to density >500 kg / m³). At the region boundary position, the gradient analysis algorithm is applied to determine the accurate boundary, and the density gradient vector g=▽ρ of adjacent voxels is calculated, and when ‖g‖>50 kg / m³ / dm, it is identified as a boundary point. In order to reduce the influence of noise, morphological processing is performed on the preliminary segmentation result, and open-close operation with a radius of 2 voxels is executed to eliminate isolated points and small regions. Finally, the spatial structure map containing region ID, boundary coordinates and density classification, i.e. the density partition map, is generated.

[0103] The material species identification and classification is based on the density partition map and the physical properties of the material. In combination with the material storage record, a density-material species control table is established: the concentrate (mainly crushed grains or concentrated material) density range is 550-650 kg / m³; the roughage (mainly straw, hay, etc.) density range is 350-450 kg / m³; and the mixed material density range is 450-550 kg / m³. For each density partition, the average density value is calculated, and the material species type is determined by the determination rule: when >520 kg / m³, it is determined as a concentrate region; when <430 kg / m³, it is determined as a roughage region; and when 430 kg / m³≤ ≤520 kg / m³, it is determined as a mixed material region. For the boundary region which is difficult to determine directly (within ±15 kg / m³ range), the spatial continuity rule is adopted: adjacent to the known concentrate region is classified as concentrate, and adjacent to the known roughage region is classified as roughage. After the classification is completed, a unique material species ID is assigned to each region: the concentrate is marked as 1, the roughage is marked as 2, and the mixed material is marked as 3, forming the material species marking data.

[0104] The partition volume calculation adopts the grid voxel method to accurately measure the volume of each density partition. The density partition map is superimposed on the fused surface contour to determine the surface boundary of each partition. Each partition is independently calculated for volume, and the region is divided into small voxel units with a side length of 5 cm. It is judged whether each voxel is located inside the region: when the voxel center point height is lower than the surface contour and located within the partition boundary, it is marked as an effective voxel. The number of effective voxels in each partition is accumulated , multiplied by the unit voxel volume (i.e. 5 cm×5 cm×5 cm), to obtain the accurate volume value For the partial voxels at the boundary, the intersection volume calculation method of hexahedron and surface triangular mesh is adopted to improve the volume calculation accuracy. The ID, spatial range, material type and volume value of each density partition are recorded to form the partition volume data.

[0105] The mass distribution correction calculation multiplies the partition volume by the corresponding bulk density coefficient. For each partition, a differentiated bulk density coefficient is applied according to its material type and compaction degree. The bulk density coefficient mapping table provides the normalized bulk density coefficients of different materials under different compaction conditions: the bulk density coefficient of fine material area ranges from 1.0 to 1.3; the bulk density coefficient of coarse material area ranges from 0.65 to 0.85; the bulk density coefficient of mixed material area ranges from 0.85 to 0.95. The same material at different depth positions uses a layered bulk density coefficient: the reference value is used for the surface layer (0-1m); the middle layer (1-3m) increases by 5-10%; the bottom layer (>3m) increases by 10-15%. The corrected mass calculation formula of each partition is: wherein is the partition volume, is the reference density of fine material (600kg / m³), is the bulk density coefficient of the material, is the depth correction coefficient. The mass distribution table containing partition ID, volume, material type, bulk density coefficient and corrected mass is generated.

[0106] The storage quantity summary and calibration accumulates the corrected mass of all partitions to calculate the total storage quantity. Let the total number of partitions be n, then the total storage quantity is The calculation result is verified by comparing the historical warehouse entry records and warehouse exit data. When the calculation deviation exceeds the allowed range ±3%, the calibration coefficient is introduced for correction: wherein is the recorded value. The calibrated storage quantity is At the same time, the distribution proportion of each type of material is calculated, and the corrected storage quantity data is finally generated, including the total storage quantity, the weight of each type of material, the distribution proportion and the volume data, and a time stamp is attached to provide the accurate inventory status at the current time. The innovation of this method lies in combining the internal structure information of multi-sensor fusion detection to realize the identification and differentiated measurement of different density areas, overcoming the limitations of the traditional single bulk density coefficient calculation method, significantly improving the measurement accuracy under mixed storage conditions, and reducing the measurement error from the original 15-20% to within 3%.

[0107] Preferably, the determination of the bulk density coefficient mapping table according to the material marking data in step S4 comprises:

[0108] extracting the material reference bulk density from the material marking data;

[0109] obtaining the material warehouse entry time, combining the material reference bulk density to calculate the stacking time factor, and obtaining the time compaction factor;

[0110] According to the material mark data and the time compaction factor, a comprehensive compaction coefficient is obtained through deep compaction correction;

[0111] The comprehensive compaction coefficient and the material reference bulk density are subjected to bulk density coefficient standardization to obtain a standardized bulk density coefficient.

[0112] In the embodiment of the present application, the bulk density coefficient mapping table determines that the material reference bulk density is first extracted from the material mark data. A standard material bulk density reference database is established to record the reference bulk density values of various materials under standard conditions (newly stored and uncompacted). The standard bulk density of concentrate (including grain powder, concentrated material, etc.) is set to 600 kg / m³; the standard bulk density of coarse material (including hay, straw, etc.) is set to 400 kg / m³; and the standard bulk density of mixed material is set to 500 kg / m³. For each region in the material mark data, the corresponding reference bulk density value is matched according to the material ID (concentrate = 1, coarse material = 2, mixed material = 3) . The system simultaneously records the physical characteristic parameters of different materials, including the bulkiness coefficient: the bulkiness coefficient of concentrate is set to 1.0 (reference value), the bulkiness coefficient of coarse material is set to 1.8 (loose structure, easy to compress), and the bulkiness coefficient of mixed material is set to 1.3 (moderate compressibility). These reference values are used for subsequent compaction effect calculation.

[0113] The time compaction factor calculation is based on the storage time of the material and the compaction law. The storage time record of the material in each region is extracted from the warehouse management system, and the time difference between the current monitoring time and the storage time is calculated (in days). The long-term stacking of the material leads to self-weight compaction, and a time compaction factor calculation formula is established: , wherein is the maximum compaction gain (concentrate , coarse material , mixed material ), is the compaction rate coefficient ( / day). This formula describes the law that the compaction effect approaches the limit value exponentially with time. When the stacking time exceeds 30 days, the compaction factor tends to be stable, the maximum compaction factor of concentrate is 1.3, the maximum compaction factor of coarse material is 1.5, and the maximum compaction factor of mixed material is 1.4. In particular, when the storage time record is missing, the medium compaction state is used by default, i.e. . The calculated time compaction factor is stored corresponding to the region ID.

[0114] The deep compaction correction considers the pressure gradient effect caused by the vertical stacking of materials. For each region in the material type marking data, a layer correction is made according to its position in the vertical direction. The vertical direction of the material pile is divided into three layers: the surface layer (0-1 m from the surface), the middle layer (1-3 m), and the bottom layer (> 3 m). Different depth layers are set with different depth correction coefficients : surface layer =1.0 (reference value); middle layer , where is the depth value (m); bottom layer , not more than 1.2. The correction coefficient increases linearly with depth, reflecting the compaction effect of the upper material gravity on the lower layer. The compression characteristics of different materials are processed differently: the calculated value is directly used for fine material area; the deep correction coefficient is enhanced to for coarse material area due to its loose structure and significant compaction effect; the correction coefficient for mixed material area is . Multiply the time compaction factor by the depth correction coefficient to get the comprehensive compaction coefficient .

[0115] The bulk density coefficient standardization converts the comprehensive compaction coefficient of each region to a standardized bulk density coefficient under the same reference. Taking the standard bulk density of fine material (600 kg / m³) as the reference value, the actual bulk density values of all materials are converted to relative coefficients. The actual bulk density value is calculated as: , where is the reference bulk density of the material, is the comprehensive compaction coefficient. The standardization bulk density coefficient calculation formula is: . For fine material area, ranges from 1.0 to 1.3; for coarse material area ranges from 0.65 to 0.85; for mixed material area ranges from 0.85 to 0.95. To eliminate the sudden change of bulk density coefficient between adjacent regions, a 10 cm wide transition zone is set at the boundary of the region, and the transition bulk density coefficient is calculated using distance weighted average in the transition zone. Finally, a complete three-dimensional bulk density coefficient mapping table is constructed, each record containing spatial coordinates (x, y, z) and corresponding standardized bulk density coefficient value. This bulk density coefficient mapping table realizes the fine expression of bulk density at different positions, different materials, and different compaction degrees in the material warehouse, providing key parameters for subsequent quality distribution correction calculation. The innovation of this method lies in combining the physical properties of materials, stacking time and spatial position in three dimensions, establishing a dynamic and adaptive bulk density coefficient calculation model, and effectively solving the metering error problem caused by traditional single bulk density coefficient.

[0116] Table 1: Stock type Stock type Stock type Stock type Stock type Stock type Stock type Newly received (<5 days) 1.00 1.03 1.08 1.13-1.20 Stock type Stock type 1.10 1.14 1.19 1.24-1.28 Stock type Stock type 1.20 1.24 1.28 1.30 Stock type Newly received (<5 days) 0.67 0.70 0.73 0.75-0.78 Stock type Stock type 0.73 0.76 0.79 0.81-0.83 Stock type Stock type 0.80 0.82 0.84 0.85 Stock type Newly received (<5 days) 0.85 0.87 0.89 0.91-0.93 Stock type Stock type 0.88 0.90 0.92 0.93-0.94 Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type Stock type 0.92 0.93 0.94 0.95

[0117] Please refer to Table 1, which is a standardized bulk density coefficient ( ) table. The following is an explanation of Table 1 (standardized bulk density coefficient ( ) table):

[0118] 1. Standardization reference: The standard bulk density of concentrate (600 kg / m³) is used as the unified reference, and the actual bulk density of all materials is converted through the standardization coefficient .

[0119] 2. Calculation method: ;

[0120] : Reference bulk density of material (concentrate 600 kg / m³, roughage 400 kg / m³, mixed material 500 kg / m³);

[0121] : Time compaction factor, which increases exponentially with the increase of stacking time and approaches the maximum value;

[0122] : Depth correction coefficient, which increases linearly with the increase of depth.

[0123] 3. Application explanation:

[0124] Actual mass calculation: ;

[0125] A 10 cm transition zone is set at the boundary of the area, and the distance weighted average smoothing bulk density coefficient change is used;

[0126] The compaction effect of concentrate is small, and the compaction effect of roughage is significant (up to 50%);

[0127] When stacked for a long time (> 30 days), the compaction effect tends to be stable, and the coefficient no longer increases significantly.

[0128] 4. Advantages: This dynamic bulk density coefficient system considers three dimensions of material characteristics, time variation and spatial position, realizes accurate measurement in material storage monitoring, and reduces the measurement error from 15-20% to less than 3%.

[0129] Therefore, from any point of view, the embodiments should be regarded as exemplary and non-limiting, the scope of the present application is defined by the appended claims rather than the above description, and therefore all changes falling within the meaning and scope of the equivalent elements of the application file are intended to be included in the present application.

[0130] The foregoing is considered as illustrative only of the principles of the application. Numerous modifications and changes will readily occur to those skilled in the art, and it is intended to embrace all such modifications and changes that fall within the scope of the application. Accordingly, the application is not to be restricted in scope to the specific embodiments disclosed herein but is to be accorded the full scope that the principles and novel features request appropriately granted.

Claims

1. A method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar, characterized in that: The following steps are involved: Step S1: Perform dual-modal synchronous scanning of the material pile using a laser-millimeter wave radar to obtain a spatiotemporally aligned point cloud, wherein the spatiotemporally aligned point cloud includes penetration echo data and surface measurement data; divide the material pile into depth layers using the penetration echo data, 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: Spatially align the layered density feature map with the fused surface contour, and calculate the support stress distribution of each spatial unit to obtain the vertical pressure distribution field and support strength matrix; identify the voxel weak areas of the support strength matrix based on the vertical pressure distribution field, and perform connectivity analysis of the weak areas to obtain weak area connectivity data; Quantitatively assess collapse risk based on weak zone connectivity data and obtain a stress risk distribution map; Step S4: identifying 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 coefficient mapping table is used to dynamically correct the partition bulk density of the material pile to obtain the corrected storage capacity data.

2. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: The dual-mode synchronous scanning of the material pile in step S1 includes: The dual radar response signals of the material pile are collected synchronously by using laser radar and millimeter wave radar; Performing coordinate system transformation on the dual radar response signals and calibrating the time stamps to obtain a calibrated response signal; The radar point cloud data of the calibrated response signal is fused to obtain the spatiotemporal aligned point cloud.

3. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: In step S1, the material pile is divided into depth layers, and the interlayer dielectric parameters of the depth layers are calculated, including: Extract the penetration echo data of the millimeter-wave radar from the time-space aligned point cloud, perform echo sequence separation, and form a multi-layer echo pulse sequence; When the time interval between two adjacent echoes in a multi-layer echo pulse sequence is greater than a preset echo time threshold, it is determined to be an independent reflection layer. The material pile is divided into multiple depth layers of varying thickness from the surface to the bottom according to the reflection layer, and the boundary coordinates of the depth layers are recorded. The interlayer dielectric parameters of the material pile are calculated based on the depth layer boundary coordinates.

4. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: In step S1, the density gradient anomaly detection based on the interlayer dielectric parameters includes: Calculate the dielectric change rate of the interlayer dielectric parameters to obtain the interlayer change rate matrix; The central area data of the material pile is selected from the inter-layer change rate matrix, and the normal change benchmark is established to obtain the stratification benchmark range; Identify the abnormal point locations of the inter-layer change rate matrix according to the layered benchmark range; Determine the void characteristics at the abnormal point location based on the interlayer dielectric parameters and obtain a void type classification table; According to the void type classification table, the abnormal area is aggregated to obtain the aggregated abnormal area; Generate density anomaly label data based on aggregated anomaly areas; The depth layer is assigned a density level value according to the density anomaly marker data to obtain a layered density feature map.

5. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: Step S2 includes: Extract the surface measurement data of the LiDAR from the spatiotemporally aligned point cloud, perform quality assessment on the LiDAR point cloud, and obtain the reliability distribution of the point cloud. Identify the data missing area of ​​the point cloud reliability distribution and obtain the missing area boundary set; According to the missing area boundary set, millimeter wave radar surface data is extracted from the layered density feature map to supplement the missing area and obtain a supplementary height dataset; Perform dual-source data fusion and splicing on the point cloud reliability distribution and supplementary height dataset to obtain a complete height point matrix; The complete height lattice is smoothed and reconstructed to obtain the fused surface contour.

6. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: In step S3, the layered density feature map is spatially registered with the fused surface contour, and the support stress distribution of each spatial unit is calculated, including: Perform spatial data registration on the layered density feature map and the fused surface contour to obtain registered spatial data; The material bin space is divided into a regular voxel grid according to the fused surface contour and the layered density feature map to form a density-labeled voxel grid; The vertical pressure transmission is calculated on the density-marked voxel grid to obtain the vertical pressure distribution field; The support strength of the density-marked voxel grid is evaluated according to the vertical pressure distribution field to obtain the support strength matrix.

7. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: In step S3, the voxel weak areas of the support strength matrix are identified based on the vertical pressure distribution field, and the connectivity analysis of the weak areas is performed, including: The support strength matrix and the vertical pressure distribution field are used to compare the support strength and pressure of the voxels, and the weak voxel locations are screened out; Establish a weakness adjacency graph for the weakness voxel locations; Perform small area filtering and merging on the weakness adjacency graph to obtain a set of valid connected areas; Calculate the regional geometric characteristics of the effective connected region set to obtain the regional geometric parameters; The overburden impact assessment is performed based on regional geometric parameters and vertical pressure distribution field to obtain weak zone connectivity data.

8. The method for accurately monitoring stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: The quantitative assessment of collapse risk based on the weak zone connectivity data in step S3 includes: The support strength matrix is ​​used to extract the peripheral support force of the weak area connectivity data and obtain the weak area support parameters; The overburden pressure accumulation calculation is performed 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 support parameters and pressure load of the weak area to obtain the graded safety factor; Calculate the standard risk index based on the graded safety factor; Divide the risk levels according to the standard risk index and obtain a graded risk assessment table; Perform spatial risk mapping on the hierarchical risk assessment table to obtain a stress risk distribution map.

9. The method for accurate monitoring of stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: Step S4 includes: Based on the stress risk distribution map, the material pile space is divided into different density level areas to obtain a density zoning map; Identify and classify materials based on the density zoning map to obtain material labeling data; Calculate partition volume data based on density partition map and fused surface contour; Perform mass distribution correction calculation on the partition volume data and bulk density coefficient mapping table to obtain a mass distribution table; The storage capacity of the mass distribution table is summarized and calibrated to obtain the corrected storage capacity data.

10. The method for accurate monitoring of stored materials based on the fusion of laser radar and millimeter wave radar according to claim 1 is characterized in that: Determining the bulk density coefficient mapping table according to the material type marking data in step S4 includes: Extracting material benchmark density from material tag data; Obtain the material storage time, calculate the stacking time factor based on the material benchmark density, and obtain the time compaction factor; According to the material type marking data and time compaction factor, the depth compaction correction is performed to obtain the comprehensive compaction coefficient; The bulk density coefficient is standardized by the comprehensive compaction coefficient and the material benchmark bulk density to obtain the standardized bulk density coefficient.

Images