Multi-sensor fusion method and system driven by non-aligned spatio-temporal data

Abstract

The invention discloses a non-aligned spatio-temporal data driven multi-sensor fusion method and system, and the method comprises the steps: building a matching measurement function through collecting radar point cloud density and visual image texture, and constructing a synchronous gradient field through employing a space depth change parameter and a visual motion contour, and determining a data synchronization standard; detecting a time sequence abnormal accumulation area based on a synchronization standard, and identifying an abnormal starting point and propagation time delay through a backtracking track to form an asynchronous sensor coordinate table; iMU angle increment and sensor response delay are collected, and sampling imbalance information is generated through time domain matching of an attitude switching signal and a data collection inflection point; performing association mapping on the asynchronous coordinate table and the sampling imbalance information to generate a feature offset, determining a matching core area, and creating a hierarchical calibration sequence; performing internal and external source difference identification analysis on the calibration sequence to obtain an offset transfer speed to form an offset coupling factor; and a fusion enhancement region is determined based on the offset coupling factor, and high-precision multi-source data fusion is realized through cascade analysis and gain coefficient coding.

Description

Technical Field

[0001] The present invention relates to the field of data processing technology, and in particular to a multi-sensor fusion method and system driven by non-aligned spatiotemporal data. Background Art

[0002] Drone multi-sensor systems integrate multiple sensors, including lidar, visual cameras, and inertial measurement units, enabling them to acquire rich environmental perception information. However, due to differences in the physical properties of each sensor, deviations in installation position, and changes in the operating environment, data synchronization errors, inconsistent response delays, and sampling frequency offsets are common among sensors. This makes it difficult to accurately align multi-source data in both the temporal and spatial domains.

[0003] Existing sensor calibration methods often rely on offline calibration or single-sensor calibration strategies, lacking in-depth analysis of the dynamic coupling relationships between sensors and unable to effectively handle real-time offset changes in complex flight environments. Traditional methods often overlook the impact of sensor asynchrony on data quality and fail to establish a dynamic model for inter-sensor offset propagation, causing calibration effectiveness to rapidly decay in dynamic environments. Existing calibration methods struggle to ensure multi-source data fusion accuracy and system stability, particularly under complex operating conditions such as high-maneuverability flight, strong interference environments, and long-duration operations. This severely restricts the intelligent perception and autonomous decision-making capabilities of drones. Summary of the Invention

[0004] The present invention discloses a multi-sensor fusion method and system driven by non-aligned spatiotemporal data, aiming to solve the technical problems of low multi-sensor data synchronization accuracy, unstable calibration effect, and poor dynamic adaptability in the existing technology. By constructing data synchronization standards, identifying asynchronous sensor distribution, and analyzing sampling imbalance characteristics, accurate modeling of the offset coupling relationship between sensors is achieved, and a hierarchical calibration sequence and cascade fusion mechanism are established. Ultimately, a high-precision and high-stability multi-source data fusion result is formed, providing reliable perception data support for UAV operations in complex environments.

[0005] The first aspect of the present invention proposes a multi-sensor fusion method driven by non-aligned spatiotemporal data, comprising the following steps: Collecting radar point cloud density and visual image texture, obtaining spatial depth variation parameters based on the point cloud density, capturing visual motion profiles through the visual image texture, and establishing a data synchronization standard through the spatial depth variation parameters and the visual motion profiles; Detecting a timing anomaly cluster area based on the data synchronization standard, determining a data loss diffusion range based on the timing anomaly cluster area, deriving frequency drift node data based on the data loss diffusion range, and forming an asynchronous sensor coordinate table based on the frequency drift node data; collecting IMU angle increments and sensor response delays, detecting an attitude switching signal using the IMU angle increments, capturing a data acquisition inflection point using the sensor response delays, and performing time domain matching between the attitude switching signal and the data acquisition inflection point to generate sampling imbalance information; Performing correlation mapping between the asynchronous sensor coordinate table and the sampling imbalance information to generate a feature offset, establishing a feature matching core area based on the feature offset, and creating a hierarchical calibration sequence based on the feature matching core area; performing internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data, performing extended path analysis on the source offset data to obtain an offset transfer velocity, and fusing the offset transfer velocity and the passive offset data to form an offset coupling factor; A fusion enhancement region is determined based on the offset coupling factor and the asynchronous sensor coordinate table, and a cascade analysis is performed on the fusion enhancement region to form a multi-source data fusion result.

[0006] A second aspect of the present invention provides a multi-sensor fusion system driven by non-aligned spatiotemporal data, comprising: a data acquisition module for collecting radar point cloud density and visual image texture, obtaining spatial depth variation parameters based on the point cloud density, capturing visual motion profiles through the visual image texture, and establishing a data synchronization standard through the spatial depth variation parameters and the visual motion profiles; an asynchronous detection module, configured to detect a timing anomaly cluster area based on the data synchronization standard, determine a data loss diffusion range based on the timing anomaly cluster area, derive frequency drift node data based on the data loss diffusion range, and form an asynchronous sensor coordinate table based on the frequency drift node data; a sampling processing module, configured to collect IMU angle increments and sensor response delays, detect an attitude switching signal using the IMU angle increments, capture a data acquisition inflection point using the sensor response delays, and perform time domain matching between the attitude switching signal and the data acquisition inflection point to generate sampling imbalance information; a feature mapping module, configured to perform correlation mapping between the asynchronous sensor coordinate table and the sampling imbalance information to generate a feature offset, establish a feature matching core region based on the feature offset, and create a hierarchical calibration sequence based on the feature matching core region; an offset analysis module for performing internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data, performing extended path analysis on the source offset data to obtain an offset transfer velocity, and fusing the offset transfer velocity and the passive offset data to form an offset coupling factor; A fusion processing module is used to determine a fusion enhancement area based on the offset coupling factor and the asynchronous sensor coordinate table, and perform cascade analysis on the fusion enhancement area to form a multi-source data fusion result.

[0007] The beneficial effects of the present invention are reflected in the following aspects: 1. The spatial correlation between the density of the radar point cloud and the texture of the visual image is established through the matching metric function, and the adaptive establishment of the dynamic data synchronization standard is realized by combining the synchronous gradient field and contour line analysis. The time series anomaly cluster detection technology is adopted to accurately locate the anomaly starting point through the data backtracking trajectory, and the complete asynchronous sensor space distribution map is constructed in combination with the propagation delay analysis, so as to achieve the improvement of the spatiotemporal synchronization accuracy of the multi-sensor system and the accurate identification of the abnormal state. 2. The collaborative analysis and processing of the IMU angle increment and the sensor response delay are adopted, and the time domain matching technology of the attitude switching signal and the data acquisition inflection point is used to accurately quantify the sampling deviation characteristics using the distribution unevenness of the spectrum jump accumulation. Combined with the vector domain analysis of the characteristic offset and the core area establishment function, a complete dynamic relationship of the offset propagation between sensors is established, and the accurate identification and effective compensation of sampling imbalance under dynamic working conditions are realized. 3. The precise separation of source offset and passive offset is realized through the internal and external source difference recognition technology, and a quantitative description of the offset transmission speed and a mathematical representation of the offset coupling factor are established. A cascade fusion architecture based on the transfer link graph is constructed, and dynamic gain coefficient allocation processing and multi-stage coding fusion operations are adopted to achieve optimal weight distribution and high-precision collaborative fusion of multi-source heterogeneous data, thereby improving the system's data processing capabilities and the reliability of fusion results in complex environments.

[0008] It should be understood that the foregoing general description and the following detailed description are exemplary and explanatory only and are not restrictive of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0009] Figure 1 It is a flow chart of a multi-sensor fusion method driven by non-aligned spatiotemporal data of the present invention.

[0010] Figure 2 This is a structural block diagram of a multi-sensor fusion system driven by non-aligned spatiotemporal data of the present invention. DETAILED DESCRIPTION

[0011] The following will be combined with the drawings in the embodiments of this application to clearly and completely describe the technical solutions in the embodiments of this application. Obviously, the embodiments described are only part of the embodiments of this application, not all of the embodiments. Based on the embodiments of this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0012] It should be noted that all directional indications in the embodiments of the present application (such as up, down, left, right, front, back, etc.) are only used to explain the relative position relationship, movement status, etc. between the various components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication will also change accordingly.

[0013] It should also be noted that when an element is referred to as being "fixed on" or "disposed on" another element, it may be directly on the other element or there may be an intermediate element. When an element is referred to as being "connected to" another element, it may be directly connected to the other element or there may be an intermediate element.

[0014] In addition, the descriptions of "first", "second", etc. in this application are for descriptive purposes only and should not be understood as indicating or implying their relative importance or implicitly indicating the number of the technical features indicated. Therefore, the features defined as "first" or "second" may explicitly or implicitly include at least one of such features. In addition, the technical solutions between the various embodiments can be combined with each other, but this must be based on the fact that they can be implemented by ordinary technicians in this field. When the combination of technical solutions is contradictory or cannot be implemented, it should be deemed that such combination of technical solutions does not exist and is not within the scope of protection required by this application.

[0015] In the following description, specific details such as specific system structures and techniques are provided for purposes of illustration rather than limitation to facilitate a thorough understanding of the embodiments of the present application. However, it will be apparent to those skilled in the art that the present application may be implemented in other embodiments without these specific details. In other cases, detailed descriptions of well-known systems, devices, circuits, and methods are omitted to avoid obscuring the description of the present application with unnecessary detail.

[0016] It should be understood that when used in the present specification and the appended claims, the term "comprising" indicates the presence of described features, integers, steps, operations, elements and / or components, but does not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components and / or collections thereof.

[0017] References to "one embodiment" or "some embodiments" in this specification mean that a particular feature, structure, or characteristic described in conjunction with that embodiment is included in one or more embodiments of the present application. Thus, phrases such as "in one embodiment," "in some embodiments," "in other embodiments," and "in other embodiments" appearing in various places in this specification do not necessarily refer to the same embodiment, but rather mean "one or more but not all embodiments," unless otherwise specifically emphasized. The terms "including," "comprising," "having," and variations thereof all mean "including but not limited to," unless otherwise specifically emphasized.

[0018] The technical solutions of the embodiments of this application are introduced below.

[0019] like Figure 1 As shown, the embodiment of the present invention provides a multi-sensor fusion method driven by non-aligned spatiotemporal data, including the following steps S110 to S160: Step S110 , collecting radar point cloud density and visual image texture, obtaining spatial depth variation parameters based on point cloud density, capturing visual motion profiles through visual image textures, and establishing data synchronization standards through spatial depth variation parameters and visual motion profiles.

[0020] Specifically, radar point cloud density and visual image texture are collected. A LiDAR sensor array deployed on a mobile platform scans the surrounding environment in real time, acquiring point cloud data distribution in three-dimensional space. The LiDAR uses pulsed laser scanning, with an appropriate scanning frequency and angular resolution to ensure effective detection range coverage. Radar point cloud density is quantified by counting the number of valid reflection points within a voxel, which is divided into a cubic grid. Point cloud density calculation uses an adaptive threshold method, dynamically adjusting the density judgment criteria based on ambient lighting and target reflectivity. Visual image texture information is simultaneously collected, using a high-resolution camera to capture visible light images of the environment. The camera is equipped with a zoom lens and has an appropriate image resolution and frame rate. Visual image texture is extracted using the gray-level co-occurrence matrix method, and texture features are calculated, including statistical parameters such as contrast, correlation, entropy, and uniformity. Texture analysis uses a multi-scale windowing technique with incrementally varying window sizes to capture texture patterns at different scales. The collected radar point cloud density and visual image texture are temporally synchronized and spatially registered to form a corresponding multimodal data set.

[0021] Spatial depth variation parameters are obtained based on point cloud density. The collected radar point cloud density data is spatially registered according to a three-dimensional coordinate system to establish a voxel-based representation of the density distribution. By calculating the density gradient between adjacent voxels, regions of spatial depth variation are identified. Spatial depth variation parameters include three main components: depth gradient amplitude, variation direction, and variation rate. The depth gradient amplitude is defined as the modulus of the density difference between adjacent voxels and reflects the intensity of the depth variation. The variation direction is determined by the azimuth angle of the density gradient vector, describing the spatial orientation of the depth variation. The variation rate is calculated by the time derivative of the time-series density data to quantify the dynamic characteristics of the depth variation. Point cloud density is preprocessed using filtering to eliminate noise interference in depth calculation. Interpolation is used to fill in missing density data to ensure the continuity and integrity of the spatial depth variation parameters. The depth variation parameters are normalized using the maximum-minimum normalization method (D_norm = (D - D_min) / (D_max - D_min)) to map the parameter values ​​to the range of 0-1.

[0022] Visual motion contours are captured through visual image textures. The captured visual image texture features are analyzed temporally, and motion regions are detected based on texture differences between consecutive frames. Optical flow methods are used to calculate pixel motion vectors, and an appropriate calculation window size is set. By analyzing the distribution pattern of optical flow vectors, sets of pixels with similar motion characteristics are identified, forming a preliminary segmentation of the motion region. Edge detection is performed on the motion region, and the region boundaries are extracted using edge detection operators to form a geometric description of the visual motion contour. Morphological filtering is applied during contour extraction to remove isolated noise points and fill contour gaps. Visual motion contours are represented using chain codes, recording the basic characteristic parameters of the contour. Contour features are normalized using a combination of length and amplitude normalization to normalize the contour feature vector to unit amplitude. A vectorized representation of the contour features is established, and key feature components are extracted through dimensionality reduction using principal component analysis. A contour quality control threshold is set to filter out invalid contours that are too small or too complex.

[0023] In some embodiments, establishing a data synchronization standard through the spatial depth change parameters and the visual motion profile includes: identifying high-fit areas and low-fit areas through the matching relationship between the spatial depth change parameters and the visual motion profile; coupling analysis of the peak value of the high-fit area and the valley value of the low-fit area to generate a synchronization potential difference; using the synchronization potential difference to construct a synchronization gradient field; and determining the data synchronization standard along the data contour line of the synchronization gradient field.

[0024] High- and low-match regions are identified by comparing the spatial depth variation parameters with the visual motion profile. A matching metric function M(x,y,t) = α·D(x,y,t) + β·C(x,y,t) is established for the spatial depth variation parameters and the visual motion profile, where M is the matching metric, D(x,y,t) is the normalized value of the depth variation parameters, C(x,y,t) is the normalized value of the motion profile features, and α and β are weighting coefficients. The weighting coefficients are determined through correlation analysis of the training data, and the sum of α and β is equal to 1 to ensure reasonable matching. A sliding window technique is used to calculate the local matching degree in the spatiotemporal domain, with the window size adaptively adjusted based on sensor resolution and motion speed. A dual-threshold segmentation method is used to divide the matching degree distribution into different regions, with a high threshold of 0.7 and a low threshold of 0.3. Regions with a matching degree above the high threshold are considered high-match regions, indicating a strong correlation between the spatial depth variation parameters and the visual motion profile. When the matching degree is lower than the low threshold, it is judged as a low-match area, indicating that there is a significant difference or inconsistency between the two data. Connectivity analysis is used to cluster the matching areas, merge spatially adjacent similar areas, and set a minimum area threshold to filter out noisy areas.

[0025] The peaks of high-fit regions and the valleys of low-fit regions are coupled and analyzed to generate synchronization potential differences. The local maximum points of matching degree are extracted from the identified high-fit regions as matching peaks. The local minimum points of matching degree are extracted from the low-fit regions as matching valleys. Each peak is paired with its adjacent valley using the shortest distance criterion. The synchronization potential difference is calculated as the difference between the peak and valley values. This parameter quantifies the degree of synchronization between spatial depth change parameters and visual motion profiles in different regions. The synchronization potential difference is spatially interpolated to construct a continuous potential difference distribution throughout the observation area. The Gaussian interpolation method is used, and the interpolation radius is adaptively adjusted according to the data sampling density. A potential difference change detection mechanism is established. When the potential difference change exceeds the preset threshold, the synchronization standard is recalculated. Outlier detection and correction are performed on the potential difference data to ensure the stability and reliability of the synchronization analysis.

[0026] Use the synchronized potential difference to construct the synchronized gradient field. Based on the spatial distribution of the synchronized potential difference, calculate the spatial gradient of the potential difference ,in is the synchronous gradient, and are the partial derivatives of the potential difference in the x and y directions, respectively. Gradients are calculated using the finite difference method, with the difference step size determined by the data sampling density. The synchronized gradient field describes the direction and intensity of change of the synchronized potential difference in space, with the gradient vector pointing in the direction of the fastest potential increase. A stability assessment mechanism for the gradient field is established, monitoring the field stability by the rate of change of the gradient amplitude. Gradient field smoothing is employed, with a 3×3 Gaussian filter used to eliminate the effects of computational noise. A directional consistency check for the gradient field is established to ensure the rationality of the gradient directions in adjacent regions. Quality control of the gradient field is performed, with a reasonable range for the gradient amplitude set and abnormal gradient values ​​filtered out.

[0027] The data synchronization standard is determined along the data contour lines of the synchronized gradient field. Potential difference contour lines are drawn in the synchronized gradient field, connecting spatial points with the same potential difference value. The key geometric parameters of the contour lines are calculated, including the contour line length L, the average curvature K, and the closedness C. The contour line selection score S_line = w1×(1 / K)+w2×L+w3×C, where w1, w2, and w3 are weight coefficients. The top 20% contour lines with the highest scores are selected as candidate boundaries. The contour lines are divided into three synchronization levels according to the potential difference value: high level corresponds to a potential difference value greater than 0.8, medium level corresponds to 0.4-0.8, and low level corresponds to less than 0.4. Quantitative indicators for determining the data synchronization standard: synchronization accuracy is defined as the average matching error between the spatial depth change parameter and the visual motion contour, synchronization delay is defined as the time required to achieve stable matching, and synchronization stability is defined as the inverse of the matching variance. Standard thresholds are set based on synchronization levels: High level requires synchronization accuracy less than 5%, synchronization latency less than 10ms, and synchronization stability greater than 0.9; Medium level requires accuracy less than 10%, latency less than 20ms, and stability greater than 0.7; Low level requires accuracy less than 20%, latency less than 50ms, and stability greater than 0.5. Ultimately, the data synchronization standards are output as a structured data format that includes contour line coordinates, synchronization level, quantitative indicator thresholds, and applicable scope.

[0028] Step S120 , detecting a timing anomaly cluster area based on a data synchronization standard, determining a data loss diffusion range through the timing anomaly cluster area, deriving frequency drift node data through the data loss diffusion range, and forming an asynchronous sensor coordinate table based on the frequency drift node data.

[0029] Specifically, the method detects clusters of timing anomalies based on data synchronization standards. Using the data synchronization standards as a detection benchmark, the real-time data stream of a drone's multi-sensor system is analyzed for timing anomalies. The data synchronization standards encompass three key parameters: synchronization accuracy, synchronization delay, and synchronization stability. Each parameter corresponds to a specific anomaly determination threshold. The data stream is continuously monitored using a sliding time window to ensure the temporal resolution of anomaly detection. Within each detection window, the deviation between the actual data synchronization level and the synchronization standard is calculated as D(t) = |S_actual(t) - S_standard|, where D(t) is the synchronization deviation at time t, S_actual(t) is the measured synchronization parameter, and S_standard is the standard synchronization parameter. When the deviation values ​​for multiple consecutive detection windows exceed the threshold, a timing anomaly is determined to exist in that time period. Clustering is used to spatially aggregate the detected anomaly time points, identifying anomalies that are close in time and space and forming anomaly clusters. Clustering is performed using the DBSCAN algorithm, with a temporal distance threshold of 5 seconds, a spatial distance threshold of 10 meters, and a minimum number of cluster points of 3. Cluster boundaries were extracted, and the convex hull algorithm was used to determine the minimum bounding box for each cluster. The center coordinates of the clusters were recorded. The anomaly intensity of each cluster was assessed, and the average and maximum anomaly values ​​within the clusters were calculated to quantify the severity of the anomaly. Spatial correlations between clusters and sensor nodes were established, and the sensor range affected by each cluster was determined.

[0030] In some embodiments, determining the data missing diffusion range through the timing anomaly cluster area includes: tracing back forward based on the timing anomaly cluster area to generate a data backtracking trajectory; identifying the starting point of the anomaly according to the data backtracking trajectory; obtaining the propagation delay through the starting point of the anomaly; and moving the starting point of the anomaly forward based on the propagation delay to determine the data missing diffusion range.

[0031] A data backtracking trajectory is generated by tracing back through time series anomaly clusters. Starting from the detected time series anomaly clusters, the historical state of the data stream is analyzed backward along the time axis to track the temporal evolution of the anomaly formation. A backward tracing method is employed, using a sliding window technique to gradually trace back, with a window size of 30 seconds and a backtracking step of 5 seconds. At each backtracking step, changes in the statistical characteristics of the data are analyzed, including key indicators such as mean shift, variance change, and distribution skew. The degree of anomaly is calculated for each backtracking point, and the intensity of the anomaly is quantified by the deviation of the current statistical parameter from the reference value. When the degree of anomaly falls below a set threshold, the moment is marked as a backtracking point. All backtracking points are connected to form a data backtracking trajectory, and the trajectory curve is smoothed using spline interpolation. The trajectory is segmented to identify key nodes where the degree of anomaly changes. Each node corresponds to a significant moment in the anomaly evolution. The resulting data backtracking trajectory fully describes the temporal path of the anomaly's evolution from a normal state to an aggregated state.

[0032] The anomaly onset is identified based on the data backtracking trajectory. Within the generated data backtracking trajectory, the precise moment at which the anomaly began is identified using the change point detection method. A cumulative sum detection method is used to calculate the cumulative anomaly value at each moment along the trajectory, reflecting the cumulative effect of the anomaly signal. A cumulative anomaly detection threshold is set at 0.5; when the cumulative anomaly value first exceeds this threshold, the corresponding moment is identified as the anomaly onset. A bisection method is used to precisely locate the onset, improving detection accuracy to the second level. A refined analysis is performed within a ±10-second time window around the candidate onset location, calculating the gradient rate of change and variance jump amplitude within this window. The reliability of the onset is confirmed using multiple validation methods, including confidence interval estimation and statistical significance testing. The detected onset is classified into three types of anomaly patterns: sudden, gradual, and oscillatory. The precise time coordinates of the anomaly onset are recorded. An uncertainty assessment of the onset is established, and the confidence level of the detection results is calculated using the Monte Carlo method. Spatial localization analysis of the onset is performed to determine its relative position within the sensor network, providing a spatial basis for propagation path analysis.

[0033] The propagation delay is obtained from the anomaly starting point. Using the identified anomaly starting point, the time delay for the anomaly to propagate from the starting point to the cluster is calculated. The propagation delay, T_delay, equals t_detect - t_start, where T_delay is the propagation delay, t_detect is the time the anomaly cluster is detected, and t_start is the time the anomaly starts. Considering the physical topology of the sensor network, the spatial distance from the starting point to the cluster center is calculated. The propagation velocity is estimated by the ratio of distance to delay, reflecting the propagation characteristics of the anomaly in the sensor network. A statistical model for propagation velocity is established, and propagation data from multiple anomaly events is collected to calculate the average propagation velocity and standard deviation. When the calculated propagation velocity exceeds a reasonable range, the historical average is used as a substitute. The propagation velocity is adjusted based on sensor type and environmental conditions to ensure the accuracy of the propagation delay calculation.

[0034] The anomaly's starting point is shifted forward based on the propagation delay to determine the data loss diffusion range. Using the calculated propagation delay, the anomaly's starting point is shifted forward on the time axis to determine the actual time when the anomaly begins to affect the system. The shifted time distance is equal to the propagation delay plus a safety margin, which is set at 20% of the propagation delay. This shift determines the actual anomaly onset time. The spatial diffusion range is determined using the shifted starting point as the center and propagation velocity information. The diffusion radius is calculated based on the propagation velocity and time span. The diffusion range is modeled using a circular region, with the center located at the spatial location of the starting point. The diffusion boundary is discretized, with boundary points evenly distributed along the circumference. The diffusion boundary is adjusted based on the actual sensor distribution to ensure that it fully covers the affected sensor nodes. The resulting data loss diffusion range defines the spatial constraints for frequency drift detection.

[0035] Frequency drift node data is derived based on the data loss diffusion range. Within the defined data loss diffusion range, the frequency characteristics of each sensor node are monitored and analyzed to identify frequency drift phenomena affected by anomalies. Frequency measurements are performed on all sensor nodes within the diffusion range, and the frequency spectrum characteristics of each node are analyzed using the FFT algorithm. The primary and secondary frequency components are extracted, and the frequency stability index is calculated. The frequency drift Δf(t) is defined as Δf(t) = f(t) - f_ref, where Δf(t) is the frequency drift, f(t) is the real-time measurement frequency, and f_ref is the standard reference frequency. The frequency drift detection threshold is set at 0.1% of the reference frequency. When the drift exceeds the threshold and lasts for more than 30 seconds, the node is considered to have frequency drift. Spatial interpolation analysis is performed on the drifting nodes, using the inverse distance weighted interpolation method to estimate the drift extent at unmonitored locations. The drift magnitude and duration are calculated for each drifting node. A frequency drift node database is established, recording the node number and drift magnitude to form a frequency drift node data set. Perform quality control on drift data, filter out pseudo-drift caused by random noise, and ensure the reliability and accuracy of the test results.

[0036] In some embodiments, forming an asynchronous sensor coordinate table based on the frequency drift node data includes: generating a synchronization deviation threshold according to the frequency drift node data; performing upper bound compression on the frequency drift node data based on the synchronization deviation threshold to generate a compressed node interval; performing a stability test on the compressed node interval to generate a stable sub-interval; and selecting the longest continuous segment from the stable sub-interval as the asynchronous sensor coordinate table.

[0037] A synchronization deviation threshold is generated based on frequency drift node data. Using the collected frequency drift node data, the drift distribution characteristics are statistically analyzed, and the frequency deviation statistical parameters of all drift nodes are calculated. The mean and standard deviation of the drift are calculated, and the skewness and kurtosis characteristics of the data distribution are analyzed. The synchronization deviation threshold is determined using the 3σ criterion, with the upper threshold equal to the mean plus three times the standard deviation, and the lower threshold equal to the mean minus three times the standard deviation. Specific deviation thresholds are set for different sensor types: the threshold coefficient for radar sensors is 1.5, the threshold coefficient for vision sensors is 1.0, and the threshold coefficient for IMU sensors is 2.0. The time-varying characteristics of the threshold are calculated, and the threshold parameters are updated using a sliding window to ensure that the threshold can adapt to dynamic changes in system status. The generated synchronization deviation threshold becomes the judgment standard for data compression processing.

[0038] Based on the synchronization deviation threshold, the frequency drift node data is compressed to an upper bound to generate a compressed node interval. Drift data exceeding the upper deviation threshold is truncated, and the compression function limits the maximum allowable value of frequency drift. The compressed data is sorted according to the drift amplitude and then grouped according to the node number. The K-means clustering algorithm is used to divide the compressed data into multiple intervals, and the number of clusters is adaptively determined based on the complexity of the data distribution. Each compressed node interval records the center coordinates, boundary range, number of nodes included, and average drift value. Statistical analysis is performed on the node data within the interval to calculate the drift mean, variance, and distribution characteristics of the interval. When the distance between adjacent intervals is less than the set threshold, they are merged to optimize the interval division results. The generated compressed node interval set constitutes the basic data structure for stability analysis.

[0039] A stability test is performed on the compressed node interval to generate stable subintervals. The generated compressed node interval is evaluated for three-dimensional stability, including temporal stability, spatial stability, and numerical stability. Temporal stability is evaluated by analyzing the time series autocorrelation of the data in the compressed node interval to detect the temporal consistency of the data. Spatial stability is evaluated by calculating the spatial distribution uniformity of the nodes in the compressed node interval, and the distribution characteristics are quantified using the coefficient of variation of the nearest neighbor distance. Numerical stability is evaluated by calculating the coefficient of variation of the drift values ​​in the compressed node interval, which reflects the stability of the numerical distribution. The overall stability score is calculated by combining the stability indicators of the three dimensions. When the overall stability score exceeds the set threshold, the corresponding compressed node interval is marked as a stable subinterval. The compressed node interval that does not meet the stability requirements is split twice, and the stability test is recursively performed until all subintervals meet the requirements.

[0040] The longest continuous segment from the stable subinterval is selected as the asynchronous sensor coordinate table. Among all identified stable subintervals, the interval combination with the best temporal and spatial continuity is found. Temporal continuity is assessed by examining the degree of temporal overlap between adjacent stable subintervals and calculating the ratio of the overlapping time to the total time. Spatial continuity is assessed by examining the spatial adjacency of adjacent stable subintervals and analyzing the distance relationship between interval boundaries. A continuous segment scoring function is established, comprehensively considering three factors: segment length, temporal continuity, and spatial continuity. The continuous segment with the highest score is selected as the optimal solution, ensuring that the selected interval combination has the best overall performance. The selected longest continuous segment is converted into an asynchronous sensor coordinate table, with fields including sensor ID, 3D coordinates, asynchronous level, timestamp, and drift magnitude. An index structure is established for the coordinate table to support fast query operations based on different fields, completing the construction of the asynchronous sensor coordinate table.

[0041] Step S130 , collecting IMU angle increments and sensor response delays, detecting attitude switching signals through the IMU angle increments, capturing data acquisition inflection points through the sensor response delays, and performing time domain matching between the attitude switching signals and the data acquisition inflection points to generate sampling imbalance information.

[0042] Specifically, IMU angle increments and sensor response delays are collected. Based on the asynchronous sensor coordinate table, an IMU angle monitoring module and response delay measurement device are deployed at the asynchronous sensor locations identified in the coordinate table. IMU angle increments are acquired using a six-axis inertial measurement unit (IMU), which includes the combined output of a three-axis gyroscope and a three-axis accelerometer. The gyroscope measures angular velocity changes along three axes, while the accelerometer measures linear acceleration changes along three axes. Angle increments are obtained using an integration method. Continuous integration yields incremental changes in three-dimensional attitude angles, including roll, pitch, and yaw angle increments. Sensor response delay data is collected synchronously, using a high-precision clock synchronization system to measure the complete delay chain from signal input to data output. Response delay measurements include three components: internal sensor processing delay, data transmission delay, and system processing delay. The delay measurement system uses known signal injection and detection to verify the accuracy of delay measurements. The acquisition system is equipped with a ring buffer to support continuous recording and real-time processing of high-speed data streams. All collected data is timestamped with high-precision timestamps to ensure accurate time synchronization between different sensor data. The IMU angle increment and sensor response delay data are time-aligned and spatially registered to form time-space corresponding data records.

[0043] Attitude switch signals are detected using IMU angle increments. Real-time detection of UAV attitude switches is achieved based on collected IMU angle increment data. Attitude switch detection criteria are defined: an attitude switch event is identified when the angle increment amplitude exceeds a preset threshold or the angular acceleration exceeds an acceleration threshold. The threshold parameters are dynamically set based on the UAV's flight mode: 2 degrees / second in hovering mode and 10 degrees / second in rapid maneuvering mode. A sliding window statistical method is used to analyze short-term trends in angle increments, with a window length of 100 milliseconds and a 50% overlap. The mean and standard deviation of the angle increments within the window are calculated, and switch detection is triggered when the standard deviation exceeds 1.5 times the mean. A Butterworth low-pass filter is used to preprocess the angle increments to eliminate high-frequency noise, with a cutoff frequency set to 20 Hz. Attitude switch signals are digitally represented, and three key parameters, namely, switch amplitude, switch rate, and duration, are recorded for each switch event. Detected attitude switch signals are filtered, using both signal amplitude and duration to filter out false detections. All valid attitude switch events are recorded chronologically to form a complete switch signal sequence.

[0044] Capture data acquisition inflection points through sensor response delay. Based on the collected sensor response delay data, identify key turning points and outliers in the data acquisition process. Analyze the statistical characteristics of the response delay data and calculate the mean, standard deviation, and coefficient of variation of the delay data. Define the data acquisition inflection point as the moment when the response delay changes significantly. Use change point detection to automatically identify the inflection point location. Change point detection uses a cumulative sum control method. The cumulative sum formula is CUSUM(t)=max(0,CUSUM(t-1)+(x(t)-μ)-k), where CUSUM(t) is the cumulative sum value at time t, x(t) is the current delay measurement, μ is the historical delay mean, and k is the drift compensation parameter. An inflection point is identified when the cumulative sum exceeds the detection threshold h, which is set to three times the standard deviation. Based on the direction and magnitude of the delay change, inflection points are classified into three types: rising, falling, and fluctuating. The impact strength of the inflection point is calculated based on the magnitude and duration of the delay change. Captured data acquisition inflection points are time-stamped, recording the precise moment and type of inflection point. The signal-to-noise ratio threshold is used to screen out high-intensity inflection points and filter out false inflection points caused by random noise.

[0045] In some embodiments, the time domain matching of the posture switching signal and the data acquisition inflection point to generate sampling imbalance information includes: establishing a motion rhythm time axis based on the posture switching signal; mapping the data acquisition inflection point to the motion rhythm time axis to form an interruption moment mark; dividing the motion rhythm time axis into an interruption front segment and an interruption back segment with the interruption moment mark as the boundary; and performing staggered superposition processing on the interruption front segment and the interruption back segment to generate sampling imbalance information.

[0046] A motion rhythm timeline is established based on attitude switching signals. Using the detected attitude switching signals as a time reference, a standardized timeline describing the drone's motion rhythm is constructed. The attitude switching signals are arranged in chronological order according to their occurrence time, forming a discrete sequence of time events. The time intervals between adjacent switching events are calculated, and their statistical distribution characteristics are analyzed. The motion rhythm timeline uses a normalized representation, mapping actual time to a standard interval of 0-1. Linear interpolation is used to form a continuous time function between switching events, ensuring the continuity of the timeline. The sampling interval of the timeline is determined based on the frequency characteristics of attitude switching, with a typical sampling frequency set at 100 Hz. The motion rhythm timeline is smoothed, using a 5-point moving average method to eliminate local time fluctuations. A time consistency check ensures the integrity of the timeline data, detecting and repairing missing time periods.

[0047] Data acquisition inflection points are mapped onto the motion rhythm timeline to form discontinuity markers. Based on the constructed motion rhythm timeline, captured data acquisition inflection points are projected onto corresponding positions on the rhythm timeline. A time coordinate transformation is used to convert inflection points in the actual time domain into markers in the rhythm time domain. A nearest neighbor search method is used to determine the best matching position of each inflection point on the rhythm timeline, with a search radius set to ±5 seconds. Mapping accuracy is calculated, and the mapping effect is quantitatively evaluated using the root mean square value of the mapping time error. Each discontinuity marker records three basic attributes: time position, inflection point type, and intensity level. The mapped discontinuity markers are sorted chronologically to form a marker sequence. The marker density distribution is calculated, and the change in the number of markers per unit time is statistically analyzed. Discontinuity markers are clustered to identify dense and sparse marker areas. A mapping error threshold is used to select high-precision markers to ensure the temporal accuracy of the markers.

[0048] The movement rhythm timeline is divided into a pre-break segment and a post-break segment using the break markers as the demarcation points. Using the identified break markers as segmentation points, the complete movement rhythm timeline is divided into multiple time segments. Each break marker divides the timeline into two parts: the pre-break segment and the post-break segment. The pre-break segment is set to a 30-second window before the marker, and the post-break segment is set to a 30-second window after the marker. If the segment length is insufficient, zero padding is used to complete it to the required length. Data integrity is checked for each time segment, and the proportion of valid data points within the segment is calculated. Each time segment is assigned a unique sequence number in chronological order. When adjacent time segments overlap, a weighted average is used to process the data in the overlapping area. High-quality time segments are selected based on data integrity, and segments with a valid data ratio below 80% are labeled as low-quality segments. Pre-break segments and post-break segments are paired to form a correspondence table between the pre- and post-break segments.

[0049] Exemplarily, the staggered superposition processing of the discontinuous front segment and the discontinuous rear segment to generate sampling imbalance information includes: converting the discontinuous front segment into a spectrum sequence; staggering and superimposing the discontinuous rear segment onto the spectrum sequence to form a spectrum difference graph; deriving a spectrum jump accumulation amount from the spectrum difference graph; and generating sampling imbalance information based on the distribution unevenness of the spectrum jump accumulation amount.

[0050] The pre-interval segment was converted into a spectral sequence. A fast Fourier transform was performed on the identified pre-interval segment time series to obtain a spectral sequence in the frequency domain. The transform length was set to 1024 points to ensure a frequency resolution of 0.1 Hz. The amplitude spectrum information was extracted from the spectral sequence, while the phase spectrum was discarded to simplify subsequent processing. The spectral sequence was normalized to its maximum value, uniformly mapping the amplitude range to the 0-1 interval. The frequency range of 0-50 Hz was selected for analysis, as this band contains the primary motion feature information. A Hanning window function was used to reduce spectral leakage effects, with the window function length being the same as the FFT length. A three-point moving average smoothing process was performed on the spectral sequence to reduce random fluctuations in the spectrum. High-quality spectral data was filtered using a signal-to-noise ratio threshold, and data with a signal-to-noise ratio below 10 dB was filtered. Statistical characteristic parameters of the spectral sequence, such as the dominant frequency, bandwidth, and total energy, were recorded.

[0051] The segment after the discontinuity is overlaid onto the spectral sequence to form a spectral difference map. The segment after the discontinuity is also converted into a spectral sequence and then subtracted from the spectral sequence before the discontinuity. The spectral difference ΔS(f) = S_after(f) - S_before(f) is calculated, where ΔS(f) is the spectral difference at frequency f, S_after(f) is the spectral amplitude at frequency f after the discontinuity, and S_before(f) is the spectral amplitude at frequency f before the discontinuity. A spectral difference vector containing the difference information at all frequency points is constructed. The spectral differences are amplitude-limited, truncating the difference values ​​to within ±1 to prevent interference from extreme values. The distribution pattern of the spectral differences is displayed in grayscale format to facilitate visual analysis of the difference features. Statistical characteristics of the spectral differences are calculated, including parameters such as the mean, standard deviation, and peak-to-peak value. The statistical significance of the spectral differences is determined using statistical tests, with a significance level of 0.05. A 3×3 median filter is used to filter the difference map to eliminate isolated noise points.

[0052] The spectrum jump accumulation is derived from the spectrum difference plot. Spectral jumps are identified in the spectrum difference plot, with a jump defined as a change in the difference between adjacent frequency points exceeding a threshold. The first-order difference method is used to calculate the rate of change of the spectrum difference and obtain the difference gradient information. A jump detection threshold is set to 0.1, and jump points are marked when the absolute value of the rate of change exceeds this threshold. The spectrum jump accumulation JCA is calculated as JCA = Σ|ΔS(f_i+1)-ΔS(f_i)|, where ΔS(f_i+1) and ΔS(f_i) are the spectrum differences at the i+1th and ith frequency points, respectively. JCA is the jump accumulation, representing the sum of the difference changes at all jump locations, covering all detected jump points. The accumulation is normalized by the total number of frequency points to obtain a jump intensity index per unit frequency. Statistical analysis is performed on the jump accumulation to calculate the distribution range and median value. Jump intensity is categorized into three levels: low, medium, and high, with cutoffs of 0.3 and 0.7, respectively. Track the change trend of jump accumulation over time and identify the growth and decay patterns of accumulation.

[0053] Sampling imbalance information is generated based on the distribution unevenness of the spectral jump accumulation. The distribution characteristics of the spectral jump accumulation are used to quantify the degree of sampling imbalance. Distribution unevenness (CV) is calculated as σ_JCA / μ_JCA, where CV is the coefficient of variation, σ_JCA is the standard deviation of the jump accumulation, and μ_JCA is the mean of the jump accumulation. This metric reflects the degree of dispersion in the jump accumulation distribution. The coefficient of variation value determines the imbalance level: CV < 0.3 indicates mild imbalance, 0.3 ≤ CV < 0.6 indicates moderate imbalance, and CV ≥ 0.6 indicates severe imbalance. Information entropy is used to quantify the complexity of the distribution; higher entropy values ​​indicate more complex imbalance patterns. A multi-dimensional imbalance description is constructed, encompassing four basic dimensions: imbalance intensity, level identifier, complexity, and timestamp. Confidence intervals are estimated for the imbalance information, and the confidence level of the results is calculated using a bootstrap method. The imbalance information is formatted into a standard data structure, containing fields such as level, intensity value, complexity indicator, and detection time. The sampling imbalance information is output for use by subsequent data processing modules, supporting the system's adaptive adjustment capabilities.

[0054] Step S140 , performing correlation mapping on the asynchronous sensor coordinate table and the sampling imbalance information to generate a feature offset, establishing a feature matching core area through the feature offset, and creating a hierarchical calibration sequence based on the feature matching core area.

[0055] Specifically, the asynchronous sensor coordinate table and sampling imbalance information are correlated and mapped to generate feature offsets. Based on the asynchronous sensor coordinate table and sampling imbalance information, spatial-temporal correlation mapping is achieved between the two types of data. The asynchronous sensor coordinate table contains key fields such as sensor ID, three-dimensional spatial coordinates, asynchrony level, and drift characteristics. The sampling imbalance information contains important attribute parameters such as imbalance type, intensity level, impact range, and duration. Spatial proximity analysis is used to match the sensor location in the coordinate table with the impact range in the imbalance information. A spatial matching radius of 50 meters is set, and relevant imbalance events are searched within this radius. For each asynchronous sensor location, all sampling imbalance events within its neighborhood are searched, and the number and type of matching events are recorded. A spatial distance weight is calculated using the formula w = exp(-d² / 2σ²), where d is the distance between the sensor and the imbalance center and σ is an attenuation parameter set to 25 meters. A comprehensive impact coefficient is calculated based on the imbalance intensity and distance weight. The impact coefficient reflects the degree of imbalance on the sensor. A feature offset is derived through weighted fusion, comprehensively considering the cumulative impact of multiple imbalance events on a single sensor. The characteristic offset consists of two components: spatial offset and performance offset. Spatial offset describes changes in sensor position, while performance offset describes changes in sensor parameters. Spatial offset consists of three coordinate components: X, Y, and Z. Performance offset consists of three parameter components: time offset, frequency offset, and amplitude offset. The calculated characteristic offset is amplitude-limited. Spatial offset is limited to ±10 meters, while performance offset is truncated to the normal range of variation for each parameter. The characteristic offsets of all sensors are organized into an offset data set to form a description of the overall offset distribution of the system.

[0056] In some embodiments, establishing the feature matching core area through the feature offset includes: converting the feature offset into an offset vector domain; locating a convergence center node in the offset vector domain; performing area expansion with the convergence center node as the starting point to form a preliminary core area; and performing boundary refinement on the preliminary core area to form a feature matching core area.

[0057] The feature offset is converted into an offset vector domain. The feature offset at each sensor location is represented as a six-dimensional vector, with the first three dimensions representing the spatial offset component and the second three dimensions representing the performance offset component. An inverse distance weighted interpolation method is used to construct a continuous vector field between sensor nodes, with an interpolation radius set to 100 meters. The three-dimensional space is divided into a grid of 20 m × 20 m × 20 m, and the interpolated offset vector value is calculated at each grid node. The offset vector domain is normalized, mapping each component to a value range of 0–1 to eliminate dimensionality differences between different dimensions. The magnitude distribution of the vector field is calculated, representing the overall magnitude of the offset. The directional distribution of the vector field is analyzed to identify the main directional characteristics of the offset. The offset vector is smoothed and filtered, using a 3×3×3 cubic kernel for convolution to reduce the influence of local noise. The gradient of the vector field is calculated, reflecting the spatial distribution of the offset variation.

[0058] Locate convergence centers within the offset vector domain. Identify convergence centers of vector flow by analyzing the flow characteristics of the offset vector domain. Calculate the divergence of the offset vector field. A negative divergence value indicates a convergence effect of vector flow at that location. Use the central difference method to calculate the divergence, with the difference step size set to the grid spacing. Define the convergence strength as the magnitude of the negative divergence value; a larger convergence strength indicates a more pronounced convergence effect. Use local extrema detection to identify local minima in the divergence field; these points correspond to potential convergence centers. Set the convergence strength threshold to -0.1 to eliminate candidate points with weak convergence effects. Merge convergence center candidates within 50 meters of each other to avoid overly dense convergence nodes. Calculate the influence strength of each convergence center, which is equal to the product of the convergence strength and the amplitude of the surrounding vectors. Select the top five locations with the largest influence strengths as primary convergence centers. Record the spatial coordinates and convergence strength values ​​of the convergence centers to provide starting parameters for region expansion.

[0059] Starting from the convergence center node, regional expansion is performed to form a preliminary core domain. Starting from each identified convergence center node, a region growing method is used to gradually expand the domain to form a continuous spatial region. Criteria for region expansion are set: when the offset vector magnitude of adjacent grid points is less than 1.2 times the current magnitude and the direction angle is less than 30 degrees, the grid point is included in the expansion region. The expansion process uses a breadth-first search strategy, starting from the convergence center and expanding outward layer by layer to ensure the spatial connectivity of the region. Expansion stop conditions are set: when grid points on the expansion boundary do not meet the expansion criteria or the expansion radius exceeds 200 meters, the region expansion is terminated. Connectivity checks are performed on the expansion region formed by each convergence center to ensure that all points in the region are connected to the convergence center through adjacency relationships. Basic geometric properties of the expansion region are calculated, including regional volume, surface area, and compactness. When multiple expansion regions overlap, region ownership is determined based on the convergence strength, with convergence centers with higher convergence strengths receiving priority in the overlapping region. The resulting preliminary core domain is regularized, isolated protrusions are removed, and small internal voids are filled. The boundary coordinates and number of grid points included in each preliminary core domain are recorded.

[0060] The preliminary core region is refined to form the feature matching core region. For the preliminary core region generated by region expansion, a boundary optimization method is used to improve the accuracy of the region boundary. The preliminary core region is processed using a morphological opening operation, first eroding and then dilating, to remove small protrusions on the boundary. A spherical core with a radius of 5 meters is set as the structuring element to ensure uniformity of the processing effect. The processed region boundary is smoothed using the Laplace smoothing method to reduce jagged edges. The number of smoothing iterations is set to 10, with a smoothing strength coefficient of 0.1 per iteration. Irregular boundary protrusions are removed using a boundary shrinkage operation, with the shrinkage distance set to the size of one grid cell. Convex hull detection is performed on the boundary, and the ratio of the convex hull area to the actual area of ​​the region is calculated. A ratio close to 1 indicates a relatively regular region shape. The minimum bounding box method is used to describe the refined region boundary, with the bounding box's orientation aligned with the region's principal axis. Geometric parameters such as the center coordinates, principal axis length, and orientation angle of the feature matching core region are calculated. The resulting feature matching core regions are numbered and sorted from highest to lowest convergence strength.

[0061] A hierarchical calibration sequence is created based on the feature matching core region. Using the established feature matching core region as a spatial reference, a hierarchical calibration strategy is designed for the drone sensor system. The calibration process is divided into three levels based on the sensor's position within the core region and the size of its offset. Level 1 calibration targets sensors located within the primary core region with an offset greater than a threshold. Calibration requires high accuracy, and calibration parameters include position correction, time synchronization, and frequency adjustment. Level 2 calibration targets sensors located within the secondary core region or with a moderate offset. It requires moderate accuracy and primarily performs position correction and time synchronization. Level 3 calibration targets sensors located at the edge of the core region or with a small offset. Calibration requires relatively low accuracy and primarily performs basic time synchronization. The calibration priority calculation formula is set as P = w1 × (1 / d) + w2 × |offset|, where d is the distance from the sensor to the center of the core region, offset is the offset size, and w1 and w2 are weighting coefficients. All sensors requiring calibration are sorted according to their priority values, with sensors with higher priority being calibrated first. A calibration time allocation strategy is established: 60 seconds for level 1 calibration, 30 seconds for level 2 calibration, and 15 seconds for level 3 calibration. Set the calibration parameter adjustment range to ±5 meters for position correction, ±1 millisecond for time synchronization accuracy, and ±100 Hz for frequency adjustment. Develop a calibration sequence execution process, including four basic steps: parameter reading, calibration calculation, parameter writing, and verification. Organize the hierarchical calibration sequence into an execution table, containing fields such as sensor ID, calibration level, priority, estimated time, and calibration parameters. Export the calibration sequence for use by the sensor calibration system, supporting automated calibration execution.

[0062] Step S150 , performing internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data, performing extended path analysis on the source offset data to obtain offset transfer velocity, and fusing the offset transfer velocity and passive offset data to form an offset coupling factor.

[0063] In some embodiments, the processing of identifying endogenous and exogenous differences on the hierarchical calibration sequence to derive source offset data and passive offset data includes: performing difference source identification detection on the hierarchical calibration sequence to derive endogenous differences and exogenous differences; performing an impact assessment on the exogenous differences based on the endogenous differences to form a difference correlation matrix; deriving a dominant difference vector based on eigenvalue decomposition of the difference correlation matrix; and performing difference classification according to the dominant difference vector to create source offset data and passive offset data.

[0064] Difference source identification and detection are performed on a hierarchical calibration sequence to derive endogenous and exogenous differences. Based on the calibration task information contained in the hierarchical calibration sequence, different sources of differences are identified. The time series of offset corrections for each calibration task is analyzed to extract basic statistical characteristics of the offsets. Key parameters such as the mean, standard deviation, amplitude, and duration of the offset series are calculated. The sources of differences are distinguished based on the temporal characteristics of the offset changes. Slow-changing and long-lasting offsets are classified as endogenous, while fast-changing and short-lasting offsets are classified as exogenous. A time threshold is set for distinguishing between endogenous and exogenous sources: offset events lasting more than 300 seconds are classified as endogenous, while those lasting less than 60 seconds are classified as exogenous. Endogenous differences primarily include internal issues such as sensor parameter drift, calibration errors, and device aging. Exogenous differences primarily include external influences such as environmental interference, electromagnetic interference, temperature changes, and vibration and shock. Each offset event in the calibration sequence is automatically classified, and the distribution of the number of endogenous and exogenous differences is calculated. The intensity distribution of endogenous and exogenous differences is calculated, with the average intensity of endogenous differences generally being lower than that of exogenous differences. Record the basic attributes of each difference event, such as the occurrence time, duration, offset magnitude, and difference type.

[0065] Based on the impact of endogenous and exogenous differences, an impact correlation matrix was constructed. The interaction between endogenous and exogenous differences was analyzed, quantifying the strength of the association between the two types of differences. The temporal correlation between endogenous and exogenous differences was calculated, and the Pearson correlation coefficient was used to measure the degree of linear association. The spatial distance between sensors was considered to moderate the impact relationship, with closer sensors exhibiting a greater impact. The impact was calculated as I = ρ × exp(-d² / 2σ²), where ρ is the temporal correlation coefficient, d is the distance between sensors, and σ is the spatial attenuation parameter, set to 50 meters. The impact values ​​between each pair of sensors were calculated to form the basic elements of the impact matrix. A impact correlation matrix was constructed, with rows corresponding to endogenous and exogenous difference events and columns corresponding to exogenous difference events. The matrix elements represent the impact values ​​between the corresponding differences. The correlation matrix was normalized, mapping the impact values ​​to a standard range of 0–1. An impact threshold of 0.3 was set, and weak correlations below the threshold were eliminated, retaining pairs with significant associations. The proportion of non-zero elements in the correlation matrix is ​​counted to reflect the density of the correlation relationship between differences. The processed difference correlation matrix is ​​stored in sparse matrix format to improve storage efficiency and calculation speed.

[0066] The dominant difference vectors (DVs) are derived from the difference correlation matrix using eigenvalue decomposition. Eigenvalue decomposition is used to analyze the main eigenmodes of the difference correlation matrix. Eigenvalue decomposition is performed as A=QΛQ^T, where A is the difference correlation matrix, Q is the eigenvector matrix, T is the transpose, and Λ is the eigenvalue diagonal matrix. Eigenvectors are sorted in descending order by eigenvalue, with larger eigenvalues ​​indicating greater importance. The eigenvectors corresponding to the first three largest eigenvalues ​​are selected as the dominant difference vectors, which reflect the main difference correlation patterns. The cumulative contribution of the first three eigenvalues ​​is calculated, which typically explains more than 80% of the total variation. The component distribution of the dominant difference vectors is analyzed; the sign and magnitude of the vector components reflect the direction and strength of the correlation between different differences. The dominant difference vectors are normalized to ensure a vector modulus of 1. Key components of the dominant difference vectors are extracted to identify the sensor locations that contribute most to the difference correlation. The dominant difference vectors are organized into a matrix format to facilitate subsequent difference classification calculations. The eigenvalue corresponding to each dominant difference vector is recorded as a quantitative indicator of its importance.

[0067] Difference classification is performed based on the dominant difference vectors to create source-type and passive-type drift data. Using the extracted dominant difference vectors as classification criteria, identified differences are categorized as source-type and passive-type drift. The difference projection score S = V^T × d is calculated, where V is the dominant difference vector matrix, T is the transpose, and d is the difference eigenvector. The difference classification is determined based on the projection score: differences with a projection score greater than 0 are classified as source-type drift, while differences with a projection score less than 0 are classified as passive-type drift. Source-type drift corresponds to a source of drift with active propagation capabilities and typically manifests as a dominant difference affecting other sensors. Passive-type drift corresponds to a responsive drift caused by other drifts and typically manifests as a follow-up response difference. The distribution of source-type and passive-type drift is statistically analyzed, with source-type drift typically being less common than passive-type drift. Key attributes of source-type drift data are recorded, including the sensor ID of the drift source, spatial location, drift intensity, and impact range. Basic information for passive-type drift data is recorded, including the ID of the affected sensor, response amplitude, response delay, and response duration. The classification results were cross-validated and the consistency of the classification was verified by using different combinations of dominant difference vectors.

[0068] Extended path analysis is performed on source-based offset data to obtain offset propagation velocity. Based on the identified source-based offset data, the propagation characteristics of offset phenomena in the sensor network are analyzed. Source-based offset data contains key information such as the spatial location of the offset source, offset intensity, and offset type. A wavefront propagation model is used to describe the process of offset diffusion from the source to surrounding sensors. The time series of offset arrival at different sensors are calculated to analyze the spatiotemporal evolution of offset propagation. The propagation velocity is calculated based on the relationship between offset arrival time and propagation distance. The offset propagation velocity v = d / Δt, where v is the propagation velocity, d is the propagation distance, and Δt is the propagation time difference. Statistical analysis of propagation velocities in different directions is performed to identify the directional characteristics of the propagation velocity. The mean and standard deviation of the propagation velocity are calculated to reflect the average characteristics and variability of the propagation process. The relationship between propagation velocity and sensor density is analyzed; high-density areas generally have faster propagation velocities. Key nodes in the propagation process are identified, as these nodes play an important relay role in offset propagation. The propagation velocity distribution of different source-based offsets is recorded to form a propagation velocity database.

[0069] The offset transfer velocity and passive offset data are combined to form an offset coupling factor. The offset transfer velocity information is combined with the passive offset data for analysis to construct a coupling model for inter-sensor offset interactions. Passive offset data contains characteristic parameters of the affected sensors, such as response delay, response amplitude, and response duration. The offset coupling factor C = (v × R) / (τ + 1) is calculated, where C is the coupling factor, v is the offset transfer velocity, R is the passive response amplitude, and τ is the response delay. The coupling factor values ​​are calculated for each sensor pair to form a coupling factor matrix. Based on the magnitude of the coupling factor, the inter-sensor interactions are classified into three levels: strong, medium, and weak. Strong coupling corresponds to sensor pairs with a coupling factor greater than 0.7, medium coupling ranges from 0.3 to 0.7, and weak coupling is less than 0.3. The spatial distribution of the coupling factor is analyzed to identify high and low coupling regions. The temporal variation of the coupling factor is calculated to track the temporal evolution of the coupling relationship. The number of sensor pairs with different coupling levels is counted to analyze the overall coupling characteristics of the system. Cluster analysis of coupling factors identifies groups of sensors with similar coupling characteristics. The offset coupling factors are organized into a standardized data format, including fields such as sensor pair ID, coupling strength, coupling level, and a timestamp. The offset coupling factors are output for use by system optimization and control modules, supporting dynamic management and performance improvement of sensor networks.

[0070] Step S160 : determining a fusion enhancement area based on the offset coupling factor and the asynchronous sensor coordinate table, and performing cascade analysis on the fusion enhancement area to form a multi-source data fusion result.

[0071] Specifically, fusion enhancement regions are determined based on the offset coupling factor and the asynchronous sensor coordinate table. The offset coupling factor contains core information such as coupling strength and coupling relationship, describing the interaction characteristics of offset phenomena between sensors. The asynchronous sensor coordinate table records key fields such as sensor ID, three-dimensional spatial coordinates, asynchronous level, and drift characteristics. Spatial matching is performed between the sensor location information in the coordinate table and the coupling node pairs in the coupling factor. Areas with dense sensors and high coupling strength generally achieve better data fusion results, and such areas are identified as candidate fusion enhancement regions. The selection criteria for fusion enhancement regions are set to require that the average coupling strength of sensors within the area is greater than 0.6 and the asynchronous level is less than 3. Sensor locations that meet these criteria are grouped using a spatial clustering method to form contiguous spatial regions. Sensor density and data quality indicators are calculated for each candidate region, with regions with high sensor density and good data quality being prioritized. Regions with an area of ​​100-500 square meters are selected as fusion enhancement regions to ensure that the region size is appropriate for data processing. The boundaries of the selected fusion enhancement regions are determined, and the minimum enclosing rectangle method is used to determine the spatial extent of the region. Record the center coordinates, boundary range, and sensor list of the fusion enhancement area to provide a spatial basis for subsequent cascade analysis. Establish an index table for the fusion enhancement area to support fast query and data access operations.

[0072] In some embodiments, the cascade analysis of the fusion enhanced area is performed to form a multi-source data fusion result, including: determining an analysis path based on detecting serial extension features of the fusion enhanced area, the serial extension features including point cloud sparse gradient, image blur degree and IMU drift rate; tracing the fusion transfer process along the analysis path to form a transfer link graph; deriving the gain coefficient of each node from the transfer link graph; and performing cascade fusion encoding on the gain coefficients to form a multi-source data fusion result.

[0073] The analysis path is determined based on the serial extension features detected in the fusion-enhanced region. Within the defined fusion-enhanced region, the main change features of the multi-source sensor data are detected. The density distribution of the radar point cloud data is analyzed, and the density change rate of the point cloud in different directions is calculated. Changes in visual image clarity are detected, and the degree of blur is assessed through image contrast and edge strength. The angular velocity and acceleration stability of the IMU sensor are monitored, and the drift rate trend is calculated. The three feature data types are normalized to eliminate dimensional differences between different data types. The spatial gradient of feature changes is calculated to identify the spatial directions with the most significant data changes. The directions with the largest gradient changes are selected as the primary analysis path, and the directions with moderate gradient changes are selected as auxiliary analysis paths. The sampling interval of the analysis path is set to 10 meters to ensure complete path coverage and analysis accuracy. The defined analysis path is spatially labeled, and the coordinates of key nodes along the path are recorded. The analysis path reflects the main patterns and spatial distribution patterns of data changes within the fusion-enhanced region.

[0074] The fusion transfer process is traced along the analysis path to form a transfer link graph. Based on the determined analysis path, the processing flow of multi-source data within the fusion enhancement area is traced. Starting from the starting point of the analysis path, transfer nodes are established according to the order of data processing. When data passes through a sensor location and undergoes fusion processing, a transfer node is set at that location. Each transfer node records basic information such as the processed data type, fusion algorithm type, and output data format. Directed connections are established between transfer nodes based on the data processing order to form a data flow path. The weight of the connecting edge reflects the importance of data transfer and is determined by the data volume and processing complexity. The complete transfer link graph is represented using a graph structure, with nodes representing data processing locations and edges representing data flow relationships. The structural characteristics of the transfer link graph are analyzed to identify the main data flow paths and key node locations. The depth and breadth of the transfer link graph are calculated to reflect the hierarchical structure and degree of parallelism of data processing. The in-degree and out-degree distribution of transfer nodes are statistically analyzed to identify data aggregation nodes and data distribution nodes.

[0075] The gain coefficient of each node is derived from the transfer link graph. By analyzing the network structure of the transfer link graph, the contribution of each node to the final fusion result is calculated. The number of connections for each node is counted. Nodes with more connections are generally more important. The average path length from a node to other nodes is calculated. Nodes with shorter paths are more important in data transmission. The distribution of nodes on the critical path is analyzed. Nodes located on the critical path are given higher weights. The node gain coefficient G_i = (C_i + P_i + K_i) / 3 is calculated, where G_i is the gain coefficient of the i-th node, C_i is the connection importance, P_i is the path importance, and K_i is the critical path importance. The gain coefficients of all nodes are normalized to ensure that the sum of the gain coefficients is equal to 1. The gain coefficients are adjusted based on the node's data processing capacity, with the gain coefficients of nodes with higher processing capabilities appropriately increased. The distribution characteristics of the gain coefficients are analyzed to identify high-importance nodes, medium-importance nodes, and low-importance nodes. The gain coefficient value of each node is recorded to provide the weight parameter for cascade fusion coding.

[0076] The gain coefficients are cascaded and fused to generate a multi-source data fusion result. The calculated gain coefficients for each node are used to perform cascade fusion encoding. Cascade fusion employs a layered processing strategy, performing data fusion sequentially according to the hierarchical structure of the transmission link graph. The first layer performs weighted fusion on the raw sensor data, with the fusion weights equal to the gain coefficients of the corresponding nodes. The second layer extracts features from the fusion results of the first layer, retaining the most important data feature components. The third layer performs a final fusion of the extracted features to generate a unified multi-source data fusion result. The weighted fusion formula R = Σ(G_i × D_i) is used, where R is the final fusion result, G_i is the gain coefficient of the i-th node, and D_i is the data vector of the i-th node. Error control is implemented during the fusion process, and the gain coefficient allocation is adjusted when the error in the fusion result exceeds a threshold. The data format of the fusion result is set, including fields such as fused data, confidence level, timestamp, and source information. The fusion result is quality assessed, ensuring its reliability through consistency checks and anomaly detection. An output interface for the fusion result is established, supporting both real-time data transmission and batch data processing. The multi-source data fusion results are stored in a standardized format to facilitate subsequent data analysis and application calls. Ultimately, the data fusion processing of the drone multi-sensor system is completed, achieving effective integration of heterogeneous sensor data.

[0077] In order to implement a non-aligned spatiotemporal data driven multi-sensor fusion method corresponding to the above method embodiment, to achieve the corresponding functions and technical effects. Figure 2 , Figure 2The following is a block diagram of a multi-sensor fusion system 200 driven by non-aligned spatiotemporal data according to an embodiment of the present application. For ease of explanation, only the parts relevant to this embodiment are shown. The multi-sensor fusion system 200 driven by non-aligned spatiotemporal data according to an embodiment of the present application includes: A data acquisition module 201 is configured to collect radar point cloud density and visual image texture, obtain spatial depth variation parameters based on the point cloud density, capture visual motion profiles through the visual image texture, and establish a data synchronization standard using the spatial depth variation parameters and the visual motion profiles; an asynchronous detection module 202 configured to detect a timing anomaly cluster based on the data synchronization standard, determine a data loss diffusion range based on the timing anomaly cluster, derive frequency drift node data based on the data loss diffusion range, and form an asynchronous sensor coordinate table based on the frequency drift node data; A sampling processing module 203 is configured to collect IMU angle increments and sensor response delays, detect an attitude switching signal using the IMU angle increments, capture a data acquisition inflection point using the sensor response delay, and perform time domain matching between the attitude switching signal and the data acquisition inflection point to generate sampling imbalance information; A feature mapping module 204 is configured to perform correlation mapping between the asynchronous sensor coordinate table and the sampling imbalance information to generate a feature offset, establish a feature matching core region based on the feature offset, and create a hierarchical calibration sequence based on the feature matching core region; An offset analysis module 205 is configured to perform internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data, perform extended path analysis on the source offset data to obtain an offset transfer velocity, and fuse the offset transfer velocity and the passive offset data to form an offset coupling factor. The fusion processing module 206 is configured to determine a fusion enhancement region based on the offset coupling factor and the asynchronous sensor coordinate table, and perform cascade analysis on the fusion enhancement region to form a multi-source data fusion result.

[0078] The aforementioned non-aligned spatiotemporal data-driven multi-sensor fusion system 200 can implement the non-aligned spatiotemporal data-driven multi-sensor fusion method of the aforementioned method embodiment. The optional options in the aforementioned method embodiment also apply to this embodiment and are not described in detail here. The remaining contents of the present application embodiment can be referenced to the contents of the aforementioned method embodiment and are not further described in this embodiment.

[0079] The above description is only a partial or preferred embodiment of the present application. Neither the text nor the drawings can limit the scope of protection of the present application. Any equivalent structural transformation made by using the contents of the present application specification and drawings under the overall concept of the present application, or direct / indirect application in other related technical fields, is included in the scope of protection of the present application.

Claims

1. A multi-sensor fusion method driven by non-aligned spatiotemporal data, characterized in that: include: Collecting radar point cloud density and visual image texture, obtaining spatial depth variation parameters based on the point cloud density, capturing visual motion profiles through the visual image texture, and establishing a data synchronization standard through the spatial depth variation parameters and the visual motion profiles; Detecting a timing anomaly cluster area based on the data synchronization standard, determining a data loss diffusion range based on the timing anomaly cluster area, deriving frequency drift node data based on the data loss diffusion range, and forming an asynchronous sensor coordinate table based on the frequency drift node data; collecting IMU angle increments and sensor response delays, detecting an attitude switching signal using the IMU angle increments, capturing a data acquisition inflection point using the sensor response delays, and performing time domain matching between the attitude switching signal and the data acquisition inflection point to generate sampling imbalance information; Performing correlation mapping between the asynchronous sensor coordinate table and the sampling imbalance information to generate a feature offset, establishing a feature matching core area based on the feature offset, and creating a hierarchical calibration sequence based on the feature matching core area; performing internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data, performing extended path analysis on the source offset data to obtain an offset transfer velocity, and fusing the offset transfer velocity and the passive offset data to form an offset coupling factor; A fusion enhancement region is determined based on the offset coupling factor and the asynchronous sensor coordinate table, and a cascade analysis is performed on the fusion enhancement region to form a multi-source data fusion result.

2. The method according to claim 1, characterized in that The establishing of a data synchronization standard using the spatial depth change parameter and the visual motion profile includes: Identifying high-fit regions and low-fit regions based on a matching relationship between the spatial depth variation parameter and the visual motion profile; Couple the peak value of the high-fit region with the valley value of the low-fit region to generate a synchronous potential difference; constructing a synchronous gradient field using the synchronous potential difference; A data synchronization criterion is determined along the data contour of the synchronization gradient field.

3. The method according to claim 1, characterized in that Determining the data missing diffusion range based on the time series abnormality clustering area includes: Generate a data backtracking trajectory based on tracing back the time series anomaly cluster area; Identify the starting point of the anomaly based on the data backtracking trajectory; Obtaining the propagation delay through the starting point of the abnormality; The starting point of the abnormality occurrence is moved forward based on the propagation delay to determine the data loss diffusion range.

4. The method according to claim 1, wherein The forming of the asynchronous sensor coordinate table based on the frequency drift node data includes: generating a synchronization deviation threshold according to the frequency drift node data; Performing upper bound compression on the frequency drift node data based on the synchronization deviation threshold to generate a compressed node interval; Performing a stability test on the compressed node interval to generate a stable subinterval; The longest continuous segment is selected from the stable subinterval as the asynchronous sensor coordinate table.

5. The method according to claim 1, wherein The step of performing time domain matching between the posture switching signal and the data acquisition inflection point to generate sampling imbalance information includes: Establishing a motion rhythm time axis based on the posture switching signal; Mapping the data collection inflection point to the motion rhythm time axis to form a discontinuity moment mark; Dividing the movement rhythm time axis into a pre-interruption segment and a post-interruption segment based on the intermittent moment mark; Performing staggered superposition processing on the discontinuous front section and the discontinuous rear section to generate sampling imbalance information.

6. The method according to claim 1, characterized in that The step of establishing a feature matching core area by using the feature offset includes: Converting the feature offset into an offset vector domain; locating a convergence center node in the offset vector domain; Taking the convergence center node as a starting point, regional expansion is performed to form a preliminary core area; Boundary refinement is performed on the preliminary core region to form a feature matching core area.

7. The method according to claim 1, characterized in that The step of performing internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data includes: performing difference source identification detection on the hierarchical calibration sequence to derive endogenous differences and exogenous differences; Based on the endogenous differences, an impact evaluation is performed on the exogenous differences to form a difference correlation matrix; deriving a dominant difference vector based on eigenvalue decomposition of the difference correlation matrix; Difference classification is performed according to the dominant difference vector to create source offset data and passive offset data.

8. The method according to claim 1, characterized in that The performing cascade analysis on the fusion enhancement area to form a multi-source data fusion result includes: Determining an analysis path based on the fused enhanced region detection serial extension features, the serial extension features including point cloud sparse gradient, image blur level, and IMU drift rate; Tracing the fusion transfer process along the analysis path to form a transfer link diagram; deriving a gain coefficient of each node from the transfer link graph; The gain coefficients are cascaded and fused to form a multi-source data fusion result.

9. The method according to claim 5, characterized in that The step of performing staggered superposition processing on the discontinuous front segment and the discontinuous rear segment to generate sampling imbalance information includes: converting the discontinuous front segment into a spectrum sequence; The discontinuous rear segment is staggered and superimposed on the spectrum sequence to form a spectrum difference map; deriving a spectrum jump accumulation amount from the spectrum difference map; Sampling imbalance information is generated based on the distribution unevenness of the spectrum hopping accumulation amount.

10. A multi-sensor fusion system driven by non-aligned spatiotemporal data, characterized in that: include: a data acquisition module for collecting radar point cloud density and visual image texture, obtaining spatial depth variation parameters based on the point cloud density, capturing visual motion profiles through the visual image texture, and establishing a data synchronization standard through the spatial depth variation parameters and the visual motion profiles; an asynchronous detection module, configured to detect a timing anomaly cluster area based on the data synchronization standard, determine a data loss diffusion range based on the timing anomaly cluster area, derive frequency drift node data based on the data loss diffusion range, and form an asynchronous sensor coordinate table based on the frequency drift node data; a sampling processing module, configured to collect IMU angle increments and sensor response delays, detect an attitude switching signal using the IMU angle increments, capture a data acquisition inflection point using the sensor response delays, and perform time domain matching between the attitude switching signal and the data acquisition inflection point to generate sampling imbalance information; a feature mapping module, configured to perform correlation mapping between the asynchronous sensor coordinate table and the sampling imbalance information to generate a feature offset, establish a feature matching core region based on the feature offset, and create a hierarchical calibration sequence based on the feature matching core region; an offset analysis module for performing internal and external source difference identification processing on the hierarchical calibration sequence to derive source offset data and passive offset data, performing extended path analysis on the source offset data to obtain an offset transfer velocity, and fusing the offset transfer velocity and the passive offset data to form an offset coupling factor; A fusion processing module is used to determine a fusion enhancement area based on the offset coupling factor and the asynchronous sensor coordinate table, and perform cascade analysis on the fusion enhancement area to form a multi-source data fusion result.

Images