Non-repeating scanning laser radar based ultra-long time domain point cloud accumulation FOD detection method

Abstract

The application discloses a non-repeated scanning laser radar-based super-long time domain point cloud accumulation FOD detection method, relates to the technical field of laser radar data processing, and comprises the following steps: dynamically setting the space voxel grid side length according to the minimum physical size of the to-be-detected foreign matter and the radar detection distance, and ensuring that the foreign matter occupies at least one voxel; establishing the relationship between the voxel activation proportion and the super-long time domain point cloud accumulation time according to the space voxel grid side length; and reversely calculating the minimum accumulation time length according to the preset false alarm rate. The non-repeated scanning laser radar-based super-long time domain point cloud accumulation FOD detection method improves the accuracy and robustness of FOD identification, and improves the reliability and real-time performance of foreign matter detection in a large range of complex environments as a whole.

Description

Technical Field

[0001] This invention relates to the field of lidar data processing technology, and specifically to a method for detecting cumulative FOD in ultra-long time-domain point clouds based on non-repeating scanning lidar. Background Technology

[0002] With the development of lidar technology, its application in airport runways, industrial parks, and other large-scale security monitoring scenarios is becoming increasingly widespread. In these applications, lidar continuously scans to acquire three-dimensional point cloud data of the environment and identifies abnormal objects in the scene based on the point cloud information, thereby achieving automatic detection of foreign objects. Due to its long-range, high-precision, and all-weather operation characteristics, lidar has become one of the important technical means for current FOD detection systems.

[0003] In practical applications, some lidar systems employ non-repeating scanning for spatial detection, meaning the scanning trajectory is not completely repeated within different scanning cycles, thus gradually covering the monitored area over a longer period. This scanning method can acquire rich spatial point cloud data over a large area, but it also results in significant differences in the number of scans at the same spatial location at different times. To improve detection stability, existing technologies typically require accumulating point cloud data over a certain time period and establishing a background distribution model based on the accumulated point cloud, so that abnormal targets can be identified through statistical features during subsequent real-time detection. However, several technical problems remain to be solved in FOD identification based on non-repeating scanning lidar. First, in large-scale monitoring scenarios, the physical size of the object to be detected may be small, while the detection range of the lidar may reach hundreds of meters. How to reasonably set the voxel grid size for spatial division so that it can cover small-sized objects without significantly increasing computational complexity due to an excessive number of voxels is a key issue in current system design. Second, when using point cloud data to establish a background model, it is usually necessary to accumulate and statistically analyze the scan data over a certain time period. However, in non-repeating scanning mode, the sampling frequency of point clouds in different spatial regions may vary significantly. If the accumulation time is too short, the background statistical features may be unstable, easily leading to false alarms; while if the accumulation time is too long, it may reduce the system's adaptability to environmental changes, thus affecting the real-time performance of foreign object detection. Therefore, how to reasonably determine the point cloud accumulation time while ensuring a controllable false negative rate to achieve stable background modeling over a long period is also an important problem that existing technologies need to solve. Furthermore, in the real-time detection phase, common methods typically determine the presence of anomalies by comparing the currently observed point cloud with the background model. However, since the number of scans generated by non-repeating scanning lidar at different spatial locations is uneven, using a uniform fixed threshold for anomaly detection may lead to significant differences in detection results between different regions, thus affecting the overall detection reliability. Therefore, how to construct a reasonable anomaly detection mechanism based on the scanning statistical features of spatial voxels in the background state to improve the accuracy of anomaly identification is also a technical challenge faced by existing technologies. Furthermore, when confirming candidate anomalies, relying solely on single-point statistical features is often insufficient to effectively distinguish between real foreign objects and complex background structures, such as ground textures or environmental noise. Therefore, it is also necessary to combine local point cloud distribution characteristics to further analyze abnormal areas in order to improve the accuracy and stability of FOD identification. Summary of the Invention

[0004] The purpose of this invention is to provide an ultra-long time-domain point cloud cumulative FOD detection method based on non-repeating scanning lidar, thereby solving the problems existing in the prior art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: an ultra-long temporal domain point cloud accumulation FOD detection method based on non-repetitive scanning lidar, comprising: S1, dynamically setting the side length of the spatial voxel grid according to the minimum physical size of the foreign object to be detected and the radar detection range, ensuring that the foreign object occupies at least one voxel; S2, establishing the relationship between the voxel activation ratio and the ultra-long temporal domain point cloud accumulation time for the spatial voxel grid side length, and calculating the minimum accumulation time in reverse according to the preset false alarm rate; S3, if the current ultra-long temporal domain point cloud accumulation time reaches the minimum accumulation time, then using the mean vector of the point cloud in the background voxel as the geometric center, calculating the covariance matrix as the discreteness and squareness features, and storing the mean vector. S4. In the real-time detection phase, obtain the point cloud coordinates generated by the non-repeating scanning lidar, index the voxels to which the point cloud coordinates belong, and read the corresponding mean vector and covariance matrix from the background distribution model; S5. Calculate the Mahalanobis distance from the point cloud coordinates to the background distribution model, and determine whether to write it into the abnormal candidate point set based on the comparison between the Mahalanobis distance and the dynamic threshold; S6. If the Mahalanobis distance is greater than the dynamic threshold, write the point cloud coordinates into the abnormal candidate point set, otherwise discard the point cloud coordinates; S7. For each point cloud coordinate in the abnormal candidate point set, perform a nearest neighbor search to obtain the local neighborhood point cloud, and calculate the real-time mean vector and real-time covariance matrix to obtain the real-time distribution model.

[0006] Preferably, step S1 includes obtaining the minimum physical size of the foreign object and the current detection distance; extracting the geometric features of the minimum physical size of the foreign object to obtain the side length of the basic bounding box; retrieving the beam cross-section diameter from a preset beam divergence model using the current detection distance; if the side length of the basic bounding box is greater than the beam cross-section diameter, obtaining the dynamic voxel reference value by multiplying the side length of the basic bounding box and the beam cross-section diameter; performing three-dimensional segmentation in the spatial coordinate system based on the dynamic voxel reference value to obtain a spatial voxel grid; and performing spatial mapping of the radar echo data through the spatial voxel grid to determine that the foreign object occupies at least one voxel.

[0007] Preferably, step S2 includes defining a discrete sampling interval in a three-dimensional coordinate system using the side length of the spatial voxel grid, calculating the distribution density of the reflected point cloud within the discrete sampling interval to obtain the voxel activation ratio; mapping the voxel activation ratio in a preset probability distribution model to determine the point cloud persistence state under different reflection intensities, and establishing the relationship between the voxel activation ratio and the ultra-long time-domain point cloud accumulation time; retrieving the corresponding confidence interval in the relationship according to the preset false alarm rate, and extracting the time-domain feature vector corresponding to the confidence interval; obtaining the minimum accumulation time by multiplying the time-domain feature vector with the sampling frequency, thus realizing the reverse calculation of the minimum accumulation time based on the preset false alarm rate.

[0008] Preferably, step S3 includes: comparing the current ultra-long temporal point cloud accumulation time with a preset minimum accumulation time; if the current ultra-long temporal point cloud accumulation time is greater than or equal to the minimum accumulation time, extracting the three-dimensional coordinate set of all discrete sampling points within the voxel grid; performing an arithmetic mean operation on the three-dimensional coordinate set to obtain the mean vector representing the geometric center of the voxel grid; using the mean vector and the three-dimensional coordinate set to perform a sum of squared differences to determine the covariance matrix reflecting the spatial distribution state of the point cloud; extracting feature vectors describing the discreteness and squareness of the point cloud based on the eigenvalue decomposition results of the covariance matrix, and mapping the feature vectors to the mean vector to obtain a background distribution model; and persistently storing the geometric center, discreteness, and squareness features of the point cloud within the background voxels by writing the background distribution model into a preset storage space.

[0009] Preferably, step S4 includes acquiring the current point cloud coordinate set generated by the non-repeating scanning lidar, determining the voxel grid position to which the point cloud coordinate set belongs by dividing the coordinate value range; extracting the mean vector and covariance matrix corresponding to the voxel grid position from the preset background distribution model as initial distribution parameters; performing matrix operations on the mean vector and covariance matrix to obtain a deviation vector group between the point cloud coordinate set and the background distribution; calculating the Mahalanobis distance value based on the deviation vector group; if the Mahalanobis distance value is less than a preset similarity threshold, then performing clustering processing on the deviation vector group to determine abnormal point cloud subsets; executing a density estimation algorithm based on the abnormal point cloud subsets to generate updated background distribution parameters; and writing the updated background distribution parameters back into the background distribution model to realize dynamic reading and adjustment of the mean vector and covariance matrix of the voxels to which the point cloud coordinates belong during the real-time detection stage.

[0010] Preferably, step S5 includes acquiring the real-time point cloud coordinates collected by the non-repeating scanning lidar and mapping them to the corresponding voxel grid; extracting the mean vector and covariance matrix corresponding to the voxel grid from a preset background distribution model; using the mean vector and covariance matrix to perform a multi-dimensional spatial difference measurement on the real-time point cloud coordinates to obtain the Mahalanobis distance; retrieving the cumulative scanning frequency of the voxel grid in the historical background state to determine a dynamic threshold; if the Mahalanobis distance exceeds the dynamic threshold, then writing the corresponding real-time point cloud coordinates into the abnormal candidate point set, wherein the dynamic threshold is positively correlated with the number of times the voxel is scanned in the background state.

[0011] Preferably, step S6 includes acquiring the real-time point cloud coordinates collected by the non-repeating scanning lidar and mapping them to the corresponding voxel grid; retrieving the historical scanning frequency of the voxel grid within a preset period to determine the dynamic threshold at the current moment; retrieving the mean vector and covariance matrix that match the voxel grid from the background distribution model; performing multi-dimensional spatial operations on the real-time point cloud coordinates using the mean vector and covariance matrix to obtain the Mahalanobis distance; if the Mahalanobis distance is greater than the dynamic threshold, then writing the real-time point cloud coordinates into the abnormal candidate point set, otherwise discarding the real-time point cloud coordinates.

[0012] Preferably, step S7 includes extracting the coordinates of a single candidate point from the set of abnormal candidate points and performing a radius search in three-dimensional space with that point as the center to obtain a local neighborhood point cloud set; inputting the local neighborhood point cloud set into the spatial moment calculation module, performing an arithmetic mean operation using the three-dimensional coordinate components of each point in the local neighborhood point cloud set to determine the real-time mean vector; performing an outer product operation on each coordinate point in the local neighborhood point cloud set after subtracting the real-time mean vector and accumulating the results, dividing the accumulated outer product matrix by the total number of points in the local neighborhood point cloud set to obtain the real-time covariance matrix; and performing a multidimensional Gaussian distribution modeling on the real-time mean vector and the real-time covariance matrix to obtain the real-time distribution model.

[0013] Preferably, the method further includes step S8: calculating the KL divergence as an anomaly score based on the real-time mean vector, the real-time covariance matrix, and the mean vector and covariance matrix of the corresponding background distribution model. Specifically, this includes obtaining the background mean vector and background covariance matrix of the background distribution model from a preset storage space; substituting the background mean vector, background covariance matrix, real-time mean vector, and real-time covariance matrix into the relative entropy analytical formula of the multidimensional Gaussian distribution; performing trace operation and determinant logarithmic operation on the background covariance matrix and the real-time covariance matrix using the relative entropy analytical formula to obtain a preliminary divergence value; using the preliminary divergence value combined with the difference term between the background mean vector and the real-time mean vector to perform a quadratic form weighted operation to determine the final Kübkelebler divergence; and performing normalization mapping processing on the Kübkelebler divergence to obtain an anomaly score reflecting the degree of deviation of the local point cloud distribution.

[0014] Preferably, the method further includes S9: if the KL divergence is greater than a preset threshold, the region corresponding to the abnormal candidate point set is determined to be a foreign object; otherwise, it is determined to be background texture. Specifically, this includes acquiring raw point cloud data from a lidar sensor and extracting candidate point sets; performing multi-dimensional feature fusion of spatial coordinates and reflection intensity on the candidate point sets to obtain a local point cloud distribution model; comparing the probability density of the local point cloud distribution model with a preset background texture model to determine the real-time Kürbeklebühler divergence; if the real-time Kürbeklebühler divergence is greater than a preset divergence threshold, the local region where the candidate point set is located is determined to be a foreign object region; if the real-time Kürbeklebühler divergence is less than or equal to the divergence threshold, the local region where the candidate point set is located is determined to be background texture.

[0015] As can be seen from the above technical solution, the present invention has the following beneficial effects: This ultra-long temporal point cloud accumulation FOD detection method based on non-repetitive scanning lidar dynamically sets the side length of the spatial voxel grid according to the minimum physical size of the foreign object to be detected and the radar detection range, ensuring that the foreign object occupies at least one voxel in the spatial division. This ensures the detectability of small-sized foreign objects while avoiding the computational burden caused by an excessive number of voxels. Simultaneously, by establishing the relationship between the voxel activation ratio and the ultra-long temporal point cloud accumulation time, and inversely calculating the minimum accumulation time based on a preset false alarm rate, it achieves reasonable long-term temporal accumulation modeling of point cloud data under non-repetitive scanning conditions, thereby improving the stability and reliability of the background model. In practice… In the time detection phase, the Mahalanobis distance from the point cloud to the background model is calculated using the mean vector and covariance matrix in the background distribution model. Anomaly point screening is performed by combining a dynamic threshold related to the number of voxel scans, which enables stable anomaly detection even when the number of scans is uneven in different spatial regions. In addition, by performing local neighborhood statistical analysis on the anomaly candidate point set, a real-time distribution model is constructed and KL divergence is used to further distinguish anomaly regions, effectively differentiating real foreign objects from background textures or environmental noise. This improves the accuracy and robustness of FOD recognition and enhances the reliability and real-time performance of foreign object detection in large-scale complex environments. Attached Figure Description

[0016] Figure 1 This is a signal transmission diagram of the present invention; Figure 2 This is a structural block diagram of the local terminal of an exemplary electronic device of the present invention; Figure 3 This is a structural block diagram of the network terminal of an exemplary electronic device of the present invention. Detailed Implementation

[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] Example 1: As Figure 1As shown, this invention provides a technical solution: an ultra-long temporal point cloud accumulation FOD detection method based on non-repetitive scanning lidar, including: S1, dynamically setting the spatial voxel grid side length according to the minimum physical size of the foreign object to be detected and the radar detection distance, ensuring that the foreign object occupies at least one voxel; S2, establishing the relationship between the voxel activation ratio and the ultra-long temporal point cloud accumulation time for the spatial voxel grid side length, and calculating the minimum accumulation time inversely based on the preset false alarm rate; S3, if the current ultra-long temporal point cloud accumulation time reaches the minimum accumulation time, then using the mean vector of the point cloud in the background voxel as the geometric center, and calculating the covariance matrix as the discreteness and squareness features, storing the mean vector and covariance matrix as the background distribution model; S4, in the real-time detection stage, acquiring the point cloud coordinates generated by the non-repetitive scanning lidar, indexing the voxel to which the point cloud coordinates belong, and reading from the background distribution model... S5. Calculate the Mahalanobis distance from the point cloud coordinates to the background distribution model. Compare the Mahalanobis distance with a dynamic threshold to determine whether to write it into the abnormal candidate point set. The dynamic threshold is positively correlated with the number of times the voxel is scanned in the background state. S6. If the Mahalanobis distance is greater than the dynamic threshold, write the point cloud coordinates into the abnormal candidate point set; otherwise, discard the point cloud coordinates. S7. For each point cloud coordinate in the abnormal candidate point set, perform a nearest neighbor search to obtain the local neighborhood point cloud. Calculate the real-time mean vector and real-time covariance matrix to obtain the real-time distribution model. S8. Calculate the KL divergence as the abnormal score based on the real-time mean vector, real-time covariance matrix, and the mean vector and covariance matrix of the corresponding background distribution model. S9. If the KL divergence is greater than a preset threshold, determine that the region corresponding to the abnormal candidate point set is a foreign object; otherwise, determine it as background texture.

[0019] This embodiment utilizes an ultra-long temporal domain point cloud accumulation and statistical distribution modeling method to achieve automatic identification of foreign objects. During continuous scanning, the scanning trajectory of a non-repeating scanning lidar will not repeatedly cover the same spatial location within a short period; through long-term accumulation, a complete spatial sample can gradually be formed. Therefore, this application first performs voxelization processing on the space based on the minimum physical size of the foreign object to be detected and the radar detection range, and dynamically sets the voxel grid side length, ensuring that the target foreign object occupies at least one voxel unit after spatial discretization, thereby guaranteeing detection sensitivity.

[0020] Furthermore, after the spatial voxel grid is established, the relationship between the proportion of voxels activated by the point cloud and time is statistically analyzed to construct a functional relationship between the voxel activation ratio and the point cloud accumulation time. The minimum point cloud accumulation time required is then derived in reverse by combining the preset false alarm rate to ensure that the background space is fully sampled.

[0021] Specifically, when the accumulated time reaches the minimum accumulated duration, statistical analysis is performed on the point cloud within each background voxel. The mean vector of the point cloud is calculated as the spatial geometric center, and the covariance matrix is ​​calculated to describe the dispersion and spatial shape characteristics of the point cloud distribution, thereby constructing a background distribution model.

[0022] Furthermore, during the real-time detection phase, the system continuously receives point cloud data output from the non-repeating scan lidar and maps the point cloud coordinates to the corresponding voxels. Subsequently, the mean vector and covariance matrix of the corresponding voxels are read from the background distribution model, and the difference between the current point cloud and the background distribution is measured based on Mahalanobis distance. Since different voxels have different background sampling times, this application introduces a dynamic threshold mechanism related to the number of scans to adapt to different statistical confidence conditions. When the Mahalanobis distance exceeds the corresponding threshold, the point cloud is written into the anomaly candidate point set.

[0023] Furthermore, a local neighborhood search is performed on the point cloud within the anomaly candidate point set. A real-time distribution model is constructed by statistically analyzing the real-time mean vector and real-time covariance matrix of the neighborhood point cloud. This real-time distribution model is then compared with the background distribution model, and the degree of distribution difference is measured by calculating the Kullback-Leibler divergence between the two. When the KL divergence exceeds a preset threshold, it indicates that the point cloud distribution in that region significantly deviates from the background distribution, thus identifying it as an alien target; otherwise, it is identified as background texture or normal environmental structure.

[0024] Example 2: S1 includes obtaining the minimum physical size of the foreign object and the current detection distance; extracting the geometric features of the minimum physical size of the foreign object to obtain the side length of the basic bounding box; using the current detection distance to retrieve the beam cross-section diameter from the preset beam divergence model; if the side length of the basic bounding box is greater than the beam cross-section diameter, then the dynamic voxel reference value is obtained by multiplying the side length of the basic bounding box and the beam cross-section diameter; the spatial voxel grid is obtained by three-dimensional subdivision in the spatial coordinate system based on the dynamic voxel reference value; and the radar echo data is spatially mapped through the spatial voxel grid to determine that the foreign object occupies at least one voxel.

[0025] In this embodiment, the voxel grid side length is set with the minimum physical size of the object to be detected and the current detection distance as input. First, the minimum spatial scale that the object must cover when spatially discretized is obtained. Then, this scale is made consistent with the spatial coverage scale of the laser beam at this distance, thereby determining the voxel side length value used for three-dimensional segmentation.

[0026] Specifically, the minimum physical size of the foreign object is determined during the system deployment phase based on the security requirements of the application scenario. The determination process includes: sorting out the types of foreign objects that must be identified in the scenario and their minimum external dimensions, selecting the minimum size that has the greatest impact on safety as the detection lower limit, and recording the lower limit as configuration parameters in the form of numerical values ​​in three directions: length, width, and height. When a unified lower limit is adopted for the scenario, the length, width, and height are taken as the same value to ensure that the minimum foreign object in any orientation meets the detection boundary conditions.

[0027] Furthermore, when obtaining the current detection range, the correspondence between the radar coordinate system and the scene coordinate system is first established. The origin of the radar coordinate system is taken as the radar ranging starting point, and the directions of the three coordinate axes are determined by the installation calibration. The calibration process includes: measuring the radar installation position, installation attitude, and relative relationship with the scene reference point, and saving this relationship as coordinate transformation parameters.

[0028] Specifically, during real-time operation, the current detection distance is taken as the representative distance of the radar to the target monitoring area. The process of determining the representative distance is as follows: after the boundary of the monitoring area is defined, the spatial distance from the radar origin to the farthest boundary point of the area is calculated as the current detection distance. The farthest boundary point is determined by the maximum coordinate range of the monitoring area. The distance is used to reflect the maximum expansion state of the radar beam in the monitoring area at this moment. When the system manages the monitoring area by partition, the distance of the farthest boundary point of each partition is taken, and the voxel side length value is set for each partition.

[0029] Furthermore, when extracting the basic bounding box side length by geometric feature extraction of the minimum physical size of the foreign object, the minimum physical size is first organized into three orthogonal direction size inputs, the input source being the configuration parameters recorded during the deployment phase.

[0030] Subsequently, a cuboid enclosure structure is constructed that can completely enclose the smallest foreign object. The three sides of this enclosure structure are aligned with three orthogonal directions. To ensure that the foreign object still occupies at least one voxel in any orientation, the side length of the basic bounding box is taken as the maximum value among the three sides of the cuboid enclosure structure. The selection process for this maximum value is as follows: the length, width, and height are compared sequentially, and the largest value is retained as the side length of the basic bounding box. The origin direction is then recorded after this side length.

[0031] When retrieving the beam cross-section diameter from the preset beam divergence model using the current detection distance, the parameter source and generation process of the beam divergence model are first established.

[0032] Specifically, the beam divergence model is determined during the equipment integration phase by combining radar factory parameters and field calibration: First, the radar's beam divergence characteristic parameters and exit port characteristic parameters are read; then, reflective targets are set at multiple known distances to collect the coverage area of ​​the echo beam and count the coverage diameter value of the beam in the lateral direction; multiple distances and their corresponding diameters are compiled into a distance-to-diameter lookup table and stored in ascending order of distance. The retrieval process during operation is as follows: Locate the two distance nodes closest to the current detection distance in the lookup table, take the diameter values ​​corresponding to these two nodes, and perform linear interpolation on the diameter values ​​according to the distance difference ratio to obtain the beam cross-section diameter corresponding to the current detection distance; before interpolation, the distance unit and diameter unit are standardized to use the same length unit to avoid distortion of voxel side lengths caused by numerical scale deviations.

[0033] Specifically, when the side length of the basic bounding box is greater than the beam cross-section diameter, the dynamic voxel reference value is calculated. The "greater than" determination here is a numerical comparison determination. The process is as follows: after unifying the side length of the basic bounding box and the beam cross-section diameter to the same length unit, the numerical values ​​of the two are directly compared. The comparison uses the same numerical precision, which is determined by the system quantization resolution. The quantization resolution is set during the deployment phase according to the minimum calibration precision of the scene coordinate system. If the side length of the basic bounding box is greater than the beam cross-section diameter, it indicates that the size of the object exceeds the lateral coverage scale of a single beam at that distance. In this case, the side length of the basic bounding box is multiplied by the beam cross-section diameter to obtain the dynamic voxel reference value. Unit consistency is achieved before multiplication, and the result is rounded after multiplication. The rounding direction is upward to ensure that the voxel side length is not lower than the spatial coverage requirement represented by the calculated result. The rounding step is the minimum grid step set during the deployment phase. This minimum grid step is jointly determined by the scene coordinate system calibration precision and the system calculation resolution and recorded as a configuration parameter.

[0034] Through this multiplication operation and rounding up, the voxel side length value changes together with the target size and beam spread, thus matching the long-range resolution capability.

[0035] Furthermore, when the side length of the basic bounding box is not greater than the beam cross section diameter, the dynamic voxel reference value uses the value of the beam cross section diameter as the basis for the voxel side length value. The process is as follows: after unifying the beam cross section diameter to the length unit used by the voxel grid, directly take the value and round it up according to the smallest grid step.

[0036] Specifically, this processing ensures that the voxel side length is not less than the beam coverage scale, avoiding the discontinuous distribution of single sampling within the voxel due to excessively small voxel side lengths, which would affect the stability of subsequent point cloud spatial mapping. At the same time, the voxel side length still satisfies the minimum external object size constraint, because when the basic bounding box side length is not greater than the beam cross-section diameter, the observable scale of the foreign object at that distance does not exceed the beam coverage scale. Taking the beam cross-section diameter as the voxel side length can ensure that the spatial unit where the foreign object is located has effective sampling support.

[0037] As a preferred embodiment, when obtaining a spatial voxel mesh by three-dimensional segmentation in a spatial coordinate system based on dynamic voxel reference values, the three-dimensional boundary range of the monitoring space is first defined. The three-dimensional boundary range is jointly determined by the geographic boundary and the height boundary of the monitoring area: the geographic boundary is obtained through scene mapping or existing facility coordinate data and converted to the scene coordinate system; the height boundary is determined by the highest possible height of the monitored object and the ground reference height, and is recorded in numerical form.

[0038] Specifically, the 3D segmentation process is as follows: The boundary span of each coordinate axis is calculated, and the span is the difference between the maximum and minimum boundary coordinates. Then, using the dynamic voxel reference value as the voxel side length, the span is divided into equally spaced segments based on the voxel side length. The number of segments is rounded up to cover the entire boundary range. The rounded-up number of segments is used to generate the voxel index range for that axis. After segmentation along the three coordinate axes, the 3D voxel mesh is generated by a Cartesian combination of the three directional indices. Each voxel corresponds to a unique 3D index number and is organized in the storage structure according to the index order.

[0039] Furthermore, when spatially mapping radar echo data using a spatial voxel grid, the measurement results of the radar echo points are first converted into spatial coordinates.

[0040] Specifically, the conversion process is as follows: The distance information measured by the radar is combined with the scanning direction information to obtain the three-dimensional coordinates of the point based on the radar coordinate system. Then, the coordinates are converted to the scene coordinate system based on the coordinate transformation parameters obtained from installation and calibration. After the coordinate transformation is completed, voxel indexing is performed on each point: In each coordinate axis direction, the point's coordinate value is subtracted from the minimum boundary coordinate of that axis to obtain the relative coordinate value. The relative coordinate value is then segmented according to the voxel side length, with the integer part of the segment number used as the axis index number. The three axis index numbers are combined to obtain the three-dimensional index number of the voxel to which the point belongs. If the coordinates of a point exceed the boundary range of any axis, the point is determined not to belong to the monitoring space and is removed from the subsequent processing chain, ensuring that the mapping result is strictly constrained by the monitoring boundary.

[0041] Furthermore, the assumption that the foreign object occupies at least one voxel is based on the constraint logic of the voxel side length value and the coverage logic of the 3D segmentation. The voxel side length value is derived from the minimum shape scale of the foreign object and the beam coverage scale, ensuring that the spatial occupancy range of the foreign object corresponds to at least one voxel index after voxelization: when the voxel side length value is not less than the minimum shape scale of the foreign object, the minimum coverage volume of the foreign object in any posture can fall into at least one voxel; when the voxel side length value is less than the minimum shape scale of the foreign object, the spatial occupancy range of the foreign object spans multiple voxel indices, thus also satisfying the requirement that at least one voxel is occupied.

[0042] At the same time, the three-dimensional segmentation covers the entire monitoring boundary range by rounding up, ensuring that when a foreign object is in the monitoring space, its spatial coordinates are mapped to the effective voxel index range, thereby forming a spatial occupancy relationship in the voxel grid that is consistent with the position of the foreign object.

[0043] Example 3: S2 includes defining discrete sampling intervals in a three-dimensional coordinate system using the side length of a spatial voxel grid, calculating the distribution density of the reflected point cloud within the discrete sampling interval to obtain the voxel activation ratio; mapping the voxel activation ratio in a preset probability distribution model to determine the point cloud persistence state under different reflection intensities, and establishing the relationship between the voxel activation ratio and the ultra-long temporal point cloud accumulation time; retrieving the corresponding confidence interval in the relationship according to the preset false alarm rate, and extracting the temporal feature vector corresponding to the confidence interval; obtaining the minimum accumulation time by multiplying the temporal feature vector with the sampling frequency, thus realizing the reverse calculation of the minimum accumulation time based on the preset false alarm rate.

[0044] In this embodiment, the determination of the minimum accumulation time revolves around the "convergence law of voxel activation ratio with accumulation time". The goal is to convert the risk constraint corresponding to the preset false alarm rate into a clear accumulation time length, so that the background sampling can reach the specified confidence level before entering the subsequent background distribution modeling process.

[0045] Specifically, the side length of the voxel grid has been derived from the minimum physical size of the foreign object and the detection distance in the previous steps, and is directly used as the scale parameter for spatial discretization here. During the system deployment phase, the side length is configured in length units and the precision bits are locked. The precision bits are jointly determined by the scene coordinate calibration accuracy and the radar measurement resolution to ensure that all subsequent volume, density, and time calculations use the same numerical precision.

[0046] Furthermore, the delineation of discrete sampling intervals is completed using the monitoring space boundary and voxel grid edge length as inputs. The monitoring space boundary is obtained through scene mapping and installation calibration, specifically: the minimum and maximum boundary coordinates in the scene coordinate system are given, and the boundary coordinates are configured numerically and remain unchanged. The delineation process is as follows: starting from the minimum boundary coordinates in each of the three coordinate axes, the segment is continuously advanced in the positive direction with the voxel grid edge length as the segment step size until the maximum boundary coordinates are covered; each advancement forms a one-dimensional segment, and the segment is assigned an index starting from 0 and incrementing. After completing the segmentation along the three axes, the three axis indices are combined to obtain the three-dimensional discrete sampling interval; each discrete sampling interval corresponds to one voxel unit, and the voxel unit is uniquely indexed in the storage structure using the three-dimensional index.

[0047] Specifically, the volume of a voxel cell is calculated once after it is defined. The calculation process is as follows: multiply the side length of the voxel mesh by itself three times to obtain the volume value, and use this volume value as the unified denominator for subsequent density normalization to avoid time fluctuations caused by repeated calculations during operation.

[0048] Furthermore, the calculation of the voxel activation ratio begins with the validity determination of the reflection point cloud. Reflection intensity, as a validity determination parameter, originates from the intensity value of the radar echo output. This intensity value undergoes zero-point drift calibration and quantization step confirmation during the deployment phase. The quantization step is directly determined and registered by the device output format. The effective reflection intensity threshold is determined upon its first occurrence as follows: Within the monitoring space, a non-physical reflection direction is selected, and echo intensity sequences for at least one cumulative duration configuration period are continuously collected. The intensity with the highest frequency in this sequence is taken as the noise principal component. Then, the upper limit of the noise intensity sequence's fluctuation is calculated. The calculation process for the upper limit of fluctuation involves first calculating the average noise intensity, then calculating the square of the difference between each noise intensity and the average, and averaging the results. The square root of this average square is then taken to obtain the fluctuation amount. Finally, the average value is added to three times the fluctuation amount to obtain the noise upper limit. The effective reflection intensity threshold is the noise upper limit plus one quantization step to ensure that pure noise echoes do not enter the valid point statistics.

[0049] After the threshold is determined, a point-by-point judgment is performed on the point cloud output in each scanning cycle: points with intensity values ​​less than the effective reflection intensity threshold are not included in the statistics, while points with intensity values ​​not less than the effective reflection intensity threshold are included in the voxel statistics.

[0050] Furthermore, the distribution density is used to quantify the spatial concentration of effective reflection points within the discrete sampling interval, and its calculation is strictly performed in the order of "spatial attribution, counting, volume normalization, and time normalization".

[0051] Specifically, spatial attribution is achieved through voxel indexing: the three-dimensional coordinates of a point are subtracted from the minimum boundary coordinates of the monitored space to obtain relative coordinates in three directions. These relative coordinates are then divided by the voxel grid side length, and the integer part is taken to obtain three direction indices. These three direction indices are combined to form the voxel cell index to which the point belongs. If any direction indice exceeds the configured range, the point is not included in the statistics. Counting is performed independently within each voxel cell: valid points entering the voxel cell are cumulatively counted, and scan cycles participating in the statistics are also cumulatively counted. The scan cycle count is based on the scan start timestamp, which is provided by the radar time source or synchronization time source and time consistency is verified during the deployment phase.

[0052] Specifically, volume normalization is performed within each voxel cell: the cumulative count of effective points in that voxel cell is divided by the voxel cell volume to obtain the number of points per unit volume. Time normalization is then performed: the number of points per unit volume is divided by the cumulative count of the scan cycle to obtain the "number of points per unit volume for a single scan," which is the distribution density used for subsequent percentage calculations.

[0053] Furthermore, the voxel activation percentage is derived from the distribution density. To avoid the lack of a reference for the percentage leading to meaningless values, a baseline density parameter needs to be introduced.

[0054] The baseline density is determined during the deployment and calibration phase. The determination process is as follows: Under the premise that the monitored space is in background condition and free of foreign objects, point cloud data covering one minimum cumulative duration configuration period is continuously collected. The distribution density sequence for each voxel unit is calculated using the aforementioned method. This sequence undergoes steady-state screening, with the density change amplitude threshold as the criterion. The density change amplitude threshold is determined upon its first occurrence as follows: The absolute values ​​of the distribution density differences between two adjacent scan periods in the background acquisition sequence are taken to form a difference sequence. The average value and fluctuation of this difference sequence are calculated. The fluctuation is calculated by first averaging the difference sequence, then averaging the squares of the differences between each difference and the average value, and finally adding the average value to three times the fluctuation to obtain the change amplitude threshold. When the difference in a certain scan period does not exceed this threshold, that period is considered to be in a steady-state segment. After the steady-state segment is determined, the average distribution density within the steady-state segment is calculated to obtain the baseline density of that voxel unit. The calculation of the voxel activation percentage is performed on a voxel-by-voxel basis during the runtime phase: the current distribution density is divided by the reference density of the voxel unit to obtain the ratio. If the ratio is greater than 1, it is taken as 1; if the ratio is less than 0, it is taken as 0; otherwise, the original value of the ratio is taken, thus obtaining an activation percentage value between 0 and 1.

[0055] Furthermore, a preset probability distribution model is used to map the "activation percentage value" to the "point cloud survival state", and further associates the "cumulative time" with the "probability of reaching the survival state".

[0056] The reflection intensity grading is determined as follows when it first appears: the minimum and maximum values ​​of the effective reflection intensity are counted in the background calibration data. The minimum value is the minimum observed value among all effective reflection intensities, and the maximum value is the maximum observed value among all effective reflection intensities. Then, the intensity width of a single grade is determined, which is 4 times the intensity quantization step. The total intensity width is obtained by subtracting the minimum value from the maximum value. The total intensity width is then divided by the intensity width of a single grade and rounded up to obtain the number of grades. Finally, starting from the minimum value, the intensity intervals of each grade are continuously divided by the intensity width of a single grade, and each grade is assigned a grade number that increases from 0. The point cloud persistence state is defined by two types: background non-converged state and background converged state. The convergence judgment threshold is determined when it first appears as follows: In the background calibration data, the activation proportion sequence is calculated for each voxel unit to form a sequence of absolute values ​​of the difference in activation proportion between two adjacent scan cycles. Then, the proportion change threshold is obtained by adding three times the fluctuation amount to the average value. When the proportion difference of multiple consecutive scan cycles does not exceed the proportion change threshold, it is judged as background converged state. The number of scan cycles for consecutive judgment is set as "the number of scan cycles corresponding to the longest fluctuation duration observed in the background calibration data plus 1" when it first appears. The longest fluctuation duration is located by counting the scan cycles.

[0057] After defining the states, the preset probability distribution model is implemented using a lookup table structure: indexed by the tier number, cumulative time scale number, and activation percentage interval number, each table entry stores the probability of being in a background convergence state and the probability of being in a background non-convergence state. The table entry probabilities are derived from the frequency statistics of the background calibration data. The statistical process involves summarizing the activation percentage of all voxel units under the same tier number and cumulative time scale number, dividing the count by the activation percentage interval, and then dividing the count by the total number of samples to obtain the probability value. Based on this lookup table structure, the relationship between the voxel activation percentage and the ultra-long temporal point cloud accumulation time is established for each cumulative time scale: starting from the accumulation start point, the cumulative time scale is incremented, and for each scale, the lookup table structure is called to obtain the background convergence probability. This probability is then registered together with the corresponding time scale, forming a relationship sequence ordered by time. This relationship sequence is maintained independently within each intensity tier.

[0058] As a preferred implementation, the preset false alarm rate is used to apply the risk indicator to the time scale selection. The preset false alarm rate is determined as follows when it first occurs: the scene manager provides the upper limit of the target number of false alarms allowed within one inspection cycle. The upper limit of the target number of false alarms is divided by the number of judgments within the inspection cycle to obtain the probability of false alarms allowed in a single judgment. The number of judgments is obtained by dividing the inspection cycle duration by the system judgment cycle duration and rounding down. The system judgment cycle duration is taken as the average interval between the start timestamps of two consecutive scans. The obtained probability of false alarms allowed in a single judgment is used as the preset false alarm rate configuration value and remains unchanged.

[0059] Furthermore, when retrieving confidence intervals based on the preset false negative rate, a time-scale determination is made in the relation sequence corresponding to each intensity level: the "probability of the background non-convergence state" is compared with the preset false negative rate. When the non-convergence probability is not greater than the preset false negative rate, the time scale is considered to meet the confidence requirement; the time scale that first meets the confidence requirement is recorded as the starting point. To avoid misjudgment of the starting point due to occasional probability fluctuations, a continuous satisfaction threshold is introduced. The continuous satisfaction threshold is taken as "the number of time scales corresponding to the longest continuous length of the non-convergence probability rebound in the background calibration data plus 1" when it first appears. The longest continuous length is obtained by counting the rebound segments of the relation sequence; the starting point is only confirmed when the number of time scales indicated by the threshold is continuously satisfied from the starting point. The endpoint of the confidence interval is taken as the end time scale of the continuous satisfaction condition segment after the starting point is confirmed. The interval is jointly determined by the starting point and the endpoint and used to extract time-domain feature information.

[0060] Furthermore, the temporal feature vector is used to carry the number of scan cycles required for each intensity level to reach the confidence requirement. Its composition and extraction process is as follows: within each intensity level, the cumulative time scale corresponding to the confirmation starting point is converted into the number of scan cycles. The conversion method is to use the cumulative scan cycle count corresponding to the time scale as the value; the number of scan cycles obtained for each intensity level is arranged in ascending order according to the level number to form a temporal feature vector and registered.

[0061] Specifically, the sampling frequency is used to convert the number of scan cycles into the actual time length. The sampling frequency is determined when it first appears as follows: the starting timestamps of no less than 300 consecutive scans are recorded, the time difference between two adjacent timestamps is calculated to form an interval sequence, the maximum and minimum values ​​in the interval sequence are removed, and the remaining intervals are averaged to obtain the single scan time length; then the unit time length is divided by the single scan time length to obtain the number of scans per unit time, and this number of scans per unit time is used as the sampling frequency configuration value.

[0062] The calculation of the minimum cumulative duration is carried out item by item according to the time domain feature vector: the number of scanning cycles of each item in the vector is multiplied by the length of a single scan to obtain the candidate value of the cumulative duration corresponding to the grade. Then, the largest candidate value is selected from all candidate values ​​as the minimum cumulative duration, so as to ensure that the background of any intensity grade reaches the confidence level required by the preset false alarm rate within the cumulative duration.

[0063] Example 4: S3 includes comparing the current ultra-long temporal point cloud accumulation time with a preset minimum accumulation time. If the current ultra-long temporal point cloud accumulation time is greater than or equal to the minimum accumulation time, the three-dimensional coordinate set of all discrete sampling points within the voxel grid is extracted. An arithmetic mean operation is performed on the three-dimensional coordinate set to obtain the mean vector representing the geometric center of the voxel grid. The mean vector and the three-dimensional coordinate set are used to perform a sum of squared differences to determine the covariance matrix reflecting the spatial distribution state of the point cloud. Based on the eigenvalue decomposition results of the covariance matrix, feature vectors describing the discreteness and squareness of the point cloud are extracted, and the feature vectors are associated and mapped with the mean vector to obtain the background distribution model. By writing the background distribution model into a preset storage space, the persistent storage of the geometric center, discreteness, and squareness features of the point cloud within the background voxels is achieved.

[0064] In this embodiment, the numerical relationship between the "current ultra-long temporal point cloud accumulation time" and the "preset minimum accumulation time" is used as the modeling trigger condition. The current ultra-long temporal point cloud accumulation time is recorded at the start of the background accumulation phase, with the start time being the timestamp of the first frame of point cloud data being received. Subsequently, at the end of each frame of point cloud processing, the current timestamp is recorded, and the current ultra-long temporal point cloud accumulation time is obtained by subtracting the start timestamp from the current timestamp. This timestamp is provided by a unified time source and maintains the same time unit.

[0065] The preset minimum accumulation time is derived from the reverse calculation result of the preceding steps and is issued as a configuration parameter during the system deployment phase. The configuration includes the time unit and numerical precision, with the numerical precision being the decimal places corresponding to the scene time synchronization precision. When a judgment is triggered, the current ultra-long temporal domain point cloud accumulation time is compared with the preset minimum accumulation time after being unified to the same time unit. If the current ultra-long temporal domain point cloud accumulation time is not less than the preset minimum accumulation time, the background distribution model establishment process is initiated; if the comparison result is less than the preset minimum accumulation time, the accumulation process is maintained and the current ultra-long temporal domain point cloud accumulation time continues to be updated.

[0066] Further, after entering the background distribution model establishment process, the first step is to form a "3D coordinate set of discrete sampling points within the voxel grid". The voxel grid is obtained by 3D subdividing the voxel grid edge length determined in the previous step in the scene coordinate system. Each voxel unit has a unique index, which is formed by combining interval numbers in three directions and represented by an integer. The extraction of the 3D coordinate set is performed one voxel unit at a time: for each frame of point cloud, the coordinate system is first transformed, converting the radar-measured coordinates to scene coordinate system coordinates; for each point after transformation, its relative displacement in the three coordinate directions is calculated, and the relative displacement is the point coordinate minus the minimum boundary coordinate of the monitoring space; then, the relative displacement is assigned to intervals according to the voxel grid edge length. The assignment process is to remove the relative position and use the voxel grid edge length and take the integer part as the interval number in that direction. The combination of the interval numbers in the three directions gives the voxel unit index to which the point belongs; then, the 3D coordinates of the point are merged into the coordinate set of the corresponding voxel unit according to the index.

[0067] Specifically, to ensure that the background distribution model is generated only from background voxels, the background voxel determination parameter is determined during the deployment phase and saved as a configuration item. The determination parameter is taken as "the point cloud persistence state of the voxel unit is in the background convergence state during the accumulation phase"; where the point cloud persistence state is obtained by the preceding probability mapping, and the determination criterion corresponding to the background convergence state and its probability output are configured together and kept consistent during the deployment phase. Voxel units determined to be background voxels retain their 3D coordinate set, while non-background voxels are not included in this modeling process. If a background voxel does not form any discrete sampling points within the accumulation window, no subsequent statistics are generated for that voxel, it is recorded as a data-free state, and the voxel index is retained for subsequent supplementary sampling and updates.

[0068] Furthermore, when performing an arithmetic mean operation on the three-dimensional coordinate set to obtain the mean vector, the mean vector consists of three components, which correspond to the center positions of the point cloud in the voxel unit in the three coordinate directions.

[0069] Specifically, the calculation process is completed independently within a single voxel unit: First, the number of points in the 3D coordinate set of the voxel unit is counted, represented by an integer and used as the divisor parameter for the mean calculation; then, the first coordinate component of all points in the set is summed point-by-point, and the summation is stored with uniform numerical precision, which is the number of decimal places corresponding to the coordinate quantization step; the sum of the first coordinate components is divided by the number of points to obtain the first mean component. The above summation and division process is repeated for the second coordinate component to obtain the second mean component, and the above process is repeated for the third coordinate component to obtain the third mean component. The three mean components are combined to form a mean vector, which is used to characterize the geometric center position of the background point cloud within the voxel unit. All subsequent discretization and squareness calculations use this mean vector as the central reference.

[0070] Furthermore, when using the mean vector and the three-dimensional coordinate set to perform the sum of squared differences to determine the covariance matrix, the covariance matrix adopts a 3-row, 3-column structure to express the discrete intensity and directional correlation of the point cloud in three directions within the voxel unit.

[0071] Specifically, the calculation process is carried out point by point and accumulated item by item within a single voxel unit: First, cumulative quantities are established for the nine positions of the covariance matrix, with the cumulative quantities initially set to 0 and using uniform numerical precision; then, each point in the three-dimensional coordinate set of the voxel unit is traversed, and the first mean component of the mean vector is subtracted from the first coordinate component of the point to obtain the first deviation component, the second mean component is subtracted from the second coordinate component to obtain the second deviation component, and the third mean component is subtracted from the third coordinate component to obtain the third deviation component. After obtaining three deviation components, calculate the product of each pair of deviation components in turn and accumulate them to the corresponding positions: the product of the first deviation component and the first deviation component is accumulated in the first row and first column; the product of the first deviation component and the second deviation component is accumulated in the first row and second column; the product of the first deviation component and the third deviation component is accumulated in the first row and third column; the product of the second deviation component and the first deviation component is accumulated in the second row and first column; the product of the second deviation component and the second deviation component is accumulated in the second row and second column; the product of the second deviation component and the third deviation component is accumulated in the second row and third column; the product of the third deviation component and the first deviation component is accumulated in the third row and first column; the product of the third deviation component and the second deviation component is accumulated in the third row and second column; the product of the third deviation component and the third deviation component is accumulated in the third row and third column.

[0072] After accumulating all points, the nine accumulated values ​​are divided by the number of points to obtain normalized results. These normalized results are arranged in 3 rows and 3 columns to form a covariance matrix. To ensure the symmetry of the covariance matrix, consistency checks are performed on the following columns after calculation: 1 row 2 column and 2 row 1 column, 1 row 3 column and 3 row 1 column, 2 row 3 column and 3 row 2 column. The check criterion is that the absolute value of the difference between the two values ​​does not exceed the variance magnitude corresponding to the coordinate quantization step. The variance magnitude is determined by squaring the coordinate quantization step, and this numerical scale is pre-calculated based on the coordinate quantization step and saved as a check parameter during the deployment phase. If the check fails, the arithmetic mean of the two values ​​is used to replace the values ​​at these two positions, ensuring that the covariance matrix meets the symmetry requirement.

[0073] As a preferred embodiment, when extracting discreteness and squareness eigenvectors based on the eigenvalue decomposition results of the covariance matrix, the covariance matrix is ​​first subjected to eigenvalue decomposition to obtain three eigenvalues ​​and their corresponding three eigenvectors.

[0074] Specifically, the decomposition process employs a symmetric real matrix diagonalization solution: the covariance matrix is ​​used as input to initialize the diagonalization iteration state. The iteration process aims to gradually reduce the off-diagonal elements until the absolute values ​​of all off-diagonal elements do not exceed the convergence threshold, at which point the process terminates. The convergence threshold is determined upon its first occurrence as follows: based on the coordinate quantization step, the variance magnitude corresponding to the coordinate quantization step is first calculated, and then this variance magnitude is used as the convergence threshold value, ensuring that the numerical stability of the decomposition result is consistent with the coordinate resolution capability. After decomposition, three eigenvalues ​​and three eigenvectors are obtained. The eigenvalues ​​are sorted according to their numerical values ​​from largest to smallest, and the eigenvectors are rearranged in the same order to ensure that the largest eigenvalue corresponds to the principal expansion direction. The discrete eigenvectors are extracted using the sorted three eigenvalues ​​as input to form three components: the first component takes the largest eigenvalue, the second component takes the middle eigenvalue, and the third component takes the smallest eigenvalue, used to directly characterize the discrete intensity of the point cloud in the three principal directions.

[0075] The squared-degree eigenvector represents morphological uniformity by the proportional relationship between eigenvalues. The extraction process involves dividing the minimum and maximum eigenvalues ​​to obtain the first proportional component, dividing the intermediate and maximum eigenvalues ​​to obtain the second proportional component, and then combining the first and second proportional components in sequence to form the squared-degree eigenvector. The numerical precision of the proportional components is taken from the decimal places of the covariance matrix to ensure consistent precision before and after the proportional calculation. After eigenvector extraction, the mean vector, covariance matrix, discrete eigenvector, and squared-degree eigenvector are associated under the same voxel index. The association uses the voxel unit index as the retrieval key, which is a unique integer key value formed by concatenating three directional interval indices. The key generation rules are determined and maintained consistently during the deployment phase to avoid key conflicts during runtime.

[0076] Furthermore, the long-term storage of the background distribution model is accomplished using a pre-configured storage space. This pre-configured storage space is configured during system initialization, with configuration parameters including storage medium identifier, storage path or storage area number, single record length, index structure type, capacity limit, and verification strategy. The single record length is determined by the number of fields and the precision of the field values. Each field must include at least a voxel index, three components of the mean vector, nine positional values ​​in the covariance matrix, three components of the discrete eigenvector, two components of the squared degree eigenvector, and a modeling completion timestamp. The storage process is executed in voxel index order: a record is generated for each background voxel, and the record content is serialized according to the field order before being submitted to the pre-configured storage space. The serialized numerical precision is the system-configured precision, and the timestamp precision is the time synchronization precision.

[0077] The index structure is updated synchronously during the saving process. The index entries are located by voxel index, and the location information includes the record start offset and record length. This is used to quickly read the corresponding background statistical information by voxel index in real time.

[0078] Specifically, after saving, a consistency check is performed. The consistency check parameters are determined as follows when they first appear: the record quantity check is taken as the number of background voxels, which is obtained by statistically analyzing the voxel grid index range and the background voxel determination results; the index range check is taken as the minimum and maximum voxel index values, which are determined by the voxel grid boundaries; the timestamp check is taken as the modeling completion timestamp being no earlier than the background accumulation start timestamp and no later than the current processing frame timestamp. After all checks are satisfied, the background distribution model enters a readable state. In the subsequent real-time detection phase, the corresponding mean vector and covariance matrix and their derived features are directly retrieved based on the voxel index, realizing the long-term storage and stable reuse of background geometric center, dispersion, and squareness information.

[0079] Example 5: S4 includes acquiring the current point cloud coordinate set generated by the non-repeating scanning lidar, determining the voxel grid position to which the point cloud coordinate set belongs by dividing the coordinate value range; extracting the mean vector and covariance matrix corresponding to the voxel grid position from the preset background distribution model as initial distribution parameters; performing matrix operations on the mean vector and covariance matrix to obtain the deviation vector group between the point cloud coordinate set and the background distribution; calculating the Mahalanobis distance value based on the deviation vector group; if the Mahalanobis distance value is less than the preset similarity threshold, clustering the deviation vector group to determine the abnormal point cloud subset; executing the density estimation algorithm based on the abnormal point cloud subset to generate updated background distribution parameters; and writing the updated background distribution parameters back into the background distribution model to realize the dynamic reading and adjustment of the mean vector and covariance matrix of the voxels to which the point cloud coordinates belong during the real-time detection stage.

[0080] In this embodiment, the real-time detection phase takes the current point cloud coordinate set output by the non-repeating scanning lidar as input. First, it determines the spatial affiliation of the point cloud coordinate set in the voxel grid. Then, it calculates the deviation vector group and Mahalanobis distance based on the mean vector and covariance matrix in the background distribution model. When the similarity constraint is met, it performs clustering separation on the deviation vector group. Subsequently, it performs density estimation based on the abnormal point cloud subset to generate updated background distribution parameters. Finally, it updates the background distribution parameters into the background distribution model, so that the mean vector and covariance matrix of the voxel in subsequent frames are continuously adjusted with real-time data.

[0081] Specifically, the voxel grid side length, as a spatial discrete scale parameter, is given by the preceding voxel setting step and remains unchanged as a system configuration parameter; the boundary range, as a monitoring space parameter, is determined by the minimum and maximum boundary coordinates of the scene coordinate system during deployment and calibration and remains unchanged; the coordinate quantization step, as a coordinate accuracy parameter, is determined and registered by the radar output resolution and the number of decimal places in the coordinate transformation link.

[0082] Furthermore, after obtaining the current point cloud coordinate set, the voxel grid position to which the point cloud coordinate set belongs is first determined by dividing the coordinate value range. The specific implementation process is as follows: For each point, its three coordinate components in the scene coordinate system are read. The scene coordinate system is determined by the installation calibration parameters, which are obtained and registered during the deployment phase by measuring the relationship between the radar installation position, installation attitude, and scene reference point. Then, boundary judgment is performed on each point. The boundary judgment parameters are the minimum and maximum boundary coordinates of the monitoring space. The judgment process is a component-by-component comparison. Points whose coordinate components are less than the corresponding minimum boundary or greater than the corresponding maximum boundary are judged as out-of-bounds points and removed from subsequent calculations in this step. For non-out-of-bounds points, the voxel index is calculated: the minimum boundary coordinate in the corresponding direction is subtracted from each coordinate component of the point to obtain the relative displacement. Then, the relative displacement is removed, and the interval number in that direction is obtained by taking the integer part of the voxel grid side length. The three interval numbers are combined to obtain the voxel grid position. This voxel grid position is represented by an integer index. The index encoding rule is determined during system initialization and remains consistent to ensure that the same voxel position corresponds to the same background record.

[0083] Furthermore, after determining the voxel grid positions, the mean vector and covariance matrix corresponding to the voxel grid positions are extracted from the preset background distribution model as initial distribution parameters. The mean vector has three components, representing the geometric center coordinates of the background point cloud of the voxel; the covariance matrix is ​​a 3x3 numerical structure, representing the dispersion and directional correlation of the background point cloud of the voxel in three directions. Both are generated during the background modeling stage and stored in the storage space.

[0084] The specific extraction process is as follows: Using the voxel grid position as the retrieval key, the background distribution model record position is located. The mean vector and the nine position values ​​of the covariance matrix are read and loaded into the runtime memory. To ensure the stability of subsequent inverse matrix calculations, a reversibility check is performed after reading. The reversibility check parameter is the lower bound of the minimum eigenvalue. The lower bound of the minimum eigenvalue is determined during the deployment phase based on the coordinate quantization step. The determination process involves multiplying the coordinate quantization step by itself to obtain the variance magnitude, and then using this variance magnitude as the lower bound value and registering it. The check process involves solving for the eigenvalues ​​of the covariance matrix and taking the minimum value. This minimum value is compared with the lower bound of the minimum eigenvalue. If the minimum value is less than the lower bound, the variance magnitude is superimposed on the three diagonal positions of the covariance matrix to ensure that the covariance matrix meets the numerical stability requirements. The superimposed covariance matrix is ​​used as the initial distribution parameter for calculation in this frame.

[0085] Furthermore, after obtaining the mean vector and covariance matrix, deviation vector sets are calculated for the point cloud coordinate set and background distribution. The deviation vector set calculation is performed point-by-point: for each point, the first component of the mean vector is subtracted from the first coordinate component to obtain the first deviation component; the second component is subtracted from the second coordinate component to obtain the second deviation component; and the third component is subtracted from the third coordinate component to obtain the third deviation component. These three deviation components form the deviation vector for that point. This subtraction operation is repeated for all points within the voxel to form a deviation vector set. Subsequently, Mahalanobis distance values ​​are calculated based on the deviation vector sets, with the Mahalanobis distance values ​​calculated point-by-point to obtain a distance sequence.

[0086] The specific calculation process is as follows: First, the covariance matrix is ​​inverted. The inversion process adopts a symmetric matrix decomposition and back-substitution solution process. In the decomposition stage, the covariance matrix is ​​converted into a decomposed structure. In the back-substitution stage, each column of the unit basis vector is solved to obtain each column of the inverse matrix. The numerical precision of the inverse matrix is ​​taken from the precision of the covariance matrix. For the deviation vector of each point, the deviation vector is first multiplied by the inverse matrix to obtain the weighted deviation vector. Then, the weighted deviation vector is multiplied by the original deviation vector to obtain the squared distance. Finally, the square root of the squared distance is taken to obtain the Mahalanobis distance value of the point. If the squared distance is negative, the negative value is set to 0 before the square root operation is performed. The threshold for negative values ​​is taken as the negative number of the variance order corresponding to the coordinate quantization step.

[0087] Furthermore, a preset similarity threshold is used to determine whether the Mahalanobis distance value meets the background similarity constraint. This threshold is determined during the deployment phase and registered as a configuration parameter. The similarity threshold is determined based on the target false alarm constraint, which is given by the scene operation procedure in the form of the maximum number of false alarms allowed within a single determination period. The single determination period is determined by the average difference between the timestamps of two adjacent point cloud frames. The average value is calculated by continuously recording the reception completion timestamps of no less than 300 point cloud frames, calculating the difference between adjacent timestamps to form an interval sequence, removing the maximum and minimum intervals, and averaging the remaining intervals to obtain the determination period duration.

[0088] The threshold calibration process is as follows: under background conditions, point cloud data covering no less than 10 minimum cumulative durations are collected, and Mahalanobis distance value sequences are calculated within each voxel after voxel mapping; the Mahalanobis distance values ​​of all voxels are merged to form a background distance set, the number of out-of-threshold points in the background distance set under different thresholds is counted, and the number of out-of-threshold points is divided by the total number of points in the background distance set to obtain the background out-of-threshold ratio; the minimum threshold that makes the background out-of-threshold ratio no greater than the allowable ratio obtained by converting the target false alarm constraint is selected as the similarity threshold.

[0089] Specifically, the conversion process for the allowable ratio is as follows: divide the maximum number of false alarms allowed in a single judgment period by the maximum number of points in a single judgment period. The maximum number of points is determined and registered by the maximum number of points of the radar in the highest frame rate and maximum echo output mode.

[0090] Furthermore, when the Mahalanobis distance value within a voxel is less than the similarity threshold, clustering is performed on the deviation vector group within that voxel to determine a subset of the outlier point cloud. The clustering process employs a density-connected separation strategy, with the input being the three deviation components of each point in the deviation vector group.

[0091] Specifically, the clustering parameters include the neighborhood radius and the minimum number of clustered points. The neighborhood radius, as a spatial proximity determination parameter, is determined during the deployment phase by the voxel grid edge length and the estimated background point spacing. The estimated background point spacing is obtained from background calibration data. The process involves calculating the nearest neighbor distance for each point within each voxel, summing the nearest neighbor distance sequence, and averaging it as the estimated background point spacing. The neighborhood radius is taken as 0.5 times the sum of the voxel grid edge length and the estimated background point spacing and is recorded to ensure that the neighborhood coverage is consistent with the voxel scale and the point density.

[0092] The minimum number of clustered points is used as a parameter to determine the validity of clustering. It is determined by the average number of points of the voxel in the background calibration stage. The average number of points is determined by counting the number of points of the voxel in each frame of the background calibration data and averaging them. The minimum number of clustered points is the average number of points multiplied by 0.1, rounded up, and recorded to ensure that the clustering results have sufficient sample support. The clustering process is as follows: First, calculate the Euclidean distance between each point in the deviation vector group and the other points, and filter the point set in the neighborhood by the neighborhood radius; when the number of points in the neighborhood of a point is not less than the minimum number of clustered points, the point is marked as a core point; starting from the core point, merge the point set in its neighborhood into the same cluster, and repeat the neighborhood expansion for newly added points in the cluster until there are no more newly added points to obtain a cluster; traverse all points until all clusters are obtained. After clusters are generated, subsets of anomalous point clouds are determined. The determination process is as follows: For each cluster, the mean Mahalanobis distance and the number of points in the cluster are calculated. The mean Mahalanobis distance is calculated by summing all Mahalanobis distance values ​​in the cluster and dividing by the number of points in the cluster. The clusters are sorted in ascending order of the mean Mahalanobis distance, and the cluster with the smallest mean is taken as the background-consistent cluster. The remaining clusters are taken as subsets of anomalous point clouds, thereby separating local deviation structures under the premise that the overall similarity constraint is met.

[0093] Furthermore, after obtaining the subset of the outlier point cloud, a density estimation algorithm is executed to generate updated background distribution parameters, which include the updated mean vector and the updated covariance matrix.

[0094] Specifically, density estimation uses the coordinates of points in the anomaly cloud subset as input. First, the density weight of each point is obtained, and then weighted statistics are performed using these density weights. The calculation process for density weights is as follows: For each point in the anomaly cloud subset, the number of neighboring points within its neighborhood radius is counted; the larger the number of neighboring points, the higher the density. The number of neighboring points is divided by the number of points in the anomaly cloud subset to obtain a normalized density value. Upper and lower limits are applied to the normalized density value, with an upper limit of 1 and a lower limit of 0.01. The lower limit is used to avoid excessively small weights leading to an excessively small weighted sum, while the upper limit is used to maintain a consistent weight scale. The pruning threshold is determined during the deployment phase based on the range of the number of neighboring points in the background calibration data. The determination process involves statistically analyzing the minimum and maximum values ​​of the number of neighboring points in the background calibration data, and mapping the normalized values ​​corresponding to the minimum and maximum values ​​to 0.01 and 1, respectively.

[0095] The updated mean vector is obtained by weighted averaging. The calculation process is as follows: weighted summation is performed on the three coordinate directions respectively. First, the coordinate component of each point in that direction is multiplied by the density weight and accumulated to obtain the weighted sum. Then, the density weights are accumulated to obtain the weighted sum. Finally, the weighted sum is divided by the weighted sum to obtain the updated mean component in that direction. The updated mean components in the three directions are combined to form the updated mean vector.

[0096] The updated covariance matrix is ​​obtained through weighted deviation accumulation. The calculation process is as follows: for each point in the outlier cloud subset, calculate its three deviation components relative to the updated mean vector, and then multiply the deviation components pairwise to obtain nine product terms; multiply each product term by the density weight of that point and accumulate it to the corresponding nine cumulative quantities; after accumulation, divide each cumulative quantity by the weight sum to obtain the normalized result, and arrange the normalized result in 3 rows and 3 columns to form the updated covariance matrix; then perform symmetry verification and minimum eigenvalue lower bound verification on the updated covariance matrix again. The symmetry verification criterion is taken as the variance level corresponding to the coordinate quantization step, and the minimum eigenvalue lower bound is taken as the aforementioned invertibility verification parameter to ensure that the updated parameters have numerical stability.

[0097] Furthermore, when updating the background distribution parameters to the background distribution model, to avoid abrupt changes in background parameters due to single-frame fluctuations, an update coefficient is introduced as a smoothing parameter. The update coefficient is determined during the deployment phase based on the environmental change rate, which is obtained by differentiating the time window of the background calibration data: In the continuous background data, adjacent windows are divided according to the window length. The mean vector and covariance matrix corresponding to each window are calculated separately. Then, the sum of the absolute values ​​of the component differences of the mean vectors between adjacent windows and the sum of the absolute values ​​of the differences at the nine positions of the covariance matrix are calculated. These differences are divided by the window duration to obtain the drift rate and the morphological rate, and the larger of the two is taken as the environmental change rate. The environmental change rate is linearly mapped to a preset minimum and maximum rate to obtain an update coefficient between 0 and 1, and the minimum and maximum rates are registered as configuration parameters.

[0098] Specifically, the update process is as follows: The three components of the mean vector are fused component by component. For each component, the original background component is multiplied by 1, the update coefficient is subtracted, and then the updated component is multiplied by the update coefficient to obtain a new component. The nine positions of the covariance matrix are fused position by position, using the same fusion method as the mean vector. After fusion, the new mean vector and the new covariance matrix are updated to the corresponding record in the background distribution model according to the voxel grid position index, and the update timestamp of the record is updated synchronously. The update timestamp is provided by a unified time source and maintains the same time unit as the point cloud timestamp. After the update, a consistency check is performed on the index structure. The consistency check includes that the voxel index has not changed, the record position has not changed, and the record field length has not changed. After the check passes, the mean vector and covariance matrix extracted from the voxel in subsequent frame processing will be the updated values.

[0099] Example 6: S5 includes acquiring the real-time point cloud coordinates collected by the non-repeating scanning lidar and mapping them to the corresponding voxel grids; extracting the mean vector and covariance matrix corresponding to the voxel grids from the preset background distribution model; using the mean vector and covariance matrix to perform a multi-dimensional spatial difference measurement on the real-time point cloud coordinates to obtain the Mahalanobis distance; retrieving the cumulative scanning frequency of the voxel grids in the historical background state to determine the dynamic threshold; if the Mahalanobis distance exceeds the dynamic threshold, the corresponding real-time point cloud coordinates are written into the abnormal candidate point set, wherein the dynamic threshold is positively correlated with the number of times the voxel is scanned in the background state.

[0100] In this embodiment, the screening of anomaly candidates for real-time point cloud coordinates is completed within the voxel grid scale. The processing chain starts from "coordinate mapping", ends at "Mahathano distance calculation" and "dynamic threshold comparison". Its core is to use the mean vector and covariance matrix provided by the background distribution model to perform multidimensional difference measurement on real-time points, and introduce a dynamic threshold that is positively correlated with the number of historical background scans, so that different voxels can maintain a consistent false alarm constraint under different levels of observation.

[0101] Among them, the voxel grid side length, as a spatial discrete scale parameter, is derived from the previous voxel setting results and registered during the system deployment phase; the minimum and maximum boundary coordinates of the monitoring space, as spatial range parameters, are derived from scene mapping and installation calibration and registered during the deployment phase; and the coordinate quantization step, as a numerical accuracy parameter, is derived from the radar output resolution and the number of decimal places in the coordinate transformation link and registered during the deployment phase.

[0102] Furthermore, when acquiring the real-time point cloud coordinates collected by the non-repeating scanning lidar and mapping them to the corresponding voxel mesh, coordinate system unification is first completed. The installation calibration parameters used for coordinate system unification are determined during the deployment phase by measuring the spatial relationship between the radar installation position, installation attitude and scene reference point, and are registered in the form of transformation parameters.

[0103] Specifically, during the runtime phase, transformation parameters are called point by point for each frame of point cloud to convert the point coordinates in the radar coordinate system to the point coordinates in the scene coordinate system. Then, boundary determination is performed for each point. The boundary determination directly uses the minimum and maximum boundary coordinates of the monitoring space. The determination process is a component-by-component comparison: if any coordinate component is less than the minimum boundary coordinate in the corresponding direction, or if any coordinate component is greater than the maximum boundary coordinate in the corresponding direction, then the point will not enter the subsequent difference measurement link in this step.

[0104] For points that pass the boundary determination, the voxel index is calculated as follows: First, the relative displacement is obtained by subtracting the minimum boundary coordinate of the corresponding direction from the coordinate components of each direction of the point; then, the relative displacement is divided by the side length of the voxel grid and the integer part is taken to obtain the interval number of each direction; finally, the three interval numbers are combined into a voxel index. The encoding rules of the voxel index are determined and registered during system initialization.

[0105] Furthermore, when extracting the mean vector and covariance matrix corresponding to the voxel grid from the preset background distribution model, the background distribution model is formed and stored in the preset storage space during the background modeling stage. The record fields include at least the voxel index, three components of the mean vector, three rows and three columns of values ​​in the covariance matrix, and the modeling timestamp. The extraction process uses the voxel index as the retrieval key, locates the corresponding record, reads the mean vector and covariance matrix, and uses them as the background distribution parameters of the voxels.

[0106] Specifically, to ensure the numerical stability of the subsequent inverse matrix solution, the invertibility of the covariance matrix is ​​checked after extraction. The invertibility check parameter is the lower bound of the minimum eigenvalue. This lower bound is determined in the deployment phase based on the coordinate quantization step. The determination process is to multiply the coordinate quantization step by itself to obtain the variance magnitude, and then register this variance magnitude as the lower bound of the minimum eigenvalue.

[0107] The verification process involves finding the eigenvalues ​​of the covariance matrix and comparing the minimum eigenvalue with the lower limit. When the minimum eigenvalue is less than the lower limit, the variance magnitudes are superimposed on the three diagonal positions of the covariance matrix to ensure that the covariance matrix meets the stability requirements for subsequent inversion. After superposition, the minimum eigenvalue is found again to complete the verification. Once the verification is passed, the distance calculation begins.

[0108] As a preferred embodiment, when using the mean vector and covariance matrix to measure the multidimensional spatial difference of real-time point cloud coordinates to obtain the Mahalanobis distance, the calculation is performed independently at each real-time point. First, the deviation vector is constructed: the first deviation component is obtained by subtracting the first component of the mean vector from the first coordinate component of the real-time point; the second deviation component is obtained by subtracting the second component of the mean vector from the second coordinate component; and the third deviation component is obtained by subtracting the third component of the mean vector from the third coordinate component. The three deviation components are combined to form the deviation vector of that point.

[0109] The covariance matrix is ​​then inverted using a numerical process of symmetric matrix decomposition and back-substitution: first, the covariance matrix is ​​decomposed into an intermediate structure that facilitates inversion; then, each column of the unit basis vector is substituted back to obtain the columns of the inverse matrix. The numerical precision of the inverse matrix is ​​consistent with the precision of the covariance matrix. After obtaining the inverse matrix, it is multiplied by the deviation vector to obtain the weighted deviation vector. Then, the weighted deviation vector is multiplied by the deviation vector to obtain the squared distance. Finally, the square root of the squared distance is performed to obtain the Mahalanobis distance. If the squared distance has a small negative value, it is set to 0 according to the numerical stability rule before the square root is performed. The threshold for this negative value is a negative number of the variance order, which is calculated stepwise by coordinate quantization and registered during the deployment phase, ensuring that this correction only covers the range of numerical rounding errors.

[0110] Furthermore, when retrieving the cumulative scan frequency of the voxel mesh under historical background conditions to determine the dynamic threshold, the criteria for determining the "historical background condition" and the statistical criteria for the "cumulative scan frequency" must first be clarified. The historical background condition is determined by the point cloud persistence status. The point cloud persistence status is established during the background calibration stage and registered during the deployment stage. The registration content includes the criteria for determining the background convergence status and the continuous satisfaction threshold required for the determination.

[0111] Specifically, the cumulative scan frequency is calculated using a voxel-based system with a scan counter and time window. The scan counter accumulates the number of times a voxel is scanned in a background convergence state. The counting rule is that if a voxel has at least one boundary-judged point within a frame and is in a background convergence state, the scan count for that voxel in that frame is incremented by 1. The time window is used to calculate the frequency. The time window length is determined during the deployment phase, using a preset minimum cumulative duration to ensure consistency between the frequency statistics and the sample size for background modeling. The frequency calculation process involves dividing the scan count within the time window by the time window length to obtain the cumulative scan frequency of the voxel in the historical background state. Simultaneously, the scan count within the time window is retained as the cumulative scan count for compensation based on the number of scans when frequency is missing or the time window is not full.

[0112] Furthermore, the dynamic threshold generation adopts a monotonic mapping structure positively correlated with the number of scans. This mapping structure is determined by background calibration data and registered as a threshold table during the deployment phase. The threshold table generation process is as follows: Point cloud data covering multiple time windows are collected under background convergence conditions. Mahalanobis distance sequences are calculated by voxel, and scan counts within the corresponding time windows are recorded simultaneously. Subsequently, the scan counts are divided into several levels according to their numerical range. The level boundaries are determined by the minimum scan count, the maximum scan count, and the level step size. During deployment, the level step size is taken as an integer multiple of the radar frame rate and registered, ensuring that the level span is consistent with the temporal resolution. For each level, the Mahalanobis distances are arranged from smallest to largest within the corresponding samples, and candidate thresholds are tested one by one. The percentage of points exceeding the candidate threshold is counted and compared with the allowable percentage corresponding to the target false alarm constraint. The smallest candidate threshold that ensures the percentage of points exceeding the threshold is not greater than the allowable percentage is selected as the threshold for that level.

[0113] Specifically, the target false alarm constraint is given by the scenario operation procedure during the deployment phase, in the form of the maximum number of false alarms allowed within a single decision period. The allowable percentage is obtained by dividing the maximum number of false alarms allowed by the upper limit of the number of points within a single decision period. The upper limit of the number of points is determined and registered by the maximum number of points per period of the radar under the configuration of the highest frame rate and the maximum echo output. After the threshold table is completed, the operation phase locates the corresponding level based on the scan count of the voxel in the current time window, and reads the threshold of that level as the dynamic threshold. As the scan count increases, it locates a higher level and reads a larger threshold, thereby satisfying the constraint that the dynamic threshold is positively correlated with the number of voxel background scans.

[0114] Furthermore, when the Mahalanobis distance exceeds the dynamic threshold, the real-time point enters the anomaly candidate set. The anomaly candidate point set is organized by voxel index buckets, with the bucket key being the voxel index to ensure that subsequent neighborhood analysis is directly carried out within the local voxel range. The specific implementation process is as follows: For each real-time point, its Mahalanobis distance is compared with the dynamic threshold corresponding to the voxel; when the Mahalanobis distance is greater than the dynamic threshold, the point's 3D coordinates, voxel index, point cloud timestamp, and Mahalanobis distance value are registered in the anomaly candidate point set. The registration action is completed in the memory structure and a consistency check is performed before the end of the same frame processing. The consistency check includes voxel index existence check and timestamp monotonicity check; when the Mahalanobis distance is not greater than the dynamic threshold, the point does not enter the anomaly candidate point set.

[0115] Example 7: S6 includes acquiring the real-time point cloud coordinates collected by the non-repeating scanning lidar and mapping them to the corresponding voxel grid; retrieving the historical scanning frequency of the voxel grid within a preset period to determine the dynamic threshold at the current moment; retrieving the mean vector and covariance matrix that match the voxel grid from the background distribution model; using the mean vector and covariance matrix to perform multi-dimensional spatial operations on the real-time point cloud coordinates to obtain the Mahalanobis distance; if the Mahalanobis distance is greater than the dynamic threshold, then the real-time point cloud coordinates are written into the abnormal candidate point set, otherwise the real-time point cloud coordinates are discarded.

[0116] In this embodiment, before the real-time point cloud enters the anomaly candidate filtering, voxel mapping, dynamic threshold generation, background distribution parameter retrieval, and Mahalanobis distance calculation are completed. Then, the value comparison between the Mahalanobis distance and the dynamic threshold determines whether the point enters the anomaly candidate point set. The voxel grid side length, as a spatial discrete scale parameter, is determined by the voxel partitioning step and registered as a numerical parameter during the deployment phase. The minimum and maximum boundary coordinates of the monitoring space, as spatial range parameters, are obtained by scene mapping and installation calibration and registered as numerical parameters during the deployment phase. The preset period, as the scanning frequency statistical window length parameter, is jointly determined by the radar frame interval and the background change rate. The background change rate is obtained by the time normalization result of the mean vector drift and covariance matrix change under the background state. The preset period is the shortest time length that can cover the stable segment of background change and is registered with a unified time source time unit.

[0117] Furthermore, when acquiring the real-time point cloud coordinates collected by the non-repeating scanning LiDAR and mapping them to the corresponding voxel mesh, coordinate system unification is first performed point by point for each frame of point cloud. The coordinate transformation parameters are derived from the installation calibration, which is obtained by measuring the spatial relationship between the radar installation position, installation attitude and scene reference coordinates. The transformation parameters are registered in the deployment phase and directly called in the runtime phase.

[0118] Specifically, after coordinate transformation, boundary determination is performed for each point. The boundary determination parameters are the minimum and maximum boundary coordinates of the monitoring space. The determination process is a component-by-component comparison: if any coordinate component is less than the minimum boundary coordinate in the corresponding direction, or if any coordinate component is greater than the maximum boundary coordinate in the corresponding direction, then the point is removed from the subsequent operation chain of this step. For points that pass the boundary determination, a voxel index is calculated. The calculation process is as follows: subtract the minimum boundary coordinate in the corresponding direction from the coordinate component of the point in each direction to obtain the relative displacement; remove the relative displacement and obtain the original value of the interval index by the voxel grid side length; take the integer part of the original value of the interval index to obtain the interval index in that direction; combine the three interval indices in each direction according to the encoding rules registered during system initialization to form the voxel index, thereby determining the voxel grid position to which the point belongs.

[0119] Furthermore, when retrieving the historical scan frequency of the voxel grid within a preset period to determine the dynamic threshold at the current moment, a sliding window counter is first maintained for each voxel dimension. The sliding window counter contains three pieces of information: the window start timestamp, the number of valid observation frames within the window, and the voxel index. The timestamp is provided by a unified time source and is timed in the time unit registered during the deployment phase. The counting rule for the number of valid observation frames within the window is: if at least one point of the voxel passes the boundary judgment and completes the voxel index positioning within a frame, then the count of the voxel is incremented by 1 for that frame.

[0120] The window update employs a sliding mechanism. The sliding step size is obtained from frame interval statistics. Frame interval statistics are obtained by continuously recording the completion timestamps of at least 300 point cloud reception frames, calculating the difference between adjacent timestamps, and averaging the difference sequence. The sliding step size is an integer multiple of this average frame interval and is registered during the deployment phase. During each window update, the window start timestamp is first moved forward to the current timestamp minus a preset period. Then, count contributions falling before the new window start timestamp are removed, based on the timestamp records of each frame. After removal, the number of valid observation frames within the window is obtained. The historical scan frequency is calculated by dividing the number of valid observation frames within the window by the preset period duration. This frequency value serves as the input for dynamic threshold retrieval.

[0121] Furthermore, the dynamic threshold is obtained by mapping historical scan frequencies. The mapping structure is generated and registered as a threshold table during the deployment phase through background calibration. When generating the threshold table, point cloud data covering multiple preset periods is first collected under background conditions. The background state is determined by the background convergence state, which is obtained from the survival state determination during the background modeling phase, and the determination criteria are registered during the deployment phase. For each voxel, the historical scan frequency is counted within each preset period window, and the Mahalanobis distance sequence is calculated for background points falling into that voxel.

[0122] The historical scanning frequencies were then divided into multiple frequency ranges, with the range boundaries determined by the minimum frequency, maximum frequency, and frequency step size: the minimum frequency was the minimum frequency value observed in the background calibration data, the maximum frequency was the maximum frequency value observed, and the frequency step size was an integer multiple of the frequency increment corresponding to the reciprocal of the frame rate, and was registered during the deployment phase. For each frequency range, the background Mahalanobis distance samples falling within that range were aggregated and sorted from smallest to largest. Candidate thresholds were tested one by one, and the over-threshold ratio was calculated. The over-threshold ratio was calculated as the number of samples exceeding the candidate threshold divided by the total number of samples in that range.

[0123] The selection of candidate thresholds is based on target false alarm constraints, which are provided by the operational procedures during the deployment phase. These constraints are given as the maximum allowed number of false alarms within a single decision period. The allowed over-threshold ratio is obtained by dividing the maximum allowed number of false alarms by the upper limit of the number of points within a single decision period. The upper limit of the number of points is determined and registered by the maximum number of points per period under the radar's highest frame rate and maximum echo output configuration. For each frequency range, the smallest candidate threshold that ensures the over-threshold ratio is not greater than the allowed over-threshold ratio is taken as the dynamic threshold for that range and registered in the threshold table. During the operational phase, the frequency range is located based on the current voxel's historical scan frequency, and the dynamic threshold for that range is read, thus obtaining the dynamic threshold for the current moment.

[0124] Furthermore, when retrieving the mean vector and covariance matrix matching the voxel grid from the background distribution model, the voxel index is used as the retrieval key to locate the background distribution model record. The background distribution model record contains the voxel index, three components of the mean vector, nine values ​​in a 3x3 covariance matrix, and a modeling timestamp. After retrieval, a numerical stability check is performed on the covariance matrix. The check parameter is the lower bound of the minimum eigenvalue, which is determined by the coordinate quantization step. The coordinate quantization step is determined by the radar coordinate output resolution and the number of decimal places retained in the coordinate transformation link, and is registered during the deployment phase. The variance magnitude is obtained by multiplying the coordinate quantization step by itself, and this variance magnitude is registered as the lower bound of the minimum eigenvalue.

[0125] The checking process involves finding the eigenvalues ​​of the covariance matrix and comparing the minimum eigenvalue with the lower limit. If the minimum eigenvalue is less than the lower limit, the variance magnitudes are superimposed on the three diagonal positions of the covariance matrix, and the minimum eigenvalue is found again to complete the verification. After the verification is passed, the Mahalanobis distance is calculated.

[0126] Furthermore, when performing multi-dimensional spatial operations on real-time point cloud coordinates to obtain Mahalanobis distance using the mean vector and covariance matrix, the calculation is performed independently for each point. First, a deviation vector is constructed. The calculation process is as follows: the first deviation component is obtained by subtracting the first component of the mean vector from the first coordinate component of the point; the second deviation component is obtained by subtracting the second component of the mean vector from the second coordinate component of the point; and the third deviation component is obtained by subtracting the third component of the mean vector from the third coordinate component of the point. These three deviation components are combined to form the deviation vector. Subsequently, the covariance matrix is ​​inverted. The inversion uses a symmetric matrix numerical decomposition and back-substitution solution process, ensuring numerical accuracy consistent with the accuracy of the stored covariance matrix.

[0127] Furthermore, after obtaining the inverse matrix, the inverse matrix is ​​multiplied by the deviation vector to obtain the weighted deviation vector. Then, the weighted deviation vector is multiplied by the deviation vector to obtain the squared distance. Finally, the square root of the squared distance is performed to obtain the Mahalanobis distance. If the squared distance has a tiny negative value, it is set to 0 according to the numerical stability rule before the square root operation is performed. The threshold for determining tiny negative values ​​is a negative number of the variance order, which is derived from the coordinate quantization step and registered during the deployment phase.

[0128] Furthermore, after setting the dynamic threshold and Mahalanobis distance, a numerical comparison is performed on each point. The comparison rule is that if the Mahalanobis distance is greater than the dynamic threshold, the point enters the anomaly candidate point set; if the Mahalanobis distance is not greater than the dynamic threshold, the point is removed from the anomaly candidate link. The anomaly candidate point set is organized using voxel index bucketing, with the bucket key being the voxel index. The record content entering the anomaly candidate point set includes the point's 3D coordinates, voxel index, point cloud timestamp, and Mahalanobis distance value. The timestamp is provided by a unified time source and kept consistent with the point cloud frame time. After recording, a consistency check is performed on the anomaly candidate point set. The consistency check includes voxel index validity check and timestamp monotonicity check. After passing the check, the anomaly candidate point set is used as input for subsequent nearest neighbor search and anomaly scoring.

[0129] Example 8: S7 includes extracting the coordinates of a single candidate point from the set of abnormal candidate points and performing a radius search in three-dimensional space with this as the center to obtain a local neighborhood point cloud set; inputting the local neighborhood point cloud set into the spatial moment calculation module, performing an arithmetic mean operation using the three-dimensional coordinate components of each point in the local neighborhood point cloud set to determine the real-time mean vector; performing an outer product operation on each coordinate point in the local neighborhood point cloud set after subtracting the real-time mean vector and accumulating the results, dividing the accumulated outer product matrix by the total number of points in the local neighborhood point cloud set to obtain the real-time covariance matrix; and performing a multidimensional Gaussian distribution modeling on the real-time mean vector and the real-time covariance matrix to obtain the real-time distribution model.

[0130] In this embodiment, after the set of abnormal candidate points is formed, a definite traversal order is established for the candidate points. Then, a local neighborhood point cloud set is constructed for each candidate point, and the real-time mean vector and real-time covariance matrix are calculated based on this set. Finally, the real-time distribution model is constructed from the real-time mean vector and the real-time covariance matrix. The candidate point traversal order is set as a dual-key order during the system deployment phase. The first key is the point cloud timestamp in ascending order, with the timestamps provided by a unified time source and using the same time unit. The second key is the voxel index in ascending order, with the voxel index obtained by interval assignment from the voxel grid side length and the monitoring space boundary and encoded as an integer.

[0131] After extracting the coordinates of a single candidate point in the above order, the coordinates of the candidate point are used as the center point of the neighborhood search. The three coordinate components of the center point are directly taken as the coordinate values ​​of the candidate point in the scene coordinate system. The scene coordinate system is determined by the installation calibration. The installation calibration parameters are measured by the relationship between the radar installation position, installation attitude and scene reference point and are registered during the deployment phase.

[0132] Furthermore, when performing a radius search in 3D space using the center point, the search radius is first determined. During the deployment phase, the search radius is determined jointly by the voxel grid edge length and the background point cloud density. The voxel grid edge length is derived from the voxel partitioning step and recorded in length units. The background point cloud density is statistically obtained from the background state. The statistical process involves counting the number of points within each voxel frame by frame during the background convergence period, accumulating a sequence of point counts within each voxel, and then calculating the arithmetic mean of this sequence to obtain the average number of points.

[0133] The search radius is calculated by multiplying the voxel grid edge length by 0.75 to obtain the first radius value. Simultaneously, a second radius value is calculated, derived from the average nearest neighbor distance in the background state. The average nearest neighbor distance is calculated by finding the point with the smallest distance within each voxel, calculating the 3D Euclidean distance between the two points and summing them into a distance sequence, then taking the arithmetic mean of this distance sequence to obtain the average nearest neighbor distance. The second radius value is then multiplied by 3. The final search radius is the larger of the first and second radius values, ensuring that the neighborhood coverage scale simultaneously satisfies the voxel spatial scale and point cloud density constraints. After determining the search radius, the square of the search radius is calculated by multiplying the search radius by itself.

[0134] Furthermore, the radius search is implemented by first establishing a spatial retrieval structure. This structure is generated after each frame of point cloud data is received. The generation process involves binning and merging point cloud coordinates according to voxel indices, and maintaining a continuously stored list of point coordinates for each voxel bin. During candidate point search, the voxel index of the center point is first located, and simultaneously, the indices of the 26 adjacent voxels are located. These 26 adjacent voxel indices are obtained by adding or subtracting 1 from the voxel index in three directions, respectively. Indexes exceeding the boundary are directly discarded.

[0135] Furthermore, the coordinates of points within the voxel to which the center point belongs and its adjacent voxels are then merged into a candidate search point set. For each candidate search point in the set, the squared 3D distance to the center point is calculated: the squares of the differences between the first and second coordinate components of the point and the center point, the squares of the differences between the third and fourth coordinate components of the point and the center point, and the squares of the differences between the third and fourth coordinate components of the point and the center point are calculated. These three squared results are then summed to obtain the squared distance. If the squared distance does not exceed the squared search radius, the point is included in the local neighborhood point cloud set. After the local neighborhood point cloud set is established, the total number of points in the set is recorded. The total number of points is represented as an integer and used as the divisor for subsequent mean and covariance normalization.

[0136] Specifically, to avoid degradation of the real-time covariance matrix due to insufficient samples, a minimum neighborhood point threshold is set. This threshold is determined during the deployment phase based on background point cloud statistics. The determination process involves performing the same radius search on each voxel under background convergence conditions, counting the number of neighborhood points obtained in each search to form a neighborhood point count sequence, and then calculating the arithmetic mean of this sequence to obtain the average neighborhood point count. The average neighborhood point count is multiplied by 0.2 and rounded up to obtain the first threshold value. A second threshold value of 12 is also set. The minimum neighborhood point threshold is the larger of the first and second threshold values, thus ensuring sufficient sample support for the real-time covariance matrix calculation.

[0137] If the total number of points in the local neighborhood point cloud set is less than the minimum neighborhood point threshold, the candidate point will not enter the real-time distribution model calculation link, and an "insufficient sample" flag will be registered in the candidate point management structure. The registration content includes the voxel index and timestamp.

[0138] As a preferred embodiment, after the local neighborhood point cloud set meets the minimum neighborhood point number threshold, the set is input into the spatial moment calculation module and the real-time mean vector is calculated.

[0139] Specifically, the real-time mean vector consists of three components. The calculation process involves summing the first coordinate component for each point in the local neighborhood point cloud set to obtain the first coordinate sum, and then dividing the first coordinate sum by the total number of points in the set to obtain the first component of the real-time mean vector. Similarly, the second coordinate component is summed and divided by the total number of points in the set to obtain the second component. The third coordinate component is summed and divided by the total number of points in the set to obtain the third component. To ensure numerical consistency, the coordinate values ​​are quantized using a step-wise precision. The coordinate quantization step is determined by the radar output resolution and the number of decimal places retained by the coordinate transformation link, and is registered during the deployment phase. The sum is kept to have the same number of decimal places during the calculation process to avoid the propagation of accumulation errors.

[0140] As a preferred implementation, the calculation of the real-time covariance matrix is ​​performed by centrally accumulating the data around the real-time mean vector. The real-time covariance matrix adopts a 3x3 structure, and the cumulative values ​​for each of the nine matrix positions are first established and set to 0. Then, for each point in the local neighborhood point cloud set, subtract the corresponding component of the real-time mean vector from each of the three coordinate components of that point to obtain the first deviation component, the second deviation component, and the third deviation component. Based on the three deviation components, construct an outer product term and accumulate it to 9 positions: the product of the first deviation component and the first deviation component is accumulated to the first row and the first column; the product of the first deviation component and the second deviation component is accumulated to the first row and the second column; the product of the first deviation component and the third deviation component is accumulated to the first row and the third column; the product of the second deviation component and the first deviation component is accumulated to the second row and the first column; the product of the second deviation component and the second deviation component is accumulated to the second row and the second column; the product of the second deviation component and the third deviation component is accumulated to the second row and the third column; the product of the third deviation component and the first deviation component is accumulated to the third row and the first column; the product of the third deviation component and the second deviation component is accumulated to the third row and the second column; the product of the third deviation component and the third deviation component is accumulated to the third row and the third column. After accumulating all points, divide the nine accumulated values ​​by the total number of points in the local neighborhood point cloud set to obtain the nine positional values ​​of the real-time covariance matrix.

[0141] Furthermore, after the real-time covariance matrix is ​​generated, symmetry consistency verification and numerical stabilization processing are performed. The tolerance threshold for symmetry consistency verification is taken as the variance level, which is calculated by the coordinate quantization step. The calculation process involves multiplying the coordinate quantization step by itself to obtain a value, which is then registered as the variance level threshold during the deployment phase. The verification process involves calculating the absolute value of the difference between the first row and second column and the second row and first column, the absolute value of the difference between the first row and third column and the third row and first column, and the absolute value of the difference between the second row and third column and the third row and second column. If the absolute value of any difference exceeds the variance level threshold, the symmetrical position is replaced with the arithmetic mean of the two values, thereby restoring the symmetrical structure.

[0142] Specifically, the numerical stabilization process uses the minimum variance lower limit threshold, which is also a threshold of the variance order of magnitude. The process involves checking the values ​​at the three diagonal positions of the real-time covariance matrix. If any diagonal position is less than the minimum variance lower limit threshold, that diagonal position is replaced with the minimum variance lower limit threshold.

[0143] Furthermore, the real-time distribution model is composed of a real-time mean vector and a real-time covariance matrix, organized according to a multidimensional Gaussian distribution modeling approach. The modeling process involves using the real-time mean vector as the distribution center parameter and the real-time covariance matrix as the distribution shape parameter, establishing associations between these two parameters and the candidate point voxel indices and candidate point timestamps. The association key is a combination of the voxel index and the timestamp, ensuring that the real-time distribution model for the same voxel within the same time segment has a unique identifier. After the real-time distribution model is generated, it enters the subsequent distribution difference calculation chain for statistical difference comparison with the background distribution model of the same voxel.

[0144] Example 9: S8 includes obtaining the background mean vector and background covariance matrix of the background distribution model from the preset storage space; substituting the background mean vector, background covariance matrix, real-time mean vector, and real-time covariance matrix into the relative entropy analytical formula of the multidimensional Gaussian distribution; performing trace operation and determinant logarithmic operation on the background covariance matrix and real-time covariance matrix using the relative entropy analytical formula to obtain the preliminary divergence value; using the preliminary divergence value combined with the difference term between the background mean vector and the real-time mean vector to perform a quadratic form weighted operation to determine the final Körbeklebler divergence; and performing normalization mapping processing on the Körbeklebler divergence to obtain an anomaly score reflecting the degree of deviation of the local point cloud distribution.

[0145] In this embodiment, two sets of statistical parameters corresponding to the same voxel index are used as inputs. One set is the background mean vector and background covariance matrix of the background distribution model, and the other set is the real-time mean vector and real-time covariance matrix of the real-time distribution model. The Körbeklebler divergence is obtained through the analytical calculation process of the relative entropy of the multidimensional Gaussian distribution, and the divergence is converted into anomaly scores through normalization mapping.

[0146] Among them, the voxel index is generated by the voxel grid mapping process as the retrieval key. The index encoding rules are registered during the system initialization phase to ensure that the same spatial voxel corresponds to the same record in the storage space and the running memory, so that the background statistics and real-time statistics can be compared under the same coordinate reference.

[0147] Furthermore, when obtaining the background mean vector and background covariance matrix from the preset storage space, the background distribution model record is first located based on the voxel index to which the candidate point belongs.

[0148] Specifically, the preset storage space is configured during system initialization. This configuration includes storage location identifiers, record field order, single record length, and index structure. The index structure uses voxel indices as keys, with each key-value pair corresponding to a record location. After locating a record, the three components of the background mean vector and the nine positions of the background covariance matrix are read sequentially. Simultaneously, the modeling timestamp of the record is read for consistency verification. The consistency verification rule is that the modeling timestamp must be no later than the current point cloud timestamp, and the voxel index must be within the voxel grid index range. If the verification fails, the voxel is not included in the divergence calculation chain for this step, avoiding misjudgments based on invalid background parameters. After reading, the background mean vector and background covariance matrix are loaded into the runtime memory as input for subsequent matrix decomposition, inversion, and determinant operations.

[0149] Furthermore, before performing relative entropy analysis, numerical stabilization is performed on the background covariance matrix and the real-time covariance matrix. The goal of this process is to ensure that the covariance matrix has a symmetric structure and is invertible. The key parameter for numerical stabilization is the minimum variance lower bound threshold, which is determined during the deployment phase based on the coordinate quantization step. The coordinate quantization step is jointly determined by the radar output resolution and the number of decimal places retained in the coordinate transformation link.

[0150] Specifically, the threshold calculation process involves multiplying the coordinate quantization step value by itself to obtain the variance magnitude, and then registering this variance magnitude as the minimum variance lower limit threshold. The tolerance threshold for symmetry consistency verification also uses the aforementioned variance magnitude. The verification process involves calculating the absolute values ​​of the differences between the three sets of symmetric positions in the covariance matrix. If the absolute value of any difference exceeds the tolerance threshold, the symmetric position is replaced with the arithmetic mean of the two values, restoring the matrix to its symmetric structure. Reversibility processing is performed after the symmetry consistency verification. The processing involves checking the values ​​at the three diagonal positions of the covariance matrix one by one. If any diagonal position is less than the minimum variance lower limit threshold, the diagonal position is replaced with the minimum variance lower limit threshold, and the symmetry consistency verification is performed again after the replacement.

[0151] Furthermore, when inputting the background mean vector, background covariance matrix, real-time mean vector, and real-time covariance matrix into the relative entropy analytical calculation process, the covariance term is first calculated to obtain the preliminary divergence value. The covariance term includes two parts: trace operation and determinant logarithm operation.

[0152] Specifically, the trace operation first calculates the inverse of the background covariance matrix. The inverse calculation employs symmetric matrix numerical decomposition and back-substitution: in the decomposition stage, the background covariance matrix is ​​decomposed into an intermediate structure that facilitates solving; in the back-substitution stage, each column of the unit basis vectors is solved to obtain the columns of the inverse matrix. The numerical precision of the inverse matrix is ​​consistent with the precision of the covariance matrix. After obtaining the inverse matrix, matrix multiplication is performed. The multiplication process is calculated item by item in a 3x3 matrix: the first row and first column of the product matrix are obtained by multiplying the three elements of the first row of the inverse matrix by the three elements of the first column of the real-time covariance matrix and summing them; the first row and second column are obtained by multiplying the three elements of the first row of the inverse matrix by the three elements of the second column of the real-time covariance matrix and summing them; the first row and third column are calculated in the same way with the third column of the real-time covariance matrix; the second and third rows are combined in corresponding rows and columns, and the above multiplication and addition process is repeated to obtain the complete product matrix. The trace operation result is obtained by summing the values ​​at the three diagonal positions of the product matrix. The summation process is to add the value at the first position to the value at the second position, and then add the value at the third position to obtain the trace operation result.

[0153] As a preferred embodiment, the determinant logarithm operation calculates the natural logarithm of the determinant for both the background covariance matrix and the real-time covariance matrix.

[0154] Specifically, to ensure numerical stability, the determinant calculation employs numerical decomposition to obtain a triangular structure: first, the covariance matrix is ​​decomposed to obtain an upper or lower triangular structure, and then the determinant is calculated using the diagonal factors of the decomposition result; the specific calculation process involves multiplying the values ​​of the three diagonal factors sequentially to obtain the determinant value. Subsequently, when calculating the natural logarithm, the natural logarithms of the three diagonal factors are taken and summed. The summation process involves adding the first logarithm to the second logarithm, and then adding it to the third logarithm to obtain the natural logarithm of the determinant. After obtaining the logarithms of the background covariance matrix and the real-time covariance matrix, the logarithmic difference is calculated according to the direction of the relative entropy analysis process. The difference is calculated as the background logarithm minus the real-time logarithm, yielding the logarithmic term of the determinant.

[0155] Specifically, the initial divergence value is composed of the trace operation result, the logarithm of the determinant, and the dimension compensation term. The dimension compensation term takes a dimension of 3, which is consistent with the three-dimensional coordinate space and is registered as a constant during the system initialization phase. The calculation process of the initial divergence value is to subtract the dimension value of 3 from the trace operation result, and then add the logarithm of the determinant to the above result to obtain the initial divergence value; this calculation process uses the same numerical precision as the covariance matrix throughout to avoid bias introduced by mixed precision.

[0156] Furthermore, after obtaining the initial divergence value, the mean difference term is calculated and a quadratic weighted operation is performed to obtain the Kuhlberg-Klebler divergence. First, the mean difference vector is calculated. The three components of the mean difference vector are the corresponding components of the real-time mean vector minus the corresponding components of the background mean vector. The subtraction is performed component by component, preserving the system configuration precision. Then, the weighted difference vector is calculated. The weights are based on the background covariance inverse matrix, which has already been obtained and reused in the trace operation. The calculation of the weighted difference vector is performed row by row using matrix-vector multiplication: the first component is obtained by multiplying the three elements of the first row of the inverse matrix by the three components of the mean difference vector and summing them; the second component is obtained by multiplying the second row of the inverse matrix by the mean difference vector in the same way; and the third component is obtained by multiplying the third row of the inverse matrix by the mean difference vector in the same way. The quadratic form value is obtained by the inner product of the mean difference vector and the weighted difference vector. The inner product is calculated by multiplying the first component of the mean difference vector by the first component of the weighted difference vector, adding the product to the product of the second component, and then adding the product to the product of the third component to obtain the quadratic form value. The Kuhlberg-Klebler divergence is calculated by adding the initial divergence value to the quadratic form value and then multiplying by 0.5 to obtain the final divergence value. The multiplication by 0.5 is achieved by dividing by 2. The divergence nonnegativity check is performed after the final divergence value is generated. The check threshold is a negative variance magnitude, which is calculated by coordinate quantization steps and registered during the deployment phase. When the final divergence value is less than this negative threshold, the divergence value is replaced with 0 to avoid negative values ​​caused by rounding, which would corrupt the divergence semantics.

[0157] Furthermore, when performing normalization mapping on the Kuhlberg-Klebler divergence to obtain outlier scores, the mapping parameters are first determined and online mapping is completed. The mapping parameters include three items: lower bound of divergence, upper bound of divergence, and upper bound of score. The lower bound of divergence is set to 0 and recorded as a constant, corresponding to a distribution without deviation.

[0158] Specifically, the upper bound of divergence is determined during the deployment phase based on background calibration data. The determination process involves acquiring multiple continuous point cloud segments under background convergence conditions, constructing a local neighborhood for each point cloud segment by voxel, and calculating a divergence sequence. The divergence sequences are then aggregated and sorted in ascending order of value. The sorting sequence position is derived from the target false alarm constraint, which is provided by the operating procedures as the maximum allowed number of false alarms within a single decision period. The allowable overshoot ratio is obtained by dividing the maximum allowed number of false alarms by the upper limit of the sample size for a single decision period. The upper limit of the sample size is determined by the maximum number of points per period under the radar's highest frame rate and maximum echo output configuration and is registered during the deployment phase. The sorting sequence position is obtained by multiplying the total number of aggregated samples by 1, subtracting the allowable overshoot ratio, and rounding down. The result is then incremented by 1 to obtain a valid sequence number. The divergence value corresponding to this sequence number is then used as the upper bound of divergence and registered. The upper bound of the score is set to 100 and registered as a constant during system initialization.

[0159] In the online mapping process, the divergence value of each candidate point is first compared with the upper bound of divergence. When the divergence value does not exceed the upper bound, it is linearly mapped to 0 to 100 according to the ratio of the divergence value to the upper bound. The linear mapping is achieved by dividing the divergence value by the upper bound to obtain the ratio value, and then multiplying the ratio value by 100 to obtain the score value. When the divergence value exceeds the upper bound, the score value is set to 100 to avoid extreme divergence compressing the dynamic range of the score. After the score is generated, it is associated with the voxel index and point cloud timestamp. The association key uses a combination of voxel index and timestamp to ensure that subsequent threshold determination and spatial aggregation processing can accurately trace back to the corresponding local neighborhood and the corresponding background model record.

[0160] Example 10: S9 includes acquiring raw point cloud data from a lidar sensor and extracting candidate point sets; performing multi-dimensional feature fusion of spatial coordinates and reflection intensity on the candidate point sets to obtain a local point cloud distribution model; comparing the probability density of the local point cloud distribution model with a preset background texture model to determine the real-time Körbeklebler divergence; if the real-time Körbeklebler divergence is greater than a preset divergence threshold, the local area where the candidate point set is located is determined to be a foreign object area; if the real-time Körbeklebler divergence is less than or equal to the divergence threshold, the local area where the candidate point set is located is determined to be background texture.

[0161] In this embodiment, the real-time determination link takes the raw point cloud data output by the lidar sensor as input. The raw point cloud data includes the three-dimensional coordinates and reflection intensity of each echo point. The three-dimensional coordinates undergo coordinate system unification processing before entering this step. The transformation parameters used for coordinate system unification are derived from the installation calibration. The installation calibration is obtained by measuring the spatial relationship between the radar installation position, installation attitude, and scene reference point, and is registered as numerical parameters during the deployment phase. The quantization step of the reflection intensity is determined by the sensor output format and registered as an intensity accuracy parameter.

[0162] Furthermore, after the raw point cloud enters the processing module, it first performs basic screening. Basic screening uses the minimum and maximum boundary coordinates of the monitoring space as spatial constraint parameters. These boundary parameters are obtained and registered by scene mapping and deployment calibration. For each point, the coordinate components are compared with the boundary parameters one by one. Points whose coordinate components fall outside the boundary are not included in the candidate extraction link.

[0163] Furthermore, the effective reflection intensity lower limit threshold is then selected. The effective reflection intensity lower limit threshold is obtained from the background state intensity noise statistics during the deployment phase: in the monitoring space, a direction without physical reflection is selected, and no less than 300 frames of intensity sequence are continuously collected. First, the arithmetic mean of the intensity sequence is calculated, then the square of the difference between each intensity value and the average value is calculated and averaged. Then, the square root of the squared average value is taken to obtain the fluctuation amount. Finally, the average value is added to 3 times the fluctuation amount to obtain the noise upper limit. The noise upper limit is increased by 1 intensity quantization step to obtain the effective reflection intensity lower limit threshold. During the operation phase, the reflection intensity of each point is compared with this threshold. Points with reflection intensity lower than the threshold are not included in the candidate extraction link, while points with reflection intensity not lower than the threshold are included in the candidate extraction link.

[0164] Furthermore, the formation of the candidate point set in this step is completed using a local aggregation approach after voxel assignment: points entering the candidate extraction link are mapped to voxel grids to obtain voxel indices. The voxel grid side length is given and registered as a spatial discrete scale parameter by the voxel partitioning step. The voxel index calculation process is to subtract the minimum boundary coordinate from the point coordinate components to obtain the relative displacement, remove the relative displacement and take the integer part of the voxel grid side length to obtain the interval number of each direction, and then combine the three interval numbers of each direction according to the encoding rules registered in the initialization stage to form the voxel index. Subsequently, the points are merged according to the voxel index to form a candidate point set in units of voxels.

[0165] Furthermore, when performing multi-dimensional feature fusion of spatial coordinates and reflection intensity for candidate point sets, the coordinates and intensity are first scaled uniformly. The coordinate normalization scale used for scale unification is the voxel grid side length, and this parameter is registered during the deployment phase. The intensity normalization scale is the background intensity dynamic range, which is calibrated according to the voxel index during the deployment phase. Point clouds are collected under background convergence conditions, and all effective intensity values ​​within a voxel are summarized after voxel merging. The minimum intensity value and the maximum intensity value are taken and subtracted to obtain the dynamic range. The minimum intensity value and the dynamic range are registered as the intensity normalization parameters of that voxel.

[0166] Specifically, the fusion features are generated point-by-point for each point in the candidate point set. The generation process involves first retrieving the background mean vector of the voxel to which the point belongs. The background mean vector is provided by the background distribution model and retrieved through voxel indexing. Then, the corresponding components of the background mean vector are subtracted from the three coordinate components of the point, and then divided by the coordinate normalization scale to obtain three coordinate normalized components. Next, the reflection intensity of the point is subtracted from the minimum background intensity value of the voxel, and then divided by the intensity normalization scale of the voxel to obtain one intensity normalized component. The three coordinate normalized components and the one intensity normalized component are combined according to the splicing order registered in the initialization phase to form a fusion feature record, thus obtaining the fusion feature set of the candidate point set. To avoid numerical amplification caused by an excessively small intensity normalization scale, a lower threshold for intensity normalization is set. During the deployment phase, the lower threshold is set to the value obtained by multiplying one intensity quantization step by 10. When the dynamic range of a voxel is less than this threshold, the intensity normalization scale is replaced with this threshold to keep the normalized intensity components within a controllable numerical range.

[0167] Furthermore, the local point cloud distribution model is obtained by statistical analysis of the fusion feature set of the candidate point set. Its parameters include the real-time mean vector and the real-time covariance matrix. The real-time mean vector contains 4 components, and the real-time covariance matrix has a 4-row, 4-column structure.

[0168] Specifically, the calculation of the real-time mean vector is performed within the candidate point set: first, the number of candidate points is counted as the total number of points; the first component of the fused feature is summed point by point and divided by the total number of points to obtain the first mean component; the second, third, and fourth components are summed and divided repeatedly to obtain the remaining mean components; the four mean components are combined to form the real-time mean vector. The calculation of the real-time covariance matrix is ​​performed by centralizing and accumulating the real-time mean vector: 16 cumulative values ​​are established for each of the 4 rows and 4 columns and set to 0; each fused feature in the candidate point set is traversed, and the corresponding components of the real-time mean vector are subtracted from the four components of the fused feature to obtain four deviation components; the four deviation components are multiplied pairwise to form 16 product terms, and the 16 product terms are accumulated to the corresponding cumulative value positions; after traversing all fused features, the 16 cumulative values ​​are divided by the total number of points to obtain the 16 position values ​​of the real-time covariance matrix.

[0169] After the covariance matrix is ​​generated, a symmetry consistency check and a minimum variance lower bound are performed. The check tolerance threshold and the minimum variance lower bound threshold use the same variance magnitude parameter. During the deployment phase, the variance magnitude parameter is jointly determined by the squared magnitude of the coordinate quantization step and the squared magnitude of the intensity quantization step: the coordinate variance magnitude is obtained by multiplying the coordinate quantization step by itself, and the intensity variance magnitude is obtained by multiplying the intensity quantization step by itself. The larger of the two values ​​is taken as the variance magnitude parameter and recorded. The symmetry consistency check process involves calculating the absolute value of the difference between each pair of symmetrical positions in the 4x4 matrix. If any absolute value of the difference is greater than the variance magnitude parameter, the symmetrical position is replaced with the arithmetic mean of the two values ​​to restore the symmetry structure. The minimum variance lower bound process involves checking the values ​​of the four diagonal positions of the real-time covariance matrix one by one. If any diagonal position is less than the variance magnitude parameter, the position is replaced with the variance magnitude parameter to avoid ill-conditioned values ​​in subsequent probability density calculations. After the above processing is completed, the real-time mean vector and the real-time covariance matrix together constitute the local point cloud distribution model.

[0170] Furthermore, the background texture model is used to compare the probability density with the local point cloud distribution model. During the deployment phase, the background texture model is established by voxel index and stored in a preset storage space, with the storage index key being the voxel index.

[0171] Specifically, the background texture model is generated using the same fusion feature caliber as the real-time side. Its parameters include the background mean vector and the background covariance matrix. The generation process involves acquiring point clouds in the background convergence state, merging the point clouds for each voxel, forming a background fusion feature set according to the same 4-component fusion rule mentioned above, and then calculating the 4-component background mean vector and the 4-row 4-column background covariance matrix respectively. Symmetry consistency check and minimum variance lower bound processing are then performed. The variance magnitude parameters used are the same as those registered in the deployment phase, so that the background texture model and the local distribution model are comparable on the same numerical scale.

[0172] Furthermore, when determining the real-time Körbeklebler divergence by comparing the probability density of the local point cloud distribution model and the background texture model, the covariance matrices of the two sets are inverted and the determinant correlation is calculated respectively.

[0173] Specifically, the inverse matrix is ​​solved using a symmetric matrix numerical decomposition and back-substitution process: first, the covariance matrix is ​​numerically decomposed to obtain an intermediate structure; then, the unit basis vectors are substituted back column by column to obtain each column of the inverse matrix, ensuring that the numerical precision of the inverse matrix is ​​consistent with the precision of the covariance matrix. The calculation of the determinant-related quantities uses the diagonal factor path after numerical decomposition: the diagonal factors of the decomposed structure are multiplied term by term to obtain the determinant value; then, the natural logarithm of each diagonal factor is taken and summed to obtain the natural logarithm of the determinant. This diagonal factor logarithmic summation method avoids numerical underflow caused by directly taking the logarithm of the product.

[0174] Furthermore, after completing the above preparations, the log probability density value of each fused feature in the candidate point set is calculated under the local point cloud distribution model and the background texture model. The calculation process is as follows: first, the model mean vector is subtracted from the fused feature to obtain the difference vector; then, the difference vector is weighted using the model inverse matrix to obtain the weighted difference vector; then, the inner product of the difference vector and the weighted difference vector is performed to obtain the weighted distance; subsequently, the weighted distance, the natural logarithm of the determinant, and the dimension constant 4 are combined to form the log probability density value. The combination process is performed on both models using the same caliber, thus obtaining two log probability density values ​​for each fused feature under the two models. The log density difference is calculated for each fused feature, which is the local model log probability density value minus the background model log probability density value; the arithmetic mean of all log density differences in the candidate point set is calculated by summing the values ​​and dividing by the number of candidate points to obtain the real-time Kuhlberg-Klebler divergence. The divergence results are subjected to nonnegativity checks, and the nonnegativity tolerance threshold is the negative of the variance order of magnitude parameter; when the divergence is less than the nonnegativity tolerance threshold, the divergence is replaced with 0.

[0175] The divergence threshold is used to map real-time Körbeklebler divergence to the determination result of foreign object regions or background texture. The divergence threshold is determined during the deployment phase based on the target false alarm constraint. The target false alarm constraint is given by the operation procedure, and is given in the form of the maximum number of false alarms allowed in a single determination period. The duration of a single determination period is obtained by continuously recording the timestamps of point cloud reception completion for no less than 300 frames, calculating the difference between adjacent timestamps, removing the maximum and minimum differences, and then averaging the results. The upper limit of the sample size in a single determination period is determined and registered by the maximum number of points in a single period under the radar's highest frame rate and maximum echo output configuration.

[0176] Furthermore, the maximum allowed number of false alarms is divided by the upper limit of the sample size to obtain the allowable out-of-bounds ratio. The threshold calibration process involves collecting multiple continuous point cloud segments under background convergence conditions, forming a background candidate point set according to the candidate extraction caliber consistent with real-time data, and obtaining a background divergence sequence according to the aforementioned divergence calculation process. The background divergence sequence is then sorted in ascending order, and the threshold index position is calculated. The threshold index position is calculated by multiplying the total number of samples by 1, subtracting the allowable out-of-bounds ratio, rounding down, and then adding 1. The divergence value corresponding to this index is taken as the divergence threshold and recorded.

[0177] Furthermore, during the runtime phase, after obtaining the real-time Kourbecklebler divergence for each candidate point set, it is numerically compared with a divergence threshold. When the real-time Kourbecklebler divergence is greater than the divergence threshold, the local region to which the candidate point set belongs is determined as a foreign object region; when the real-time Kourbecklebler divergence is not greater than the divergence threshold, the local region to which the candidate point set belongs is determined as background texture. The determination results are associated with voxel indices and point cloud timestamps to map the determination results back to specific spatial locations and specific time segments during the output phase.

[0178] Example 11: Figure 2 This is a structural block diagram of the local terminal of an exemplary electronic device (machine) of the present invention; as shown... Figure 2 As shown, the electronic device of the present invention includes a processor 11, a memory 12, a storage space 13 for storing program code, and program code 14 for executing the method steps according to the present invention. The program code 14 for executing the method steps according to the present invention is used to execute the above-described control logic.

[0179] Figure 3 This is a structural block diagram of the network end of an exemplary electronic device of the present invention; as shown below. Figure 3 As shown, the present invention also provides an electronic device (machine), which may include at least one processor 210, at least one memory 230 communicatively connected to the processor, and a communication bus 240 and a communication interface 220 connecting different system components (including the memory 230 and the processor 210). The processor 210, the memory 230 and the communication interface 220 are connected through the communication bus 240 and communicate with each other. The communication interface 220 is used for data interaction with external devices. The memory 230 stores a machine-executable program that can be executed by the processor, and the processor 210 can execute the above-mentioned control logic by calling the machine-executable program.

[0180] Communication bus 240 represents one or more of several bus architectures, including a memory bus or memory controller, peripheral bus, graphics acceleration port, processor, or local bus using any of the various bus architectures. Examples of these architectures include, but are not limited to, Industry Standard Architecture (ISA) buses, Micro Channel Architecture (MAC) buses, Enhanced ISA buses, Video Electronics Standards Association (VESA) local buses, and Peripheral Component Interconnect (PCI) buses.

[0181] Electronic devices typically include a variety of computer system readable media, which can be any available media that can be accessed by the electronic device, including volatile and non-volatile media, and removable and non-removable media.

[0182] Memory 230 may include computer system readable media in the form of volatile memory, such as random access memory (RAM) and / or cache memory. The electronic device may further include other removable / non-removable, volatile / non-volatile computer system storage media. Memory 230 may include at least one program product having a set (e.g., at least one) of program modules configured to perform the control logic described above.

[0183] A program / utility having a set (at least one) of program modules can be stored in memory 230. Such program modules include, but are not limited to, an operating system, one or more applications, other program modules, and program data. Each or some combination of these examples may include an implementation of a network environment.

[0184] Machine-executable programs for performing this invention can be written in one or more programming languages ​​or a combination thereof. These programming languages ​​include object-oriented programming languages ​​such as Java, C++, and Python, and may also include specialized engineering languages ​​such as R. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network—including a local area network (LAN) or a wide area network (WAN)—or can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0185] The present invention also discloses a storage medium on which a machine-executable program as described above is stored.

[0186] The aforementioned storage medium may be any combination of one or more computer-readable media. Computer-readable media may be, for example, computer-readable signal media or computer-readable storage media. Computer-readable storage media include, but are not limited to, electrical, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatuses, or devices, or any combination thereof. More specific examples of computer-readable storage media (a non-exhaustive list) include: electrical connections having one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM), or flash memory, optical fiber, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof. In this invention, a computer-readable storage medium may be, for example, any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0187] Computer-readable signal media may include data signals propagated in baseband or as part of a carrier wave, carrying computer-readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. Computer-readable signal media may also be any computer-readable medium other than computer-readable storage media, capable of sending, propagating, or transmitting programs for use by or in connection with an instruction execution system, apparatus, or device.

[0188] Program code contained on a computer-readable medium may be transmitted using any suitable medium, including but not limited to wireless, wire, optical fiber, RF, etc., or any suitable combination thereof.

[0189] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for detecting cumulative FOD in ultra-long temporal point clouds based on non-repeating scanning lidar, characterized in that, include: S1. Based on the minimum physical size of the foreign object to be detected and the radar detection range, dynamically set the side length of the spatial voxel grid to ensure that the foreign object occupies at least one voxel. S2. For the side length of the spatial voxel mesh, establish the relationship between the voxel activation ratio and the accumulation time of the ultra-long temporal point cloud, and calculate the minimum accumulation time in reverse according to the preset false alarm rate. S3. If the current ultra-long time domain point cloud accumulation time reaches the minimum accumulation time, then the mean vector of the point cloud in the background voxel is used as the geometric center, and the covariance matrix is ​​calculated as the discreteness and squareness features. The mean vector and covariance matrix are stored as the background distribution model. S4. In the real-time detection phase, obtain the point cloud coordinates generated by the non-repeating scanning lidar, index the voxels to which the point cloud coordinates belong, and read the corresponding mean vector and covariance matrix from the background distribution model. S5. Calculate the Mahalanobis distance from the point cloud coordinates to the background distribution model, and determine whether to write it into the abnormal candidate point set based on the comparison between the Mahalanobis distance and the dynamic threshold. S6. If the Mahalanobis distance is greater than the dynamic threshold, the point cloud coordinates are written into the abnormal candidate point set; otherwise, the point cloud coordinates are discarded. S7. For each point cloud coordinate in the abnormal candidate point set, perform a nearest neighbor search to obtain the local neighborhood point cloud, and calculate the real-time mean vector and real-time covariance matrix to obtain the real-time distribution model. S8. Based on the real-time mean vector, real-time covariance matrix, and the mean vector and covariance matrix of the corresponding background distribution model, calculate the KL divergence as the anomaly score, specifically including: Obtain the background mean vector and background covariance matrix of the background distribution model from the preset storage space; Substitute the background mean vector, background covariance matrix, real-time mean vector, and real-time covariance matrix into the analytical formula for the relative entropy of the multidimensional Gaussian distribution; Preliminary divergence values ​​are obtained by performing trace operations and determinant logarithmic operations on the background covariance matrix and the real-time covariance matrix using the relative entropy analytical formula. The final Körbeklebler divergence is determined by using the preliminary divergence value combined with the difference between the background mean vector and the real-time mean vector through a quadratic weighted operation. Normalization mapping of the Körbeklebler divergence yields anomaly scores that reflect the degree of deviation in the local point cloud distribution. S9. If the KL divergence is greater than a preset threshold, the region corresponding to the abnormal candidate point set is determined to be a foreign object; otherwise, it is determined to be background texture. Specifically, this includes: Acquire raw point cloud data from lidar sensors and extract candidate point sets; A local point cloud distribution model is obtained by fusing multidimensional features of spatial coordinates and reflection intensity for the candidate point set. The real-time Körbeklebler divergence is determined by comparing the probability density of a local point cloud distribution model with a preset background texture model. If the real-time Körbeklebler divergence is greater than the preset divergence threshold, the local region where the candidate point set is located is determined to be a foreign object region. If the real-time Körbeklebler divergence is less than or equal to the divergence threshold, then the local region where the candidate point set is located is determined to be the background texture.

2. The method for detecting cumulative FOD in ultra-long time-domain point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S1 includes: Obtain the minimum physical size of the foreign object and the current detection distance, and extract the basic bounding box side length based on the minimum physical size of the foreign object using geometric features. The beam cross-section diameter is retrieved from the preset beam divergence model using the current detection range. If the side length of the base bounding box is greater than the beam cross-section diameter, the dynamic voxel reference value is obtained by multiplying the side length of the base bounding box and the beam cross-section diameter. Based on the dynamic voxel reference value, a spatial voxel grid is obtained by three-dimensional subdivision in the spatial coordinate system. The radar echo data is then spatially mapped using the spatial voxel grid to determine that the foreign object occupies at least one voxel.

3. The method for detecting cumulative FOD in ultra-long time-domain point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S2 includes: Discrete sampling intervals are defined in a three-dimensional coordinate system by using the side length of the spatial voxel mesh, and the distribution density of the reflection point cloud within the discrete sampling interval is statistically analyzed to calculate the voxel activation ratio. By mapping the voxel activation ratio in a pre-defined probability distribution model, the point cloud persistence state under different reflection intensities is determined, and the relationship between the voxel activation ratio and the ultra-long temporal point cloud accumulation time is established. Based on the preset false alarm rate, retrieve the corresponding confidence interval in the relation and extract the time-domain feature vector corresponding to the confidence interval; The minimum cumulative duration is obtained by multiplying the time-domain feature vector with the sampling frequency, thus realizing the reverse calculation of the minimum cumulative duration based on the preset false alarm rate.

4. The method for detecting FOD (Foreign Object Depth) in ultra-long time-domain point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S3 includes: By comparing the current ultra-long time domain point cloud accumulation time with the preset minimum accumulation time, if the current ultra-long time domain point cloud accumulation time is greater than or equal to the minimum accumulation time, the three-dimensional coordinate set of all discrete sampling points in the voxel grid is extracted. The arithmetic mean of the three-dimensional coordinate set is calculated to obtain the mean vector representing the geometric center of the voxel mesh. The covariance matrix reflecting the spatial distribution of point clouds is determined by performing a sum of squared differences operation on the mean vector and the three-dimensional coordinate set. Based on the eigenvalue decomposition results of the covariance matrix, feature vectors describing the discreteness and squareness of the point cloud are extracted, and the feature vectors are associated and mapped with the mean vector to obtain the background distribution model; By writing the background distribution model into a preset storage space, the geometric center, discreteness, and squareness features of the point cloud within the background voxel are persistently stored.

5. The method for detecting cumulative FOD in ultra-long time-domain point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S4 includes: Obtain the current point cloud coordinate set generated by the non-repeating scanning lidar, and determine the voxel grid position to which the point cloud coordinate set belongs by dividing the coordinate value range; The mean vector and covariance matrix corresponding to the voxel grid positions are extracted from the preset background distribution model as initial distribution parameters; Matrix operations are performed on the mean vector and covariance matrix to obtain the deviation vector set between the point cloud coordinate set and the background distribution; The Mahalanobis distance value is calculated based on the deviation vector group. If the Mahalanobis distance value is less than the preset similarity threshold, the deviation vector group is clustered to determine the subset of the outlier point cloud. A density estimation algorithm is performed on a subset of the outlier point cloud to generate updated background distribution parameters. The updated background distribution parameters are written back into the background distribution model to enable dynamic reading and adjustment of the mean vector and covariance matrix of the voxels to which the point cloud coordinates belong during the real-time detection phase.

6. The method for detecting FOD (Foreign Object Depth) in ultra-long time-domain point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S5 includes: Obtain the real-time point cloud coordinates acquired by the non-repeating scanning lidar and map them to the corresponding voxel grid; Extract the mean vector and covariance matrix of the voxel grid from the preset background distribution model; Mahalanobis distance is obtained by using the mean vector and covariance matrix to measure the multidimensional spatial difference of real-time point cloud coordinates. The cumulative scan frequency of the voxel mesh under historical background conditions is retrieved to determine the dynamic threshold; If the Mahalanobis distance exceeds the dynamic threshold, the corresponding real-time point cloud coordinates are written into the abnormal candidate point set, where the dynamic threshold is positively correlated with the number of times the voxel is scanned in the background state.

7. The method for detecting FOD (Foreign Object Depth) in ultra-long temporal point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S6 includes: Obtain the real-time point cloud coordinates acquired by the non-repeating scanning lidar and map them to the corresponding voxel grid; The historical scan frequency of the voxel grid within a preset period is retrieved to determine the dynamic threshold at the current moment; Retrieve the mean vector and covariance matrix that match the voxel grid from the background distribution model; Mahalanobis distance is obtained by performing multidimensional spatial calculations on real-time point cloud coordinates using the mean vector and covariance matrix. If the Mahalanobis distance is greater than the dynamic threshold, the real-time point cloud coordinates are written into the abnormal candidate point set; otherwise, the real-time point cloud coordinates are discarded.

8. The method for detecting cumulative FOD in ultra-long time-domain point clouds based on non-repeating scanning lidar according to claim 1, characterized in that: S7 includes: Extract the coordinates of a single candidate point from the set of abnormal candidate points and perform a radius search in three-dimensional space with that point as the center to obtain a local neighborhood point cloud set. The local neighborhood point cloud set is input into the spatial moment calculation module, and the arithmetic mean is calculated using the three-dimensional coordinate components of each point in the local neighborhood point cloud set to determine the real-time mean vector. For each coordinate point in the local neighborhood point cloud set, the real-time mean vector is subtracted, and the outer product operation is performed and accumulated. The accumulated outer product matrix is ​​divided by the total number of points in the local neighborhood point cloud set to obtain the real-time covariance matrix. The real-time mean vector and the real-time covariance matrix are modeled using a multidimensional Gaussian distribution to obtain the real-time distribution model.

Images