An unmanned aerial vehicle indoor coal yard autonomous coal piling method

Abstract

The application discloses an unmanned aerial vehicle indoor coal yard autonomous coal yard method, and belongs to the field of intelligent fuel management of thermal power plants. The steps are as follows: a multiple discrimination robust extended Kalman filter algorithm fusing digital twin priori and signal characteristics is adopted to realize autonomous flight; point clouds are collected and double echo physical difference pre-filtering and motion distortion compensation based on Lie algebra tangent space cubic Hermite manifold interpolation are performed; a ground station performs adaptive dust filtering, fixed structure segmentation of multi-scale feature fusion, global surface smoothing reconstruction of coal pile point clouds based on a moving least square method and incremental block volume calculation based on flight trajectory driving; and the quality of the point clouds is evaluated and automatic rescan is performed. The application solves the problem of false range error caused by metal multipath interference, significantly improves the point cloud fidelity and segmentation purity under high dust working conditions, eliminates residual high-frequency noise and surface undulation, and meets the precision and efficiency requirements of industrial coal yarding.

Description

Technical Field

[0001] This invention belongs to the field of intelligent fuel management technology for thermal power plants, specifically relating to a method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs). Background Technology

[0002] Dynamic and accurate monitoring of coal inventory in coal yards of thermal power plants is a core aspect of fuel management throughout the entire process, directly impacting the plant's fuel cost control and production scheduling efficiency. Currently, fully enclosed coal yards are gradually replacing traditional open-air coal yards, but the unique conditions of enclosed spaces present a series of intractable technical challenges to measuring coal inventory.

[0003] The enclosed space completely blocks GPS signals, rendering traditional GPS-dependent drone navigation solutions utterly ineffective. Furthermore, the numerous steel columns and beams within the interior cause severe multipath effects and non-line-of-sight (NLOS) errors in UWB signals. Existing robust extended Kalman filters (AEKF) primarily rely on the chi-square test of the innovation sequence residuals or the received signal strength index (RSSI) for NLOS determination. However, when drones fly close to large metal structures during inspections, strong total internal reflection leads to extremely high RSSI and temporary stability of the residuals, resulting in fatal "false line-of-sight" misjudgments in traditional algorithms. This causes the robust fusion model to fail, ultimately leading to drone positioning drift or even crashes. The safety and accuracy issues of manual coal inventory are becoming increasingly prominent. Workers venturing into the high-dust, complex structures of coal yards not only face high work intensity and extremely low efficiency but also significant safety hazards such as coal pile collapses and dust inhalation. Measurement errors are also greatly influenced by subjective human factors such as operator experience and location selection, making them completely uncontrollable. While the track-mounted laser scanner installed on the top of the coal yard can achieve localized automated measurement, it has inherent drawbacks such as high installation and maintenance costs and insufficient equipment flexibility. It cannot adapt to the frequent shape changes of the coal pile caused by feeding and unloading, and it also has special requirements for the structure of the coal yard top, resulting in poor versatility.

[0004] Indoor coal yards often experience high dust concentrations and uneven spatiotemporal distribution. Typical coal dust particles with diameters of 10-100 μm exhibit strong Mie scattering and absorption of 905 nm wavelength laser light, resulting in a large amount of outlier noise and false points mixed into the collected point cloud data, severely compromising the integrity and accuracy of the point cloud. While the Low Intensity Dynamic Radius Outlier Removal (LIDROR) algorithm proposed by the academic community can alleviate dust interference to some extent, it relies on fixed neighborhood search parameters and cannot adapt to the drastic fluctuations in dust concentration caused by feeding and reclaiming operations in coal yards. It generally suffers from problems such as "over-filtering," which mistakenly deletes effective feature points at the edge of the coal pile, or "under-filtering," which leaves a large number of noise points. Traditional linear interpolation or spatiotemporal synchronization methods based on kinematic Taylor expansion and quadratic extrapolation cannot simultaneously satisfy the attitude boundary constraints at both ends of the time interval. Under the high-frequency attitude jitter of UAVs, point cloud "ghosting" and integral divergence distortion are prone to occur. While existing Delaunay partitioning-driven volume calculation methods (DTVC) have improved the accuracy of traditional voxel slicing methods, they are still essentially global partitioning with extremely high time complexity. Executing global DTVC under the limited memory of airborne edge computing platforms can easily lead to memory overflow and crashes. If a crude physical block division is adopted, the breakage of normal vectors and overlapping areas at the seams of sub-blocks will cause serious volume redundancy errors, making it difficult to meet the stringent standards of high accuracy for industrial coal inventory.

[0005] Currently, drones equipped with lidar are used for coal inventory in open-air coal yards. However, these solutions heavily rely on GPS positioning and lack specialized algorithms designed for indoor high-dust and multipath metal environments, making them unsuitable for direct application in enclosed coal yard scenarios. While some research has explored UWB-based indoor drone positioning technology, it largely focuses on achieving positioning functionality, failing to address NLOS interference suppression and false line-of-sight (FLOS) misjudgments, and lacking deep integration with 3D lidar coal inventory systems. Furthermore, it lacks an integrated solution for point cloud noise reduction, feature extraction, and volume calculation in dusty environments. Existing indoor coal inventory technologies generally employ generic point cloud processing algorithms, failing to consider the loose, irregular, and rough surface characteristics of coal piles and the dynamic changes in indoor dust environments. This results in poor-quality point clouds and low volume calculation accuracy, failing to meet the technical requirements for industrial-grade coal inventory. Summary of the Invention

[0006] To address the problems existing in the prior art, this invention proposes an autonomous coal inventory method for indoor coal yards using unmanned aerial vehicles (UAVs). The method is reasonably designed, overcomes the shortcomings of the prior art, and has good results.

[0007] A method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs) includes the following steps: Step 1: Load the digital twin model of the coal yard, calibrate the UWB base station, and automatically generate a collision-free scanning route covering the coal storage area based on the digital twin model; Step 2: During flight, a triple robust extended Kalman filter algorithm based on digital twin prior is used to tightly couple UWB and IMU data, calculate the high-precision pose of the UAV in real time, and control the UAV to fly autonomously along the scanned route. Step 3: Control the 3D lidar to acquire point cloud data of the coal pile, perform point cloud preprocessing, the preprocessing includes dust noise filtering based on dual-echo physical difference and motion distortion compensation based on Lie algebra tangent space interpolation, and then download the preprocessed point cloud data to the ground station. Step 4: After receiving the point cloud data at the ground station, perform adaptive dust filtering, fixed structure-coal pile point cloud segmentation based on multi-scale features and digital twin priors, and volume calculation using an incremental block division and boundary fusion method driven by flight trajectory to obtain the coal storage amount of the coal pile. Step 5: Evaluate the quality of the reconstructed point cloud. If it does not meet the preset standards, automatically generate a rescanning route and control the UAV to perform rescanning. Step 6: After the point cloud quality meets the standards, generate the coal inventory results and standardized reports, and complete the data storage and archiving.

[0008] Furthermore, step 1 specifically includes: Load a digital twin model of the coal yard area with a resolution of ≤5cm. The model includes the three-dimensional coordinates and geometric properties of the coal yard boundary, ground elevation benchmark, and all fixed steel structures. The global coordinates of four UWB ground base stations were calibrated using a total station; Configure the core parameters for the coal inventory task, including preset point cloud density, relative flight altitude, lidar scanning frequency, and coal pile density; Based on the digital twin model, a gridded reciprocating parallel coverage scanning track is automatically generated. The distance between adjacent lines is set to 70% of the effective scanning width of the lidar, with a 1m safety distance reserved. The flight path spacing is dynamically adjusted based on the point cloud density collected in real time by the lidar, with the spacing reduced in sparse areas and increased in dense areas.

[0009] Further, step 2 specifically includes: defining the first System state vector at time t ,in Let be the position vector of the UAV in the navigation coordinate system. Let V be the velocity vector of the UAV in the navigation coordinate system. Let be the attitude quaternion of the drone. The zero bias vector of the accelerometer. This is the zero bias vector of the gyroscope; Based on prior state estimation Perform triple NLOS interference discrimination: First level: Calculate the... New information vector of a UWB base station : ; in, For the first UWB base stations Distance observations at time [time] For the first The observation matrix corresponding to each base station; No. The new information covariance of each UWB base station for: ; in, Let be the prior state covariance matrix. For the first Initial observation noise covariance of each base station; Construct the first Residual discriminant factor for each base station for: ; Perform a chi-square distribution significance test: if If the confidence threshold is exceeded, the observation of the UWB base station is determined to be abnormal; Second stage: Extracting the first Actual received signal strength of each UWB base station Calculate its relationship with the theoretical maximum received strength. Ratio deviation: ; The third step: Extracting prior state estimates Positional components in As predicted coordinates for drones ,connect With the Coordinates of UWB base stations Generate spatial rays; utilize a memory-resident digital twin model and an octree acceleration structure to perform ray tracing and collision detection. Geometric occlusion probability of a UWB base station for: ; Based on the above triple discrimination results, an exponential dynamic observation noise covariance matrix is ​​constructed: ; in, For the first Time of the first The updated observation noise covariance matrix of each UWB base station For the first Time of the first The initial observation noise covariance matrix of each UWB base station , , These are the weighting coefficients for the residual discriminant factor, RSSI deviation, and geometric occlusion probability, respectively. Take the maximum value; Calculate Kalman gain : ; in, Let be the covariance matrix of the prior state estimate. This is a global observation matrix, composed of all UWB base stations. It is pieced together. Indicates transpose; Update system state vector From the updated optimal state vector Extracting positional components and attitude components This is the output high-precision pose. The global innovation vector (composed of all base stations) (assembled from pieces)

[0010] Further, step 3 specifically includes: using the Mie scattering characteristics of coal dust to remove suspended dust particles, with the following discrimination condition: ; in, This represents the physical distance difference between the two echoes. , These are the measured distances of the first and second echoes, respectively. The wavelength of the laser; , These are the reflection intensities of the first and second echoes, respectively. If the discrimination criteria are met, the current measuring point is determined to be a suspended dust point and is directly removed; For any frame of lidar scan time interval Obtain the attitude rotation matrix and angular velocity at the start and end times. and Time span Define the measurement point time. The corresponding normalization time is ; By using a logarithmic mapping, the final attitude of the UAV is mapped to the starting tangent space, and the relative rotation vector is calculated: ; in, It is a relative rotation vector; For Li Qun up to its Lie algebra Logarithmic mapping; It is the inverse of the starting attitude rotation matrix; The final attitude rotation matrix; based on the relative rotation vector. Perform cubic Hermite interpolation in the tangent space to obtain the interpolated tangent space rotation vector. ; Generate the precise absolute rotation matrix of the measurement point through exponential mapping: ; in, For a moment The original instantaneous rotation matrix of the UAV relative to the navigation coordinate system; Lie algebra To Liqun Exponential mapping; for The antisymmetric matrix; Will Converted to roll angle Pitch angle Yaw angle The attitude rotation matrix is ​​obtained by representing the Euler angles. ; Pose interpolation compensation is performed point-by-point for each measurement point within the frame to control the time synchronization deviation within 0.3ms; The rotation matrix is ​​corrected in real time based on the UWB positioning deviation: ; in, It is the identity matrix; This is the attitude error compensation matrix; Perform point cloud coordinate calculation to convert polar coordinates in the radar coordinate system to global coordinates in the navigation coordinate system: ; in, , , The global coordinates of the measurement point in the navigation coordinate system; This is the attitude rotation matrix from the body coordinate system to the navigation coordinate system. , , These are the roll angle, pitch angle, and yaw angle of the drone, respectively. The rotation matrix from the radar coordinate system to the body coordinate system is obtained through ground-based three-dimensional calibration. This is the translation vector from the origin of the radar coordinate system to the origin of the body coordinate system; , , These are the distance to the measuring point, the horizontal azimuth angle, and the vertical azimuth angle measured by the lidar, respectively. , , The coordinates of the UAV in the navigation coordinate system; The preprocessed point cloud data is sorted and packaged according to timestamps and then transmitted to the ground station in real time via a 5GHz wireless data transmission module, while the original data is backed up locally.

[0011] Furthermore, step 4 includes the following sub-steps: Step 4.1: Adaptive dust filtering is performed using the MEA-LIDROR algorithm, which integrates Mie scattering multi-echo physical difference with real-time concentration sensing feedback to achieve precise noise reduction in extreme dust environments; Step 4.2: Perform intelligent segmentation of fixed structure and coal pile point cloud by multi-scale feature fusion. Based on the ground elevation benchmark of the coal yard, perform initial elevation segmentation. Calculate the normal vector coherence index by eigenvalue decomposition of the covariance matrix in the multi-scale neighborhood to screen candidate points of fixed structure and coal pile. Use the region growing algorithm to perform spatial connectivity analysis to remove residual coal pile point cloud. Use the digital twin model to perform prior verification of the segmentation results and automatically adjust the segmentation threshold to obtain a pure set of original coal pile point clouds. Step 4.3: The moving least squares (MLS) method is used to perform global surface smoothing reconstruction on the original point cloud of the clean coal pile. By fitting the local polynomial surface and projecting the point cloud, residual high-frequency noise and surface undulations are eliminated to generate a smooth and continuous coal pile surface point cloud. Step 4.4: Incremental volume calculation is performed using the IP-DTVC method. An orthogonal sub-block sequence is dynamically generated along the UAV flight trajectory, and the overlapping area of ​​adjacent sub-blocks is retained. Incremental constrained Delaunay triangulation is performed in each sub-block. Dynamic spatial elevation fusion based on projection axis distance is performed in the overlapping area to eliminate boundary splicing errors. The volume of the geometry formed by each triangular facet and the ground reference plane is integrated, and the total volume of the coal pile is obtained by summing them up. The coal storage quantity is calculated by combining the coal pile density.

[0012] Further, step 4.1 specifically includes: Based on the real-time environmental dust concentration C(t) collected by the airborne dust sensor, the collected point cloud is divided into three layers according to the degree of dust interference. Combined with the vertical channel of the three-dimensional lidar, the point cloud of each layer and different channels is subjected to noise reduction processing separately. Construct a KD-Tree, search for K nearest neighbors for each measurement point in each layer, and calculate the average distance from the measurement point to the K nearest neighbors. ; Calculate the dynamic search radius: ; in, For the first Dynamic search radius for each measuring point; For the first Detection distance of each measuring point; This refers to the radar angular resolution. This refers to the multiplier factor; For the first The reflection intensity at each measuring point; This represents the maximum reflection intensity of the point cloud at this layer. Combining real-time concentration data from airborne dust sensors , construct the first Dynamic distance detection convergence threshold for each measuring point: ; in, For the first The filtering threshold for each measurement point; The global average distance; The global standard deviation; For the first The local density of each measuring point within this layer; This represents the average global density of the point cloud at this layer. , This is the concentration-sensitive control coefficient; Let be the global real-time dust concentration at time t; Eliminate Outlier noise points; for points at the edge of a coal pile, if the curvature of the point... Greater than the curvature threshold ,when This should also be retained; Anisotropic diffusion completion based on Gaussian process regression is performed on the denoised point cloud to fill the tiny holes caused by noise removal, while preserving the real texture of the coal pile surface.

[0013] Further, step 4.2 specifically includes: Based on the coal yard ground elevation benchmark By setting the elevation threshold range for fixed structures, most of the point clouds of fixed structures are initially separated. Within neighborhoods at three scales of 5cm, 10cm, and 20cm, the normal vectors and curvatures of the point cloud are calculated using eigenvalue decomposition of the covariance matrix, and the normal vector coherence index is calculated. ; in, For the first The coherence index of the normal vector of each measurement point ranges from 0 to 1; For the first Each measurement point corresponds to three eigenvalues ​​of the neighborhood variance matrix; If the first Each measuring point at the corresponding scale If the value is ≥ 0.8, it is considered a candidate point for a fixed structure on a plane or cylindrical surface. The measuring points were determined to be candidate points for the coal pile; A region growing algorithm is used, with smoothness as the growth criterion. Connectivity analysis is performed on the separated candidate points of fixed structure. A threshold for the size of the connected region is set. Connected regions smaller than the threshold are judged as residual coal pile point clouds and are removed. The segmented fixed-structure point cloud is compared with the preset fixed-structure coordinates in the digital twin model of the coal yard to calculate the degree of overlap. ,like It automatically adjusts the elevation threshold range and connected component size threshold to re-divide the data.

[0014] Furthermore, step 4.4 specifically includes: A dynamic orthogonal sub-block sequence is generated along the global flight tangent direction, with the sub-block physical width... Inverse adjustment based on local point cloud density: ; in, This is the actual width of the sub-block; The base width; This represents the global average point cloud density. This represents the local point cloud density. The number of vertices in the Delaunay triangulation operation for each sub-block is controlled to be 200-300, and adjacent sub-blocks retain a 5% overlap. ; For the edge sub-blocks of the coal pile, a constrained Delaunay triangulation is adopted, using the edge feature points as constraint points to force the feature lines to not cross, so that the triangulation result fits the shape of the coal pile edge; for the main sub-blocks of the coal pile, an unconstrained Delaunay fast triangulation algorithm is adopted, which only follows the geometric criterion of the empty circumcircle, in order to improve the processing speed. In overlapping areas Perform dynamic spatial elevation fusion based on projection axis distance: ; in, The merged elevation value; , The coordinates of two adjacent sub-blocks A and B along the flight path are respectively... Elevation value at the location; , These are overlapping regions. The distance from the inner test point to the boundaries of sub-blocks A and B; Integrating the volume of the pentahedron formed by each triangular facet and the coal yard ground reference plane, the volume corresponding to a single triangular facet is: ; in, For the first The volume corresponding to each triangular facet; For the first The area of ​​each triangular facet; For the first The average elevation of the three vertices of a triangular facet; Used as a benchmark for the ground elevation of the coal yard; For the depression area of ​​the coal pile, a layered integration method is adopted to increase the number of integration layers and improve the calculation accuracy of the depression area; Total volume of coal pile The sum of the volumes of all triangular facets, combined with the coal bulk density. Calculate the amount of coal in stock. .

[0015] Furthermore, step 5 specifically includes the following sub-steps: Step 5.1: Construct point cloud quality evaluation indicators from three dimensions: point cloud density, noise rate, and integrity. The formula for calculating point cloud density is: ; in, For point cloud density, the requirement is... points / cm²; This represents the number of point clouds; The area of ​​the scanned region; The formula for calculating noise rate is: ; in, For noise rate, the requirement is... ; This represents the number of noise points; The formula for calculating integrity is: ; in, For integrity, it is required ; For effective coverage area; The area of ​​the scanned region; Step 5.2: If the point cloud density, noise rate and integrity of a certain area do not meet the above requirements, it is determined that the point cloud quality of the area is substandard. The ground station automatically generates a rescanning route for the area and controls the UAV to return and rescan. Step 5.3: After the supplementary scanning is completed, repeat steps 2 to 5 until all areas meet the point cloud density requirement. Points / cm², Noise Rate Completeness The industrial-grade coal quality standard.

[0016] The beneficial technical effects of this invention are as follows: This invention employs a multi-fusion robust positioning technology based on digital twin priors, which eliminates the "false line-of-sight" positioning drift phenomenon from a physical geometric perspective, enabling centimeter-level high-reliability positioning and fully autonomous flight of UAVs in an indoor enclosed coal yard.

[0017] By combining Lie algebra-cut space cubic Hermite manifold interpolation with a multi-echo concentration prior adaptive denoising algorithm, the motion distortion of point clouds caused by high-frequency jitter of UAVs and noise interference in high-dust environments are significantly suppressed, greatly improving the accuracy and fidelity of point cloud acquisition.

[0018] This invention integrates multi-scale geometric features with digital twin prior information to achieve high-precision intelligent segmentation of coal pile point clouds and fixed structure point clouds, with a segmentation accuracy of over 99.5%, providing a clean and reliable data foundation for volume calculation.

[0019] This invention adopts a volume calculation architecture that combines streaming incremental block division with boundary elevation weighting, which breaks through the computing power and memory bottlenecks of edge computing platforms, eliminates block splicing errors, improves processing efficiency by more than 70% compared with traditional methods, and controls the relative error of volume measurement within ±0.2%.

[0020] This invention, combined with a dustproof design and obstacle avoidance system for drones, exhibits excellent environmental adaptability and can be directly applied to the automated transformation and upgrading of various enclosed coal yards. Attached Figure Description

[0021] Figure 1 This is a flowchart of an indoor coal yard autonomous coal inventory method using unmanned aerial vehicles (UAVs) according to the present invention.

[0022] Figure 2 This is a flowchart of the real-time high-precision pose calculation of the UAV in this invention.

[0023] Figure 3 This is a flowchart of the point cloud data preprocessing process in this invention.

[0024] Figure 4 This is a hardware structure diagram of the UAV coal storage system based on UWB and 3D laser in this invention.

[0025] Figure 5 This is a comparative experimental image of the DT-Tri-Fusion AEKF algorithm in an indoor coal yard.

[0026] Figure 6 This is a comparison image of point cloud distortion compensation based on Lie algebra manifold interpolation in this invention; Among them, (a) is the ghosting of cross-sectional point cloud caused by jitter without correction; (b) is the fidelity point cloud cross-section after compensation based on Lie algebra manifold interpolation.

[0027] Figure 7 This is a comparison chart of the filtering performance of MEA-LIDROR and traditional LIDROR under different dust concentration conditions in this invention.

[0028] Figure 8 This is a simulation comparison image before and after MLS smoothing reconstruction of coal pile point cloud in this invention; Among them, (a) is the simulated original coal pile point cloud superimposed with residual dust noise and sampled discrete fluctuations, and (b) is the smooth continuous coal pile surface reconstructed by the MLS algorithm of this invention. Detailed Implementation

[0029] The specific embodiments of the present invention will be further described below with reference to specific examples: A method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs), such as Figure 1 As shown, it includes the following steps: Step 1: Load the digital twin model of the coal yard, calibrate the UWB base station, and automatically generate a collision-free scanning route covering the coal storage area based on the digital twin model; This step aims to complete the preliminary preparations for the coal inventory task. Based on the digital twin model of the coal yard, it achieves autonomous flight path planning without blind spots or collisions, providing a fundamental guarantee for subsequent data collection and overcoming the shortcomings of traditional manual route planning, which is inefficient and prone to missing blind spots. Specifically, it includes: Step 1.1 After the ground station software is started, load the digital twin model of the coal storage area with a resolution of ≤5cm. The model includes the precise three-dimensional coordinates and geometric properties of the coal yard boundary, ground elevation benchmark, and all fixed steel structures (columns, beams, conveyor belts, and grids). Step 1.2: Use a total station to calibrate the global coordinates of the four UWB ground base stations, with a calibration error of ≤2cm. Establish a unified navigation coordinate system and verify the communication link between the UWB base stations and the airborne tags to ensure signal coverage of the entire coal yard. Step 1.3: Configure the core parameters for the coal inventory task; including preset point cloud density ≥ 5 points / cm², relative flight altitude 5~8m, default lidar scanning frequency 10Hz, and coal pile density measured on-site according to MT / T739-1997 standard; Step 1.4: Based on the digital twin boundary model, the ground station automatically generates a gridded reciprocating parallel coverage scanning track. The spacing between adjacent routes is set to 70% of the effective scanning width of the lidar to ensure uniform coverage of the entire coal yard without blind spots. At the same time, it automatically avoids all known fixed structures and reserves a 1m safety distance. A new adaptive route adjustment function is added to dynamically adjust the route spacing according to the point cloud density collected by the lidar in real time. The spacing is reduced in sparse point cloud areas and increased in dense point cloud areas. Step 1.5: Simultaneously send the flight mission and digital twin model to the airborne edge computing microprocessor unit, complete the communication link test between the UAV, all airborne sensors and the ground station, and confirm that each module is working properly.

[0030] Step 2: During flight, a triple robust extended Kalman filter algorithm that integrates digital twin geometric ray tracing, novelty residual chi-square test and RSSI intensity discrimination is used to tightly couple UWB and IMU data, calculate the high-precision pose of the UAV in real time, and control the UAV to fly autonomously along the scanned route. This step aims to solve the problem of "false line-of-sight" misjudgment caused by metallic multipath interference in indoor GPS-denied environments. It employs the digital twin triple-fusion robust extended Kalman filter algorithm (DT-Tri-Fusion AEKF) designed in this invention to achieve centimeter-level high-reliability positioning, providing a precise pose reference for UAV autonomous flight and point cloud coordinate calculation. The algorithm uses a unified system state vector as its core computational carrier, and through a closed-loop process of "prior state prediction → triple NLOS interference discrimination → dynamic weight penalty → optimal state update," it eradicates the "false line-of-sight" defect of traditional algorithms from the physical geometric dimension. Figure 2 As shown, it specifically includes: Step 2.1: Define the system's core state vector; This algorithm encapsulates all physical quantities that need to be estimated in real time during the UAV's movement into a 16-dimensional system state vector, which serves as the basic data carrier for the entire process. All subsequent positioning calculations, flight control, and point cloud solutions directly call this vector or its components. ; in, For the first The system state vector at time t is the core object of the EKF algorithm's recursive estimation, and the final high-precision pose is directly extracted from this vector. This is the position vector of the UAV in the navigation coordinate system, in meters. It will be used as the starting point for geometric ray collision detection in step 2.2, the core data for pose output in step 2.4, and the global reference for point cloud coordinate calculation in step 3.5. This is the velocity vector of the UAV in the navigation coordinate system, in m / s. It is subsequently used to predict the position of the next moment in the state prediction stage, supporting the short-term calculation of the IMU when the UWB signal is lost. This is the attitude quaternion for the UAV, used to represent the rotational relationship between the body coordinate system and the navigation coordinate system. It subsequently serves as the core data for attitude output in step 2.4, the input for point cloud distortion compensation in step 3.4, and the rotational reference for point cloud coordinate calculation in step 3.5. This is the zero bias vector of the accelerometer, in m / s², which is subsequently used for online calibration of IMU measurement errors and is the key to achieving stable navigation without UWB signal in step 2.4; This is the zero bias vector of the gyroscope, in rad / s, which will be used later for online calibration of IMU attitude measurement error to avoid attitude angle divergence over time; The vector is initialized before the drone takes off: Initialize to the takeoff point coordinates calibrated by UWB. Initialize to [0,0,0] Initialize to the initial attitude after IMU calibration. , Initialize to IMU offline calibration values.

[0031] Step 2.2: Estimation based on prior state Perform triple NLOS interference discrimination; Based on the system state vector defined in step 2.1, a progressive NLOS discrimination logic is constructed from three dimensions: signal statistics, energy characteristics, and physical geometry. This logic accurately identifies UWB signal reflection interference caused by metal structures. All discrimination results are directly input into the dynamic observation noise covariance matrix calculation in section 2.3.

[0032] First step: Chi-square test of new information residuals; Based on the Prior state estimation at time 1 (from the previous moment) (Derived recursively from IMU data), the difference between UWB observations and predictions is calculated to reflect the degree of inconsistency between observations and predictions, i.e., the innovation vector. : ; in, For the first UWB base stations Distance observation at time (unit: m). For the prediction based on prior state Each base station ranging value, The expression is: ; in, The position vector (3D) of the UAV in the navigation coordinate system; For the first Global coordinates of a UWB base station (3D, total station calibration); UWB ranging Gaussian noise (mean 0, variance 0) ); The observation matrix is ​​used to extract... Positional components in The expression is: ; Only the first 3 elements are non-zero (corresponding to the partial derivatives of the position components), while the last 13 elements are 0 (velocity, attitude, and zero bias do not directly affect the ranging observation). No. The new information covariance of each UWB base station for: ; in, Let be the prior state covariance matrix. For the first Initial observation noise covariance of each base station; Construct the first Residual discriminant factor for each base station for: ; Perform a chi-square distribution significance test: if If the confidence threshold is exceeded, the observation of the UWB base station is determined to be abnormal.

[0033] Second step: RSSI energy ratio deviation test; Extracting the first from the hardware driver layer Actual received signal strength of each UWB base station Calculate its relationship with the theoretical maximum received strength. Ratio deviation: ; like If the value is too high, it indicates that the signal is attenuated by reflection.

[0034] The third layer: Digital twin geometric ray collision detection; Extracting prior state estimates Positional components in As predicted coordinates for drones ,connect With the Coordinates of UWB base stations Generate spatial rays; utilize a memory-resident high-precision digital twin model of the coal yard and an octree acceleration structure to perform ray tracing and collision detection. Geometric occlusion probability of a UWB base station For the first Whether the direct signal path between a UWB base station and the drone's current predicted location is physically blocked by fixed steel structures (columns, beams, conveyor belts, grid structures) in the coal yard is denoted as: ; This method directly determines whether a signal is blocked at the physical level, completely solving the defect of traditional algorithms where "false line-of-sight is caused by artificially high RSSI".

[0035] Step 2.3: Based on the above triple discrimination results, construct an exponential dynamic observation noise covariance matrix, apply weight penalties to the interfered UWB base stations, and achieve adaptive response to different levels of interference. ; in, For the first Time of the first The updated observation noise covariance matrix of each UWB base station corrects the relationship between the base station's observations and the system state vector. Update weights; For the first Time of the first The initial observation noise covariance matrix of each UWB base station , , These are the weighting coefficients for the residual discriminant factor, RSSI deviation, and geometric occlusion probability, respectively, with typical values ​​of 0.5, 0.3, and 1.2. The largest value reflects the highest priority of physical and geometric discrimination. When a base station is identified as NLOS (Normally Inaccessible Signal), NLOS occurs when the signal is blocked by a steel structure and reaches the drone after one or more reflections. Increase exponentially, this base station The update weight approaches 0, thus achieving automatic shielding of interference signals.

[0036] Step 2.4: EKF optimal state update and IMU short-term estimation; Based on the dynamic observation noise covariance matrix obtained in step 2.3, EKF observation update is performed to obtain the optimal system state estimate and output the high-precision pose of the UAV: Calculate Kalman gain : ; in, Let be the covariance matrix of the prior state estimate. The global observation matrix (composed of all UWB base stations) (assembled) Indicates transpose; Update system state vector , The global information vector is composed of all base stations. Concatenated from the updated optimal state vector Extracting positional components and attitude components As the output high-precision pose, the positioning frequency is 20Hz, the static error is ≤3cm, and the dynamic error is ≤5cm.

[0037] Step 2.5: Constant relative height control based on real-time pose; Call the system state vector output in step 2.4 height component By combining the real-time distance from the top of the coal pile measured by the downward-looking lidar, the edge computing microprocessor unit dynamically adjusts the drone's flight altitude, always maintaining a constant relative altitude of 5-8m, ensuring the optimal lidar scanning angle and avoiding point cloud blind spots.

[0038] Step 2.6: Adopt a combined strategy of "global obstacle avoidance + local obstacle avoidance". Global obstacle avoidance is based on a digital twin model to avoid all known fixed structures in advance; local obstacle avoidance uses a front-facing binocular camera to acquire a depth map, converts it into a local cost map, uses the VFH+ algorithm to calculate the collision-free direction in real time, generates a detour path, and returns to the original route after bypassing the obstacle, thus avoiding point cloud acquisition gaps.

[0039] To verify the positioning robustness and anti-interference performance of the DT-Tri-Fusion AEKF algorithm proposed in this invention under strong multipath interference, a flight comparison experiment was conducted. Traditional AEKF is affected by "false line-of-sight" interference in strong reflection zones, with Z-axis drift exceeding 1.5m. This invention eliminates misjudgments through triple NLOS discrimination, achieving a static positioning error ≤3cm and a dynamic error ≤5cm, enabling centimeter-level robust flight even in GPS-denied environments. Figure 5 This image shows a comparative indoor coal yard flight experiment of the DT-Tri-Fusion AEKF algorithm of this invention. The red curve represents the Z-axis positioning trajectory of the traditional AEKF algorithm, the blue curve represents the Z-axis positioning trajectory of the algorithm of this invention, and the green dashed line represents the theoretical flight altitude. The comparative experiment verifies the positioning robustness of this algorithm under strong metal multipath interference.

[0040] Step 3: Control the 3D lidar to acquire point cloud data of the coal pile, perform point cloud preprocessing, including dust and noise filtering based on dual-echo physical difference and motion distortion compensation based on Lie algebra tangent space interpolation, and then download the preprocessed point cloud data to the ground station. This step aims to achieve high-precision acquisition and real-time preprocessing of point cloud data, solving the problem of point cloud motion distortion caused by high-frequency attitude jitter of UAVs, while also removing most dust and noise in advance to reduce the computational burden on the ground station. Figure 3 As shown, it specifically includes: Step 3.1: The three-dimensional lidar performs a 360° omnidirectional scan at a set frequency. The flight control processing unit outputs a PWM synchronization pulse that is consistent with the radar scanning cycle through the auxiliary channel. The rising edge triggers the GPIO interrupt of the edge computing microprocessor unit. The interrupt response time is ≤0.1ms, realizing hardware-level time alignment between IMU sampling and laser echo, with a timestamp accuracy of 1μs.

[0041] Step 3.2: The edge computing microprocessor unit receives UDP data packets from the 3D LiDAR via the network port, parses them using the MSOP data packet format, and extracts the distance d, horizontal azimuth angle θ, vertical azimuth angle α, and reflection intensity I for each measurement point. At the same time, it extracts the first echo and second echo data in the dual echo mode.

[0042] Step 3.3: Perform airborne physical echo distance differential pre-filtering to remove suspended dust points using the Mie scattering characteristics of coal dust. The discrimination criteria are as follows: ; in, The physical distance difference between the two echoes is expressed in meters (m). , These are the measured distances of the first and second echoes, respectively. The wavelength of the laser is 905 nm. , These are the reflection intensities of the first and second echoes, respectively. If the discrimination criteria are met, the current measuring point is determined to be a suspended dust point and is directly removed in the hardware resolution layer, which can remove 30%-50% of the noise in advance; Step 3.4: Point cloud distortion compensation is performed using a method based on Lie algebra-cutting space cubic Hermite manifold interpolation, abandoning the Taylor expansion extrapolation method which easily leads to integral divergence, and ensuring the pose and velocity boundary constraints at both ends of the time interval. First, for any given frame of LiDAR scanning time interval Obtain the attitude rotation matrix and angular velocity at the start and end times. and Time span Define the measurement point time. Corresponding normalized time : ; in, This is the normalized time, with a value ranging from 0 to 1; By using a logarithmic mapping, the final attitude of the UAV is mapped to the starting tangent space, and the relative rotation vector is calculated: ; in, This is a relative rotation vector, in rad. For Li Qun up to its Lie algebra Logarithmic mapping; It is the inverse of the starting attitude rotation matrix; The final attitude rotation matrix; Based on relative rotation vector Perform cubic Hermite interpolation of the tangent space to ensure the velocities at both ends. Continuity: ; in, For normalized time The corresponding tangent space rotation vector; These are the angular velocities at the start and end times, respectively, in rad / s; Perform an exponential mapping to map the tangent space vector back to the global manifold space; generate the exact absolute rotation matrix for the measurement points: ; in, For a moment The corresponding attitude rotation matrix; Lie algebra To Liqun Exponential mapping; for The antisymmetric matrix; To adapt to subsequent point cloud coordinate calculation requirements, Attitude rotation matrix converted to aeronautical ZYX regular Euler angle representation ,in For roll angle, For pitch angle, Yaw angle; set up Matrix elements are ( For line numbers, (For column number), the inverse Euler angle calculation formula is as follows: ; Pose interpolation compensation was performed point-by-point on 12 data blocks and 192 measurement points within the frame, controlling the time synchronization deviation to within 0.3ms, thus completely eliminating point cloud distortion caused by UAV motion. To verify the effectiveness of the point cloud distortion compensation method based on Lie algebra-cut space cubic Hermite manifold interpolation described in step 3.4 of this invention, a comparative experiment on airborne laser point cloud motion distortion compensation was conducted, and the results are as follows. Figure 6 As shown.

[0043] Figure 6 (a) shows that, in the uncorrected state, the attitude and velocity boundary constraints at both ends of the time interval cannot be satisfied simultaneously. High-frequency attitude jitter of the UAV causes the cross-sectional point cloud to deviate from the theoretical physical ridge line, resulting in severe spatial multi-layer ghosting and misalignment, causing distortion in subsequent elevation projection and volume integral. The effect after processing with the method of this invention is as follows: Figure 6 As shown in (b), this invention maps the endpoint pose to the starting point tangent space using a logarithmic mapping to calculate the relative rotation vector, and performs cubic Hermite interpolation in the tangent space to ensure velocity continuity, thus achieving constrained convergence of the algorithm. Experimental results show that this invention controls the time synchronization deviation to within 0.3ms, significantly suppresses spatial ghosting and tortuosity, and the compensated point cloud cross-section closely matches the theoretical ridge line, ensuring the macroscopic geometric fidelity of the point cloud data from the source.

[0044] Step 3.5: Add an attitude error compensation term to correct the rotation matrix in real time based on the UWB positioning deviation, further improving the coordinate calculation accuracy. ; in, It is the identity matrix; The attitude error compensation matrix is ​​calculated in real time from the UWB positioning deviation. Step 3.6: Perform point cloud coordinate calculation, converting the polar coordinates of the radar coordinate system to global coordinates in the navigation coordinate system: ; in, , , These are the global coordinates of the measurement point in the navigation coordinate system, in meters (m). This is the attitude rotation matrix from the body coordinate system to the navigation coordinate system. , , These are the roll angle, pitch angle, and yaw angle of the UAV, respectively, in rad. The rotation matrix from the radar coordinate system to the body coordinate system is obtained through ground-based three-dimensional calibration. This is the translation vector from the origin of the radar coordinate system to the origin of the body coordinate system; , , are respectively the target point distance, horizontal azimuth and vertical azimuth measured by the LiDAR; , , are the position coordinates of the UAV in the navigation coordinate system; Sort and package the preprocessed point cloud data by time stamp, transmit it to the ground station in real time through a 5GHz wireless data transmission module, and back up the original data locally at the same time to prevent data loss caused by transmission interruption.

[0045] Step 4: After receiving the point cloud data at the ground station, perform adaptive dust filtering, fixed structure-coal pile point cloud segmentation based on multi-scale features and digital twin priors, global surface smoothing reconstruction of coal pile point cloud based on moving least squares method, and adopt the volume calculation of incremental blocking and boundary fusion method driven by flight trajectory to obtain the coal storage volume of the coal pile; This step aims to perform in-depth processing on the point cloud data transmitted from the airborne end, complete dust noise reduction, fixed structure segmentation, global surface smoothing reconstruction of coal pile point cloud and volume calculation, and solve the problems of much point cloud noise in high-dust environments, difficult segmentation of semi-buried components, poor stability of volume calculation caused by residual high-frequency noise and surface fluctuation, and memory overflow in massive point cloud processing. It specifically includes: Step 4.1, implement the MEA-LIDROR adaptive dust filtering algorithm, which fuses Mie scattering multi-echo physical difference and real-time concentration perception feedback to achieve accurate noise reduction in extreme dust environments: First perform point cloud layering: according to the environmental dust concentration C(t) collected in real time by the airborne dust sensor, divide the collected point cloud into three layers according to elevation: high-concentration layer (C(t)≥150mg / m³, with strong dust interference, mostly distributed in the upper part of the coal yard), medium-concentration layer (50mg / m³<C(t)<150mg / m³, the main area of the coal pile), low-concentration layer (C(t)≤50mg / m³, with weak dust interference, mostly distributed in the lower part of the coal yard), combined with the vertical channel of the three-dimensional LiDAR (-15°~15°), perform noise reduction processing separately for the point clouds of different channels in each layer; Construct a KD-Tree, search K nearest neighbor points for each point in each layer, and calculate the average distance from this point to the K nearest neighbor points ; Calculate the dynamic search radius: ; wherein, is the -th measurement point's dynamic search radius; is the -th measurement point's detection distance; is the angular resolution of the radar; is the magnification coefficient; is the The reflection intensity at each measuring point; This represents the maximum reflection intensity of the point cloud at this layer. Combining real-time concentration data from airborne dust sensors , construct the first Dynamic distance detection convergence threshold for each measuring point: ; in, For the first The filtering threshold for each measurement point; The global average distance; The global standard deviation; For the first The local density of each measuring point within this layer; This represents the average global density of the point cloud at this layer. , This is the concentration-sensitive control coefficient; typical values ​​are 0.8 and 0.02, respectively. Let be the global real-time dust concentration at time t; Perform noise removal and feature protection, and remove noise. Outlier noise points; for coal pile edge points (locally low density, high curvature), an additional curvature judgment condition is added, if the curvature of the point... Greater than the curvature threshold ,even though Slightly larger ( (Time) is also retained; curvature The calculation method is as follows: based on the first The covariance matrix is ​​constructed using the K nearest neighbors of each measurement point. Eigenvalue decomposition is then performed on the covariance matrix to obtain three eigenvalues. Then the curvature of the measuring point The calculation formula is: curvature The value ranges from 0 to 1; a larger value indicates that the surface at that point is steeper or closer to the edge. (Curvature threshold) Typical values ​​range from 0.1 to 0.15; Finally, post-processing optimization is performed, using anisotropic diffusion completion based on Gaussian process regression to fill in the tiny holes caused by noise removal, while preserving the true texture of the coal pile surface. To verify the dynamic environment adaptability of the MEA-LIDROR adaptive dust filtering algorithm of this invention, comparative tests were conducted simulating time-varying dust conditions caused by coal yard loading and unloading operations. Experiments show that, as Figure 7The performance comparison chart of the two filtering algorithms under different dust concentration conditions shows that the traditional LIDROR algorithm has a noise removal rate of <60% in the high dust range and an edge retention rate of only 69% in the low dust range. The MEA-LIDROR algorithm of this invention relies on the dynamic threshold and curvature discrimination mechanism to achieve a noise removal rate of ≥96.5% and an edge feature retention rate of ≥98.5% throughout the entire time period, balancing noise reduction accuracy and feature fidelity.

[0046] Step 4.2: Perform intelligent segmentation of fixed structure and coal pile point clouds using multi-scale feature fusion to solve the problem of difficult segmentation of semi-buried columns and complex supports; perform initial elevation segmentation based on the coal yard ground elevation benchmark; calculate the normal vector coherence index through covariance matrix eigenvalue decomposition within the multi-scale neighborhood to screen candidate points for fixed structures and coal piles; use a region growing algorithm for spatial connectivity analysis to remove residual coal pile point clouds; and use a digital twin model to perform prior verification of the segmentation results and automatically adjust the segmentation threshold; specifically including: First, initial elevation segmentation is performed: based on the coal yard ground elevation datum. Set the elevation threshold range for the fixed structure ( or , (Based on the maximum theoretical height of the coal pile), most of the fixed structure point cloud was initially separated; Multi-scale feature extraction is performed: In neighborhoods at three scales (5cm, 10cm, and 20cm), the normal vectors and curvatures of the point cloud are calculated using eigenvalue decomposition of the covariance matrix, and the normal vector coherence index is calculated. ; in, For the first The coherence index of the normal vector of each measurement point ranges from 0 to 1; For the first Each measurement point corresponds to three eigenvalues ​​of the neighborhood variance matrix; Then perform normal vector coherence screening: if the first... Each measuring point at the corresponding scale If it is, then it is determined to be a candidate point for a fixed structure on a plane or cylindrical surface. The points were determined as candidate points for the coal pile; Next, spatial connectivity screening is performed: a region growing algorithm is used, with smoothness as the growth criterion, to perform connectivity analysis on the separated fixed structure candidate points, and a threshold for the size of the connected domain is set. Connected domains smaller than the threshold are judged as residual coal pile point clouds and are removed. Finally, a priori verification of the digital twin is performed: the segmented fixed structure point cloud is compared with the preset fixed structure coordinates in the coal yard digital twin model, and the degree of overlap is calculated. For each fixed-structure point cloud obtained through segmentation, matching points with a Euclidean distance of less than 5 cm are searched in the preset fixed-structure point set of the digital twin model. The overlap ratio γ is defined as the ratio of the number of matching points to the total number of fixed-structure point clouds obtained through segmentation, i.e.: ; In the formula: For the number of matching points, This represents the total number of fixed-structure point clouds obtained through segmentation.

[0047] like It automatically adjusts the elevation threshold range and connected domain size threshold to re-segment, and automatically removes interference points that appear above the conveyor belt or outside the wall.

[0048] Step 4.3: After completing the segmentation in Step 4.2, obtain a clean set of original point clouds of the coal pile. To eliminate high-frequency noise residual in high-dust environments and minor surface undulations caused by LiDAR sampling, and to improve the realism of 3D visualization rendering and the stability of volume calculation, global surface smoothing reconstruction is performed: [Specifically targeting...] any point in Using KD-Tree to search for its smooth radius The local neighborhood point set is used. A local bivariate orthogonal polynomial surface is fitted using the moving least squares method with this neighborhood point set. Its core lies in minimizing the weighted sum of squared errors. : ; in, For the current target point to be smoothed, For the first in the neighborhood One sample point, For sample points The actual three-dimensional coordinates, The Euclidean distance between the target point and the sample point. The Gaussian weighting function decays with Euclidean distance. Local fitting surface function This represents the number of neighborhood points. After completing the local surface fitting, the points... Projected along its normal vector direction onto the fitted surface Above, obtain the smoothed new coordinates. By traversing all points, a point cloud reconstructing the coal pile with continuous, smooth surface features is finally generated. This smooth point cloud not only boasts extremely high visual rendering fidelity but also significantly reduces local abrupt errors during subsequent triangulation. Moving least squares (MLS) is employed to reconstruct the surface of the segmented pure coal pile point cloud. By fitting a locally continuous polynomial surface and projecting the discrete point cloud onto the fitted surface, the high-frequency micro-undulations and random noise of the original point cloud are effectively eliminated. To verify the algorithm's performance, simulation experiments were conducted based on real coal yard 3D scanning benchmark data: 0-5cm high-frequency sampling noise and random surface undulations were superimposed on a standard coal pile surface model to simulate the characteristics of the original point cloud collected by lidar in a high-dust environment. The algorithm generates a continuous, smooth, and high-fidelity coal pile surface model while fully preserving the overall outline, slope, depressions, and protrusions of the coal pile, providing data support for subsequent high-precision volume calculations. Figure 8 (a) is a simulated original coal pile point cloud with superimposed dust residual noise and sampling discrete fluctuations. Directly using it for volume calculation will produce a random error of more than ±0.4%. Figure 8 (b) is the smooth and continuous coal pile surface reconstructed by the MLS algorithm of this invention. While fully preserving the macroscopic geometric contour and edge features of the coal pile, the algorithm effectively removes residual high-frequency noise and eliminates small surface undulations. The surface continuity of the reconstructed point cloud is improved by 91%, providing a high-fidelity data foundation for high-precision coal inventory.

[0049] Step 4.4: Receive the smoothed reconstructed point cloud output from Step 4.3 The IP-DTVC method is used for incremental volume calculation to solve the problems of global Delaunay triangulation memory overflow and block stitching boundary error. An orthogonal sub-block sequence is dynamically generated along the UAV flight trajectory, preserving overlapping areas between adjacent sub-blocks. Incremental constrained Delaunay triangulation is performed within each sub-block, and dynamic spatial elevation fusion based on projection axis distance is performed within overlapping areas to eliminate boundary stitching errors. Volume integration is performed on the geometry formed by each triangular facet and the ground reference plane, and the total volume of the coal pile is accumulated. The coal storage quantity is then calculated in conjunction with the coal pile density. Specifically, this includes: First, trajectory-driven dynamic segmentation is performed. The system subscribes to the UAV's high-frequency pose and generates a sequence of dynamically orthogonal sub-blocks along the global flight tangent direction, with the sub-blocks having a physical width of [missing information]. Inverse adjustment based on local point cloud density: ; in, This represents the actual width of the sub-block, in meters (m). The base width; This represents the global average point cloud density. This represents the local point cloud density. The number of vertices in the Delaunay triangulation operation for each sub-block is controlled to be 200-300, enabling "calculation on demand" and completely solving the memory overflow problem; adjacent sub-blocks retain a 5% overlap area. ; Incremental constrained Delaunay subdivision is performed. For the edge sub-blocks of the coal pile, constrained Delaunay triangulation is used, and edge feature points including the coal pile-ground boundary, abrupt slope changes, and local depressions / protrusions are used as constraint points. The feature lines are forced to be non-crossable so that the subdivision result fits the shape of the coal pile edge. For the main sub-blocks of the coal pile, an unconstrained Delaunay fast subdivision algorithm is used. It is suitable for pure coal pile interior regions with no internal boundaries and uniform and continuous point cloud distribution. It only follows the geometric criterion of the empty circumscribed circle, which improves the processing speed. Subsequently, boundary distance-weighted elevation fusion was performed to eliminate volume overlap and step effects at sub-block boundaries in overlapping areas. Perform dynamic spatial elevation fusion based on projection axis distance: ; in, The elevation values ​​are the merged values, in meters (m). , The coordinates of two adjacent sub-blocks A and B along the flight path are respectively... Elevation value at the location; , These are overlapping regions. The distance from the inner test point to the boundaries of sub-blocks A and B; Integrating the volume of the pentahedron formed by each triangular facet and the coal yard ground reference plane, the volume corresponding to a single triangular facet is: ; in, For the first The volume corresponding to each triangular facet; For the first The area of ​​each triangular facet; For the first The average elevation value of the three vertices of a triangular facet; Used as a benchmark for the ground elevation of the coal yard; For the depression area of ​​the coal pile, a layered integration method is adopted to increase the number of integration layers and improve the calculation accuracy of the depression area; Total volume of coal pile The sum of the volumes of all triangular facets, combined with the coal bulk density determined according to MT / T739-1997. Calculate the coal inventory: ; in, This refers to the amount of coal in stock, expressed in tons (t). This refers to the bulk density of coal, expressed in t / m³. This represents the total volume of the coal pile, expressed in m³.

[0050] Step 5: Evaluate the quality of the reconstructed point cloud. If it does not meet the preset standards, automatically generate a supplementary scanning route and control the UAV to perform supplementary scanning.

[0051] This step aims to quantitatively evaluate point cloud quality and automatically perform supplementary scanning for areas that do not meet the requirements, ensuring that the point cloud data for the entire coal yard meets industrial-grade coal inventory standards. Specifically, it includes: Step 5.1: Construct point cloud quality evaluation indicators from three dimensions: point cloud density, noise rate, and integrity. The formula for calculating point cloud density is: ; in, For point cloud density, the requirement is... points / cm²; This represents the number of point clouds; The area of ​​the scanned region; The formula for calculating noise rate is: ; in, For noise rate, the requirement is... ; This represents the number of noise points; The formula for calculating integrity is: ; in, For integrity, it is required ; For effective coverage area; The area of ​​the scanned region; Step 5.2: If the point cloud density, noise rate and integrity of a certain area do not meet the above requirements, it is determined that the point cloud quality of the area is substandard. The ground station automatically generates a rescanning route for the area and controls the UAV to return and rescan. Step 5.3: After the supplementary scanning is completed, repeat steps 2 to 5 until all areas meet the point cloud density requirement. Points / cm², Noise Rate Completeness The industrial-grade coal quality standard.

[0052] Step 6: After the point cloud quality meets the standards, generate the coal inventory results and standardized reports, and complete the data storage and archiving; This step aims to visualize and standardize the coal inventory results, providing data support for power plant fuel management.

[0053] Step 6.1: Generate a 3D point cloud model of the entire coal yard and a contour map of the coal pile, supporting viewing from any angle, generating profile maps, measuring distances, and statistically analyzing regional selections.

[0054] Step 6.2: Output the core data of the coal inventory, including the total amount of coal stored, the total volume of the coal pile, the measurement time, the relative error of the volume measurement (≤±0.2%), and the distribution of coal in each area.

[0055] Step 6.3: Automatically generate a coal yard inventory report conforming to the MT / T739-1997 standard, including point cloud screenshots, historical data comparisons, error analysis, etc., and supports exporting to PDF / Excel format.

[0056] Step 6.4: Store the coal inventory data in the local database to support historical data query and trend analysis, providing data support for power plant fuel scheduling, cost accounting and production planning.

[0057] To implement the above method, this embodiment provides an indoor coal yard UAV coal inventory system based on UWB positioning and 3D laser scanning. The system hardware consists of two main units: an airborne unit and a ground unit. Each module communicates through a standardized interface to complete real-time data acquisition, transmission, and processing, such as... Figure 4 As shown.

[0058] The airborne hardware unit uses an edge computing microprocessor unit as its core computing node, integrating a flight control processing unit, a 3D LiDAR, a UWB airborne tag, an IMU inertial array, an obstacle avoidance module, and a power supply module. The flight control processing unit communicates with the edge computing microprocessor unit via an airborne flight control communication protocol and includes an added attitude pre-compensation module to proactively correct UAV attitude fluctuations, providing a stable attitude reference for point cloud coordinate calculation. The IMU inertial array is integrated within the flight control processing unit, outputting acceleration, angular velocity, and attitude angle data in real time via a communication interface to ensure positioning and attitude detection accuracy. The edge computing microprocessor unit adopts an industrial-grade configuration and is equipped with a custom point cloud preprocessing module, enabling real-time point cloud caching, preliminary noise reduction, and format conversion, effectively reducing the data processing burden on the ground station.

[0059] The drone platform is adapted to the complex working conditions of indoor coal yards. The fuselage is dustproof and equipped with a three-axis vibration-damping gimbal-mounted fixed 3D LiDAR to reduce the impact of flight vibration on the accuracy of point cloud acquisition. The 3D LiDAR adopts an industrial-grade configuration, possessing 360° omnidirectional scanning capability, supporting dual-echo mode, and is suitable for high-dust industrial environments. It outputs ranging data via Ethernet port using the UDP protocol, meeting the high-precision point cloud acquisition requirements of coal piles.

[0060] The UWB positioning module employs a deployment architecture of four fixed base stations and one airborne tag. The four fixed base stations are symmetrically positioned on top of the steel structure columns at the four corners of the enclosed coal yard, with an installation height 1.5-2 meters higher than the maximum designed stacking height of the coal yard, and the base station antennas facing the central area of ​​the coal yard. The global coordinates of each base station are precisely calibrated using a total station, with a calibration error ≤2cm, and the coordinate data is synchronously stored in the digital twin model. The positioning module is equipped with a dedicated UWB positioning chip, employing a hybrid positioning system of bidirectional ranging and time difference of arrival, while integrating real-time RSSI signal acquisition and multipath energy analysis functions. It can identify the signal multipath effect and non-line-of-sight interference caused by the steel structure of the coal yard, effectively suppressing positioning deviations caused by metal reflections. The airborne UWB tag communicates and links with the airborne edge computing microprocessor unit through a serial port, and integrates IMU inertial navigation data for tightly coupled calculation, achieving centimeter-level high-precision positioning output in indoor environments without GPS.

[0061] The obstacle avoidance module adopts a combination structure of binocular camera and downward-looking single-line LiDAR, and is connected to the edge computing microprocessor unit through USB interface to realize obstacle detection and terrain following, avoiding interruption or distortion of point cloud acquisition caused by drone collision.

[0062] The ground hardware unit is based on an industrial computer and integrates wireless data transmission modules, display terminals and storage devices. It runs self-developed ground station software to realize the integrated functions of overall system control, intelligent point cloud processing, data analysis and report generation, and complete data interaction and task management with the airborne terminal.

[0063] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs), characterized in that, Includes the following steps: Step 1: Load the digital twin model of the coal yard, calibrate the UWB base station, and automatically generate a collision-free scanning route covering the coal storage area based on the digital twin model; Step 2: During flight, a triple robust extended Kalman filter algorithm based on digital twin prior is used to tightly couple UWB and IMU data, calculate the high-precision pose of the UAV in real time, and control the UAV to fly autonomously along the scanned route. Step 3: Control the 3D lidar to acquire point cloud data of the coal pile, perform point cloud preprocessing, the preprocessing includes dust noise filtering based on dual-echo physical difference and motion distortion compensation based on Lie algebra tangent space interpolation, and then download the preprocessed point cloud data to the ground station. Step 4: After receiving the point cloud data at the ground station, perform adaptive dust filtering, fixed structure-coal pile point cloud segmentation based on multi-scale features and digital twin priors, and volume calculation using an incremental block division and boundary fusion method driven by flight trajectory to obtain the coal storage amount of the coal pile. Step 5: Evaluate the quality of the reconstructed point cloud. If it does not meet the preset standards, automatically generate a rescanning route and control the UAV to perform rescanning. Step 6: After the point cloud quality meets the standards, generate the coal inventory results and standardized reports, and complete the data storage and archiving; Step 4 includes the following sub-steps: Step 4.1: Adaptive dust filtering is performed using the MEA-LIDROR algorithm, which integrates Mie scattering multi-echo physical difference with real-time concentration sensing feedback to achieve precise noise reduction in extreme dust environments; Step 4.2: Perform intelligent segmentation of fixed structure and coal pile point cloud by multi-scale feature fusion. Based on the ground elevation benchmark of the coal yard, perform initial elevation segmentation. Calculate the normal vector coherence index by eigenvalue decomposition of the covariance matrix in the multi-scale neighborhood to screen candidate points of fixed structure and coal pile. Use the region growing algorithm to perform spatial connectivity analysis to remove residual coal pile point cloud. Use the digital twin model to perform prior verification of the segmentation results and automatically adjust the segmentation threshold to obtain a pure set of original coal pile point clouds. Step 4.3: The moving least squares (MLS) method is used to perform global surface smoothing reconstruction on the original point cloud of the clean coal pile. By fitting the local polynomial surface and projecting the point cloud, residual high-frequency noise and surface undulations are eliminated to generate a smooth and continuous coal pile surface point cloud. Step 4.4: Incremental volume calculation is performed using the IP-DTVC method. An orthogonal sub-block sequence is dynamically generated along the UAV flight trajectory, and the overlapping area of ​​adjacent sub-blocks is retained. Incremental constrained Delaunay triangulation is performed in each sub-block. Dynamic spatial elevation fusion based on projection axis distance is performed in the overlapping area to eliminate boundary splicing errors. The volume of the geometry formed by each triangular facet and the ground reference plane is integrated, and the total volume of the coal pile is obtained by summing them up. The coal storage quantity is calculated by combining the coal pile density.

2. The method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, Step 1 specifically includes: Load a digital twin model of the coal yard area with a resolution of ≤5cm. The model includes the three-dimensional coordinates and geometric properties of the coal yard boundary, ground elevation benchmark, and all fixed steel structures. The global coordinates of four UWB ground base stations were calibrated using a total station; Configure the core parameters for the coal inventory task, including preset point cloud density, relative flight altitude, lidar scanning frequency, and coal pile density; Based on the digital twin model, a gridded reciprocating parallel coverage scanning track is automatically generated. The distance between adjacent lines is set to 70% of the effective scanning width of the lidar, with a 1m safety distance reserved. The flight path spacing is dynamically adjusted based on the point cloud density collected in real time by the lidar, with the spacing reduced in sparse areas and increased in dense areas.

3. The method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, Step 2 specifically includes: defining the first System state vector at time t ,in Let be the position vector of the UAV in the navigation coordinate system. Let V be the velocity vector of the UAV in the navigation coordinate system. Let be the attitude quaternion of the drone. The zero bias vector of the accelerometer. This is the zero bias vector of the gyroscope; Based on prior state estimation Perform triple NLOS interference discrimination: First level: Calculate the... New information vector of a UWB base station : ； in, For the first UWB base stations Distance observations at time [time] For the first The observation matrix corresponding to each base station; No. The new information covariance of each UWB base station for: ； in, Let be the prior state covariance matrix. For the first Initial observation noise covariance of each base station; Construct the first Residual discriminant factor for each base station for: ； Perform a chi-square distribution significance test: if If the confidence threshold is exceeded, the observation of the UWB base station is determined to be abnormal; Second stage: Extracting the first Actual received signal strength of each UWB base station Calculate its relationship with the theoretical maximum received strength. Ratio deviation: ； The third step: Extracting prior state estimates Positional components in As predicted coordinates for drones ,connect With the Coordinates of UWB base stations Generate spatial rays; utilize a memory-resident digital twin model and an octree acceleration structure to perform ray tracing and collision detection. Geometric occlusion probability of a UWB base station for: ； Based on the above triple discrimination results, an exponential dynamic observation noise covariance matrix is ​​constructed: ； in, For the first Time of the first The updated observation noise covariance matrix of each UWB base station For the first Time of the first The initial observation noise covariance matrix of each UWB base station , , These are the weighting coefficients for the residual discriminant factor, RSSI deviation, and geometric occlusion probability, respectively. Take the maximum value; Calculate Kalman gain : ； in, Let be the covariance matrix of the prior state estimate. This is a global observation matrix, composed of all UWB base stations. It is pieced together. Indicates transpose; Update system state vector , The global information vector is composed of all base stations. Concatenated from the updated optimal state vector Extracting positional components and attitude components As the output, it is a high-precision pose.

4. The method for autonomous coal inventory in an indoor coal yard using a drone according to claim 1, characterized in that, Step 3 specifically includes: using the Mie scattering characteristics of coal dust to remove suspended dust particles, with the following criteria: ； in, This represents the physical distance difference between the two echoes. , These are the measured distances of the first and second echoes, respectively. The wavelength of the laser; , These are the reflection intensities of the first and second echoes, respectively. If the discrimination criteria are met, the current measuring point is determined to be a suspended dust point and is directly removed; For any frame of lidar scan time interval Obtain the attitude rotation matrix and angular velocity at the start and end times. and Time span Define the measurement point time. The corresponding normalization time is By using a logarithmic mapping, the final attitude of the UAV is mapped to the starting tangent space, and the relative rotation vector is calculated: ； in, It is a relative rotation vector; For Li Qun up to its Lie algebra Logarithmic mapping; It is the inverse of the starting attitude rotation matrix; The final attitude rotation matrix; Based on relative rotation vector Perform cubic Hermite interpolation in the tangent space to obtain the interpolated tangent space rotation vector. ; Generate the precise absolute rotation matrix of the measurement point through exponential mapping: ； in, For a moment The original instantaneous rotation matrix of the UAV relative to the navigation coordinate system; Lie algebra To Liqun Exponential mapping; for The antisymmetric matrix; Will Converted to roll angle Pitch angle Yaw angle The attitude rotation matrix is ​​obtained by representing the Euler angles. ; Pose interpolation compensation is performed point-by-point for each measurement point within the frame to control the time synchronization deviation within 0.3ms; The rotation matrix is ​​corrected in real time based on the UWB positioning deviation: ； in, It is the identity matrix; This is the attitude error compensation matrix; Perform point cloud coordinate calculation to convert polar coordinates in the radar coordinate system to global coordinates in the navigation coordinate system: ； in, , , The global coordinates of the measurement point in the navigation coordinate system; This is the attitude rotation matrix from the body coordinate system to the navigation coordinate system. , , These are the roll angle, pitch angle, and yaw angle of the drone, respectively. The rotation matrix from the radar coordinate system to the body coordinate system is obtained through ground-based three-dimensional calibration. This is the translation vector from the origin of the radar coordinate system to the origin of the body coordinate system; , , These are the distance to the measuring point, the horizontal azimuth angle, and the vertical azimuth angle measured by the lidar, respectively. , , The coordinates of the UAV in the navigation coordinate system; The preprocessed point cloud data is sorted and packaged according to timestamps and then transmitted to the ground station in real time via a 5GHz wireless data transmission module, while the original data is backed up locally.

5. The method for autonomous coal inventory in an indoor coal yard using a drone according to claim 1, characterized in that, Step 4.1 specifically includes: Based on the real-time environmental dust concentration C(t) collected by the airborne dust sensor, the collected point cloud is divided into three layers according to the degree of dust interference. Combined with the vertical channel of the three-dimensional lidar, the point cloud of each layer and different channels is subjected to noise reduction processing separately. Construct a KD-Tree, search for K nearest neighbors for each measurement point in each layer, and calculate the average distance from the measurement point to the K nearest neighbors. ; Calculate the dynamic search radius: ； in, For the first Dynamic search radius for each measuring point; For the first Detection distance of each measuring point; This refers to the radar angular resolution. This refers to the multiplier factor; For the first The reflection intensity at each measuring point; This represents the maximum reflection intensity of the point cloud at this layer. Combining real-time concentration data from airborne dust sensors , construct the first Dynamic distance detection convergence threshold for each measuring point: ； in, For the first The filtering threshold for each measurement point; The global average distance; The global standard deviation; For the first The local density of each measuring point within this layer; This represents the average global density of the point cloud at this layer. , This is the concentration-sensitive control coefficient; Let be the global real-time dust concentration at time t; Eliminate Outlier noise points; for points at the edge of a coal pile, if the curvature of the point... Greater than the curvature threshold ,when This should also be retained; Anisotropic diffusion completion based on Gaussian process regression is performed on the denoised point cloud to fill the tiny holes caused by noise removal, while preserving the real texture of the coal pile surface.

6. The method for autonomous coal inventory in an indoor coal yard using a drone according to claim 1, characterized in that, Step 4.2 specifically includes: Based on the coal yard ground elevation benchmark By setting the elevation threshold range for fixed structures, most of the point clouds of fixed structures are initially separated. Within neighborhoods at three scales of 5cm, 10cm, and 20cm, the normal vectors and curvatures of the point cloud are calculated using eigenvalue decomposition of the covariance matrix, and the normal vector coherence index is calculated. ； in, For the first The coherence index of the normal vector of each measurement point ranges from 0 to 1; For the first Each measurement point corresponds to three eigenvalues ​​of the neighborhood variance matrix; If the first Each measuring point at the corresponding scale If it is, then it is determined to be a candidate point for a fixed structure on a plane or cylindrical surface. The measuring points were determined to be candidate points for the coal pile; A region growing algorithm is used, with smoothness as the growth criterion. Connectivity analysis is performed on the separated candidate points of fixed structure. A threshold for the size of the connected region is set. Connected regions smaller than the threshold are judged as residual coal pile point clouds and are removed. The segmented fixed-structure point cloud is compared with the preset fixed-structure coordinates in the digital twin model of the coal yard to calculate the degree of overlap. ,like It automatically adjusts the elevation threshold range and connected component size threshold to re-divide the data.

7. The method for autonomous coal inventory in an indoor coal yard using unmanned aerial vehicles (UAVs) according to claim 1, characterized in that, Step 4.4 specifically includes: A dynamic orthogonal sub-block sequence is generated along the global flight tangent direction, with the sub-block physical width... Inverse adjustment based on local point cloud density: ； in, This is the actual width of the sub-block; The base width; This represents the global average point cloud density. This represents the local point cloud density. The number of vertices in the Delaunay triangulation operation for each sub-block is controlled to be 200-300, and adjacent sub-blocks retain a 5% overlap. ; For the edge sub-blocks of the coal pile, a constrained Delaunay triangulation is adopted, using the edge feature points as constraint points to force the feature lines to not cross, so that the triangulation result fits the shape of the coal pile edge; for the main sub-blocks of the coal pile, an unconstrained Delaunay fast triangulation algorithm is adopted, which only follows the geometric criterion of the empty circumcircle, in order to improve the processing speed. In overlapping areas Perform dynamic spatial elevation fusion based on projection axis distance: ； in, The merged elevation value; , The coordinates of two adjacent sub-blocks A and B along the flight path are respectively... Elevation value at the location; , These are overlapping regions. The distance from the inner test point to the boundaries of sub-blocks A and B; Integrating the volume of the pentahedron formed by each triangular facet and the coal yard ground reference plane, the volume corresponding to a single triangular facet is: ； in, For the first The volume corresponding to each triangular facet; For the first The area of ​​each triangular facet; For the first The average elevation of the three vertices of a triangular facet; Used as a benchmark for the ground elevation of the coal yard; For the depression area of ​​the coal pile, a layered integration method is adopted to increase the number of integration layers and improve the calculation accuracy of the depression area; Total volume of coal pile The sum of the volumes of all triangular facets, combined with the coal bulk density. Calculate the amount of coal in stock. .

8. The method for autonomous coal inventory in an indoor coal yard using a drone according to claim 1, characterized in that, Step 5 specifically includes the following sub-steps: Step 5.1: Construct point cloud quality evaluation indicators from three dimensions: point cloud density, noise rate, and integrity. The formula for calculating point cloud density is: ； in, For point cloud density, the requirement is... points / cm²; This represents the number of point clouds; The area of ​​the scanned region; The formula for calculating noise rate is: ； in, For noise rate, the requirement is... ; This represents the number of noise points; The formula for calculating integrity is: ； in, For integrity, it is required ; For effective coverage area; The area of ​​the scanned region; Step 5.2: If the point cloud density, noise rate and integrity of a certain area do not meet the above requirements, it is determined that the point cloud quality of the area is substandard. The ground station automatically generates a rescanning route for the area and controls the UAV to return and rescan. Step 5.3: After the supplementary scanning is completed, repeat steps S2 to S4, and then proceed to step 5.1 for quality evaluation until all areas meet the point cloud density requirements. Points / cm², Noise Rate Completeness The industrial-grade coal quality standard.

Images