Digital Twin-based UAV Safety Evaluation Method

Abstract

The embodiments of the present invention provide a method for drone safety evaluation based on digital twins. The steps include constructing a high-precision digital twin, generating a dynamic airflow disturbance field, optimizing flight control parameters through reinforcement learning, evaluating the risk probability, and importing the parameters into a physical drone through a hardware-in-the-loop system for verification. On the one hand, this method can cover most scenarios, with the scenario coverage rate in the research reaching more than 90%. At the same time, it is easy to reproduce scenarios under extreme conditions. On the other hand, the test cycle is significantly shortened. On the third hand, it is easy to automatically adjust and optimize the flight control system, and the efficiency is greatly improved.

Description

Technical Field

[0001] The present invention relates to the technical field of UAV safety, and particularly to a UAV safety evaluation method based on digital twin in a complex airflow environment. Background Art

[0002] A complex airflow environment is an important factor affecting UAV safety. Generally, for the UAV safety evaluation method in a complex airflow environment, the methods are to use field measurement and / or wind tunnel test, or to use special equipment given in documents such as CN109297673A, CN207197779U, CN110346109A, CN220472930U, etc. for detection or evaluation.

[0003] However, all of the above methods face the following three problems:

[0004] 1. Insufficient scene coverage. At the same time, it is difficult to reproduce scenes under extreme conditions. For example, it is difficult to generate a transient wind field with a speed exceeding 15 m / s.

[0005] 2. Long test cycle. For example, a single wind shear scene verification requires 3 - 5 days.

[0006] 3. Parameter tuning depends on experience: flight control engineers need to manually adjust more than 200 parameters.

[0007] Therefore, a new technical solution is needed to solve the above technical problems. Summary of the Invention

[0008] To this end, the present invention provides a UAV safety evaluation method based on digital twin to solve the above technical problems.

[0009] A UAV safety evaluation method based on digital twin includes the following steps:

[0010] Step S100: Construct a high-precision digital twin of the UAV, including an aerodynamic model, a flight control system mirror image, and a sensor noise model, where the aerodynamic model meets the simulation accuracy with Reynolds number Re≥1×10^5.

[0011] Step S200: Generate a dynamic airflow disturbance field in a virtual environment, and the disturbance field includes:

[0012] A gust model based on the Davenport spectrum, with the wind speed fluctuation range , as the reference wind speed;

[0013] Adopt the k-ω SST turbulence model to simulate the wind tunnel effect of the building complex, with a spatial resolution ≤0.1 m.

[0014] Step S300: Optimize the flight control parameters through a reinforcement learning algorithm, and perform Monte Carlo tests in the digital twin environment. When the attitude angle deviation δ satisfies: Trigger parameter adjustment, where is the pitch angle deviation, is the roll angle deviation;

[0015] Step S400: Establish a risk probability model to evaluate the probability of dangerous working conditions occurring:

[0016] ;

[0017] Among them, : extreme value distribution shape parameter; s: extreme event threshold; u: location parameter; σ: scale parameter;

[0018] At the same time, a preset threshold is preset, and after step S400, it includes comparing the probability of dangerous working conditions occurring with the preset threshold ; If ≤ ， Then execute step S500; If > and the current iteration number < , then return to step S300 and adjust the reinforcement learning reward function based on historical optimization data; If ≥ still does not satisfy ≤ , terminate the test and trigger a safety warning;

[0019] Step S500: Import the optimized flight control parameters into the physical unmanned aerial vehicle through the hardware-in-the-loop system to complete the safety verification.

[0020] Among them, in step S200, generating a dynamic airflow disturbance field in the virtual environment includes the following steps:

[0021] Step S2010: Construct an intelligent disturbance generation framework;

[0022] Step S2020: Basic wind field modeling, where the basic wind field is a gust model based on the Davenport spectrum;

[0023] Step S2030: Turbulence and building wake modeling;

[0024] Step S2040: Data-driven intelligent disturbance enhancement, including the following steps:

[0025] Construct a physical constraint generative adversarial network (PC-GAN):

[0026] Among them,

[0027] Generator G(T, Mt; θG): Adopts the U-Net structure to generate a high-resolution perturbation field ΔU;

[0028] Discriminator D(U; θD): Based on PatchGAN, verifies the Navier-Stokes equation residual:

[0029] ;

[0030] Total loss function:

[0031] ;

[0032] Among them, in the formula: T: Topographic elevation map; M t : Real-time meteorological data; θG : Trainable weight parameters of the generator U-Net; U: Composite three-dimensional wind speed field including the generated perturbation field and the basic wind field; θD : Discriminator network parameters; : Three-dimensional differential operator; : Kinematic viscosity; : Pressure gradient; : Material derivative; : Turbulence energy spectrum; : Kolmogorov theoretical energy spectrum; : Expected value; ΔU: High-resolution perturbation field output by the generator; : Adversarial loss weight; : Physical constraint weight; : Energy spectrum matching weight;

[0033] Step S2050: Real-time perturbation field fusion and injection.

[0034] Among them, in step S200, the flight control parameters are optimized through a reinforcement learning algorithm, and a Monte Carlo test is performed in the digital twin environment, which is carried out through the following steps:

[0035] Step S3010: Define the reinforcement learning state-action space;

[0036] Step S3020: Design a multi-objective reward function;

[0037] Step S3030: Implement the policy optimization algorithm;

[0038] Step S3040: Perform a Monte Carlo test, including the following steps:

[0039] Test scenario generation;

[0040] Randomly sample the wind shear parameter: wind speed gradient , duration ;

[0041] Among them, the wind speed gradient is uniformly distributed, and the duration is normally distributed;

[0042] Perturbation injection position: Randomly set on the UAV flight path = 20 perturbation areas, and the spacing follows a Poisson distribution;

[0043] Among them, the test execution and termination conditions:

[0044] Single test duration: = 180 s;

[0045] Termination condition: a. The attitude angle deviation δ ≥ 30°, lasting for more than 2 s;

[0046] b. Altitude loss ∆h ≥ 50 m;

[0047] c. The number of times the servo is saturated N ≥ 10;

[0048] The described termination condition means that the test stops when any one of a, b, or c is satisfied.

[0049] Step S3050: Adaptive sampling optimization;

[0050] Step S3060: Parameter adjustment trigger and verification.

[0051] Among them, the parameter adjustment trigger and verification in step S3060 are carried out according to the following rules:

[0052] Trigger logic: Real-time monitor the composite attitude deviation:

[0053] ;

[0054] When and it lasts for t ≥ 0.5 s, start parameter adjustment:

[0055] , ;

[0056] In the formula: : Saturation function; : Adjustment coefficient;

[0057] Online verification: Immediately execute N = 5 fast Monte Carlo tests after adjustment ( = 10 s);

[0058] Verification index: After adjustment: ≤ 25°; Overshoot: ≤ 15%;

[0059] If the verification is passed, update the flight control parameter library; otherwise, roll back and explore new strategies.

[0060] Among them, in step S3020, the multi-objective reward function is designed as: Composite reward function: , where

[0061] Stability reward:

[0062] ;

[0063] Trajectory tracking reward:

[0064] ;

[0065] Energy efficiency reward:

[0066] ;

[0067] In the formula:

[0068] : Time window length; : Angle deviation; : Number of trajectory points; , : Desired and actual position vectors; : Motor power; : Maximum motor power; : Reward weight coefficient, W1 = 0.6, W2 = 0.3, W3 = 0.1.

[0069] Among them, step S3050: Adaptive sampling optimization is carried out as follows: Dynamically adjust the sampling weight based on the test results:

[0070] ;

[0071] In the formula: : Operating condition parameter vector that triggers termination; : Parameters including wind speed gradient, duration, number of disturbance areas, etc., with a dimension of 3; : Initial weight; : Gaussian kernel bandwidth, = 0.1.

[0072] Among them, the policy optimization algorithm in step S3030 is the Proximal Policy Optimization (PPO) algorithm.

[0073] Among them, in step S500, the optimized flight control parameters are imported into the physical UAV through the hardware-in-the-loop system to complete the safety verification, which specifically includes the following steps:

[0074] Step S510: Perform CRC check and ECC encryption on the optimized flight control parameters to generate a data packet with a digital signature;

[0075] Step S520: Inject it into the flight control system through a dual-redundant CAN FD bus with a delay of ≤8 ms, and synchronously adopt the PTP clock protocol;

[0076] Step S530: Activate the parameters in segments and verify the stability margin in real time. If the servo saturation limit is exceeded, trigger a rollback;

[0077] Step S540: Execute step response, frequency sweep, and Monte Carlo stress tests, and calculate the comprehensive success rate;

[0078] Step S550: When a Level 3 anomaly is detected, roll back to the safe parameter version within 50 ms and start the forced landing procedure; the Level 3 anomaly is the case where the attitude angle > 45°.

[0079] Among them, in step S100, injecting noise into the digital twin is also included, and injecting noise into the digital twin includes the following steps:

[0080] Step S1110: Identify noise parameters based on Allan variance analysis:

[0081] Step S1120: Generate noise using a second-order autoregressive model;

[0082] Step S1130: Superimpose the noise signal on the flight control system mirror.

[0083] Among them, the noise is calibrated in the following manner:

[0084] Step S1210: Collect raw sensor data;

[0085] Step S1220: Quantitatively analyze the noise characteristics;

[0086] Step S1230: Calibrate the noise model parameters;

[0087] Step S1240: Verify the injection of the calibrated noise model into the digital twin. The KL divergence of the probability distribution between the digital twin and the output of the physical sensor is ≤0.1, and the time-domain correlation error is ≤15%.

[0088] Beneficial effects: The embodiment of the present invention provides a method for evaluating the safety of an unmanned aerial vehicle (UAV) based on digital twin. The steps include constructing a high-precision digital twin, generating a dynamic airflow disturbance field, optimizing flight control parameters through reinforcement learning, evaluating the risk probability, and verifying the parameters by importing them into the physical UAV through a hardware-in-the-loop system. On the one hand, this method can cover most scenarios, and the scenario coverage rate in the research reaches more than 90%. At the same time, it can easily reproduce scenarios under extreme conditions. On the other hand, the test cycle is greatly shortened. On the third hand, it is easy to automatically adjust and optimize the flight control system, and the efficiency is greatly improved. Description of the Drawings

[0089] Figure 1 It is a schematic flow chart of the method for evaluating the safety of an unmanned aerial vehicle based on digital twin according to the present invention. Detailed Embodiments

[0090] The embodiment of the present invention provides a method for evaluating the safety of an unmanned aerial vehicle based on digital twin, which is used to evaluate the safety of the unmanned aerial vehicle, especially to evaluate the safety of coping with wind shear in the wind field. At the same time, the unmanned aerial vehicle can be in the form of a helicopter, a multi-rotor unmanned aerial vehicle, or a fixed-wing form, which is not limited here.

[0091] As Figure 1 shown, the method for evaluating the safety of an unmanned aerial vehicle based on digital twin includes the following steps.

[0092] Step S100: Construct a high-precision digital twin of the unmanned aerial vehicle, including an aerodynamic model, a flight control system mirror image, and a sensor noise model, where the aerodynamic model meets the simulation accuracy of Reynolds number Re≥1×10^5.

[0093] Specifically, the aerodynamic model solves the unsteady Reynolds-averaged Navier-Stokes equation and is closed by using the k-ω SST turbulence model.

[0094] The unsteady Reynolds-averaged Navier-Stokes equation:

[0095] (1);

[0096] In the formula:

[0097] : Air density, taking the value under standard atmospheric conditions (such as 1.225 kg / m³ at sea level);

[0098] : Average velocity vector, Cartesian coordinate system components ( u , v , w ), simulation solution result;

[0099] : Aerodynamic viscosity, calculated by the Sutherland formula, 1.789×10 -5 Pa @ 15°C;

[0100] : Static pressure, standard atmospheric pressure at sea level 101325 Pa;

[0101] : Reynolds stress tensor, closed and solved by the k-ω SST turbulence model;

[0102] : Body force;

[0103] In addition, the grid resolution has the following requirements:

[0104] The first layer of the near-wall grid satisfies y + ≤ 1, and the grid size in the wake region Δx ≤ 0.05c (c is the chord length of the airfoil);

[0105] The sliding mesh method is used in the rotor region, and the number of overlapping grids at the interface between the rotating domain and the non-rotating domain ≥ 3 layers.

[0106] Dynamic aerodynamic effect modeling: The dynamic stall model uses the Leishman-Beddoes model:

[0107] (2);

[0108] In the formula:

[0109] : Dynamic normal force coefficient;

[0110] : Quasi-steady normal force coefficient;

[0111] : Normal force loss caused by separation;

[0112] : Time;

[0113] : Starting time of airflow separation;

[0114] : Dynamic stall characteristic time constant, = 6.0.

[0115] Furthermore, it is necessary to verify the aerodynamic model and make the error of the lift coefficient CL ≤ 3%, and the root mean square error RMSE of the pressure distribution Cp curve ≤ 0.015.

[0116] Specifically, the aerodynamic model can be verified through the following steps.

[0117] Step S1010: Define the reference operating condition.

[0118] Define the standard operating condition set:

[0119] Angle of attack range: ∈[-5°, 15°], with an interval of 1°;

[0120] Reynolds number: Re = 1.2×10 5 ;

[0121] Mach number: Ma = 0.15.

[0122] Step S1020: High-precision simulation calculation.

[0123] Solve the unsteady Reynolds-averaged Navier-Stokes equation (1) using the k-ω SST turbulence model;

[0124] Among them, the boundary conditions are:

[0125] Inlet: Velocity inlet, turbulence intensity I = 0.1%;

[0126] Outlet: Pressure outlet, static pressure p atm = 101325 Pa;

[0127] Wall: No-slip condition, the first layer of grid near the wall y + ≤1.

[0128] Step S1030: Verification of dynamic stall effect.

[0129] Specifically, the dynamic stall effect is verified through the following steps:

[0130] Step S1031: Perform pitch oscillation simulation, angle of attack change:

[0131] .

[0132] Step S1032: Compare the prediction of the Leishman-Beddoes model with the simulation results:

[0133] (3);

[0134] In the formula:

[0135] : Lift coefficient;

[0136] : Relative deviation of the maximum lift coefficient;

[0137] : Maximum value of the simulated lift coefficient

[0138] : The maximum value of the lift coefficient predicted by the Leishman-Beddoes model;

[0139] 5%: The maximum allowable relative error threshold.

[0140] Step S1040: Wind tunnel experiment comparison and verification.

[0141] Specifically, compare the simulated pressure coefficient with the measured data in the wind tunnel:

[0142] (4);

[0143] In the formula:

[0144] : Root mean square error;

[0145] : Pressure coefficient;

[0146] : Simulated pressure coefficient;

[0147] : Experimental pressure coefficient;

[0148] : Normalized chordwise position;

[0149] N: Number of measurement points;

[0150] : The maximum allowable RMSE threshold.

[0151] And, lift coefficient error verification:

[0152] (5);

[0153] In the formula: is the angle of attack;

[0154] : Experimental lift coefficient; feng

[0155] : Experimental lift coefficient at the angle of attack of.

[0156] Furthermore, the sensor noise model includes IMU noise, airspeed meter noise, and GPS positioning error.

[0157] Among them, the IMU noise includes:

[0158] Angular velocity random walk: ;

[0159] Accelerometer bias instability:

[0160] , where ;

[0161] The airspeed meter noise mentioned above includes:

[0162] Power spectral density of wind speed measurement noise:

[0163] ( ) (6);

[0164] In the formula:

[0165] : Power spectral density of wind speed measurement noise;

[0166] : Standard deviation of airspeed noise;

[0167] : Frequency;

[0168] : Cut-off frequency.

[0169] The GPS positioning error mentioned above includes:

[0170] Horizontal positioning error ellipse:

[0171] (7);

[0172] In the formula:

[0173] : Standard deviation of latitude and longitude positioning errors;

[0174] : Horizontal dilution of precision;

[0175] : Standard deviation of user range error.

[0176] Furthermore, the following steps are adopted to inject noise into the digital twin:

[0177] Step S1110: Noise parameter identification based on Allan variance analysis:

[0178] (Angle random walk coefficient) (8);

[0179] In the formula:

[0180] : Angle random walk coefficient;

[0181] : Fitting slope of the Allan variance curve in a specific time period.

[0182] Step S1120: Generate correlated noise using a second-order autoregressive model.

[0183] Step S1130: Superimpose the noise signal on the flight control system mirror image:

[0184] (9);

[0185] In the formula:

[0186] : Measured angular velocity value with noise;

[0187] : True angular velocity;

[0188] : Gaussian white noise;

[0189] : Time-varying term of angular velocity bias.

[0190] In addition, the sensor noise parameters need to be calibrated, and the calibration is carried out according to the following steps.

[0191] Step S1210: Sensor raw data acquisition.

[0192] Fix the UAV in a windless environment and continuously collect the raw outputs of each sensor:

[0193] IMU: Angular velocity ω x , ω y , ω z , acceleration a x , a y , a z ;

[0194] Air speed indicator: Indicated airspeed Uind, static pressure p s ;

[0195] GPS: Latitude and longitude (λ, φ), altitude h, sampling frequency ≥ 100Hz.

[0196] Step S1220: Noise characteristic quantization analysis.

[0197] Use Allan variance analysis to quantitatively analyze the noise characteristics.

[0198] Taking the gyroscope as an example, calculate the Allan variance at different correlation times τ:

[0199] (10);

[0200] In the formula:

[0201] : Allan variance;

[0202] : Time;

[0203] : Total number of sampling points;

[0204] : Number of interval points;

[0205] : Discretized angle sequence, which can be obtained by integrating the angular velocity of the gyroscope.

[0206] Extract noise parameters:

[0207] Angle random walk:

[0208] (11)

[0209] In the formula:

[0210] : Standard deviation of Allan variance at τ = 1s.

[0211] Bias instability:

[0212] (12)

[0213] In the formula:

[0214] : Bias instability;

[0215] : Standard deviation of Allan variance at = 100s.

[0216] Perform power spectral density estimation:

[0217] Calculate SU(f) using the Welch method to verify whether it conforms to the model:

[0218] (13);

[0219] In the formula:

[0220] : Measured airspeed noise power spectral density;

[0221] : White noise base power, corresponding to = 0.5m / s.

[0222] Step S1230: Calibrate noise model parameters.

[0223] Establish the IMU angular velocity noise model:

[0224] ;

[0225] In the formula:

[0226] : True angular velocity;

[0227] : Bias term;

[0228] : Standard deviation of Allan variance at τ = 100 s;

[0229] : Standard deviation of angular velocity white noise, .

[0230] Among them, the bias term b(t): a first-order Gaussian Markov process:

[0231] ;

[0232] : Time constant of the bias process;

[0233] : Gaussian white noise driving term;

[0234] GPS positioning error model parameter extraction:

[0235] Standard deviation of horizontal positioning error:

[0236] , ≤2.4 m;

[0237] : GPS horizontal dilution of precision;

[0238] Elevation error: ≤1.5× .

[0239] Step S1230: Digital twin noise injection verification.

[0240] Inject the calibrated noise model into the digital twin. The KL divergence of the probability distributions of the digital twin and the physical sensor output is ≤0.1, and the time-domain correlation error (maximum deviation of the autocorrelation function) is ≤15%.

[0241] Step S200: Generate a dynamic airflow disturbance field in the virtual environment. The disturbance field includes:

[0242] A gust model based on the Davenport spectrum, with the wind speed fluctuation range , is the reference wind speed;

[0243] The k-ω SST turbulence model is used to simulate the wind tunnel effect of the building complex, and the spatial resolution ≤ 0.1 m.

[0244] Specifically, it includes the following steps:

[0245] Step S2010: Construct an intelligent perturbation generation framework. Define a perturbation field generation architecture based on the fusion of physical constraints and data-driven, satisfying:

[0246] Input conditions: terrain elevation map T(x, y), real-time meteorological data Mt, and UAV pose qt;

[0247] Output form: three-dimensional wind speed field U(x, y, z, t) = [u, v, w]T.

[0248] Step S2020: Basic wind field modeling.

[0249] Among them, the gust model based on the Davenport spectrum:

[0250] ;

[0251] Among them, the spectral density function:

[0252] (16);

[0253] Among them, the boundary layer wind shear profile:

[0254] ;

[0255] In formulas (15), (16), and (17):

[0256] : the reference wind speed;

[0257] : Davenport spectrum constant;

[0258] : friction velocity;

[0259] : von Kármán constant, ;

[0260] : height above the ground;

[0261] : surface roughness length;

[0262] : total number of frequency components, usually N = 100;

[0263] : Frequency resolution,;

[0264] : Random phase angles uniformly distributed in [0, 2 π )

[0265] : For energy normalization, so that the variance of the synthesized wind speed is consistent with the spectral density integral;

[0266] Step S2030: Turbulence and building wake modeling.

[0267] Use the k-ω SST turbulence model to solve the Reynolds-averaged equations:

[0268] ;

[0269] Among them, the turbulent viscosity: ;

[0270] In addition, the building wake effect is generated:

[0271] ;

[0272] Among them, r = distance from the building, σ = 0.1L;

[0273] In equations (18) and (19):

[0274] : Kinematic viscosity;

[0275] : Turbulent viscosity;

[0276] : Horizontal distance from the building center, unit: m

[0277] : Turbulent kinetic energy;

[0278] : Specific dissipation rate;

[0279] : Free-stream velocity;

[0280] : Building drag coefficient, ;

[0281] : Wake diffusion coefficient; = 0.1L, where L is the building characteristic length.

[0282] Step S2040: Data-driven intelligent perturbation enhancement.

[0283] Specifically, construct a Physical Constraint Generative Adversarial Network (PC-GAN):

[0284] Among them,

[0285] Generator G(T, Mt; θG): Adopt a U-Net structure to generate a high-resolution perturbation field ΔU;

[0286] Discriminator D(U; θD): Based on PatchGAN, verify the Navier-Stokes equation residuals:

[0287] ;

[0288] Total loss function:

[0289] ;

[0290] In Eqs. (20) and (21):

[0291] T: Topographic elevation map. Specifically, the topographic elevation map T is a Digital Elevation Model (DEM) with a grid resolution ≤ 10 m and an elevation accuracy ≤ 0.5 m to cover the topographic data of the UAV test area; M t : Real-time meteorological data, which usually can include real-time wind speed, pressure gradient, air temperature, relative humidity, etc.; θG : Trainable weight parameters of the generator U-Net. Specifically, θG is the weight matrix of the convolutional layer, which can be optimized by minimizing the total loss function; U: Composite three-dimensional wind speed field containing the generated perturbation field and the basic wind field; θD : Discriminator network parameters. Specifically, the parameter θD can be updated through adversarial training, and the objective function is to maximize the classification accuracy of the real perturbation field and the generated perturbation field;

[0292] : Three-dimensional differential operator;

[0293] : Kinematic viscosity;

[0294] : Pressure gradient;

[0295] : Material derivative;

[0296] : Turbulent energy spectrum;

[0297] : Kolmogorov theoretical energy spectrum;

[0298] : Expected value, approximated through batch processing (batch size = 128);

[0299] ΔU: High-resolution perturbation field output by the generator (spatial resolution ≤ 0.1 m), superimposed on the basic wind field to form the composite wind field U;

[0300] : Adversarial loss weight; : Physical constraint weight; : Energy spectrum matching weight; = 1, = 10, = 0.1; Specifically, the discriminator D is of the PatchGAN architecture, outputting local discrimination results of 70×70, with a convolutional kernel size of 4×4 and the number of channels [64, 128, 256, 512].

[0301] Step S2050: Real-time perturbation field fusion and injection.

[0302] Assimilate real-time data through Ensemble Kalman Filtering (EnKF):

[0303] ;

[0304] In the formula:

[0305] : Forecast field;

[0306] : Analysis field;

[0307] : Kalman gain matrix;

[0308] : Observation operator;

[0309] : Analysis field;

[0310] : Sensor observation vector.

[0311] Generate the final composite perturbation field:

[0312] (23);

[0313] In the formula:

[0314] : Basic wind field, generated by the Davenport spectrum in formula (15);

[0315] : Turbulence component, solved by the k-ω SST model in formula (18);

[0316] : The perturbation increment generated by GAN, which is the output of PC-GAN in Equation (21);

[0317] : The EnKF assimilation correction term, in Equation (22) .

[0318] Step S300: Optimize the flight control parameters through a reinforcement learning algorithm, and perform a Monte Carlo test in the digital twin environment. When the attitude angle deviation δ satisfies: Trigger parameter adjustment when, where is the pitch angle deviation, is the roll angle deviation.

[0319] Specifically, the following steps are used to optimize the flight control parameters through a reinforcement learning algorithm and perform a Monte Carlo test in the digital twin environment.

[0320] Step S3010: Define the reinforcement learning state-action space.

[0321] Among them, the state space S:

[0322] The dynamic state of the UAV:

[0323] (24);

[0324] In the formula:

[0325] : Pitch, roll, and yaw angle deviations;

[0326] : Body angular rates (roll, pitch, yaw);

[0327] : Airspeed;

[0328] : Altitude;

[0329] : Mean of historical attitude deviations.

[0330] Environmental state: Current wind shear intensity, Turbulence integral scale Lt.

[0331] The action space A:

[0332] Flight control parameter adjustment amount:

[0333] (25);

[0334] PID parameter adjustment range:

[0335] (26);

[0336] Filter time constant:

[0337] ;

[0338] Constraint conditions:

[0339] Stability constraint: Phase margin ≥45°, Gain margin ≥6 dB.

[0340] Physical constraint: Servo deflection rate .

[0341] In equations (25), (26), and (27):

[0342] : Pitch-axis proportional coefficient adjustment amount. Specifically, ;

[0343] : Roll-axis integral coefficient adjustment amount. Specifically, ;

[0344] : Yaw-axis differential coefficient adjustment amount. Specifically, ;

[0345] : Nominal PID parameters;

[0346] : Filter time constant adjustment amount.

[0347] Step S3020: Design a multi-objective reward function.

[0348] Define a composite reward function: .

[0349] Stability reward:

[0350] ;

[0351] Trajectory tracking reward:

[0352] ;

[0353] Energy efficiency reward:

[0354] ;

[0355] In equations (28), (29), and (30):

[0356] : Time window length;

[0357] : Angle deviation;

[0358] : Number of trajectory points;

[0359] : Expected and actual position vectors;

[0360] : Maximum motor power;

[0361] : Motor power;

[0362] : Reward weight coefficient, .

[0363] Step S3030: Implementation of the policy optimization algorithm.

[0364] Algorithm selection: Proximal Policy Optimization (PPO) algorithm is adopted, and its update formula:

[0365] (31);

[0366] In the formula:

[0367] : Parametric distribution of the policy network, Gaussian distribution;

[0368] : Take the expectation of the state-action distribution at time step t;

[0369] : Historical parameters of the policy network. Preferably, it is synchronously updated every 5 iterations;

[0370] : Generalized Advantage Estimation;

[0371] : Policy clip threshold, ;

[0372] and : Action is the rotational speed command (4D vector) of four motors; State includes the attitude angle, angular velocity, position, and wind speed information of the UAV (12D vector).

[0373] Network architecture:

[0374] Actor network: 3-layer MLP (256-128-64), outputting the mean and variance of the Gaussian distribution;

[0375] Critic network: 2-layer LSTM (128 units) + fully connected layer, outputting the estimated state value.

[0376] Step S3040: Perform Monte Carlo testing.

[0377] Test scenario generation:

[0378] Randomly sample wind shear parameters: wind speed gradient , duration ;

[0379] Among them, the wind speed gradient is uniformly distributed, and the specific meaning of the wind speed gradient is that for every 100-meter increase in height, the wind speed changes by 5 - 15 m / s, and the duration is normally distributed, where 5 2 is the variance;

[0380] Perturbation injection location: Randomly set N disturb = 20 perturbation areas on the UAV flight path, and the spacing follows a Poisson distribution.

[0381] Among them, the test execution and termination conditions:

[0382] Single test duration: = 180 s;

[0383] Termination condition: a. The attitude angle deviation δ ≥ 30°, lasting for more than 2 s;

[0384] b. The height loss ∆h ≥ 50 m;

[0385] c. The number of times the servo is saturated ≥ 10;

[0386] The described termination condition means that the test stops when any one of a, b, or c is satisfied.

[0387] Among them, is the cumulative number of times the servo command ≥ 85%.

[0388] Step S3050: Adaptive sampling optimization.

[0389] Dynamically adjust the sampling weight based on the test results:

[0390] ;

[0391] In the formula:

[0392] : The working condition parameter vector that triggers termination;

[0393] : Includes parameters such as wind speed gradient, duration, and the number of perturbation areas, with a dimension of 3;

[0394] : Initial weight. Preferably, the initial weight = 1 / N s, N where s = 1000 is the total number of working condition samples;

[0395] : Gaussian kernel bandwidth, = 0.1.

[0396] Step S3060: Parameter adjustment trigger and verification.

[0397] Trigger logic:

[0398] Real-time monitor the composite attitude deviation:

[0399] (33);

[0400] When δ ≥ 30 ∘ and it lasts for t ≥ 0.5 s, start parameter adjustment:

[0401] (34);

[0402] In the formula:

[0403] : Saturation function;

[0404] : Adjustment coefficient, ,

[0405] : Pitch angle deviation;

[0406] : Roll angle deviation.

[0407] Online verification:

[0408] Immediately after adjustment, perform N = 5 times of fast Monte Carlo tests (T short = 10 s)

[0409] Verification indicators:

[0410] After adjustment: δmax ≤ 25°;

[0411] Overshoot: Mp ≤ 15%;

[0412] If the verification passes, update the flight control parameter library; otherwise, roll back and explore new strategies.

[0413] Step S400: Establish a risk probability model to evaluate the probability of dangerous working conditions occurring:

[0414] ;

[0415] In the formula:

[0416] : Shape parameter of extreme value distribution;

[0417] : Extreme event threshold, specifically, the composite attitude deviation threshold for triggering flight control failure;

[0418] : Location parameter;

[0419] : Scale parameter;

[0420] Specifically, the shape parameter 、location parameter 、scale parameter are determined by fitting the generalized extreme value distribution with historical flight data (including failure cases).

[0421] Furthermore, there is a preset threshold , and after step S400, it includes comparing the occurrence probability of dangerous working conditions with the preset threshold ;

[0422] If ≤ ， then execute step S500; If > and the current iteration number < , then return to step S300 and adjust the reinforcement learning reward function based on historical optimization data;

[0423] If ≥ , still does not meet ≤ , terminate the test and trigger a safety warning.

[0424] Among them, the risk threshold:

[0425] ;

[0426] Among them, = 1 / λ is the mean time between failures.

[0427] In the formula:

[0428] λ: System failure rate;

[0429] : Mission duration.

[0430] Further, > and the current iteration number < , then return to step S300 and adjust the reinforcement learning reward function based on historical optimization data in the following manner:

[0431] Parameter adjustment direction guidance:

[0432] Generate a targeted reinforcement learning reward function according to the type of dangerous working conditions (such as pitch out of control, roll oscillation):

[0433] (36);

[0434] In the formula:

[0435] ;

[0436] : The optimal action of the previous iteration. Specifically, it can be the mean of the action sequences with the top 5% of the reward values selected from the historical optimization database;

[0437] : Regularization coefficient. Specifically, Prevent overfitting.

[0438] Monte Carlo test scenario focus:

[0439] Based on the risk contribution analysis, increase the sampling weight for high-probability dangerous scenarios:

[0440] (37);

[0441] In the formula:

[0442] : .

[0443] In addition, the maximum number of iterations is carried out in the following manner: Set K max = 5 - 10 to prevent infinite loops.

[0444] Further, if ≥ still does not meet ≤ when, also lock the flight control parameter write protection and prohibit automatic update.

[0445] Step S500: Import the optimized flight control parameters into the physical unmanned aerial vehicle through the hardware-in-the-loop system to complete the safety verification.

[0446] Specifically, import the optimized flight control parameters into the physical unmanned aerial vehicle through the hardware-in-the-loop system through the following steps.

[0447] Step S510: Perform CRC check and ECC encryption on the optimized flight control parameters to generate a data packet with a digital signature.

[0448] Specifically, it also includes:

[0449] Preprocessing of flight control parameters: Pack the optimized PID parameters into a standardized data structure.

[0450] Calculate the checksum: Use the CRC-16-CCITT algorithm,

[0451] ;

[0452] In the formula:

[0453] : Byte sequence of the data packet;

[0454] : Polynomial variable;

[0455] : Number of bytes of the data packet.

[0456] Furthermore, the data packet format:

[0457] Structure: Frame header (2 bytes) | PID parameter (20 bytes) | Timestamp (4 bytes) | Checksum (2 bytes);

[0458] Byte order: Big-endian (high-order byte first).

[0459] Calculation steps:

[0460] Initialize the CRC register to 0xFFFF;

[0461] XOR byte by byte and shift right, reflecting the remainder;

[0462] Take the complement of the final CRC value (0xFFFF ^ CRC).

[0463] Dynamic encrypted transmission:

[0464] Use session key negotiation based on Elliptic Curve Cryptography (ECC):

[0465] ;

[0466] In the formula:

[0467] : Private key of the Hardware-in-the-Loop (HIL) system;

[0468] : Public key of the unmanned aerial vehicle;

[0469] : Session key.

[0470] Parameter encryption: Use the AES-GCM mode with an additional authentication tag (MAC).

[0471] Step S520: Inject into the flight control system through a dual-redundant CAN FD bus with a delay of ≤8 ms, and synchronize using the PTP clock protocol.

[0472] Among them, the physical layer configuration:

[0473] Interface type: Dual-redundant CAN FD bus, baud rate 5 Mbps;

[0474] Synchronization mechanism: Based on the IEEE 1588 PTP protocol, the clock synchronization error is ≤1 μs.

[0475] Communication protocol stack: The application layer protocol uses the drone-specific standard DDS-XRCE:

[0476] Real-time guarantee:

[0477] Set the QoS policy: The transmission priority is Critical, and the end-to-end delay is ≤8 ms;

[0478] The hardware interrupt response time is ≤5 μs.

[0479] Step S530: Segmentally activate the parameters and verify the stability margin in real time. If the servo saturation limit is exceeded, trigger a rollback.

[0480] Among them, the segmented injection strategy: Divide the parameters into a base group (Base) and an optimized group (Optimized), and use the gray release mechanism:

[0481] (40);

[0482] In the formula:

[0483] : Warm-up time, preferably = 2 s, used to stabilize the initial parameters;

[0484] : Transition time, ;

[0485] : Mixing coefficient, linearly decaying from 1 to 0.

[0486] Activation condition monitoring:

[0487] Online verification of the flight control system stability margin: Gain margin Gm ≥ 6 dB, phase margin ϕm ≥ 45°.

[0488] Servo travel saturation monitoring: If any servo ≥85%, and it lasts for more than 0.2 s continuously, a rollback is triggered.

[0489] Step S540: Perform step response, frequency sweep, and Monte Carlo pressure tests, and calculate the comprehensive success rate.

[0490] Among them, the step response test:

[0491] Apply a pitch-axis step command θcmd = 10°, and verify that the rise time tr ≤ 1.5 s, the overshoot Mp ≤ 15%, and the steady-state error ess ≤ 0.5°.

[0492] Among them, the frequency sweep test:

[0493] Input a sweep signal , ;

[0494] It is required that the amplitude attenuation ≤ -3 dB, the frequency bandwidth ≥ 5 Hz, and when the frequency is 10 Hz, the phase lag does not exceed 90 degrees.

[0495] Among them, the Monte Carlo pressure test:

[0496] Inject 20 groups of randomly combined wind shear scenarios, and the success rate threshold:

[0497] ;

[0498] In the formula:

[0499] : The number of times the attitude deviation meets the standard in the test;

[0500] : The total number of tests.

[0501] Step S550: When a Level 3 anomaly is detected, roll back to the safe parameter version within 50 ms and start the forced landing procedure.

[0502] Multi-level anomaly detection:

[0503] Level 1 (mild): Record the log, and the warning code is 0x1A;

[0504] Level 2 (moderate): Freeze the current parameters and start the auxiliary controller;

[0505] Level 3 (severe: such as the attitude angle > 45°): Perform a full parameter rollback to the safe version and trigger an emergency forced landing.

[0506] Rollback protocol:

[0507] Rollback to the Golden version within ≤ 50 ms;

[0508] Verify the authenticity of historical parameters using a hash chain:

[0509] ;

[0510] Where:

[0511] : The hash value of the nth parameter version;

[0512] : The (n - 1)th parameter set;

[0513] : The data concatenation operation, representing concatenation by big - endian bytes;

[0514] In addition, the initial hash value is the root hash at secure startup (hard - coded in the flight control firmware).

[0515] The following further illustrates the present invention through specific embodiments: Embodiment 1

[0516] Test object: A six - axis drone, with a weight of 4.2 kg, a wheelbase of 650 mm, and a maximum hovering wind speed of 15 m / s.

[0517] Wind field configuration: Building wake model: A rectangular building with a height of 80 m, using the formula:

[0518] ,

[0519] Superimposed downwind shear: Superimposed downwind shear: Gradient ∇U = 10 m / s / 50 m, lasting for 20 s.

[0520] Optimization process: Adjust parameters using reinforcement learning: .

[0521] The Monte Carlo test coverage rate is increased to 91%.

[0522] Test results:

[0523] The maximum attitude deviation is reduced from 34° to 19°;

[0524] The RMSE of the trajectory tracking error at the building corner is reduced from 3.5 m to 1.2 m. Embodiment 2

[0525] A fixed - wing drone, with a wingspan of 2.8 m, a stall speed of 13 m / s, and a cruise speed of 25 m / s;

[0526] Disturbance setting: Microburst model: Microburst model: Sinking velocity w = -12 m / s, influence radius 200 m;

[0527] Adopt the dynamic stall model: .

[0528] Flight control optimization: The PPO algorithm adjusts the parameters of the angle of attack limiter: α max is increased from 18° to 24°.

[0529] Introduce the energy reward term R energy , reducing the overshoot risk of the power system.

[0530] Test results:

[0531] The recovery time after stall is shortened from 8.2 s to 3.1 s.

[0532] The airspeed fluctuation during the gust penetration phase is reduced by 42%.

[0533] The above are only the embodiments of the present invention, and do not limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made by using the content of the specification and drawings of the present invention, or directly or indirectly applied in other related technical fields, shall be included in the patent protection scope of the present invention by the same token.

Claims

1. A method for evaluating the safety of drones based on digital twins, characterized in that, It includes the following steps: Step S100: Construct a high-precision digital twin of the drone, including an aerodynamic model, a flight control system mirror image, and a sensor noise model, where the aerodynamic model meets the simulation accuracy with Reynolds number Re≥1×10^5; Step S200: Generate a dynamic airflow disturbance field in the virtual environment, and the disturbance field includes: Gust model based on Davenport spectrum, wind speed fluctuation range , is the reference wind speed; Simulate the wind tunnel effect of the building complex using the k-ω SST turbulence model, with a spatial resolution ≤0.1m; Step S300: Optimize the flight control parameters through the reinforcement learning algorithm, perform Monte Carlo tests in the digital twin environment, and trigger parameter adjustment when the attitude angle deviation δ satisfies: wherein, is the pitch angle deviation, is the roll angle deviation; Step S400: Establish a risk probability model to evaluate the probability of dangerous working conditions: ； Among them, : extreme value distribution shape parameter; s: extreme event threshold; u: location parameter; σ: scale parameter; Meanwhile, a preset threshold is set , and after step S400, it includes comparing the probability of occurrence of a dangerous working condition with the preset threshold ; if ≤ ， then execute step S500; if > and the current iteration number < , then return to step S300 and adjust the reinforcement learning reward function based on historical optimization data; if ≥ still does not satisfy ≤ , terminate the test and trigger a safety warning; Step S500: Import the optimized flight control parameters into the physical drone through the hardware-in-the-loop system to complete the safety verification.

2. The drone safety evaluation method according to claim 1, wherein, In step S200, generating a dynamic airflow disturbance field in the virtual environment includes the following steps: Step S2010: Construct an intelligent disturbance generation framework; Step S2020: Perform basic wind field modeling, where the basic wind field is a gust model based on the Davenport spectrum; Step S2030: Model turbulence and building flow around; Step S2040: Data-driven intelligent disturbance enhancement, including the following steps: Construct a physical constraint generative adversarial network: Among them, Generator G(T, Mt; θG): Adopting a U-Net structure, it generates a high-resolution perturbation field ; Discriminator D(U;θD): Based on PatchGAN, verify the Navier-Stokes equation residual: ； Total loss function: ； Among them, in the formula: T: topographic elevation map; M t : real-time meteorological data; θG : trainable weight parameters of the generator U-Net; U: composite three-dimensional wind speed field containing the generated perturbation field and the basic wind field; θD : discriminator network parameters; : three-dimensional differential operator; : kinematic viscosity; : pressure gradient; : material derivative; : turbulence energy spectrum; : Kolmogorov theoretical energy spectrum; : expected value; : high-resolution perturbation field output by the generator; : adversarial loss weight; : physical constraint weight; : energy spectrum matching weight; Step S2050: Real-time disturbance field fusion and injection.

3. The drone safety evaluation method according to claim 2, wherein In step S200, optimize the flight control parameters through a reinforcement learning algorithm and perform Monte Carlo tests in the digital twin environment, which are carried out through the following steps: Step S3010: Define the reinforcement learning state-action space; Step S3020: Design a multi-objective reward function; Step S3030: Implement the policy optimization algorithm; Step S3040: Perform Monte Carlo tests, including the following steps: Generate test scenarios; Randomly sample wind shear parameters: wind speed gradient , duration ; Among them, the wind speed gradient is uniformly distributed, and the duration is normally distributed; Perturbation injection location: Randomly set on the UAV flight path = 20 perturbation areas, and the spacing follows a Poisson distribution; Among them, the test execution and termination conditions: Single test duration: = 180 s; Termination condition: a. The attitude angle deviation δ≥30° lasts for more than 2s; b. The height loss ∆h≥50m; c. Number of times the steering gear is saturated ≥10; The above termination condition means that the test stops as long as any one of a, b, or c is satisfied; Step S3050: Adaptive sampling optimization; Step S3060: Parameter adjustment triggering and verification.

4. The drone safety evaluation method according to claim 3, wherein In step S3060, the parameter adjustment triggering and verification are carried out according to the following rules: Trigger logic: Real-time monitor the composite attitude deviation: When and it lasts for t≥0.5s, start parameter adjustment: ， ； In the formula: : saturation function; : adjustment coefficient; Online verification: Immediately after adjustment, perform N = 5 quick Monte Carlo tests, where, ; Verification index: After adjustment: ≤ 25°; Overshoot: ≤ 15%; If the verification is passed, update the flight control parameter library; otherwise, roll back and explore new strategies.

5. The drone safety evaluation method according to claim 3, wherein, In step S3020, the multi-objective reward function is designed as: Composite reward function: , where Stability reward: ； Among them, ; Trajectory tracking reward: ； Energy efficiency reward: ； In the formula: : Time window length; : Angle deviation; : Number of trajectory points; , : Expected and actual position vectors; : Motor power; : Maximum motor power; : Reward weight coefficients, W1 = 0.6, W2 = 0.3, W3 = 0.

1.

6. The drone safety evaluation method according to claim 3, wherein, Step S3050: Adaptive sampling optimization is carried out in the following way: Dynamically adjust the sampling weight based on the test results: ； Among them, : The working condition parameter vector that triggers termination; : It includes parameters such as wind speed gradient, duration, and the number of disturbance areas, with a dimension of 3; : Initial weight; : Gaussian kernel bandwidth, = 0.

1.

7. The drone safety evaluation method according to claim 3, wherein The policy optimization algorithm in step S3030 is the proximal policy optimization algorithm.

8. The drone safety evaluation method according to claim 3, wherein In step S500, import the optimized flight control parameters into the physical drone through the hardware-in-the-loop system to complete the safety verification, which specifically includes the following steps: Step S510: Perform CRC check and ECC encryption on the optimized flight control parameters to generate a data packet carrying a digital signature; Step S520: Inject it into the flight control system through a dual-redundancy CAN FD bus with a delay ≤8ms, and synchronously adopt the PTP clock protocol; Step S530: Activate parameters in segments and verify the stability margin in real time. If the saturation of the servo exceeds the limit, trigger a rollback; Step S540: Perform step response, frequency sweep, and Monte Carlo stress tests, and calculate the comprehensive success rate; Step S550: When a Level 3 anomaly is detected, roll back to the safe parameter version within 50 ms and start the emergency landing procedure; the Level 3 anomaly is a situation where the attitude angle > 45°.

9. The drone safety evaluation method according to claim 6, wherein In step S100, it also includes injecting noise into the digital twin, and injecting noise into the digital twin includes the following steps: Step S1110: Identify noise parameters based on Allan variance analysis; Step S1120: Generate noise using a second-order autoregressive model; Step S1130: Superimpose the noise signal on the flight control system mirror image.

10. The drone safety evaluation method according to claim 7, characterized in that, The noise is calibrated in the following manner: Step S1210: Collect raw sensor data; Step S1220: Quantitatively analyze the noise characteristics; Step S1230: Calibrate the noise model parameters; Step S1240: Verify the injection of digital twin noise. Moreover, inject the calibrated noise model into the digital twin, and the KL divergence of the probability distributions between the digital twin and the output of the physical sensor ≤ 0.1, and the time-domain correlation error ≤ 15%.

Images