A Method and System for Anomaly Detection of UAV Atmospheric Data Based on Spectral Analysis

Abstract

This invention discloses a method and system for anomaly detection in UAV atmospheric data based on graph analysis, belonging to the field of data detection technology. The method involves real-time acquisition of atmospheric observation data, time synchronization processing, and spatial mapping processing to generate an atmospheric data sequence. Time-frequency analysis is performed on the atmospheric data sequence, introducing a time-frequency distribution matrix driven by dual-condition vectors to decompose the atmospheric data sequence into multi-scale perturbation vectors, generating multi-scale perturbation vectors. These multi-scale perturbation vectors are output to a preset anomaly graph detection model. Adaptive learning of the graph structure and a multi-head attention mechanism are used to propagate and fuse node features, identifying single and compound anomalies, inferring the root causes of the anomalies, and outputting the anomaly detection results and root causes. Based on the anomaly detection results and root causes, targeted corrections are made to different observed anomalies to obtain corrected atmospheric data, thereby improving the quality and reliability of UAV atmospheric observation data.

Description

Technical Field

[0001] This invention relates to the field of data detection technology, specifically to a method and system for detecting anomalies in UAV atmospheric data based on spectral analysis. Background Technology

[0002] Unmanned aerial vehicles (UAVs) are now widely used in environmental monitoring and meteorological observation. However, existing methods for detecting anomalies in UAV atmospheric data still have shortcomings. First, during UAV flight, the collected atmospheric data is easily affected by the combined effects of multiple factors, such as aircraft dynamic disturbances and flight status interference. Existing disturbance analysis methods mostly use single-scale or fixed models for decomposition, making it difficult to effectively distinguish disturbance components from different sources. Second, UAV atmospheric data anomalies are often not single anomalies, but rather complex anomalies formed by the superposition of multiple factors. Existing anomaly detection methods mostly rely on single thresholds and simple diagnostic models for judgment, which can only identify the overall abnormal state and cannot further distinguish the components of the anomaly or locate multiple root causes. This leads to a lack of specificity in subsequent dynamic corrections and makes it difficult to guarantee correction accuracy. Summary of the Invention

[0003] To address the problems existing in the background technology, this invention discloses a method and system for detecting anomalies in UAV atmospheric data based on spectral analysis, which improves the quality and reliability of UAV atmospheric observation data.

[0004] The technical solution to achieve the objective of this invention is as follows:

[0005] On the one hand, this invention provides a method for detecting anomalies in UAV atmospheric data based on spectral analysis, including the following steps:

[0006] Acquire atmospheric observation data collected in real time during the flight of the UAV, perform time synchronization processing on the atmospheric observation data, and map it uniformly to the ground reference coordinate system to generate an atmospheric data sequence;

[0007] Time-frequency analysis is performed on the atmospheric data sequence to construct dynamic condition vectors and motion condition vectors. Based on the two types of condition vectors, the corresponding time-frequency separation parameters are calculated. A time-frequency distribution matrix driven by dual condition vectors is introduced to perform multi-scale perturbation decomposition on the atmospheric data sequence to generate multi-scale perturbation vectors. The multi-scale perturbation vectors include real atmospheric signals, motion coupling perturbations, body dynamic perturbations, and atmospheric turbulence perturbations.

[0008] The multi-scale perturbation vector is mapped to a preset anomaly map detection model, which includes an observation layer map and an anomaly layer map. The map structure adaptive learning and multi-head attention mechanism are used to propagate and fuse node features, calculate the anomaly probability of anomaly layer nodes, identify single anomalies and compound anomalies, and infer the root cause of anomalies based on the anomaly propagation path, and output the anomaly detection result and the root cause of anomalies.

[0009] Based on the anomaly detection results and root causes, targeted corrections are made to different observational anomalies to obtain corrected atmospheric observation data.

[0010] Specifically, a multi-source atmospheric sensor and a multi-source flight status sensor are integrated at a fixed position on the UAV fuselage to simultaneously collect atmospheric observation data and flight status data. The atmospheric observation data includes at least atmospheric parameters such as temperature, humidity, air pressure, wind speed, wind direction, and particulate matter concentration. The flight status data includes attitude parameters, navigation parameters, and power system parameters. The attitude parameters include three-axis attitude angles and corresponding angular velocities, the navigation parameters include the UAV's flight altitude, speed, and acceleration, and the power system parameters include information such as motor speed and throttle output. Before the UAV takes off, all types of sensors undergo unified initialization processing, including zero-point calibration and communication status detection, to ensure the consistency of measurement benchmarks and the stability of the data link. During the UAV's monitoring mission, various real-time atmospheric observation data and flight status data are simultaneously collected based on the preset sampling frequency of each sensor.

[0011] Furthermore, a reference time axis synchronization mechanism is employed to perform time-series alignment processing on the multi-source acquired data. Specifically, sensor data with higher sampling frequencies and better data stability are selected from the multi-source atmospheric sensors as the reference time axis. After obtaining the reference time axis, time calibration processing is performed on the other sensor data. For sensor data with sampling frequencies lower than the reference time axis, linear interpolation is used for time resampling to generate matching data values ​​at the corresponding time nodes of the reference time axis. For data with sampling frequencies higher than the reference time axis, synchronous downsampling processing is used to ensure consistency with the reference time axis time series. In addition, abnormal jumps or lost data points that occur during the acquisition process are repaired or removed by interpolation to ensure the continuity and stability of the time series. After the above time synchronization processing, each reference timestamp corresponds to a complete set of atmospheric observation data and flight status data, resulting in time-calibrated atmospheric observation data.

[0012] Furthermore, to eliminate the coupled attitude effects caused by changes in aircraft attitude, a unified spatial ground reference coordinate system is established, and attitude compensation and spatial mapping processing are performed on atmospheric observation data, specifically including:

[0013] Establish a body coordinate system and a ground reference coordinate system. The body coordinate system takes the center of the UAV body as the origin, with the X-axis pointing forward, the Y-axis pointing to the right side of the UAV body, and the Z-axis perpendicular to the body and upward. The ground reference coordinate system takes the projection point of the UAV's instantaneous center of gravity on the ground as the origin, with the X-axis pointing due north as located by the GPS module, the Y-axis pointing due east, and the Z-axis perpendicular to the ground and upward.

[0014] During flight, the three-axis attitude angles of the UAV are acquired in real time, converted into quaternion representations, and the attitude transformation matrix of the body coordinate system relative to the ground reference coordinate system is calculated in real time using a quaternion attitude calculation algorithm. The attitude calculation frequency is kept consistent with the IMU's sampling frequency to update the attitude transformation matrix in real time. ;

[0015] Using attitude transformation matrix The vector atmospheric data in the body coordinate system is rotated and mapped to the ground reference coordinate system to obtain atmospheric data under a unified spatial reference. The wind direction angle is recalculated in the ground reference coordinate system. The wind direction angle is recalculated with the X-axis of the ground reference coordinate system, i.e., due north, as the reference. The wind direction angle is defined as zero degrees with due north and increasing in a clockwise direction, so as to ensure that the wind direction information can accurately reflect the real atmospheric flow direction and eliminate the directional deviation introduced by the UAV attitude change.

[0016] After the above time synchronization and spatial attitude compensation processes, atmospheric observation data under time synchronization and a unified spatial reference coordinate system can be obtained, thus forming an atmospheric data sequence with consistent spatiotemporal semantics.

[0017] Furthermore, the atmospheric data sequence is subjected to low-pass filtering to remove high-frequency noise generated during sensor sampling. Let the atmospheric data collected by the i-th sensor at time t be... The preprocessed atmospheric data sequence was obtained. ;

[0018] Decoupling of multi-source coupled disturbances in the observed signals is performed. Considering the characteristics of low-altitude flight of the UAV and the complex operating conditions, the observed signals are mainly subjected to the superposition and coupling of three types of disturbances. The atmospheric data collected by the airborne sensors can be decomposed into four independent signals, represented as follows:

[0019] ,

[0020] The above decomposition formula satisfies the energy conservation constraint, where, Indicates that the i-th sensor of the drone is in Atmospheric data collected at all times This is a true atmospheric signal. For motion-coupled disturbances, This refers to the dynamic disturbance of the machine body. This is due to atmospheric turbulence disturbance;

[0021] Before performing disturbance decoupling, it is necessary to consider the impact of the differences in the installation positions of multiple sensors on the measurement results. Sensors at different positions are subject to different degrees of interference from the rotor airflow. The influence of sensor position is mainly reflected in the dynamic disturbance of the airframe. In this embodiment, during the system initialization phase, the installation positions of all atmospheric sensors are calibrated, and their three-dimensional position coordinates in the airframe coordinate system are recorded. Define the location influence factor This is used to quantify the sensitivity of the sensor installation location to dynamic disturbances of the aircraft; as an optional implementation, a position influence factor... It can be represented as:

[0022] ,

[0023] in, The three-dimensional position coordinates of the aerodynamic disturbance source of the UAV rotor are given. It is a tiny constant. The larger the value, the closer the sensor is to the source of dynamic disturbance of the machine body, and the stronger it is affected by dynamic disturbance of the machine body;

[0024] Furthermore, the atmospheric data undergoes time-domain trend feature analysis and frequency-domain statistical feature analysis, respectively. Then, by constructing a time-frequency distribution matrix driven by dual operating conditions, multi-scale perturbation decomposition is performed on the atmospheric data sequence to generate multi-scale perturbation vectors, specifically including:

[0025] The atmospheric data sequence is subjected to time-domain trend feature analysis. Specifically, an adaptive sliding window regression model is used to fit the atmospheric data sequence to low frequencies. A sliding time window L is set, and local weighted regression is used within each sliding time window L to obtain the low-frequency trend curve and output the low-frequency trend component. This is used to reflect the low-frequency, slowly varying characteristics of real atmospheric changes; subsequently, the rate of change of the signal is calculated using first-order difference on the residual signal after trend removal. ;

[0026] Further frequency domain feature analysis was performed on the atmospheric data sequence. Specifically, a fast Fourier transform was performed on the atmospheric data sequence to obtain the signal spectrum. Peak frequency sets were extracted through peak frequency detection and energy density statistics. and the corresponding narrowband energy percentage It is used to identify narrowband periodic frequency regions in the observed signal. When the proportion of narrowband energy exceeds the preset proportion threshold, it indicates that there is a periodic structure in the signal. This periodic structure is usually related to rotor aerodynamic disturbance or airframe vibration.

[0027] Based on the aforementioned time-domain trend characteristics and frequency-domain statistical characteristics, the atmospheric data sequence is initially coarsely segmented, including the low-frequency trend component. It mainly corresponds to real atmospheric signals, while the narrow-band periodic frequency region corresponds to the dynamic disturbance of the aircraft.

[0028] A dual-condition modulation function is constructed by introducing flight state information to dynamically modulate the time-frequency separation parameters of subsequent multi-scale disturbance decomposition, specifically including:

[0029] Based on the real-time acquired power system parameters, the dynamic parameters causing body vibration are extracted. To make parameters of different dimensions comparable, each parameter is normalized to obtain a dynamic condition vector. :

[0030] ,

[0031] in, This refers to the motor speed. This is the rate of change of motor speed, calculated using the motor speed at adjacent sampling times, reflecting how quickly the motor speed changes. Indicates the throttle output signal. The throttle output change rate is calculated using the throttle output signals at adjacent sampling times and is used to characterize the dynamic changes in power output. To assess the vibration intensity of a multi-rotor UAV, the inconsistent rotational speeds of different motors can also cause structural vibration. By calculating the speed deviations between multiple motors, the overall speed dispersion can be further calculated as an indicator of the vibration intensity. When the value is large, it indicates a significant difference in rotational speed between the rotors, which may cause strong vibrations in the airframe structure.

[0032] Based on dynamic operating condition vector A dynamic disturbance modulation function (DMC) for the airframe was constructed using a combination of a multiple linear regression model and a piecewise adaptive method. The weight coefficients of each component were obtained by calibration from historical flight data. Based on the DMC, the dynamic intensity index was calculated. ;

[0033] Based on real-time acquired flight state parameters, parameters related to the UAV's motion state are extracted, and these parameters are normalized to obtain a motion state vector. :

[0034] ,

[0035] in, This represents the flight speed of the drone at time t. Indicates flight acceleration. Indicates pitch angle, Indicates the roll angle. The rate of ascent is calculated from the rate of change of altitude by the altitude sensor over continuous time.

[0036] Based on motion condition vectors A motion disturbance modulation function (MDM) is constructed using a multiple linear regression model combined with a piecewise adaptive method. The weight coefficients of each parameter are obtained by calibration based on historical flight data. When a parameter of a certain operating condition exceeds its preset stability threshold, the weight coefficient of the corresponding parameter is automatically multiplied to adapt to extreme flight conditions. The motion intensity index is obtained based on the motion disturbance modulation function. ;

[0037] After obtaining the dual-condition modulation function, a short-time Fourier transform is performed on the atmospheric data sequence to obtain the time-frequency distribution matrix of the observed signal. The time-frequency distribution matrix is ​​then marked as a multi-channel time-frequency distribution matrix according to the channels, wherein the number of channels is consistent with the number of sensors. Based on the spatial layout relationship of each sensor, a position weighting factor is assigned to each channel.

[0038] According to motor speed Predicting the dominant frequency of dynamic disturbances in the machine body Combined with peak frequency set and the corresponding narrowband energy percentage Correction is performed on the predicted dominant frequency of the body dynamics disturbance, specifically, in the peak frequency set. Finding and predicting the dominant frequency of dynamic disturbances in the organism closest peak frequency At the same time, determine the peak frequency. Corresponding narrowband energy percentage Does it exceed a preset percentage threshold? If the condition is met, it indicates that the peak frequency is [value missing]. If it exhibits obvious periodic structural characteristics, then the peak frequency will be... As the main frequency of dynamic disturbance of the body Conversely, it indicates that the peak frequency is... It does not have obvious periodic structural characteristics, so the dominant frequency of the predicted dynamic disturbance of the organism is maintained. As the main frequency of dynamic disturbance of the body ;

[0039] Extracting the spectral bandwidth near the dominant frequency of the body dynamic disturbance ,in For bandwidth, determined by dynamic strength index Dynamic adjustment:

[0040] ,

[0041] in Based on bandwidth, The adjustment coefficient is determined by the main frequency and spectral bandwidth of the body dynamic disturbance. The time-frequency region of the body dynamic disturbance is located in the time-frequency distribution matrix. An inverse short-time Fourier transform is performed on this time-frequency region to obtain the body dynamic disturbance signal, which is then combined with the position influence factors of each sensor. Obtain the corresponding body dynamic disturbance signal ;

[0042] After stripping away the dynamic disturbance signal from the machine body, the signal change rate is used. Exercise intensity index Identify abrupt broadband energy regions in the time-frequency distribution matrix; specifically, calculate motion intensity indices using motion perturbation modulation functions. Based on statistical analysis of historical flight data, the range thresholds for the motion intensity index are determined, and a first intensity threshold is set. Second intensity threshold ,when When the current flight motion state is determined to be at a low intensity level, At that time, the current flight motion state is determined to be of medium intensity level. When the current flight motion state is determined to be of high intensity, corresponding target mutation detection thresholds and target spectral bandwidth thresholds are preset for different flight motion state levels.

[0043] Based on this, the rate of change of the signal is utilized Detect the abrupt change time interval in the signal, when When the target mutation detection threshold corresponding to the current intensity level is greater than the threshold value, a mutation behavior is determined to exist in that time period. The spectral bandwidth of the motion coupling disturbance is calculated within the current mutation time interval. When the spectral bandwidth is greater than the target spectral bandwidth threshold value corresponding to the current intensity level, the current time-frequency region is determined to be a motion coupling disturbance region. A time-frequency mask is constructed for the identified time-frequency region, and the time-frequency mask is applied to the original time-frequency distribution matrix to obtain the time-frequency region of the motion coupling disturbance. Subsequently, an inverse short-time Fourier transform is performed on this time-frequency region to obtain the motion coupling disturbance signal. ;

[0044] After separating the dynamic disturbance and motion-coupled disturbance signals from the body dynamics, the atmospheric turbulence disturbance signal is further separated by calculating the spectral bandwidth and spectral entropy characteristics of the remaining time-frequency distribution matrix. Specifically, the low-frequency trend component obtained from the aforementioned time-domain analysis is used. As an initial estimate of the real atmospheric signal, the remaining signal is calculated relative to... The residual is used to locate the time-frequency region of the residual in the time-frequency distribution matrix, and the effective spectral bandwidth and spectral entropy characteristics of the time-frequency region are calculated to characterize the coverage range and uniformity of energy distribution of the signal energy in the frequency domain, respectively. When the effective spectral bandwidth is greater than the effective spectral bandwidth threshold and the spectral entropy value is greater than the spectral entropy threshold, it indicates that the signal exhibits a high spectral entropy distribution over a wide frequency range. The signal is then determined to have turbulent characteristics and is identified as an atmospheric turbulence disturbance region. An inverse short-time Fourier transform is performed on this time-frequency region to obtain the atmospheric turbulence disturbance signal. ;

[0045] Further, smoothing constraints and low-rank constraints are introduced to optimize the structure of the remaining signal in both the time and channel dimensions. This suppresses high-frequency noise or random disturbances remaining during the disturbance separation process, while maintaining the continuous variation trend of the real atmospheric signal and the correlation structure between multiple parameters, thus obtaining the real atmospheric signal. .

[0046] Furthermore, after completing the multi-scale perturbation decoupling, the multi-scale perturbation vector corresponding to each observation parameter can be obtained, including the real atmospheric component, the motion-coupled perturbation component, the airframe dynamic perturbation component, and the atmospheric turbulence perturbation component. These multi-scale perturbation components are input into a pre-constructed anomaly map detection model to achieve real-time anomaly detection, anomaly type determination, and anomaly root cause localization for UAV atmospheric data. The anomaly map detection model includes an observation layer map and an anomaly layer map. The observation layer map is used to characterize the correlation between the atmospheric parameters collected by the UAV and their perturbation components, while the anomaly layer map is used to characterize the anomaly type and anomaly root cause, specifically including:

[0047] The observation layer map uses various parameters from atmospheric observation data as nodes and decoupled multi-component signals as node features. It captures the dynamic correlations between atmospheric parameters through graph convolution and multi-head attention mechanisms. Specifically, the observation layer map... Represented as ,in This represents the set of observation layer nodes, directly corresponding to atmospheric observation parameters. Each node includes a multi-dimensional feature vector. The decoupled signals—the real atmospheric signal, the motion-coupled disturbance signal, the body dynamic disturbance signal, and the atmospheric turbulence disturbance signal—are represented as follows: , To establish the connecting edges between nodes in the observation layer, edges are established based on meteorological and physical relationships. The initial weights of the edges are initialized using historical atmospheric observation data statistics, resulting in the observation layer adjacency matrix. ;

[0048] To integrate the correlation information between observation layer nodes, a graph convolution feature propagation mechanism is introduced into the observation layer map. Through the graph convolution propagation mechanism, each observation layer node can integrate the information of its neighboring nodes, thereby reflecting the correlation of changes between atmospheric data parameters.

[0049] Considering the time-varying physical correlations between different atmospheric parameters, dynamic modeling of the influence relationships between nodes is required. A multi-head attention mechanism is introduced to achieve adaptive aggregation of node features. For each node in the observation layer map, the attention weights with neighboring nodes are first calculated. Then, a linear transformation of the node features is performed using the feature mapping matrix B, and the attention score between any node z and its neighboring node j is calculated. The attention scores of all neighboring nodes are normalized using the Softmax function to obtain normalized attention weights. Finally, new node features are obtained through feature aggregation.

[0050] To improve the ability to identify complex anomalies, K independent attention heads are set up. Each attention head independently calculates the attention weights and aggregated features between nodes. The output features of multiple attention heads are then concatenated to obtain the final node features.

[0051] The anomaly layer map Represented as ,in This is a set of anomaly layer nodes, covering various anomaly types, including but not limited to: temperature anomalies, air pressure anomalies, wind speed anomalies, humidity anomalies, dynamic system anomalies, motion state anomalies, and atmospheric turbulence anomalies. The edges connecting nodes in the anomaly layer represent the causal propagation logic and coupling relationships between anomaly types. The weights of the initial anomaly layer graph edges can be initialized based on historical flight anomaly data statistics to obtain the anomaly layer adjacency matrix. Then, based on graph convolutional network adaptive optimization, the node features of the anomaly layer nodes are obtained by mapping the node features output from the observation layer graph. The high-dimensional features of the observation layer nodes are passed to the anomaly layer through a learnable mapping function, realizing the correlation modeling between observation features and anomaly types.

[0052] Furthermore, historical evolution pattern constraints are introduced into the anomaly layer graph reasoning. Specifically, historical time series data is stored for each anomaly layer node, and it is decomposed into trend components, periodic components, and random residuals using a time series decomposition method. Based on the trend components and periodic components, the normal characteristic value of the current time node is predicted by a time series prediction model. Then, the deviation between the current observation value and the predicted value is calculated. If the deviation exceeds a preset deviation threshold, the node is considered to have an abnormal trend. Based on this, the original anomaly probability of the anomaly layer node is corrected to obtain the corrected anomaly probability.

[0053] After processing through graph reasoning and historical evolution pattern constraints, anomaly detection is achieved by calculating the activation level of nodes in the anomaly layer. The detection function is:

[0054] ,

[0055] in, Features of abnormal layer nodes Indicates the probability of an anomaly, if If the value exceeds the preset abnormal threshold, the abnormal layer node is determined to be activated. The abnormal layer node has an abnormal state at the current time and is marked as an abnormal activated node. Then, the linkage relationship between abnormal activated nodes is analyzed. If at least one of the neighboring nodes of an abnormal activated node is an abnormal activated node, it is determined that the abnormal activated node participates in and forms a compound abnormality.

[0056] Furthermore, the propagation paths between anomalously activated nodes are analyzed to infer the root cause of the anomaly, specifically including:

[0057] Based on the adjacency matrix of the anomaly layer Based on the node activation state, extract the directed path of anomaly propagation to obtain the anomaly propagation path;

[0058] For each abnormal propagation path, calculate its total weight. The larger the total weight, the higher the propagation probability of the abnormal propagation path.

[0059] An abnormal activation node without a preceding abnormal activation node in the propagation path is taken as a candidate root cause node. The root cause contribution of the candidate root cause node is calculated by combining the total weight of the abnormal propagation path and the abnormal probability corresponding to the abnormal activation node.

[0060] The candidate node with the highest contribution is selected as the core root cause node to complete the anomaly root cause localization.

[0061] Furthermore, the abnormal observation data undergoes dynamic correction processing. This dynamic correction is based on the anomaly type, propagation path, and root cause obtained from anomaly detection. The abnormal observation data is adaptively corrected. If a single anomaly is detected, correction is performed according to the anomaly type corresponding to the anomaly activation node. If a compound anomaly is identified, its propagation path and root cause are obtained, and targeted correction is performed based on the root cause, including:

[0062] When the abnormal root cause node is the abnormal dynamic disturbance of the body, the disturbance component attenuation correction method is used to attenuate and correct the dynamic disturbance component of the body.

[0063] When the abnormal root cause node is an abnormal motion state, the adaptive filtering and smoothing correction method is used to filter and smooth the motion disturbance component to correct it.

[0064] When the root cause node of the anomaly is an environmental anomaly such as turbulence anomaly or temperature anomaly, it reflects the real atmospheric environmental change. It is marked and output as effective atmospheric observation information. Specifically, the anomaly type, anomaly occurrence time and anomaly propagation path corresponding to the anomaly activation node are recorded and marked as a real atmospheric anomaly event.

[0065] After correcting each disturbance component, the corrected atmospheric observation data is obtained and output as a meteorological observation result to the subsequent meteorological analysis module.

[0066] On the other hand, the present invention provides a remote backtracking and fault diagnosis system for flight parameters of long-endurance unmanned aerial vehicles, including a data processing module, a disturbance decoupling module, an anomaly detection module and a dynamic correction module;

[0067] The data processing module acquires atmospheric observation data collected in real time during the flight of the UAV, performs time synchronization processing on the atmospheric observation data, and maps it to the ground reference coordinate system to generate an atmospheric data sequence.

[0068] The disturbance decoupling module performs time-domain and frequency-domain feature analysis on the atmospheric data sequence, introduces a time-frequency distribution matrix driven by dual-condition vectors to perform multi-scale disturbance decomposition on the atmospheric data sequence, and generates multi-scale disturbance vectors, which include real atmospheric signals, motion coupling disturbances, body dynamic disturbances and atmospheric turbulence disturbances.

[0069] The anomaly detection module maps the multi-scale perturbation vector to a preset anomaly graph detection model. The anomaly graph detection model includes an observation layer graph and an anomaly layer graph. It uses graph structure adaptive learning and multi-head attention mechanism to propagate and fuse node features, calculates the anomaly probability of anomaly layer nodes, identifies single anomalies and compound anomalies, and infers the root cause of anomalies based on the anomaly propagation path and anomaly activation nodes, and outputs the anomaly detection result and the root cause of anomalies.

[0070] The dynamic correction module performs targeted corrections on different observation anomalies based on the anomaly detection results and root causes, thereby obtaining corrected atmospheric observation data.

[0071] Compared with the prior art, the significant advantages of this invention are:

[0072] 1. Based on the UAV attitude angle information, a transformation relationship between the airframe coordinate system and the ground reference coordinate system is constructed to realize unified spatial coordinate mapping and attitude compensation of atmospheric observation data. At the same time, a dual-condition modulation function is constructed based on the UAV's dynamic state and flight state to dynamically modulate the time-frequency separation parameters of multi-scale disturbance decomposition. A time-frequency distribution matrix driven by dual-condition vector is introduced to perform multi-scale disturbance decomposition of atmospheric data, effectively identifying disturbance components from different sources in complex flight environments.

[0073] 2. Construct an anomaly detection model that includes observation layer graphs and anomaly layer graphs. In graph inference, node features are dynamically fused through graph convolutional propagation and multi-head attention mechanism. At the same time, the anomaly probability is corrected by combining historical evolution pattern constraints, thereby realizing the identification of single anomalies and compound anomalies and root cause localization, identifying the source of anomalies, and realizing corresponding dynamic correction. Attached Figure Description

[0074] Figure 1 This is a flowchart of a method for detecting anomalies in UAV atmospheric data based on spectral analysis.

[0075] Figure 2 The flowchart for generating multi-scale perturbation vectors in this invention is shown below;

[0076] Figure 3 This is a flowchart of the anomaly detection process in this invention;

[0077] Figure 4 This is a flowchart of an anomaly detection system for UAV atmospheric data based on graph analysis. Detailed Implementation

[0078] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0079] Example 1

[0080] like Figure 1 As shown, this invention discloses a method for detecting anomalies in UAV atmospheric data based on spectral analysis, comprising the following steps:

[0081] Acquire atmospheric observation data collected in real time during the flight of the UAV, perform time synchronization processing on the atmospheric observation data, and map it uniformly to the ground reference coordinate system to generate an atmospheric data sequence;

[0082] Time-frequency analysis is performed on the atmospheric data sequence to construct dynamic condition vectors and motion condition vectors. Based on the two types of condition vectors, the corresponding time-frequency separation parameters are calculated. A time-frequency distribution matrix driven by dual condition vectors is introduced to perform multi-scale perturbation decomposition on the atmospheric data sequence to generate multi-scale perturbation vectors. The multi-scale perturbation vectors include real atmospheric signals, motion coupling perturbations, body dynamic perturbations, and atmospheric turbulence perturbations.

[0083] The multi-scale perturbation vector is mapped to a preset anomaly map detection model, which includes an observation layer map and an anomaly layer map. The map structure adaptive learning and multi-head attention mechanism are used to propagate and fuse node features, calculate the anomaly probability of anomaly layer nodes, identify single anomalies and compound anomalies, and infer the root cause of anomalies based on the anomaly propagation path, and output the anomaly detection result and the root cause of anomalies.

[0084] Based on the anomaly detection results and root causes, targeted corrections are made to different observational anomalies to obtain corrected atmospheric observation data.

[0085] In this embodiment, a multi-source atmospheric sensor and a multi-source flight status sensor are integrated at a fixed position on the UAV fuselage to simultaneously collect atmospheric observation data and flight status data. The atmospheric observation data includes at least atmospheric parameters such as temperature, humidity, air pressure, wind speed, wind direction, and particulate matter concentration. The flight status data includes attitude parameters, navigation parameters, and power system parameters. The attitude parameters include three-axis attitude angles and corresponding angular velocities. The navigation parameters include the UAV's flight altitude, speed, and acceleration. The power system parameters include information such as motor speed and throttle output. In a specific implementation, the multi-source atmospheric sensor and the multi-source flight status sensor can be devices such as digital meteorological sensors, inertial measurement units (IMUs), GPS positioning modules, and fuselage accelerometers. For example, a miniature wind speed and direction sensor installed at the front of the fuselage can collect incoming wind speed and direction information, a temperature and humidity sensor in the middle of the fuselage can collect environmental temperature and humidity parameters, and the UAV's flight attitude and spatial position data can be obtained synchronously through the IMU and GPS modules.

[0086] Before the drone takes off, all types of sensors undergo unified initialization, including zero-point calibration and communication status detection, to ensure the consistency of measurement benchmarks and the stability of data links. During the drone's monitoring mission, various real-time atmospheric observation data and flight status data are collected synchronously based on the preset sampling frequency of each sensor.

[0087] Because different sensors have different sampling frequencies—for example, temperature and humidity sensors typically have a sampling frequency of 1Hz-10Hz, while barometric pressure sensors typically have a sampling frequency of 5Hz-20Hz—multi-source sensor data exhibits inconsistencies in sampling frequency and time asynchrony on the time axis. To ensure the temporal consistency of multi-source acquired data, this embodiment employs a reference time axis synchronization mechanism to perform time alignment processing on the multi-source acquired data. Specifically, sensor data with a higher sampling frequency and better data stability is selected from the multi-source atmospheric sensors as the reference time axis. For example, in wind field monitoring tasks, wind speed or barometric pressure sensor data can be selected as the reference time axis, and their timestamps are used as a unified time reference. After obtaining the parameters... After referencing the time axis, time calibration processing is performed on other sensor data. For sensor data with sampling frequencies lower than the reference time axis, linear interpolation is used for time resampling to generate matching data values ​​at the corresponding time nodes on the reference time axis. For data with sampling frequencies higher than the reference time axis, synchronous downsampling is used to ensure consistency with the reference time axis time series. In addition, abnormal jumps or lost data points that occur during the acquisition process are repaired or removed by interpolation to ensure the continuity and stability of the time series. After the above time synchronization processing, each reference timestamp corresponds to a complete set of atmospheric observation data and flight status data, resulting in time-calibrated atmospheric observation data.

[0088] Because the attitude of a drone continuously changes during flight, directly analyzing the raw vector data would introduce coupling errors between the drone's attitude changes and the direction of atmospheric flow, leading to spatial distortions in vector parameters such as wind direction and speed. To eliminate the coupled attitude effects caused by differences in sensor installation locations and drone attitude changes, this embodiment establishes a unified spatial reference coordinate system and performs attitude compensation and spatial mapping processing on the atmospheric observation data, specifically including:

[0089] Establish a body coordinate system and a ground reference coordinate system. The body coordinate system takes the center of the UAV body as the origin, with the X-axis pointing forward, the Y-axis pointing to the right side of the UAV body, and the Z-axis perpendicular to the body and upward. The ground reference coordinate system takes the projection point of the UAV's instantaneous center of gravity on the ground as the origin, with the X-axis pointing due north as located by the GPS module, the Y-axis pointing due east, and the Z-axis perpendicular to the ground and upward.

[0090] During flight, the three-axis attitude angles of the UAV, including roll angle, are acquired in real time. Pitch angle Yaw angle The attitude transformation matrix of the body coordinate system relative to the ground reference coordinate system is converted into a quaternion representation and then calculated in real time using a quaternion attitude calculation algorithm. The attitude calculation frequency is kept consistent with the IMU's sampling frequency to update the attitude transformation matrix in real time. ;

[0091] Using attitude transformation matrix By rotating and mapping the vector atmospheric data in the body coordinate system to the ground reference coordinate system, the directional deviation caused by attitude changes is eliminated, and atmospheric data under a unified spatial reference is obtained. For example, the wind speed vector measured in the body coordinate system is converted into a three-dimensional wind field vector in the ground reference coordinate system.

[0092] The wind direction angle is recalculated in the ground reference coordinate system. The wind direction angle is recalculated with the X-axis of the ground reference coordinate system, i.e., due north, as the reference. The wind direction angle is defined as zero degrees with due north and increasing in a clockwise direction, so as to ensure that the wind direction information can accurately reflect the real atmospheric flow direction and eliminate the direction deviation introduced by the UAV attitude change.

[0093] After the above time synchronization and spatial attitude compensation processes, atmospheric observation data under a unified time axis and a unified spatial reference coordinate system can be obtained, thus forming an atmospheric data sequence with consistent spatiotemporal semantics.

[0094] Furthermore, the spatiotemporally consistent atmospheric data sequence is subjected to low-pass filtering to remove high-frequency noise generated during sensor sampling. Let the atmospheric data collected by the i-th sensor at time t be... The preprocessed atmospheric data sequence was obtained. ;

[0095] To improve the accuracy and stability of atmospheric data collected by UAVs during flight, it is necessary to decouple the multi-source coupled disturbances in the observation signals. When a UAV performs atmospheric monitoring tasks, the atmospheric data collected by its onboard sensors not only contains real atmospheric change information but is also affected by various factors such as changes in flight motion state, power system operation, and external airflow environment, resulting in superimposed disturbances in the observation signals. Through systematic analysis of the UAV's flight mechanism and aerodynamic environment, the observation signals can be represented as a superposition of real atmospheric signals and multiple types of disturbance signals. In this embodiment, considering the characteristics of UAV low-altitude flight and complex operating scenarios, the observation signals are mainly affected by the superposition and coupling of three types of disturbances. The atmospheric data collected by the onboard sensors can be decomposed into four independent signals, represented as follows:

[0096] ,

[0097] The above decomposition formula satisfies the energy conservation constraint, where, Indicates that the i-th sensor of the drone is in Atmospheric data collected at all times This is a true atmospheric signal. For motion-coupled disturbances, This refers to the dynamic disturbance of the machine body. The disturbances are categorized into atmospheric turbulence disturbances and dynamic disturbances. Real atmospheric signals originate from changes in natural atmospheric conditions, such as the gradual distribution of temperature, air pressure, or wind speed in space. These signals typically exhibit relatively gentle changes and strong spatial continuity, displaying low-frequency, slowly varying characteristics in atmospheric data. Dynamic disturbances primarily arise from mechanical vibrations and rotor aerodynamic interference generated during the operation of the power system, such as motor rotational vibrations, rotor downwash, and propeller periodic oscillations. These disturbances typically manifest as relatively stable, narrow-frequency periodic oscillations, with frequencies closely related to the rate of change of motor speed and throttle. Motion-coupled disturbances are mainly caused by the UAV's flight operation status. Changes in motion state, such as acceleration, turning, or climbing, can cause changes in the local airflow inlet angle of the sensor, mainly manifested as broadband disturbances with abrupt changes. Atmospheric turbulence disturbances mainly originate from random aerodynamic phenomena in the external atmospheric environment, such as gusts, wind shear, and vortices, and are usually manifested as broadband non-periodic random fluctuations. Based on the above analysis, the atmospheric data sequence can be divided into a superposition structure of real atmospheric signals and three types of disturbance signals. Since the frequencies of these three types of disturbances overlap, single-domain decoupling methods are difficult to achieve accurate separation, and under extreme conditions, the degree of disturbance coupling will be further aggravated, resulting in a significant decrease in decoupling accuracy.

[0098] Furthermore, before performing disturbance decoupling, it is necessary to consider the impact of the differences in the installation positions of multiple source sensors on the measurement results. Sensors at different positions are affected by the rotor airflow to varying degrees. For example, sensors near the rotor area are more susceptible to the influence of the rotor downwash airflow, while sensors located in the nose area are closer to a free-flow state. The influence of sensor position is mainly reflected in the dynamic disturbance of the airframe. In this embodiment, during the system initialization phase, the installation positions of all atmospheric sensors are calibrated, and their three-dimensional position coordinates in the airframe coordinate system are recorded. Define the location influence factor This is used to quantify the sensitivity of the sensor installation location to dynamic disturbances of the aircraft; as an optional implementation, a position influence factor... It can be represented as:

[0099] ,

[0100] in, The three-dimensional position coordinates of the aerodynamic disturbance source of the UAV rotor are given. It is a tiny constant used to avoid the denominator being zero; The larger the value, the closer the sensor is to the source of dynamic disturbance of the machine body, and the stronger it is affected by dynamic disturbance of the machine body;

[0101] like Figure 2As shown, time-domain trend feature analysis and frequency-domain statistical feature analysis are performed on the atmospheric data respectively. Then, by constructing a time-frequency distribution matrix driven by dual operating conditions, multi-scale perturbation decomposition is performed on the atmospheric data sequence to generate multi-scale perturbation vectors, specifically including:

[0102] The atmospheric data sequence is subjected to time-domain trend feature analysis. Specifically, an adaptive sliding window regression model is used to perform low-frequency fitting on the original atmospheric data sequence. A sliding time window L is set, and local weighted regression is used within each sliding time window L to obtain the low-frequency trend curve and output the low-frequency trend component. This is used to reflect the low-frequency, slowly varying characteristics of real atmospheric changes; subsequently, the rate of change of the signal is calculated using first-order difference on the residual signal after trend removal. ;

[0103] Further frequency domain statistical feature analysis was performed on the atmospheric data sequence. Specifically, the signal spectrum was obtained by performing a fast Fourier transform on the atmospheric data sequence.

[0104] ,

[0105] in, To determine the frequency band energy distribution, a set of peak frequencies is extracted through spectral peak detection and energy density statistics. and the corresponding narrowband energy percentage It can identify narrowband periodic frequency regions in the signal. When the proportion of narrowband energy exceeds the preset proportion threshold, it indicates that there is a periodic structure in the signal. This periodic structure is usually related to rotor aerodynamic disturbance or airframe vibration. The preset proportion threshold is determined by statistical analysis of UAV model, sensor type and historical flight data, and can also be adaptively adjusted according to actual flight conditions.

[0106] Based on the aforementioned time-domain trend characteristics and frequency-domain statistical characteristics, the atmospheric data sequence is initially coarsely segmented, including the low-frequency trend component. It mainly corresponds to real atmospheric signals, while the narrow-band periodic frequency region corresponds to the dynamic disturbance of the aircraft.

[0107] Further, flight state information is introduced to construct a dual-condition modulation function, which dynamically modulates the time-frequency separation parameters of subsequent multi-scale disturbance decomposition, specifically including:

[0108] The UAV power system is modeled. Specifically, the airframe dynamic disturbances mainly consider the thrust fluctuations caused by rotor periodic disturbances and speed changes. For example, when the propeller rotates at a certain speed, its blades continuously cut the airflow, creating periodic aerodynamic pulsations around the airframe. In addition, changes in motor speed cause transient thrust fluctuations, resulting in airframe vibration. The power system's output thrust is also related to the throttle control signal. Based on the above analysis and real-time acquired power system parameters, the dynamic parameters causing airframe vibration are extracted. To make parameters of different dimensions comparable, each parameter is normalized to obtain the dynamic operating condition vector. :

[0109] ,

[0110] in, This refers to the motor speed. This is the rate of change of motor speed, calculated using the motor speed at adjacent sampling times, reflecting how quickly the motor speed changes. Indicates the throttle output signal. The throttle output change rate is calculated using the throttle output signals at adjacent sampling times and is used to characterize the dynamic changes in power output. To assess the vibration intensity of a multi-rotor UAV, the inconsistent rotational speeds of different motors can also cause structural vibration. By calculating the speed deviations between multiple motors, the overall speed dispersion can be further calculated as an indicator of the vibration intensity. When the value is large, it indicates a significant difference in rotational speed between the rotors, which may cause strong vibrations in the airframe structure.

[0111] Dynamic operating condition vector As an optional implementation method, a body dynamics disturbance modulation function is constructed using a multiple linear regression model combined with a dynamic weight adaptive strategy, expressed as:

[0112] ,

[0113] Among them, the weighting coefficient Calibrated using the least squares method based on historical flight data, and with the actual physical response of the power intensity index as the optimization objective, the output of the airframe dynamic disturbance modulation function is made consistent with the actual airframe vibration and wing disturbance intensity. Simultaneously, a dynamic threshold triggering mechanism is employed: when a certain power condition vector parameter exceeds its preset stability threshold, the weight coefficient of the corresponding parameter is automatically multiplied. For example, the weight coefficients, after fitting with historical data, are respectively assigned the following values: [Values ​​for the motor speed term are missing from the original text]. The weight of the speed change rate term is 0.35. The weight of the throttle output term is 0.25. The weight of the throttle change rate term is 0.20. The weight of the body vibration intensity term is 0.15. The weight is 0.05, and the sum of all weight coefficients is 1. When the motor speed exceeds the preset stable threshold, the weight of the motor speed term is increased. The power intensity index is calculated based on the body dynamic disturbance modulation function. ;

[0114] To model the flight state of the UAV and characterize the formation mechanism of motion-coupled disturbances, specifically, during UAV flight, in complex conditions such as acceleration, turning, or climbing, the airflow direction and local airflow velocity of the onboard sensors change, leading to changes in the measured airflow velocity. Based on real-time acquired flight state parameters, parameters related to the UAV's motion state are extracted, and these parameters are normalized to obtain the motion state vector. :

[0115] ,

[0116] in, This represents the flight speed of the drone at time t. Indicates flight acceleration. Indicates pitch angle, Indicates the roll angle. The rate of ascent is calculated from the rate of change of altitude by the altitude sensor over continuous time.

[0117] Based on motion condition vectors As an optional implementation method, a motion perturbation modulation function is constructed using a multiple linear regression model combined with a dynamic weight adaptive strategy, expressed as follows:

[0118] ,

[0119] Among them, the weighting coefficient Based on historical flight data calibrated using the least squares method, the optimization objective is the actual physical response to motion disturbance intensity. A threshold-triggered strategy is employed: when a motion condition vector parameter exceeds its preset stability threshold, the weight coefficient of the corresponding parameter is automatically multiplied, while the weight coefficients of other terms are correspondingly reduced. This enhances the representation capability of extreme flight conditions. For example, the fitted values ​​for each weight coefficient are: weight of flight speed term... The weight of the flight acceleration term is 0.30. The weight of the pitch angle term is 0.35. The weight of the roll angle term is 0.15. The weight of the climb rate item is 0.15. The value is 0.05, and the sum of the weighting coefficients is 1; the motion intensity index is obtained based on the motion perturbation modulation function. ;

[0120] After obtaining the dual-condition modulation function, a short-time Fourier transform is performed on the atmospheric data sequence to obtain the time-frequency distribution matrix of the observed signal. This time-frequency distribution matrix is ​​then labeled as a multi-channel time-frequency distribution matrix according to the number of channels. The number of sensors matches the number of channels, with each atmospheric observation sensor corresponding to an independent channel, as shown below:

[0121] ,

[0122] in, The short-time Fourier transform result of the i-th channel is a two-dimensional function of time t and frequency f, representing the time-frequency energy distribution of the signal at time t and frequency f. For sliding window functions, such as the Hanning window, This represents the amount of shift of the sliding window function on the time axis. The atmospheric data sequence of the i-th channel, As a continuous-time variable, it represents the time-dimensional sampling points of the signal. For the complex exponential kernel of the Fourier transform;

[0123] To achieve precise separation of dynamic disturbances in the machine body by incorporating dynamic operating condition information, based on motor speed... Predicting the dominant frequency of dynamic disturbances in the machine body :

[0124] ,

[0125] Where tk is the proportionality coefficient, used to characterize the proportional relationship between the motor speed and the main frequency of the machine body's dynamic disturbance. It is calibrated through fitting a large amount of power system test data. For example, the proportionality coefficient tk is 0.02Hz / (r / s), combined with the peak frequency set. and the corresponding narrowband energy percentage Correction is performed on the predicted dominant frequency of the body's dynamic disturbances, specifically, within the peak frequency set. Finding and predicting the dominant frequency of dynamic disturbances in the organism closest peak frequency Simultaneously determine the narrowband energy percentage corresponding to the peak frequency. Does it exceed a preset percentage threshold? If the condition is met, it indicates that the peak frequency is [value missing]. If it exhibits obvious periodic structural characteristics, then the peak frequency will be... As the main frequency of dynamic disturbance of the body Conversely, it indicates that the peak frequency is... It does not have obvious periodic structure characteristics, so the dominant frequency of the predicted dynamic disturbance of the machine body is maintained. As the main frequency of dynamic disturbance of the body And in the main frequency of the body dynamic disturbance Nearby spectrum bandwidth extraction ,in For bandwidth, determined by dynamic strength index Dynamic adjustment:

[0126] ,

[0127] in Based on bandwidth, The adjustment coefficient is determined by the main frequency and spectral bandwidth of the body dynamic disturbance. The time-frequency region of the body dynamic disturbance is located in the time-frequency distribution matrix. An inverse short-time Fourier transform is performed on this time-frequency region to obtain the body dynamic disturbance signal, which is then combined with the position influence factors of each sensor. Obtain the body dynamic disturbance signals corresponding to each channel. ;

[0128] After stripping away the dynamic disturbance signal from the machine body, the signal change rate is used. Exercise intensity index Identify abrupt broadband energy regions in the time-frequency distribution matrix; specifically, calculate motion intensity indices using motion perturbation modulation functions. Based on statistical analysis of historical flight data, the range thresholds for the motion intensity index are determined, and a first intensity threshold is set. Second intensity threshold ,when When the current flight motion state is determined to be at a low intensity level, At that time, the current flight motion state is determined to be of medium intensity level. When the current flight motion state is determined to be of high intensity, corresponding target mutation detection thresholds and target spectral bandwidth thresholds are preset for different flight motion state levels.

[0129] Based on this, the rate of change of the signal is utilized Detect the abrupt change time interval in the signal, when When the target mutation detection threshold corresponding to the current intensity level is greater than the threshold value, a mutation behavior is determined to exist in that time period. The spectral bandwidth of the motion coupling disturbance is calculated within the current mutation time interval. When the spectral bandwidth is greater than the target spectral bandwidth threshold value corresponding to the current intensity level, the current time-frequency region is determined to be a motion coupling disturbance region. A time-frequency mask is constructed for the identified time-frequency region, and the time-frequency mask is applied to the time-frequency distribution matrix to obtain the time-frequency region of the motion coupling disturbance. Subsequently, an inverse short-time Fourier transform is performed on this time-frequency region to obtain the motion coupling disturbance signal. ;

[0130] After separating the dynamic disturbance and motion-coupled disturbance components of the aircraft, the remaining signal mainly contains real atmospheric signals and atmospheric turbulence disturbances. Real atmospheric parameters change relatively smoothly over short time scales, while atmospheric turbulence disturbances typically exhibit a relatively dispersed broad-spectrum energy distribution in the frequency domain. Therefore, this embodiment separates atmospheric turbulence disturbances by calculating the spectral bandwidth and spectral entropy characteristics of the remaining time-frequency distribution matrix. Specifically, the low-frequency trend component obtained from the aforementioned time-domain analysis is used. As an initial estimate of the real atmospheric signal, the remaining signal is calculated relative to... The residual is used to locate the time-frequency region of the residual in the time-frequency distribution matrix. The effective spectral bandwidth and spectral entropy characteristics of this time-frequency region are calculated to characterize the coverage range and uniformity of the signal energy distribution in the frequency domain, respectively. For example, based on historical flight data, the effective spectral bandwidth threshold is 15Hz and the spectral entropy threshold is 0.7. When the effective spectral bandwidth is greater than the effective spectral bandwidth threshold and the spectral entropy value is greater than the spectral entropy threshold, it indicates that the signal exhibits a high spectral entropy distribution over a wide frequency range. The signal is then determined to have turbulent characteristics and is identified as an atmospheric turbulence disturbance region. An inverse short-time Fourier transform is performed on this region to obtain the atmospheric turbulence disturbance signal. ;

[0131] After separating the components of body dynamics disturbance, motion coupling disturbance, and atmospheric turbulence disturbance, the remaining signal mainly consists of the real atmospheric signal and a small amount of residual noise. Since the real atmospheric environment typically exhibits continuity and smoothness, and there is physical correlation between multi-channel atmospheric observation parameters, smoothing constraints and low-rank constraints are further introduced to optimize the structure of the remaining signal. This suppresses high-frequency noise or random disturbances remaining from the disturbance separation process, while maintaining the continuous variation trend of the real atmospheric signal and the multi-channel correlation structure. As an optional implementation method, the smoothing constraint uses first-order difference. The regularization term constrains the continuous temporal variation characteristics of the real atmospheric signal, and corresponds to the penalty coefficient. Setting it to 0.01, the low-rank constraint uses a nuclear norm regularization term to constrain multi-channel correlation structures, with a corresponding penalty coefficient. By setting the value to 0.005 and solving the optimization problem with the aforementioned double regularization term, the remaining signal is structurally optimized to obtain the true atmospheric signal. .

[0132] like Figure 3As shown, after completing the multi-scale perturbation decoupling, the multi-scale perturbation vector corresponding to each observation parameter can be obtained, containing four independent signals: the true atmospheric component, the motion-coupled perturbation component, the airframe dynamic perturbation component, and the atmospheric turbulence perturbation component. To detect anomalies and identify anomaly sources in atmospheric observation data in complex flight environments, this embodiment constructs a dynamic weighted anomaly map detection model to achieve real-time anomaly detection, anomaly type determination, and anomaly root cause localization of UAV atmospheric data. The anomaly map detection model includes an observation layer map and an anomaly layer map. The observation layer map is used to characterize the correlation between the atmospheric parameters collected by the UAV and their perturbation components, while the anomaly layer map is used to characterize the anomaly type and anomaly root cause, specifically including:

[0133] The observation layer map uses various parameters of atmospheric data as nodes and the decoupled multi-component signals as node features. It captures the dynamic correlations between atmospheric parameters through graph convolution and multi-head attention mechanisms. Specifically, the observation layer map... Represented as ,in This represents the set of observation layer nodes, directly corresponding to atmospheric observation parameters. Each node includes a multi-dimensional feature vector, which consists of the decoupled real atmospheric signal, motion-coupled disturbance signal, body dynamic disturbance signal, and atmospheric turbulence disturbance signal. Represented as: , To establish the connecting edges between observation layer nodes, edges between nodes are established based on meteorological physical relationships. The initial weights of the edges are initialized using historical atmospheric observation data statistics to obtain the observation layer adjacency matrix. ;

[0134] To integrate the correlation information between observation layer nodes, this embodiment introduces a graph convolution feature propagation mechanism in the observation layer graph. The update formula for node features is as follows:

[0135] ,

[0136] in, This represents the feature matrix of the nodes in the l-th layer. Represents the learned weight matrix. For activation function, The normalized adjacency matrix is ​​obtained by processing the normalized adjacency matrix. With node feature matrix The product of the two elements is used to apply an activation function, and combined with the graph convolution kernel, node feature aggregation and hierarchical propagation of the observation layer map are achieved to capture the graph structure correlation features of atmospheric data. The calculation is expressed as:

[0137] ,

[0138] in for diagonal matrix, The identity matrix and the adjacency matrix are... Based on the meteorological and physical correlations and disturbance propagation logic among the atmospheric observation parameters of UAVs, and through the graph convolution propagation mechanism described above, each observation layer node can fuse information from its neighboring nodes, thereby reflecting the changing correlations between atmospheric data parameters.

[0139] Furthermore, considering the time-varying physical correlations between different atmospheric parameters, dynamic modeling of the influence relationships between nodes is required. This embodiment introduces a multi-head attention mechanism to achieve adaptive aggregation of node features. The specific process includes:

[0140] For each node in the observation layer graph, the node features are first linearly transformed using the feature mapping matrix B, and the attention score between any node z and its neighboring node j is calculated. A higher score indicates a stronger correlation between the two nodes, expressed as:

[0141] ,

[0142] in The attention weight vector, obtained through backpropagation iterative learning, is used to fit the nonlinear relationships between node features. and Let z be the feature vectors of node z and its neighbor node j, respectively. For feature splicing operations, The activation function is used, and the attention scores of all neighboring nodes are normalized using the Softmax function to obtain the normalized attention weights of node z and neighboring node j. Finally, the new node features under single-head attention are obtained through weighted aggregation.

[0143] To improve the ability to identify complex anomalies, K independent attention heads are set up. Each attention head independently calculates the attention weights and aggregated features between nodes. Then, the output features of multiple attention heads are concatenated to obtain the final node features.

[0144] ,

[0145] in For the final node features, This represents the feature aggregation result for the k-th attention head.

[0146] Through the multi-head attention mechanism, node association information at different scales can be captured simultaneously, making it easier to identify the linkage between multiple nodes. For example, in actual flight observation scenarios, if the wind speed node and the wind direction node both show abnormal fluctuations, while the air pressure node changes little, the attention mechanism will automatically increase the weight between the wind speed node and the wind direction node, thereby identifying potential composite anomalies.

[0147] The anomaly layer map Represented as ,in This is a set of anomaly layer nodes, covering various anomaly types, including but not limited to: temperature anomalies, air pressure anomalies, wind speed anomalies, humidity anomalies, dynamic system anomalies, motion state anomalies, and atmospheric turbulence anomalies. The edges connecting nodes in the anomaly layer represent the causal propagation logic and coupling relationships between anomaly types. The weights of the initial anomaly layer graph edges can be initialized based on historical flight anomaly data statistics to obtain the anomaly layer adjacency matrix. The node features of the anomaly layer nodes are obtained by feature mapping output from the observation layer graph. The high-dimensional features of the observation layer nodes are passed to the anomaly layer through a learnable mapping function, thereby realizing the correlation modeling between observation features and anomaly types.

[0148] Furthermore, to improve the reliability of anomaly detection, this embodiment introduces historical evolution pattern constraints in the anomaly layer graph reasoning. Specifically, historical time series data is stored for each anomaly layer node, and it is decomposed into trend components, periodic components, and random residuals using a time series decomposition method. Based on the trend components and periodic components, the normal characteristic value of the current time node is predicted by a time series prediction model. Then, the deviation between the current observed value and the predicted value is calculated. If the deviation exceeds a preset deviation threshold, the anomaly layer node is considered to have an abnormal trend. On this basis, the original anomaly probability of the anomaly layer node is corrected to obtain the corrected anomaly probability. The preset deviation threshold is a safe deviation upper limit obtained based on a large amount of historical data under normal flight conditions. By introducing the constraints of historical evolution patterns, misjudgments caused by short-term disturbances can be avoided.

[0149] Optionally, the time-series prediction model adopts an autoregressive integral moving average model as its basic architecture, adapting to the linear time-series evolution law of the anomaly layer nodes. The model is expressed as follows: Where p is the order of the autoregressive term, d is the order of the difference term, and q is the order of the moving average term. For example, p=2, d=1, and q=1. The historical normal time series feature data of the abnormal layer nodes are used as the training set, and the mean square error between the predicted value and the observed value is used as the loss function. The model parameters are iteratively optimized by the gradient descent algorithm. The initial learning rate is set to 0.001, and the iteration is repeated for 500 rounds until the model converges.

[0150] After processing through graph reasoning and historical evolution pattern constraints, anomaly detection is achieved by calculating the activation level of nodes in the anomaly layer. The detection function is:

[0151] ,

[0152] in, Features of abnormal layer nodes Indicates the probability of an anomaly, if If the value exceeds the preset abnormal threshold, the abnormal layer node is determined to be activated and marked as an abnormal activated node. For example, the preset abnormal threshold ranges from 0.8 to 0, and can be adaptively adjusted according to the false judgment rate requirements of the actual application scenario. Then, the linkage relationship between abnormal activated nodes is analyzed. If at least one abnormal activated node exists among the neighboring nodes of an abnormal activated node, it is determined that the abnormal activated node participates in and forms a compound abnormality.

[0153] Specifically, the training process of the anomaly map detection model includes: acquiring historical observation data from UAVs, including atmospheric parameters such as temperature, air pressure, humidity, and wind speed, as well as corresponding flight state data; organizing the historical observation data in chronological order and segmenting it according to a preset time window length to generate a training sample sequence; performing multi-source perturbation decoupling to decompose the signal and obtain multi-scale perturbation signal data; constructing observation layer map training samples to obtain an initial training dataset; constructing anomaly label data, where anomaly labels are derived from anomaly events in historical flight records, such as abnormal flight attitude records and abnormal environmental records; generating anomaly label vectors for each training sample based on historical annotation information; the training objective is to make the anomaly probability predicted by the model consistent with the true label; using the cross-entropy loss function as the optimization objective; and dividing the dataset into training, validation, and test sets in a 7:2:1 ratio, where the training set is used for... Model parameter learning is performed using a validation set for hyperparameter tuning and early stopping control, while a test set is used to evaluate the final model performance. During training, mini-batch gradient descent is used to update model parameters. Training samples are divided into multiple batches of 64 according to a preset batch size, and each batch is input into the model for forward computation. The loss function is calculated based on the true anomaly labels, and the model parameters are updated using the gradient backpropagation algorithm. The Adam adaptive optimization algorithm is used to iterate the parameters. The initial learning rate is set to 0.001, and learning rate decay optimization is used with a decay rate of 0.8. When the validation set loss value no longer decreases after 5 consecutive training rounds, the learning rate is decayed according to the decay rate. The total number of training rounds is 200. If the validation set loss value does not decrease for 10 consecutive rounds, the early stopping mechanism is triggered. After training terminates, the optimal model parameters are saved, resulting in the trained anomaly detection model.

[0154] Furthermore, in order to locate the source of the anomaly, this embodiment also analyzes the propagation path between the anomaly-activated nodes, specifically including:

[0155] Based on the adjacency matrix of the anomaly layer Based on the activation state of the abnormal layer nodes, the directed path of abnormal propagation is extracted to obtain the abnormal propagation path;

[0156] For each abnormal propagation path, calculate its total weight. The larger the total weight, the higher the propagation probability of the abnormal propagation path.

[0157] An abnormally activated node without a preceding abnormally activated node in the propagation path is selected as a candidate root cause node. The root cause contribution of the candidate root cause node is calculated by combining the total weight of the abnormal propagation path and the abnormal probability of the abnormally activated node.

[0158] The candidate node with the highest contribution is selected as the root cause node to complete the anomaly root cause localization.

[0159] Furthermore, the abnormal observation data is dynamically corrected. The dynamic correction is based on the anomaly type, anomaly propagation path and anomaly root cause obtained from anomaly detection. Combined with the actual engineering needs of UAV atmospheric observation data, the abnormal observation data is adaptively corrected to ensure that the corrected data can both suppress false anomaly interference and retain the true atmospheric environment change characteristics.

[0160] Based on the anomaly detection results, single anomalies and compound anomalies are distinguished, and a differentiated correction strategy is adopted. If a single anomaly is detected, correction is performed according to the anomaly type corresponding to the anomaly activation node. If a compound anomaly is identified, its anomaly propagation path and root cause are obtained, and targeted correction is performed based on the root cause, including:

[0161] When the root cause of the abnormality is the abnormality of the machine body dynamic disturbance, it indicates that the power system, such as the motor speed or throttle output, has generated abnormal fluctuations. For this type of abnormality, the disturbance component attenuation correction method is used to attenuate and correct the machine body dynamic disturbance component.

[0162] When the root cause of the anomaly is an abnormal motion state, such as when a drone suddenly accelerates or changes its attitude drastically, it causes inertial interference to the atmospheric observation sensor, which in turn causes false fluctuations in the atmospheric observation data. To address this anomaly, an adaptive filtering and smoothing correction method is used to filter and smooth the motion disturbance component for correction.

[0163] When the root cause node of the anomaly is an environmental anomaly such as turbulence anomaly or temperature anomaly, the anomaly is not caused by the UAV's own system or motion disturbance, but reflects the real atmospheric environment change. It is marked and output as effective atmospheric observation information. Specifically, the anomaly type, anomaly occurrence time and anomaly propagation path corresponding to the anomaly activation node are recorded and marked as a real atmospheric anomaly event.

[0164] After correcting each disturbance component, the corrected atmospheric observation data is obtained and output as a meteorological observation result to the subsequent meteorological analysis module.

[0165] like Figure 4 As shown, this embodiment also provides an anomaly detection system for UAV atmospheric data based on spectral analysis, including a data processing module, a disturbance decoupling module, an anomaly detection module, and a dynamic correction module;

[0166] The data processing module acquires atmospheric observation data collected in real time during the flight of the UAV, performs time synchronization processing on the atmospheric observation data, and maps it to the ground reference coordinate system to generate an atmospheric data sequence.

[0167] The disturbance decoupling module performs time-domain and frequency-domain feature analysis on the atmospheric data sequence, introduces a time-frequency distribution matrix driven by dual-condition vectors to perform multi-scale disturbance decomposition on the atmospheric data sequence, and generates multi-scale disturbance vectors, which include real atmospheric signals, motion coupling disturbances, body dynamic disturbances and atmospheric turbulence disturbances.

[0168] The anomaly detection module maps the multi-scale perturbation vector to a preset anomaly graph detection model. The anomaly graph detection model includes an observation layer graph and an anomaly layer graph. It uses graph structure adaptive learning and multi-head attention mechanism to propagate and fuse node features, calculates the anomaly probability of anomaly layer nodes, identifies single anomalies and compound anomalies, and infers the root cause of anomalies based on the anomaly propagation path and anomaly activation nodes, and outputs the anomaly detection result and the root cause of anomalies.

[0169] The dynamic correction module performs targeted corrections on different observation anomalies based on the anomaly detection results and root causes, thereby obtaining corrected atmospheric observation data.

[0170] Example 2

[0171] This embodiment describes the overall process of the technical solution of the present invention in the specific application scenario of UAV low-altitude meteorological inspection. In this embodiment, the UAV performs the task of monitoring the urban low-altitude atmospheric environment. The predetermined flight task is as follows: take off automatically from the ground base station, climb to a cruising altitude of 100 meters, cruise along the preset route at a flight speed of 10 m / s, complete the continuous observation of urban temperature, humidity, air pressure, wind speed, wind direction and atmospheric turbulence intensity, etc., and finally return and land smoothly.

[0172] The UAV is equipped with multi-source atmospheric sensors and a flight status acquisition unit to simultaneously collect atmospheric observation data and UAV status data. The atmospheric observation data includes temperature, humidity, air pressure, wind speed, and wind direction; the flight status data includes three-axis attitude angles and corresponding angular velocities, flight altitude, flight speed, three-axis acceleration, motor speed, throttle output signal, and airframe vibration intensity. Before takeoff, the system completes sensor zero-point calibration and communication link testing to ensure consistent data acquisition benchmarks and stable transmission.

[0173] After the drone takes off, each sensor collects data in real time at different frequencies. The system selects the wind speed and direction data with the highest sampling frequency as the reference time axis, and uses linear interpolation to complete time-series resampling of the temperature, humidity and air pressure data collected at low frequencies. It also performs interpolation to repair abnormal jumps and missing data points, thus obtaining time-aligned atmospheric observation data.

[0174] Simultaneously, a body coordinate system and a ground reference coordinate system are established. The pitch, roll, and yaw angles of the UAV are acquired in real time, converted into quaternions, and the attitude transformation matrix is ​​calculated. The wind speed and direction vectors in the body coordinate system are rotated and mapped to the ground reference coordinate system, with true north as the reference coordinate system. The wind direction angle is recalculated clockwise to eliminate the observation bias caused by the drone's attitude deflection, and finally, a spatiotemporally consistent atmospheric data sequence is generated.

[0175] An adaptive sliding window low-frequency fitting was performed on the atmospheric data sequence to obtain the low-frequency trend component, and the signal change rate was calculated by first-order difference. A fast Fourier transform was performed on the atmospheric data sequence to extract the peak frequency set and narrowband energy ratio statistics. The narrowband energy ratio threshold was preset to 30%, and the preliminary analysis of time-domain and frequency-domain characteristics was completed.

[0176] Real-time acquisition of UAV power system parameters is used to construct a normalized power condition vector, including motor speed, motor speed change rate, throttle output, throttle output change rate, and airframe vibration intensity. An airframe dynamic disturbance modulation function is constructed using a multiple linear regression model to calculate the power intensity index. Simultaneously, flight speed, acceleration, pitch angle, roll angle, and climb rate are acquired to construct a motion condition vector and motion disturbance modulation function, thereby obtaining the motion intensity index.

[0177] A short-time Fourier transform is performed on the atmospheric data sequence to generate a multi-channel time-frequency distribution matrix consistent with the number of sensors. Position influence factors are assigned based on the relative positions of the sensors and the rotor. The dominant frequency of the airframe dynamic disturbance is predicted based on the motor speed. The dominant frequency of the airframe dynamic disturbance is corrected by combining the spectral peak set and narrowband energy percentage statistics. The spectral bandwidth is dynamically adjusted by the power intensity index, and a base bandwidth is set. adjustment coefficient The body dynamic disturbance signal is obtained by locating the narrow-band periodic time-frequency region through the main frequency of the body dynamic disturbance and the frequency spectrum bandwidth, and then by inverse short-time Fourier transform.

[0178] Based on statistical analysis of historical flight data, the interval thresholds for the motion intensity index are determined, and a first intensity threshold is set. Second intensity threshold The flight status is divided into three levels: low, medium, and high. Target mutation detection thresholds and target spectral bandwidth thresholds are set for each level, as detailed below:

[0179] Low intensity level: Target mutation threshold is 0.5, target spectral bandwidth threshold is 8Hz;

[0180] Medium intensity level: Target mutation threshold is 1.0, target spectral bandwidth threshold is 12Hz;

[0181] High intensity level: Target mutation threshold is 1.5, target spectral bandwidth threshold is 18Hz;

[0182] The signal change rate and motion intensity indices are used to identify abrupt change time intervals, and further determine the time-frequency region of motion-coupled disturbances. A time-frequency mask is constructed and the motion-coupled disturbance signal is extracted. For the remaining time-frequency region, the effective spectral bandwidth and spectral entropy value are calculated. Based on historical flight data, the effective spectral bandwidth threshold is 15Hz and the spectral entropy threshold is 0.7. When the effective spectral bandwidth is greater than the effective spectral bandwidth threshold and the spectral entropy value is greater than the spectral entropy threshold, it indicates that the signal exhibits a high spectral entropy distribution over a wide frequency range. This signal is then determined to have turbulent characteristics and is identified as an atmospheric turbulence disturbance region, from which the atmospheric turbulence disturbance signal is extracted. Finally, smoothing constraints and low-rank constraints are introduced to optimize the remaining signal. The smoothing constraint uses first-order difference. Regularization term, corresponding to the penalty coefficient Setting it to 0.01, the low-rank constraint uses a nuclear norm regularization term, with a corresponding penalty coefficient. By setting it to 0.005, and solving the optimization problem with the above double regularization term, the remaining signal is structurally optimized to obtain the real atmospheric signal.

[0183] Multi-channel, multi-scale perturbation components are input into the anomaly map detection model. The observation layer map uses atmospheric observation parameters such as temperature, air pressure, wind speed, and humidity as nodes and four types of perturbation components as node features. An initial adjacency matrix is ​​constructed based on meteorological and physical correlations. Node features are propagated through graph convolution, and a multi-head attention mechanism is used to adaptively learn the dynamic correlations between nodes to capture the linkage features between multiple parameters.

[0184] The anomaly layer map uses temperature anomalies, air pressure anomalies, wind speed anomalies, dynamic system anomalies, motion state anomalies, and atmospheric turbulence anomalies as nodes. Connection edges are constructed based on the causal relationships of the anomalies, and the initial weights are obtained from historical anomaly data statistics. The node features of the anomaly layer are obtained by transforming the output features of the observation layer through a learnable mapping function and introducing historical evolution mode constraints. The ARIMA time series prediction model is used to predict the node feature sequence. The deviation between the observed value and the predicted value is compared with the preset deviation threshold to complete the anomaly probability correction and reduce the misjudgment of short-term disturbances.

[0185] The probability of anomalies in the anomaly layer nodes is calculated and compared with a preset anomaly threshold to determine the node activation state. The preset anomaly threshold is set to 0.85. When there is only a single abnormal activation node, it is determined as a single anomaly. When there are multiple abnormal activation nodes, it is determined as a compound anomaly. Anomaly propagation directed path is extracted based on the anomaly layer adjacency matrix and node activation state. The total weight of each path is calculated. Anomaly activation nodes without preceding abnormal activation nodes are selected as candidate root causes. The root cause contribution is calculated by combining the total weight of the anomaly propagation path and the anomaly probability. The node with the highest contribution is selected as the final anomaly root cause.

[0186] In this embodiment, the wind speed data of the UAV continued to fluctuate abnormally during the cruise phase. After model detection, the root cause of the abnormality was found to be abnormal dynamic disturbance of the air body, which corresponded to abnormal fluctuation of the speed of motor No. 3. This caused the vibration of the air body to intensify and was transmitted to the atmospheric observation sensor, forming a complex anomaly.

[0187] Targeted corrections are performed based on the root causes of anomalies: For anomalies in body dynamic disturbances, attenuation coefficients are constructed based on the intensity of the anomalies, and linear attenuation corrections are applied to the body dynamic disturbance components. Normal dynamic disturbance components are retained, and anomalous fluctuations are eliminated. After correction, the real atmospheric signal, the corrected body dynamic disturbance signal, the motion coupling disturbance signal, and the atmospheric turbulence disturbance signal are reconstructed by energy conservation to obtain high-quality corrected atmospheric observation data.

[0188] For real-world environmental anomalies such as atmospheric turbulence and sudden temperature changes, the system does not perform filtering corrections. Instead, it marks the anomaly type, occurrence time, and propagation path, and outputs them as real atmospheric anomalies along with the corrected data to the ground meteorological analysis module. This provides accurate and reliable data support for urban low-altitude wind field assessment, micro-meteorological inversion, and atmospheric environment monitoring.

[0189] After the mission is completed, the ground terminal can view the complete anomaly detection record, root cause localization results, and corrected data, realizing fully automated processing of UAV atmospheric observation from data acquisition, disturbance decoupling, anomaly identification, root cause localization to dynamic correction, which greatly improves the accuracy and reliability of atmospheric observation data under complex flight conditions.

[0190] The above description is merely a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiments. All technical solutions falling within the scope of the present invention's concept are within the scope of protection of the present invention. It should be noted that for those skilled in the art, any improvements and modifications made without departing from the principles of the present invention should also be considered within the scope of protection of the present invention.

Claims

1. A method for detecting anomalies in UAV atmospheric data based on spectral analysis, characterized in that, Includes the following steps: Acquire atmospheric observation data collected in real time during the flight of the UAV, perform time synchronization processing on the atmospheric observation data, and map it uniformly to the ground reference coordinate system to generate an atmospheric data sequence; Time-frequency analysis is performed on the atmospheric data sequence to construct dynamic condition vectors and motion condition vectors. Based on the two types of condition vectors, the corresponding time-frequency separation parameters are calculated. A time-frequency distribution matrix driven by dual condition vectors is introduced to perform multi-scale perturbation decomposition on the atmospheric data sequence to generate multi-scale perturbation vectors. The multi-scale perturbation vectors include real atmospheric signals, motion coupling perturbations, body dynamic perturbations, and atmospheric turbulence perturbations. The multi-scale perturbation vector is mapped to a preset anomaly map detection model, which includes an observation layer map and an anomaly layer map. The map structure adaptive learning and multi-head attention mechanism are used to propagate and fuse node features, calculate the anomaly probability of anomaly layer nodes, identify single anomalies and compound anomalies, and infer the root cause of anomalies based on the anomaly propagation path and anomaly activation nodes, and output the anomaly detection result and the root cause of anomalies. Based on the anomaly detection results and root causes, targeted corrections are made to different observational anomalies to obtain corrected atmospheric observation data.

2. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 1, characterized in that, The generation of atmospheric data sequences includes: The atmospheric observation data is time-aligned using a reference time axis synchronization mechanism to remove outliers and obtain time-calibrated atmospheric observation data. Spatial mapping processing is performed on the time-calibrated atmospheric observation data to establish the body coordinate system and the ground reference coordinate system. The three-axis attitude angles of the UAV are acquired in real time and converted into quaternion representations. The attitude transformation matrix of the body coordinate system relative to the ground reference coordinate system is calculated in real time using the quaternion attitude calculation algorithm. By using an attitude transformation matrix, the vector atmospheric observation data in the body coordinate system is rotated and mapped to the ground reference coordinate system. The wind direction angle is then recalculated in the ground reference coordinate system to obtain an atmospheric data sequence with consistent spatiotemporal semantics.

3. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 1, characterized in that, The generation of the multi-scale perturbation vector includes: Time-domain trend analysis was performed on the atmospheric data sequence to obtain the low-frequency trend component and signal change rate. Frequency-domain statistical analysis was performed to obtain the peak frequency set and narrowband energy proportion statistical characteristics. Dual-condition modulation functions are constructed by combining the UAV's power condition information and flight motion condition information; A time-frequency distribution matrix is ​​constructed based on the short-time Fourier transform. The time-frequency distribution matrix is ​​then labeled as a multi-channel time-frequency distribution matrix according to the channels, wherein the number of channels is consistent with the number of sensors. A position weight factor is assigned to each channel according to the spatial layout relationship of each sensor. In the time-frequency domain of each channel, the time-frequency regions corresponding to different disturbance types are identified, and inverse Fourier transforms are performed on the time-frequency regions of each disturbance to obtain the real atmospheric signal, motion-coupled disturbance signal, body dynamic disturbance signal and atmospheric turbulence disturbance signal of multiple channels.

4. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 3, characterized in that, The dual-mode modulation function includes: Real-time acquisition of UAV power system parameters, construction of power condition vector, the power system parameters include motor speed, motor speed conversion rate, throttle output signal, throttle output change rate and airframe vibration intensity, construction of airframe dynamic disturbance modulation function based on power condition vector, used to calculate power intensity index; The flight status parameters of the UAV are acquired in real time, and a motion condition vector is constructed. The flight status parameters include flight speed, flight acceleration, pitch angle, roll angle and climb rate. A motion disturbance modulation function is constructed based on the motion condition vector to calculate the motion intensity index.

5. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 4, characterized in that, Extraction of dynamic disturbance signals from the body includes: The dominant frequency of the body dynamic disturbance is predicted based on the motor speed, and the narrow-band periodic time-frequency region corresponding to the predicted dominant frequency of the body dynamic disturbance is located in the time-frequency distribution matrix. The dominant frequency of the body dynamic disturbance is determined by combining the peak frequency set and the statistical characteristics of the narrow-band energy ratio. The machine body dynamic disturbance bandwidth is calculated based on the dynamic strength index. The time-frequency region corresponding to the machine body dynamic disturbance is determined based on the main frequency of the machine body dynamic disturbance and the bandwidth of the machine body dynamic disturbance. An inverse short-time Fourier transform is performed on the time-frequency region to obtain the machine body dynamic disturbance signal.

6. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 5, characterized in that, The extraction of the motion-coupled disturbance signal and the atmospheric turbulence disturbance signal includes: The target mutation detection threshold is determined based on the signal change rate and motion intensity index. The mutation broadband energy region in the time-frequency distribution matrix is ​​identified. The spectral bandwidth of the running coupling disturbance is determined according to the motion intensity index. A time-frequency mask is constructed for the mutation broadband energy region. The motion coupling disturbance signal is obtained through inverse short-time Fourier transform. The spectral bandwidth and spectral entropy characteristics are calculated for the remaining time-frequency distribution matrix. When the effective spectral bandwidth is greater than the effective spectral bandwidth threshold and the spectral entropy value is greater than the spectral entropy threshold, the corresponding time-frequency region is determined to be an atmospheric turbulence disturbance region, and an inverse short-time Fourier transform is performed to obtain the atmospheric turbulence disturbance signal. By introducing smoothing constraints and low-rank constraints, the remaining signal is structurally optimized to obtain the true atmospheric signal.

7. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 6, characterized in that, The anomaly detection model includes an observation layer map and an anomaly layer map: The observation layer map uses each atmospheric observation parameter as a node, and uses the decoupled real atmospheric component, motion coupling disturbance component, body dynamic disturbance component and atmospheric turbulence disturbance component as node features. The connection edges between nodes are established according to the meteorological and physical correlation, and the dynamic correlation between nodes is adaptively learned through graph convolutional network and multi-head attention mechanism. The anomaly layer graph uses anomaly types as nodes, establishes node connection edges based on the causal propagation relationship between anomaly types, and receives high-dimensional features output by the observation layer graph through a mapping function to characterize anomaly types and their correlations.

8. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 7, characterized in that, The inferred root cause of the anomaly includes: The abnormal probability of the abnormal layer node is calculated based on graph adaptive learning. When the abnormal probability exceeds the preset abnormal threshold, the corresponding abnormal layer node is determined to be an abnormal activation node. Analyze the linkage relationship of abnormal activation nodes based on the connection relationship between abnormal layer nodes to identify single abnormalities and compound abnormalities; Based on the adjacency matrix of the anomaly layer and the activation state of the anomaly layer nodes, the directed path of anomaly propagation is extracted, the anomaly propagation path is obtained, and the total weight of each anomaly propagation path is calculated. An abnormal activation node without a preceding abnormal activation node in the abnormal propagation path is taken as a candidate root cause node. The root cause contribution of the candidate root cause node is calculated based on the abnormal probability of the candidate root cause node and the total weight of the corresponding abnormal propagation path, and the node with the highest contribution is selected as the abnormal root cause node.

9. The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in claim 1, characterized in that, Based on the anomaly detection results and root causes, targeted corrections are made to different observed anomalies, including: When the anomaly is an abnormality of the body dynamic disturbance, the body dynamic disturbance component is attenuated and corrected. When the anomaly is a motion state anomaly, the motion coupling disturbance component is filtered and smoothed for correction. When the anomaly is an environmental anomaly, the corresponding anomaly observation data will be marked as a real atmospheric anomaly event and output. After correcting the disturbance components, the atmospheric data is reconstructed and the corrected atmospheric observation data is output.

10. A UAV atmospheric data anomaly detection system based on spectral analysis, characterized in that, The method for detecting anomalies in UAV atmospheric data based on spectral analysis as described in any one of claims 1-9 includes: The data processing module acquires atmospheric observation data collected in real time during the flight of the UAV, performs time synchronization processing on the atmospheric observation data, and maps it to the ground reference coordinate system to generate an atmospheric data sequence. The disturbance decoupling module performs time-domain and frequency-domain feature analysis on the atmospheric data sequence, introduces a time-frequency distribution matrix driven by dual-condition vectors to perform multi-scale disturbance decomposition on the atmospheric data sequence, and generates multi-scale disturbance vectors. The anomaly detection module maps the multi-scale perturbation vector to a preset anomaly graph detection model, which includes an observation layer graph and an anomaly layer graph. It uses graph structure adaptive learning and multi-head attention mechanism to propagate and fuse node features, identify single anomalies and compound anomalies, infer the root cause of anomalies, and output anomaly detection results and anomaly root causes. The dynamic correction module performs targeted corrections on different observation anomalies based on the anomaly detection results and root causes, thereby obtaining corrected atmospheric observation data.

Images