Data processing method, apparatus, and medium for aircraft atmospheric turbulence warning

Abstract

The present application relates to the technical field of computer data processing and artificial intelligence, in particular to a data processing method, device and medium for aircraft atmospheric turbulence early warning, the method comprising: offline training a first scale meteorological risk assessment model and a second scale turbulence disturbance model; in the process of aircraft flight, real-time acquisition of time series flight state data and corresponding first scale meteorological data, and extraction of first scale meteorological features and second scale turbulence features; based on the ensemble Kalman filtering algorithm, the output of the first scale meteorological risk assessment model is taken as the background trend item, and based on the key indicator factor determined by the second scale turbulence disturbance model, a nonlinear correction item is generated from real-time data, and the two are fused to generate a prediction result of atmospheric turbulence intensity; generating a turbulence early warning signal according to the prediction result. Through multi-scale information fusion, the present application realizes accurate and early warning of atmospheric turbulence, effectively improving flight safety.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of computer data processing and artificial intelligence, and particularly relates to a data processing method, device and medium for aircraft atmospheric turbulence early warning. BACKGROUND

[0002] In the field of aviation flight safety, atmospheric turbulence (especially the clear air turbulence which is difficult to accurately detect) is the main factor causing aircraft jolt and threatening flight safety and comfort. From the perspective of aerodynamics, jolt is essentially the sharp and irregular change of the flow field pressure distribution on the aerodynamic interface of the aircraft when it passes through a turbulence field composed of vortices of different scales, which leads to the destruction of the dynamic balance of the lift, drag and control moment acting on the aircraft, and finally results in the violent oscillation of key flight state parameters such as flight altitude, airspeed and attitude angle. The strength of jolt is determined by the energy of atmospheric turbulence (the recommended quantitative indicator is eddy dissipation rate (EDR)) and the flight speed, wing load and other characteristics of the aircraft.

[0003] At present, the main flight jolt early warning technology mainly relies on the following two technical routes, but they all have significant limitations:

[0004] 1. Early warning method based on large-scale numerical weather prediction model

[0005] This method relies on the numerical prediction products published by global or regional meteorological centers. Its fundamental defect is:

[0006] Insufficient spatial and temporal resolution: the calculation grid of the model is usually sparse (usually tens of kilometers), and the time update frequency is low (once every several hours). The characteristics of such coarse-grained data make it difficult to effectively capture and characterize local, rapidly growing and decaying small and medium scale turbulence systems (such as clear air turbulence, mountain waves, etc.) on the flight route.

[0007] Low matching degree of data and object: the model output is a macroscopic estimate of the atmosphere, not a precise perception of a specific flight, a specific flight path, and a specific moment. Therefore, its early warning results have inherent uncertainty and lag, and false alarms and missed reports occur frequently.

[0008] 2. Onboard early warning method based on special detection equipment

[0009] This type of method is represented by airborne laser radar (LIDAR) technology. Its working principle is advanced and can detect aerosols in front of the flight path to indirectly retrieve turbulence. However, its application is greatly limited:

[0010] Economic cost and deployment feasibility: the laser radar sensor itself is expensive, and its installation, calibration and maintenance require high cost and involve complex aircraft modification, which is difficult to achieve large-scale and popular deployment in the large number of existing and future commercial aircraft fleets.

[0011] Technical limitations: its detection capability in clean air (such as clear sky conditions) is better, but its performance under certain weather conditions will be degraded.

[0012] In summary, the existing technology system is trapped in the dilemma of "macro data is not fine, and fine detection is not economical" in dealing with the problem of aircraft turbulence warning. The deep technical contradiction lies in the lack of a core processing method that can extract and utilize the implicit value in the existing popular and low-cost data sources through efficient computing models and information fusion technology.

[0013] The "uncertainty" of the current warning result is essentially caused by the heterogeneity of the data source, the inefficiency of data processing, and the insufficient adaptability of the warning model. Therefore, the industry urgently needs a revolutionary technical solution that does not rely on a single high-cost special device, but through innovative data processing algorithms and fusion models, realizes the deep mining and intelligent analysis of multi-source heterogeneous information, so as to achieve the "low cost, high precision, high timeliness" warning goal. SUMMARY

[0014] To solve the above technical problems, the technical scheme adopted by the present application is:

[0015] According to the first aspect of the present application, a data processing method for aircraft atmospheric turbulence warning is provided, which comprises the following steps:

[0016] S100, based on the target historical meteorological data and the target time series flight state data, offline constructing a first scale meteorological risk assessment model and a second scale turbulence disturbance model; wherein the spatial coverage range of the first scale is larger than that of the second scale, and the construction process of the second scale turbulence disturbance model includes determining at least one key indicator factor.

[0017] S200, in the process of aircraft flight, real-time acquiring time series flight state data and corresponding first scale meteorological data, and extracting first scale meteorological features and second scale turbulence features based on the acquired data.

[0018] S300, generating a background trend term based on the first scale meteorological feature through the first scale meteorological risk assessment model, and generating a nonlinear correction term from the second scale turbulence feature based on at least one key indicator factor determined in the construction process of the second scale turbulence disturbance model; using an ensemble Kalman filtering algorithm to fuse the background trend term and the nonlinear correction term to generate a prediction result of atmospheric turbulence intensity.

[0019] S400, generating a turbulence early warning signal according to the prediction result.

[0020] According to the second aspect of the application, an electronic device is provided, comprising a processor and a memory; the processor is used to execute the steps of the method of the first aspect of the application by calling the program or instruction stored in the memory.

[0021] According to the third aspect of the application, a computer readable storage medium is provided, which stores a program or instruction, and the program or instruction is used to execute the steps of the method of the first aspect of the application.

[0022] The application has at least the following beneficial effects:

[0023] 1. Precise prediction of vertical wind speed change is realized

[0024] The application realizes high-precision prediction of the aircraft encountering vertical wind speed change in a future period of time by constructing a double-scale model offline and dynamically fusing large-scale background trend and small-scale bumping precursor correction signal in the flight process combined with real-time observation data. The prediction accuracy is further improved by introducing a bias correction model.

[0025] 2. Quantitative bumping intensity early warning is provided

[0026] Based on the predicted vertical wind speed sequence, the running standard deviation is calculated by using the sliding window, and the bumping intensity of each point on the future flight path is quantitatively calculated by using the vortex dissipation rate estimation formula. By comparing with the preset bumping level threshold, the specific time and intensity can be predicted before the bumping actually occurs, which provides effective early warning for flight and significantly reduces the unexpected bumping risk.

[0027] 3. It has wide scene adaptability

[0028] By fusing multi-source and multi-scale information, combined with the ensemble Kalman filtering (EnKF) framework which can handle nonlinear processes and quantify uncertainty, the application can provide accurate early warning and effective early warning time for different flight scenes (including persistent bumping, sudden strong bumping), different bumping intensity and atmospheric background conditions, showing excellent adaptability and robustness.

[0029] It is to be understood that the details set forth herein do not limit the scope of the embodiments of the application to the specific embodiments described. Rather, the scope of the embodiments of the application is to be defined by the appended claims. BRIEF DESCRIPTION OF DRAWINGS

[0030] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following will briefly introduce the drawings needed in the embodiments description. Obviously, the drawings in the following description are only some embodiments of the present application, and for those skilled in the art, other drawings can also be obtained without creative labor.

[0031] Figure 1 A flow chart of a data processing method for aircraft atmospheric turbulence early warning provided by the embodiments of the present application. DETAILED DESCRIPTION

[0032] The technical solutions in the embodiments of the present application will be described clearly and completely in the following with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only some embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0033] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this application belongs. The terminology used in the description of the application herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the application. As used in this description, the singular forms "a", "an" and "the" include plural references unless the context clearly dictates otherwise. The term "and / or" includes any and all combinations of one or more of the associated listed items.

[0034] It should be noted that some of the example embodiments are described as processes that are depicted as a flow diagram or a flowchart. Although a flow diagram or a flowchart can describe processes as a sequential process, many of the steps can be performed in parallel, concurrently or simultaneously. In addition, the order of the steps can be re-arranged. A process can be terminated when its operations are completed, but could also occur under some other condition, such as in response to a user command. The processes might also correspond in whole or in part to the procedures created by the application, a function, a procedure, a subroutine, a subprogram, etc. When the processes of the application are implemented partially or fully in software, the software can comprise one or more instructions that are resident in memory and / or storage and that, when executed by one or more processors, cause the procedures described herein to be executed.

[0035] The embodiments of the present application provide a data processing method for aircraft atmospheric turbulence early warning. The method aims to realize accurate and early warning of atmospheric turbulence which endangers flight safety through advanced data processing technology, so as to provide sufficient decision and response time for pilots and effectively avoid or reduce aircraft jolt.

[0036] In the present application, atmospheric turbulence is a kind of fluid motion state caused by different scale vortices and irregular motion in airflow, and the physical properties (such as wind speed, temperature) change randomly and violently. In the context of aviation application, it specifically refers to the disturbed airflow that causes irregular shaking, shaking, height and attitude change of the aircraft, and is the direct physical cause of flight turbulence.

[0037] As Figure 1 shown, the data processing method for aircraft atmospheric turbulence early warning provided by the embodiment of the present application comprises the following steps:

[0038] S100, based on the target historical meteorological data and the target time series flight state data, offline constructing a first scale meteorological risk assessment model and a second scale turbulence disturbance model; wherein the spatial coverage range of the first scale is larger than that of the second scale.

[0039] In the embodiment of the present application, the historical meteorological data is used to reflect the state of the historical large-scale atmospheric background field, and provides an analysis basis for constructing the first scale meteorological risk assessment model. In a specific implementation, the ERA5 reanalysis data set published by the European Centre for Medium-Range Weather Forecasts is used, which contains gridded data of multiple meteorological elements such as wind field, temperature field and geopotential height field in the global range.

[0040] The time series flight state data is a flight parameter sequence from an onboard data bus system with a sampling rate not less than 1Hz. Specifically, the embodiment uses the quick access recorder data recorded by the commercial aircraft in actual operation as the data source, which records the spatial position parameters, attitude parameters and maneuvering parameters in the flight process with high time series resolution, including latitude, longitude, barometric height, vertical acceleration and airspeed information. These parameters provide both real labels for identifying turbulence events and input data for calculating vortex dissipation rate.

[0041] To further improve the representativeness and pertinence of the model training data, the embodiment performs spatial clustering analysis on the turbulence events recorded in the quick access recorder (QAR) data to filter out the most representative training samples, thereby obtaining the target historical meteorological data and the target time series flight state data for constructing the model. Specifically, the following steps are included: S10, extracting the occurrence position coordinates of all turbulence events recorded in the historical flight state data to form a spatial point set to be clustered.

[0042] In the embodiment of the present application, the occurrence position coordinates can include latitude and longitude coordinates.

[0043] S11, according to the geographical distribution characteristics of the spatial point set, presetting the initial neighborhood radius and minimum sample number parameters of the DBSCAN (density-based noise application spatial clustering) algorithm.

[0044] Specifically, by analyzing the k-distance distribution diagram of all coordinate points in the spatial point set (usually k takes the initial value of the minimum sample number), the distance value corresponding to the inflection point or elbow position in the distance distribution curve is set as the initial neighborhood radius; the minimum sample number is initialized and selected in the interval [3, 10] according to the total number of spatial point sets and the clustering sensitivity requirement combined with the empirical value.

[0045] S12, based on the initial neighborhood radius and the minimum sample number parameters set, the DBSCAN analysis is performed on the coordinate of the occurrence position of the turbulence event, and a plurality of clustering clusters are obtained. The DBSCAN algorithm divides the points with a mutual distance within the initial neighborhood radius and a sample number reaching the minimum sample number requirement into the same clustering cluster by measuring the geographical distance between the points, thereby realizing the spatial clustering based on density and obtaining a plurality of clustering clusters.

[0046] S13, the turbulence event density of each clustering cluster is calculated, and the geographical range represented by the clustering cluster with the highest turbulence event density is identified as the target training area; the target training area represents the core area of high incidence of atmospheric turbulence in space.

[0047] In the embodiment of the application, the turbulence event density is calculated by the number of turbulence events per unit area, and the specific formula is: turbulence event density = total number of turbulence events in the clustering cluster / geographical area covered by the clustering cluster, wherein the geographical area is obtained by calculating the convex hull area of the clustering cluster. The application adopts a mature method in computational geometry to calculate the convex hull area, specifically: for each clustering cluster composed of longitude and latitude coordinate points, the convex hull of the clustering cluster is calculated to define the geographical boundary of the clustering cluster, and the convex hull is a minimum convex polygon containing all points in the clustering cluster; subsequently, based on the WGS84 earth ellipsoid model, the actual projection area of the convex hull on the earth's surface is calculated by numerical integration method.

[0048] S14, from the historical meteorological database and the time series flight state data set, the data completely corresponding to the space-time range of the target training area is extracted, and the data is used as the target historical meteorological data and the target time series flight state data for subsequent model training and construction.

[0049] In the embodiment of the present application, the first scale refers to a scale of a meteorological field that can be effectively distinguished and output by an existing numerical weather prediction model, and the horizontal spatial range thereof is generally more than several hundred to two thousand kilometers, and the time scale is several hours to several days. The first scale is directly output by the numerical weather prediction model, and describes a macro atmospheric circulation background such as a pressure system, a front, a jet stream, and the like, which is an environmental field of turbulence generation and development. In contrast, the second scale or small scale refers to an aircraft encountering scale, and the spatial range thereof is generally less than a model grid scale (within several tens of kilometers), and the time scale is a minute level, which is expressed as a local vortex that rapidly appears and disappears in the macro background field. One of the innovations of the present application is to dynamically correct the first scale background field by the second scale information. Taking the ECMWF ERA5 reanalysis data used in the embodiment as an example, the grid resolution of the ECMWF ERA5 reanalysis data is about 0.25°*0.25° (about 31 kilometers), which defines the upper limit of the first scale in the present application; and the disturbance of the second scale, that is, occurs in the grid and cannot be directly analyzed by the model.

[0050] Further, in S100, the process of constructing the first scale meteorological risk assessment model includes:

[0051] S101, spatiotemporally match the target historical meteorological data and the target time-series flight state data, calculate a plurality of turbulence prediction indexes based on the matched data, and select turbulence prediction indexes whose prediction performance representation values reach a preset threshold to constitute a basic feature set.

[0052] In the embodiment of the present application, in order to realize accurate association of the aircraft measured data and the macro meteorological background, the following process is used to spatiotemporally match each turbulence event:

[0053] Step 1, determine a spatiotemporal reference point: take one turbulence event recorded in the target time-series flight state data as a processing unit, and take the occurrence time (t0) and the spatial position (longitude λ0, latitude φ0, and height h0) of the turbulence event as the reference of this matching.

[0054] Step 2, time dimension matching:

[0055] In order to ensure that the used meteorological data can effectively represent the macro background field at the time of turbulence event occurrence, first, according to the time resolution (Δt, such as 6 hours) of the historical meteorological data (such as ECMWF ERA5), the maximum allowed time tolerance window is determined as ±Δt / 2 (that is, [t0-Δt / 2, t0+Δt / 2]).

[0056] Select the meteorological data three-dimensional grid field closest to the turbulence event occurrence time t0 as the basis data source matched with the turbulence event. To ensure data timeliness, only the meteorological data field with an absolute value of the difference between the meteorological data timestamp and the turbulence event occurrence time less than or equal to Δt / 2 can be used as the basis data source for the event matching. For example, for 6-hourly updated data, the latest one of the data three-dimensional grid fields within 3 hours of t0 is selected.

[0057] Step 3, spatial dimension matching: In the meteorological data three-dimensional grid field selected in step 2, spatial interpolation is performed for the reference position (λ0, φ0, h0) to obtain the meteorological element value at the precise point:

[0058] Horizontal interpolation: According to the coordinates of (λ0, φ0), locate the grid cell in which it is located in the Earth coordinate system. Using the bilinear interpolation method, the horizontal interpolation result at the position of (λ0, φ0) is calculated using the meteorological element values of the four corner points of the grid of the cell.

[0059] Vertical interpolation: According to h0 (convertible to pressure value P0), locate the adjacent two standard pressure layers in the vertical coordinate of the meteorological field. Using the linear interpolation or logarithmic pressure interpolation method, the vertical interpolation result at the height of h0 is calculated using the meteorological element values of the two standard layers.

[0060] Step 4, generate a matching vector: Through the above steps, a unique feature vector containing multiple basic meteorological elements (such as U, V, T, P obtained after interpolation) is generated for each turbulence event.

[0061] Steps 1 to 4 are executed in a loop until all turbulence events are processed, thereby completing the spatio-temporal matching of the entire data set.

[0062] After the target historical meteorological data and the target time series flight data are accurately matched in space and time, multiple turbulence prediction indices related to clear air turbulence are calculated based on the matched basic meteorological elements (such as wind field U / V, temperature T, pressure P, etc.). The turbulence prediction indices include but are not limited to:

[0063] Richardson number Ri: It is used to diagnose the turbulence tendency by calculating the atmospheric stability and wind shear, and its calculation formula is: where g is the acceleration of gravity, θ is the potential temperature, U and V are the zonal (eastward) and meridional (northward) wind speeds, respectively, and z is the height. The smaller the Ri number, the stronger the dynamic instability, which is more conducive to the development of turbulence.

[0064] Total deformation index D: It is used to identify strong flow regions by calculating the deformation rate of the horizontal wind field, and its calculation formula is: , U1 and V1 are the components of the wind vector in a local orthogonal coordinate system tangent to the Earth's surface, pointing to the positive east (x-axis) and positive north (y-axis) directions, respectively. In practical calculations, the wind field components are provided directly in the orthogonal coordinate system by numerical weather prediction grid data or converted by projection. x and y are two orthogonal axes of the local orthogonal coordinate system. The larger the D value, the stronger the horizontal wind shear, and the greater the probability of turbulence.

[0065] Euler number: the percentage of grid points in the calculation area with Richardson number Ri less than the critical value 0.25, used to identify local unstable areas where turbulence may occur.

[0066] Kolmogorov scale spectrum index: based on the theory of turbulent energy dissipation, the turbulent kinetic energy dissipation rate is estimated.

[0067] In the embodiments of the present application, the area under the receiver operating characteristic curve (AUC) is used as a performance characterization value to evaluate the prediction ability of each turbulence prediction index for real turbulence events. The closer the AUC value is to 1, the stronger the prediction ability of the index.

[0068] The preset threshold is determined by any of the following ways:

[0069] a) statistical significance criterion: calculate the statistical significance of the AUC value of each index relative to the no-prediction skill baseline value (AUC=0.5). Set the preset threshold to the lowest level that is significantly different from the baseline value, for example, AUC>0.55, or AUC>0.60 and the significance level (p-value) <0.05.

[0070] b) ranking selection criterion: according to the preset number of features required, all indices are sorted in descending order of AUC value, and the top N indices (e.g. the top 10) are selected; at this time, the preset threshold is dynamically determined as the AUC value corresponding to the last index in the sorted list.

[0071] S102, pre-processing the basic feature set to obtain a pre-processed feature time series; the pre-processing includes outlier processing, data standardization and data transformation.

[0072] In the embodiments of the present application, the outlier processing adopts the identification and truncation method based on the interquartile range (IQR). Specifically, for each type of turbulence prediction index, calculate the first quartile (Q1) and the third quartile (Q3), and define the normal value range as [Q1-1.5×IQR, Q3+1.5×IQR]; the data points outside this range are truncated to the upper and lower boundaries of the range.

[0073] Data standardization is to perform Z-score standardization on the feature data after the completion of outlier processing, so as to eliminate the dimensional difference between different prediction indexes. The standardization formula is: (a-μ) / σ, wherein a is the original feature value, μ is the mean value of the feature in all samples, and σ is the standard deviation.

[0074] Data transformation is to apply power transformation to the standardized data, aiming to make its distribution closer to the normal distribution to meet the assumption of the subsequent model on the data distribution. Specifically, Box-Cox transformation or Yeo-Johnson transformation is adopted, and the optimal transformation parameter lambda (λ) is automatically determined by maximum likelihood estimation.

[0075] S103, signal decomposition and screening are performed on the preprocessed feature time series to extract key dynamic modes related to the generation and decay mechanism of turbulence.

[0076] The generation and decay mechanism of turbulence is a functional and physical general description, which specifically refers to the sum of physical processes that cause the generation, maintenance and decay of turbulence in the atmosphere. These processes are usually associated with different scale weather phenomena, such as:

[0077] Generation mechanism: such as Kelvin-Helmholtz instability (triggered by wind shear), terrain lifting, convection activity, etc.

[0078] Maintenance and dissipation mechanism: such as energy cascade of turbulence, interaction with background wind field, and viscous dissipation, etc.

[0079] The key dynamic mode described in the application aims to extract those signal components with predictive significance from the historical sequence of prediction indexes, which are synchronized or ahead of the time evolution of the above-mentioned physical processes.

[0080] Further, S103 specifically includes:

[0081] S1031, for the time series of each turbulence prediction index in the basic feature set, a wavelet decomposition technique is used for multi-scale signal decomposition to extract candidate dynamic modes.

[0082] Specifically, a suitable wavelet basis function (such as Daubechies-4 wavelet) and decomposition level (for example, 5 layers) are selected to decompose each original signal into multiple intrinsic mode function (IMF) components on different frequency subbands. These IMF components constitute an initial candidate dynamic mode set, which is a signal component representing different time scale fluctuation characteristics separated from the original signal.

[0083] S1032, a multi-index quantitative evaluation system is constructed to screen key dynamic modes from the candidate dynamic modes.

[0084] The multi-index quantitative evaluation system is performed on the model training set, that is, for each time series of the candidate dynamic modality, the following four indicators between the time series and the true turbulence event label are calculated:

[0085] Pearson correlation coefficient: used to evaluate linear correlation.

[0086] Spearman rank correlation coefficient: used to evaluate monotonic correlation.

[0087] Mutual information: used to evaluate linear and nonlinear general dependence.

[0088] AUC: used to evaluate the discriminant ability of turbulence events.

[0089] Subsequently, the weighted average or arithmetic mean of the four indicator values of each candidate modality is calculated as the comprehensive score of the candidate modality. Finally, all candidate modalities are ranked from high to low according to the comprehensive score, and the top K candidate modalities are selected as the final key dynamic modalities for subsequent model training. Wherein K is a positive integer preset according to the model complexity and performance requirements, for example, K=2.

[0090] S104, based on the selected key dynamic modalities, an integrated classification model is trained to form the first scale meteorological risk assessment model.

[0091] The integrated classification model specifically refers to a machine learning model composed of multiple base classifiers (Gaussian process classifiers in this embodiment) according to a preset weight combination. The final output of the integrated classification model is the weighted combination of the outputs of each base classifier, which is used to improve the accuracy and stability of the prediction.

[0092] Further, S104 specifically includes:

[0093] S1041, base classifier training: for each key dynamic modality selected, the time series data of the key dynamic modality is used as an independent input feature, and the corresponding true turbulence event label is used as an output, and a base classifier, i.e., a Gaussian process classifier (GPC), is trained. The Gaussian process classifier uses a radial basis function as a covariance function and optimizes the hyperparameters by maximizing the marginal likelihood function. Each trained GPC can independently output a preliminary prediction value of the turbulence occurrence probability.

[0094] S1042, integrated weight calculation: according to the comprehensive score of each key dynamic mode obtained in step S103, the weight of the key dynamic mode in the integrated classification model is calculated by using the entropy weight method. The specific calculation process is as follows: first, the comprehensive score of each key dynamic mode is normalized; then, according to the principle of entropy weight method, the information entropy of each key dynamic mode score is calculated, and the difference coefficient of the key dynamic mode is further calculated according to the information entropy; finally, the objective weight of each mode is determined according to the difference coefficient. Under this mechanism, the key dynamic mode with stronger correlation with the turbulence event and higher comprehensive score is given a larger weight.

[0095] S1043, model integration: all the trained Gaussian process classifiers in step S1041, the weight calculated in step S1042, and a weighted average fusion algorithm are associated to define the first scale meteorological risk assessment model.

[0096] When the first scale meteorological risk assessment model receives new input data, it works in the following way: all GPC-based classifiers are called in parallel, and the probability values output by them are weighted and averaged according to the weights, and finally an integrated turbulence occurrence probability value, i.e. integrated prediction probability, is output, which represents the turbulence risk level under the macro atmospheric background field. This weighted average method is a mainstream and efficient integration strategy. It reflects the relative importance of the information carried by different key dynamic modes (objectively determined by entropy weight method according to the comprehensive score) through weights, so as to fuse the judgments of multiple weak classifiers into a stronger and more accurate classifier.

[0097] Further, in S100, the process of constructing the second scale turbulence disturbance model includes:

[0098] S110, identifying turbulence events from the target time series flight state data, and dividing the process data of the turbulence events into three stages of pre-event, event and post-event based on the peak time of vertical acceleration.

[0099] Based on the vertical acceleration sequence in the target time series flight state data, the effective turbulence event is identified. When the absolute value of vertical acceleration of a data point continuously exceeds a preset threshold (such as ±0.15g), it is marked as the beginning of a potential event. From the starting point, the vertical acceleration is continuously recorded until it falls below the threshold and stabilizes for a period of time (such as 30 seconds), which is marked as the end of the event. A complete turbulence event must meet the minimum duration requirement (such as 3 seconds) to exclude transient interference. The start time, end time and precise time t c .

[0100] Pre-event stage process data is all flight state data (such as vertical acceleration, airspeed, attitude angle, etc.) in a time window of [t c – Δt pre , t c ]. Wherein Δt pre is a preset forward time length (for example, 5 minutes), and the data set of this stage independently represents the early process of the evolution of the atmospheric background state to the disturbance state before the turbulence occurs.

[0101] Event stage process data is all flight state data in a time window of [t c – Δt in , t c + Δt in ]. Wherein Δt in is a preset core event time length (for example, 1 minute), and the data set of this stage independently represents the strong disturbance process of the core influence of turbulence on the aircraft.

[0102] Post-event stage process data is all flight state data in a time window of [t c , t c + ΔT post ]. Wherein ΔT post is a preset subsequent time length (for example, 3 minutes), and the data set of this stage independently represents the process of the decay of turbulence influence and the gradual recovery of the aircraft and atmospheric state to stability.

[0103] S120, respectively extracting the time domain statistical features and the frequency domain statistical features of the data of each stage, and averaging the same features of the data of each stage to construct average feature templates representing the typical state of each stage.

[0104] Wherein, the time domain statistical features and the frequency domain statistical features of the data of each stage are obtained in the following manner:

[0105] S1201, time domain statistical feature acquisition:

[0106] For the vertical acceleration sequence of each stage, the following statistical quantities are directly calculated from the time sequence data thereof:

[0107] Mean: representing the average level of the vertical acceleration of the stage.

[0108] Variance: representing the dispersion degree of the data around the mean.

[0109] Skewness: representing the asymmetry of the data distribution.

[0110] Kurtosis: representing the sharpness of the data distribution.

[0111] Root mean square value: representing the overall energy level of the vertical acceleration of the stage.

[0112] Threshold exceedance ratio: the number of data points with absolute value exceeding a preset threshold (e.g. ±0.1g) in a sequence, divided by the total length of the sequence, representing the proportion of strong disturbance events.

[0113] S1202, time-domain statistical feature acquisition:

[0114] For each phase of the vertical acceleration sequence a z (t) is subjected to a fast Fourier transform. Specifically, it includes:

[0115] (1) Spectrum estimation

[0116] Hanning window is applied to the data to reduce spectral leakage.

[0117] Perform a fast Fourier transform (FFT) to convert the signal from the time domain to the frequency domain, obtaining a complex spectrum A z (f).

[0118] Calculate the power spectral density estimate: P(f) = |A z (f)| 2 / (N×f s ), where N is the number of data points, and f s is the sampling frequency.

[0119] (2) Frequency domain feature calculation

[0120] Based on the above power spectrum P(f), within the preset effective frequency band [f min , f max ] (for example, f min =0.1Hz, f max =10.0Hz), the following feature quantities (in the following integrals, f is the frequency variable) are calculated:

[0121] Peak frequency: f peak =argmax f (P(f)), that is, the frequency corresponding to the maximum power value of the power spectrum within the effective frequency band.

[0122] Mid-frequency power: calculate the power integral within the preset mid-frequency band [f low , f high ] (for example, f low =0.5Hz, f high =2.0Hz, which is the most sensitive frequency band of the aircraft to turbulence response).

[0123] Spectral barycenter: , representing the center frequency of the power spectrum energy distribution within the effective frequency band.

[0124] Spectral width: characterizing the degree of spread of the power spectrum within the effective frequency band around the center of gravity.

[0125] In actual discrete computation, the above integral is realized by summing the discrete frequency points and their corresponding power spectrum values within the effective frequency band.

[0126] Through S1201 and S1202, a feature vector containing the above time-domain and frequency-domain statistics is calculated for each phase of each turbulence event.

[0127] The same features of the phase data are averaged to construct average feature templates representing the typical state of each phase, including:

[0128] The same features (e.g., "pre-event peak frequency" of all events) of the pre-event phase data of all turbulence events are taken as the arithmetic mean to obtain the average pre-event peak frequency.

[0129] Similarly, the arithmetic mean of all other time-domain and frequency-domain features of this phase is calculated.

[0130] The average values of all features of the pre-event phase data are combined into a vector, i.e., the pre-event average feature template is constructed.

[0131] Similarly, the event average feature template and the post-event average feature template are constructed.

[0132] Finally, the three average feature templates of pre-event, event, and post-event jointly constitute a complete and standardized second-scale turbulence disturbance model. This model statistically describes the typical state evolution of turbulence from inception, occurrence to dissipation in the complete life cycle.

[0133] S130, at least one key indicator factor is determined from the time-domain and frequency-domain statistics by comparing the pre-event average feature template with the event average feature template.

[0134] This step aims to select the most indicative signal of turbulence occurrence from a large number of features to refine the second-scale turbulence disturbance model. It includes the following sub-steps:

[0135] S1301, quantitatively compare the numerical difference of each time-domain and frequency-domain statistic between the pre-event average feature template and the event average feature template.

[0136] Specifically, for each feature F (such as peak frequency), the relative difference DF is calculated to quantify the relative change of the feature from the pre-event to the event phase. The calculation formula of DF is: DF = |F m -Ff | / [(F m +F f ) / 2], where F m Let F be the value of feature F in the average feature template of the event. f This is the value of feature F in the average feature template before the event.

[0137] This formula calculates the relative change in characteristic values ​​between two stages, effectively eliminating the influence of the dimensions and orders of magnitude of different characteristics, thus making the differences between different characteristics comparable.

[0138] S1302, based on the calculated differences of each feature and their physical meaning, determine the key indicator factors.

[0139] First, sort all features in descending order of their relative dissimilarity (DF).

[0140] Subsequently, combining knowledge of fluid mechanics and aerospace dynamics, features with clear physical significance are selected from the top-ranked features as key indicator factors. In one embodiment of the present invention, the selected key indicator factors include:

[0141] Peak frequency: Since it directly reflects the dominant frequency of the aircraft response excited by turbulent vortices, its significant forward or backward shift is direct evidence of changes in the energy structure of the flow field.

[0142] Mid-frequency power: Because it characterizes the disturbance energy in the most sensitive frequency band of an aircraft, its sharp increase during an event is a core indicator of turbulence intensity.

[0143] Threshold exceedance ratio: Because it intuitively quantifies the cumulative time of strong acceleration disturbances, it is a direct correlation factor with the perceived intensity of turbulence for pilots and passengers.

[0144] S140, together with the pre-event average feature template, the mid-event average feature template, and the post-event average feature template, and the determined key indicator factors, constitutes the second-scale turbulence disturbance model.

[0145] In this embodiment of the invention, the second-scale turbulence perturbation model is a structured knowledge base, consisting of the following two parts:

[0146] The three average characteristic templates—the pre-event average characteristic template, the mid-event average characteristic template, and the post-event average characteristic template—together define the standard reference state for the turbulence lifetime.

[0147] The key indicator factors identified from the template, along with their calculation methods, pinpoint the core signals that need to be focused on and tracked in real-time early warning.

[0148] S200, during the flight of the aircraft, real-time acquisition of time-series flight state data and corresponding first-scale meteorological data, and extraction of first-scale meteorological features and second-scale turbulence features based on the acquired data.

[0149] Two types of data are acquired in real time through an airborne system and a data link:

[0150] Time-series flight state data: high sampling rate (not less than 1 Hz) flight parameters are acquired in real time from the aircraft bus (such as ARINC429), including but not limited to vertical acceleration, pitch angle, roll angle, airspeed, and latitude and longitude coordinates.

[0151] First-scale meteorological data: through the aircraft communication addressing and reporting system or other data links, the latest gridded numerical weather prediction product data covering the current and future route is acquired.

[0152] Based on the acquired real-time first-scale meteorological data, physical quantities for representing macro-background risk are calculated. The first-scale meteorological features are obtained by querying numerical prediction products or performing physical formula calculation, including:

[0153] Vortex dissipation rate: directly read from the turbulence diagnostic product provided by the meteorological data, or calculated according to the turbulence kinetic energy parameterization scheme.

[0154] Vertical wind speed shear: the difference between the horizontal wind speed vectors of different pressure layers (such as 250 hPa and 300 hPa) is calculated.

[0155] Temperature advection: calculated according to the wind field and temperature field gradient by the formula -U h ·▽T, where U h is the horizontal wind vector and ▽T is the temperature gradient.

[0156] Based on the real-time acquired time-series flight state data, indices for representing the local disturbance state in which the aircraft is currently located are calculated. The second-scale turbulence features are obtained by statistical or signal processing on the flight parameter sequence within a time window, including:

[0157] Vertical acceleration variance: the variance of the vertical acceleration sequence within a sliding time window is calculated.

[0158] Aircraft attitude change rate: the sum or root mean square value of the absolute values of the pitch angle and roll angle change rates within a sliding time window is calculated.

[0159] Airspeed fluctuation: the standard deviation or detrended fluctuation amplitude of the airspeed sequence within a sliding time window is calculated.

[0160] S300, generating a background trend term based on the first scale meteorological feature through the first scale meteorological risk assessment model, and generating a nonlinear correction term from the second scale turbulence feature based on at least one key indicator factor determined in the construction process of the second scale turbulence disturbance model; using the ensemble Kalman filtering algorithm, fusing the background trend term and the nonlinear correction term to generate a prediction result of atmospheric turbulence intensity.

[0161] Further, generating a background trend term based on the first scale meteorological feature through the first scale meteorological risk assessment model, specifically comprising:

[0162] S301, inputting the first scale meteorological feature into the first scale meteorological risk assessment model to obtain a meteorological risk degree, and converting the meteorological risk degree into the background trend term.

[0163] Specifically, the meteorological risk degree is determined by inputting the first scale meteorological feature acquired in real time into the first scale meteorological risk assessment model, the first scale meteorological risk assessment model processing the input feature and outputting a scalar value as the meteorological risk degree. The scalar value is a continuous value between 0 and 1, representing the probability P of turbulence occurrence under the current and future short-term (for example, future 0-30 minutes) macroscopic meteorological background. LS Or risk score. The higher the value, the more conducive the macroscopic condition is to the generation and development of turbulence.

[0164] In order to be integrated with the subsequent ensemble Kalman filtering framework based on physics, it is necessary to convert the dimensionless probability prediction into a physically meaningful driving quantity. The present application converts the probability into a physical quantity with the dimension of vertical wind speed change rate (for example, m / s 2 ) through a preset and calibrated conversion relationship, that is, the background trend term Trend LS . The conversion relationship can be a lookup table obtained based on historical data statistical analysis, or a simple linear or nonlinear mapping function, for example: Trend LS = α × P LS + β, where α and β are coefficients determined by data fitting.

[0165] After this conversion, the physical meaning of the background trend term is clear: it is a predictive vertical wind speed evolution tendency dominated by macroscopic meteorological conditions. Its value is positive, indicating that the macroscopic background field tends to enhance the vertical wind speed and is conducive to the development of turbulence; its value is negative, indicating that it tends to suppress the vertical wind speed.

[0166] Further, the at least one key indicator determined in the process of constructing the second scale turbulence disturbance model is used to generate a non-linear correction term from the second scale turbulence characteristics, specifically including:

[0167] S302, from the second scale turbulence characteristics, the characteristic value of the key indicator determined by the second scale turbulence disturbance model is calculated.

[0168] Specifically, from the second scale turbulence characteristics extracted in real time in step S200, the characteristic value corresponding to the key indicator is selected.

[0169] The characteristic value of the key indicator is calculated within a preset sliding time window. The length of the preset sliding time window is 10-30 seconds. Preferably, the length of the preset sliding time window is 30 seconds, so as to provide timely turbulence intensity estimation under the premise of ensuring statistical reliability. The key indicator includes but is not limited to at least one of the following:

[0170] Spectrum peak factor: its value is the frequency value corresponding to the maximum amplitude in the power spectrum after performing spectrum analysis (such as fast Fourier transform, FFT) on the vertical acceleration data in the preset sliding time window, and the unit is usually hertz (Hz).

[0171] Mid-frequency band energy factor: its value is the power spectrum density integral value of the vertical acceleration data in the preset sliding time window in the preset mid-frequency band (for example, 0.5-2.0 Hz), which represents the turbulence energy intensity in this frequency band.

[0172] Threshold exceeding factor: its value is the percentage of the number of sample points whose vertical acceleration exceeds the preset intensity threshold (for example, ±0.1g) in the total number of sample points in the preset sliding time window, which directly reflects the intensity of strong disturbance events.

[0173] S303, the characteristic value is converted into a correction amount by a preset conversion function to obtain the non-linear correction term.

[0174] Specifically, this step converts the characteristic value of the key indicator calculated in S302 into a unified, dimensionless deviation degree as the direct input of the correction amount. This process includes:

[0175] (1) Calculate the single-factor deviation degree: for each key indicator, compare its characteristic value with a reference value representing a steady flight state obtained from the second scale turbulence disturbance model, calculate the relative change rate of the key indicator as the deviation degree δ r .

[0176] For example, for the peak frequency factor, its deviation δ peak The calculation formula is:

[0177] δ peak = (F peak-real-time - F peak_baseline ) / F peak_baseline .

[0178] Where F peak-real-time is the real-time peak frequency, which refers to the frequency corresponding to the maximum power value identified from the power spectrum calculated by Fast Fourier Transform (FFT) based on the vertical acceleration data in a preset real-time sliding time window (e.g., the past 30 seconds) during the current flight process. This parameter is a scalar that dynamically updates over time. F peak_baseline is the reference peak frequency, which is the peak frequency value representing the typical state before turbulence occurs, directly obtained from the pre-event average feature template of the second scale turbulence disturbance model. This parameter is a constant that remains unchanged after being determined during the model training phase.

[0179] Similarly, the deviation of the mid-frequency band energy factor, threshold overshoot factor, and other factors can be calculated.

[0180] (2) Generate a comprehensive correction precursor signal: weight and sum all the deviation of the key indicator factors according to their preset weights in the second scale turbulence disturbance model to generate a scalar form of precursor signal strength value S precursor . The formula can be expressed as: where w r is the weight of the rth key indicator factor, ΔFactor r is the change rate of the rth key indicator factor, and r takes values from 1 to Q, where Q is the total number of key indicator factors.

[0181] (3) Input the precursor signal strength value S precursor and the current vertical wind speed fluctuation V fluctuation in the real-time second scale turbulence feature into a preset conversion function g() to generate a non-linear correction term Adjust SS : Adjust SS = g(S precursor , V fluctuation ). V fluctuation is directly represented by the vertical acceleration standard deviation in the real-time second scale turbulence feature or obtained by conversion.

[0182] The preset conversion function is configured to have the following characteristics:

[0183] When S precursorWhen S precursor is below a set activation threshold, the output is zero or a tiny amount close to zero, indicating no significant correction is needed; the set activation threshold is determined based on the statistical distribution of S precursor in historical smooth flight data, for example, set between the 90th and 99th percentile of the statistical distribution, corresponding to a value range of about 0.2 to 0.5.

[0184] When S precursor is greater than or equal to the set activation threshold, the output is a monotonic function of S precursor and V fluctuation , and the growth relationship of the monotonic function is nonlinear to simulate the synergistic amplification effect of the two parameters on the turbulence intensity. A typical implementation is that the output is positively correlated with the product or nonlinear combination of the two parameters.

[0185] The specific form of the preset conversion function can be a piecewise linear function, an exponential function, a power function, or a lightweight neural network calibrated by historical data, etc.

[0186] The direction of action of the nonlinear correction term is consistent with that of the background trend term. That is, Adjust SS is always a non-negative value, and its role is to unidirectionally enhance the evolution tendency indicated by the background trend term (whether the trend is rising or falling), thereby achieving a physically self-consistent fine-tuning adjustment.

[0187] Through the above process, the local disturbance information perceived by the aircraft in real time is converted into a dynamic, conditionally activated nonlinear correction quantity, achieving intelligent and physically self-consistent fine-tuning adjustment of the macro trend.

[0188] Further, when generating the nonlinear correction term, a cross-scale coupling factor (denoted as CF) is introduced to modulate the correction strength;

[0189] The cross-scale coupling factor is calculated from the first-scale meteorological feature and the second-scale turbulence feature, and is used to quantify the nonlinear interaction strength between the macro background field and the local disturbance.

[0190] The calculation method of the cross-scale coupling factor CF is: CF = |▽T| • (ΔF peak / F peak_baseline ).

[0191] Where |▽T| is the absolute value of the temperature advection in the first-scale feature, which is used to represent the instability energy of the background field. ΔF peak is the difference between the real-time peak frequency and the reference peak frequency.

[0192] Temperature advection is the key dynamic factor driving the development of weather systems, and the absolute value of temperature advection reflects the potential instability energy of the background field. The change of peak frequency directly reflects the change of the scale of turbulent eddies. When the instability energy of the background field is high (|▽T| is large) and the turbulent eddies are evolving towards a more dangerous scale (ΔF peak ), the product of the two (i.e. coupling factor CF) will significantly amplify the nonlinear correction term, which accurately simulates the core atmospheric physics process that a favorable background field will amplify local disturbances.

[0193] The cross-scale coupling factor CF is positively correlated with the nonlinear correction term. The nonlinear correction term Adjust SS is calculated by the following formula:

[0194] Adjust SS =CF•g(S precursor , V fluctuation ).

[0195] Further, the background trend term and the nonlinear correction term are fused by using the ensemble Kalman filter algorithm to generate the final prediction result of the atmospheric turbulence intensity, which can be realized through the following embodiments.

[0196] (Embodiment 1)

[0197] In this embodiment, the background trend term and the nonlinear correction term are fused by using the ensemble Kalman filter algorithm to generate the final prediction result of the atmospheric turbulence intensity, specifically including:

[0198] S3041, filter system initialization:

[0199] Initialize the ensemble Kalman filter system, and set the number of ensemble members as N (for example, N=50).

[0200] The initial state x0 of each ensemble member is set according to the initial observation or the climatological distribution, representing the possible estimate of the initial vertical wind speed state of the system.

[0201] S3042, state propagation (prediction step):

[0202] At each time t, the state of each ensemble member is forward predicted according to the following state propagation equation to obtain the propagated state PropagatedState:

[0203] PropagatedState=x t ·α persist +(Trend LS +Adjust SS )·(1-α persist );

[0204] Wherein, PropagatedState represents the prior estimate of the vertical wind speed state at future times after fusing the background trend term and the nonlinear correction term, x t α represents the state of the set members at the current time t; persist It is a preset persistent weighting factor used to balance the inertia of the current state with the driving force of new prediction information.

[0205] In this embodiment of the invention, the persistence weighting factor α persist The value of is not fixed, but dynamically configured based on the stability of the macroscopic atmospheric background field. Specifically, the dynamic configuration is as follows:

[0206] When numerical weather prediction products show that the atmospheric stratification is in a stable state, it is α. persist Assign values ​​within the first numerical range to enhance the persistence of the current state; when the state is displayed as unstable, the value is α. persist Values ​​within a second numerical interval are assigned to give higher weight to new information such as background trends and correction terms. The lower limit of the first numerical interval is greater than the upper limit of the second numerical interval. The first numerical interval is [0.45, 0.55], and the second numerical interval is [0.2, 0.35]. The stability of the atmospheric stratification is determined using the Richardson number or convective available potential energy in the numerical weather prediction product.

[0207] Furthermore, the α persist The specific values ​​were determined through cross-validation using historical data (sample size ≥ 1000), and the optimization objective was to minimize the root mean square error of the ensemble Kalman filter state prediction.

[0208] α persist From a hyperparameter that required manual tuning, it has been transformed into an intelligent switch with clear physical meaning and adaptive adjustment. This is no longer a simple parameter optimization, but a methodological innovation that deeply embeds domain knowledge (meteorology) into the core of the algorithm, significantly improving the model's adaptability and forecast accuracy.

[0209] This equation of state is expressed through α persist It dynamically balances the inertial continuation of the state with multi-scale driving terms (Trend) LS +Adjust SS The update function of ). PropagatedState represents the set of prior estimates of the vertical wind speed state at future time moments after fusing multi-scale information.

[0210] S3043, Observation assimilation (update step):

[0211] When there is new observation data (e.g. vertical wind speed from real-time vertical acceleration inversion of aircraft), the ensemble Kalman filter performs a standard update step.

[0212] This step compares the prior estimate of PropagatedState with real observations, and according to the uncertainty of the two, calculates the Kalman gain, modifies the prior estimate, and obtains the final posterior state estimate x{t+1}.

[0213] S3044, final forecast generation:

[0214] The updated posterior state estimate (i.e. the predicted future vertical wind speed) is taken as the final output.

[0215] (Example 2)

[0216] This embodiment is based on Example 1, and adds a dynamic disturbance generation process to more accurately represent the nonlinear growth of system uncertainty, so as to better capture the sudden change characteristics near the turbulence. In this embodiment, the ensemble Kalman filter algorithm is used to fuse the background trend term and the nonlinear correction term to generate the final prediction result of the atmospheric turbulence intensity, which specifically includes:

[0217] S304a, filter system initialization:

[0218] Same as Example 1.

[0219] S304b, dynamic disturbance generation:

[0220] According to the matching result of the second scale turbulence disturbance model and the real-time characteristics, the statistical characteristics of the random disturbance term η to be applied to the state propagation process are dynamically generated and adjusted. The specific generation process is as follows:

[0221] (1) Trigger and input: This process is driven by the precursor signal strength S precursor calculated in step S302. precursor As a core input parameter, it reflects the potential risk of turbulence occurrence in real time.

[0222] (2) Disturbance variance calculation:

[0223] The variance σ 2 of the random disturbance is a monotonically increasing function of S precursor , and this relationship is realized through a preset mapping function h(): σ 2 =h(S precursor ).

[0224] The mapping function h() can be a linear function (such as σ 2 =k×S precursor, k is a coefficient greater than 0), or a nonlinear function (such as an exponential function) to ensure that the variance can be quickly amplified in the high-risk interval.

[0225] (3) Adjustment of the distribution form of the disturbance:

[0226] The probability distribution form of the random disturbance adopts a mixed model of Gaussian distribution and Laplace distribution.

[0227] The mixing weight is S precursor Dynamic control. Specifically, a mixing coefficient λ1 is defined, whose value is a function of S precursor , i.e. λ1=f(S precursor ), and f() is monotonically increasing, with a value range of [0, 1].

[0228] Subsequently, according to the calculated target total variance σ 2 , the variances of the two components in the mixed distribution are allocated respectively:

[0229] The variance of the Gaussian distribution component is: σ 2 G =k G •σ 2 ;

[0230] The variance of the Laplace distribution component is: σ 2 L =k L •σ 2 ;

[0231] Where k G and k L are preset variance allocation coefficients, and satisfy k G +k L =1.

[0232] The final disturbance distribution p(η) is defined as:

[0233] p(η)=(1-λ1)×N(η|0, σ 2 G )+λ1×L(η|0, b)

[0234] Where N(η|0, σ 2 G ) is a Gaussian distribution with mean 0 and variance σ² G .

[0235] L(η|0, b) is a Laplace distribution with scale parameter b=σ L / 2 1 / 2 .

[0236] σ² G and σ² LTotal variance σ calculated from step (2) 2 Proportionally distribute.

[0237] In the embodiments of the present application, when the value of S precursor is less than or equal to a first threshold, λ1 approaches 0, and the distribution p(η) is closer to Gaussian, resulting in mild and smooth perturbations. The first threshold is determined by statistics (e.g., taking the 95th percentile of the distribution) of S precursor values in the historical data for "no turbulence" or "clear sky" samples.

[0238] When the value of S precursor is greater than or equal to a second threshold, λ1 approaches 1, and the distribution p(η) is more biased towards Laplace with heavy-tailed characteristics, thus being able to generate extreme perturbation values with larger amplitude and more explosive with higher probability. The second threshold is determined by statistics (e.g., taking the 5th percentile of the distribution) of S precursor values in the historical data within a certain time window (e.g., 5 minutes before the event) before "moderate and above turbulence" events.

[0239] When the value of S precursor is between the first threshold and the second threshold, it is in a transition state. In this state, the mixing coefficient λ1 can be smoothly changed between 0 and 1 according to a pre-set transition function (e.g., linear interpolation).

[0240] (4) Perturbation term generation:

[0241] At each time, according to the calculated variance σ 2 and the mixing coefficient λ1, a random sample is taken from the above mixed distribution p(η).

[0242] For each set member, a perturbation value is independently sampled, forming a perturbation vector of the same size as the set, which is used in the state propagation equation of S304c.

[0243] Through the above process, the random perturbation term η is generated dynamically and adaptively, which can be used as an accurate uncertainty controller, significantly improving the prediction ability of the ensemble Kalman filter for turbulence, a nonlinear sudden change phenomenon.

[0244] S304c, state propagation (prediction step):

[0245] At each time t, the state of each ensemble member is forward predicted according to the following state propagation equation:

[0246] PropagatedState = x t · α persist + (Trend LS + AdjustSS + η) · (1 - α persist ); wherein, η is a random disturbance term with dynamic statistical characteristics generated from S304b step.

[0247] S304d, observation assimilation (update step):

[0248] Same as embodiment 1.

[0249] S304e, final forecast generation:

[0250] Same as embodiment 1.

[0251] S400, generating a turbulence warning signal according to the prediction result.

[0252] Further, S400 specifically includes:

[0253] S401, calculating the eddy dissipation rate at the future track point based on the vertical wind speed state obtained by the set Kalman filter prediction.

[0254] Based on the future vertical wind speed sequence w(t) obtained by the prediction of S303 step, the running standard deviation σ w of the vertical wind speed within a set time window (such as 10 seconds) is calculated. Then, the eddy dissipation rate is calculated according to the internationally recognized algorithm based on the aircraft response. The physical symbol of the eddy dissipation rate is , and its unit is m 2 / s 3 . In order to establish a more direct connection with the turbulence intensity in aviation applications, the International Civil Aviation Organization (ICAO) and the aviation meteorology field usually use the cubic root of the eddy dissipation rate, that is, . As a standard turbulence intensity measurement index, , the unit of which is m 2 / 3 / s. In the present application, the term EDR specifically refers to this physical quantity. The EDR calculation formula is:

[0255] ;

[0256] wherein, σ w is the running standard deviation of the vertical wind speed calculated within a 10-second sliding window, Va is the true airspeed, ω1 and ω2 are the cutoff frequencies, and ω1 is 0.15 Hz and ω2 is 2 Hz, C is the model constant, and is usually taken as 1.05.

[0257] S402, comparing the eddy dissipation rate with a plurality of preset turbulence level thresholds.

[0258] The preset multiple turbulence level thresholds are determined by referencing the turbulence intensity classification standards in the International Civil Aviation Organization (ICAO) Flight Weather Information (Doc 8896) and by systematically calibrating the EDR calculation model of this invention based on data from over 1000 historical measured turbulence events. The specific thresholds are as follows:

[0259] Mild bumps: EDR < 0.15m 2 / 3 / s, corresponding phenomenon: slight shaking of the aircraft, with no obvious change in altitude / airspeed;

[0260] Moderate bumps: 0.15m 2 / 3 / s≤EDR<0.45m 2 / 3 / s, corresponding phenomenon: aircraft swaying noticeably, altitude / airspeed fluctuations ±30m / ±10kt);

[0261] Severe turbulence: EDR ≥ 0.45m 2 / 3 / s, corresponding phenomenon: the aircraft shakes violently, and the altitude / airspeed fluctuates significantly (>±30 meters / >±10 knots), posing a risk of injury to personnel.

[0262] S403 generates and outputs corresponding graded turbulence early warning signals based on the comparison results.

[0263] Based on the comparison results, corresponding graded turbulence warning signals (such as "mild", "moderate", "severe") are generated and output, along with recommended avoidance suggestions.

[0264] Furthermore, the method also includes the following steps:

[0265] S500, based on the turbulence warning signal, generates a spatial distribution map of turbulence risk downstream of the flight path, providing pilots with intuitive route planning assistance information. Specifically, it includes the following sub-steps:

[0266] S501, based on the aircraft's precise current location, numerical weather prediction 3D wind field data, and the location information of turbulence warning signals, simulates the transport and diffusion path of turbulence signals along the downstream route over a preset time period using a Lagrange particle diffusion model, outputting the 3D spatial trajectory data of all virtual particles evolving over time. Specifically, this includes:

[0267] Particle Release: Using the latitude, longitude, and altitude of the turbulence warning signal as the source point, a set of N (N≥1000) virtual particles is released at the initial moment to ensure sufficient statistical representativeness of the diffusion path.

[0268] Advection transport: At each time step Δt (e.g., 30 seconds), based on the three-dimensional wind field (U, V, W) provided by numerical weather prediction, the wind speed at the current position of each particle is calculated using bilinear interpolation, and the coordinates of the particle are updated accordingly: new position = original position + wind speed × Δt.

[0269] Turbulent diffusion: To simulate turbulent mixing at the subgrid scale, a random perturbation (Δx, Δy, Δz) is added to the three-dimensional coordinates of each particle after the advection calculation at each time step. This perturbation follows a normal distribution with a mean of zero and a variance equal to 2g0×Δt, where g0 is the turbulent diffusion coefficient, the value of which is obtained by querying a lookup table determined by the Richardson number (Ri), and the Ri value is provided by numerical weather prediction products.

[0270] S502, based on the particle trajectory data output by S501, the intensity value of the turbulence warning signal, and preset attenuation parameters, assigns the turbulence warning signal as the initial source strength to the particles, calculates the intensity attenuation of the particles during the diffusion process, and generates a structured turbulence risk area dataset accordingly. Specifically, it includes:

[0271] Intensity binding and decay: Each particle in S501 is assigned an initial turbulence intensity I0, equal to the warning EDR value. This intensity decays exponentially with simulation time t: I(t) = I0•exp(-t / τ), where I(t) is the turbulence intensity at time t, in m³ / s. 2 / 3 / s, where I0 is the initial turbulence intensity value in meters. 2 / 3 / s, where t is the simulation time after particle release in seconds, and τ is the decay time constant used to control the rate of turbulence intensity decay. When the simulation time t equals τ, the turbulence intensity I(t) will decay to 1 / e (approximately 36.8%) of the initial turbulence intensity I0. The value of τ needs to be determined by calibration using observational data from historical turbulence events (e.g., by least-squares fitting of a large amount of historical turbulence event observational data, set to 600 seconds). The larger the value of τ, the slower the turbulence intensity decays, and the longer the prediction validity period. exp() is the natural exponential function.

[0272] Grid-based statistics: The particle activity region is divided into regular latitude and longitude grids (e.g., 0.1° × 0.1°). For each grid, the entire simulation duration (e.g., 30 minutes) is traversed to find all particles that have passed through that grid, and the maximum value of the turbulence intensity when they are within that grid is recorded. This maximum value is used as the risk value for that grid.

[0273] Classification: According to the same ICAO standard-based turbulence level threshold as in S400, the risk value of each grid is classified into "mild", "moderate", "severe" levels, and a structured dataset containing the geographical boundaries of the grid and the corresponding risk level code is output.

[0274] S503, based on the turbulence risk area dataset output by S502 and the electronic navigation map data, superimposes the turbulence risk area dataset on the electronic navigation map for visual rendering. Specifically, it includes:

[0275] The risk level code output by S502 is mapped to a preset visual style. For example, "mild" is mapped to semi-transparent green, "moderate" is mapped to semi-transparent yellow, and "severe" is mapped to semi-transparent red.

[0276] On the onboard electronic navigation map system, the risk area dataset is rendered as a new graphic layer and superimposed on the underlying map at the pixel level.

[0277] The risk layer is dynamically updated at a preset refresh rate (e.g., every 1 minute) to ensure that the displayed risk information is consistent with the real-time simulation results, providing continuous downstream route risk situational awareness for pilots.

[0278] The technical effect of S500 is to upgrade single-point, abstract early warning information to intuitive and operable route planning tools. It solves the core pain point in the prior art that "knows there is turbulence, but does not know the impact range", and provides forward-looking situational awareness and decision support for pilots.

[0279] The embodiment of the application also provides an electronic device, comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, and the instructions are configured to execute the method described in the embodiment of the application.

[0280] The embodiment of the application also provides a computer-readable storage medium storing computer executable instructions, and the computer executable instructions are used to execute the method described in the embodiment of the application.

[0281] It should be understood that the steps shown above can be reordered, added or deleted using various forms of flow. For example, the steps described in the present application can be executed in parallel, sequentially or in different orders, as long as the desired results of the technical solutions disclosed in the present application can be achieved, and this document does not limit it here.

[0282] The above detailed description does not limit the scope of the application. Various modifications, combinations, sub-combinations and alternatives can be made to the detailed description. Any modification, equivalent replacement and improvement etc. made within the spirit and principle of the application shall be included in the scope of the application.

Claims

1. A data processing method for aircraft atmospheric turbulence warning, characterized in that, The method comprises the following steps: S100, based on target historical meteorological data and target time-series flight state data, offline constructing a first scale meteorological risk assessment model and a second scale turbulence disturbance model; wherein the spatial coverage of the first scale is larger than that of the second scale, and the construction process of the second scale turbulence disturbance model comprises determining at least one key indicator factor; S200, in the process of aircraft flight, real-time acquiring time-series flight state data and corresponding first scale meteorological data, and extracting first scale meteorological features and second scale turbulence features based on the acquired data; S300, based on the first scale meteorological features, generating a background trend term through the first scale meteorological risk assessment model, and generating a nonlinear correction term from the second scale turbulence features based on at least one key indicator factor determined in the construction process of the second scale turbulence disturbance model; using an ensemble Kalman filtering algorithm to fuse the background trend term and the nonlinear correction term to generate a prediction result of atmospheric turbulence intensity; S400, generating a turbulence warning signal according to the prediction result; The method comprises the following steps: S100, based on target historical meteorological data and target time-series flight state data, offline constructing a first scale meteorological risk assessment model and a second scale turbulence disturbance model; wherein the spatial coverage of the first scale is larger than that of the second scale, and the construction process of the second scale turbulence disturbance model comprises determining at least one key indicator factor; The method comprises the following steps: From the second scale turbulence features, the characteristic values of the key indicator factors determined by the second scale turbulence disturbance model are calculated; the key indicator factors include at least one of the following: a spectral peak factor, a mid-frequency band energy factor, and a threshold exceeding factor; The characteristic values are converted into correction amounts by a preset conversion function to obtain the nonlinear correction term, which comprises: (1) calculating a single-factor deviation degree; (2) generating a comprehensive correction precursor signal: weighting and summing all key indicator factor deviation degrees according to the preset weights of the key indicator factors in the second scale turbulence disturbance model to generate a scalar form of precursor signal intensity value; (3) inputting the precursor signal intensity value and the current vertical wind speed fluctuation in real-time second scale turbulence features into a preset conversion function to generate a nonlinear correction term; S400 specifically comprises: Based on the vertical wind speed state obtained by the ensemble Kalman filtering prediction, the vortex dissipation rate at the future track point is calculated; The vortex dissipation rate is compared with a plurality of preset turbulence level thresholds; According to the comparison result, the corresponding graded turbulence warning signal is generated and output.

2. The method of claim 1, wherein, The target historical meteorological data and the target time-series flight state data are acquired by the following steps: S10, extracting the occurrence position coordinates of all turbulence events recorded in the historical flight state data to form a spatial point set to be clustered; S11, preset initial neighborhood radius and minimum sample number parameters of a density-based noise application spatial clustering algorithm according to geographical distribution characteristics of the spatial point set; S12, perform density-based noise application spatial clustering analysis on the occurrence position coordinates of the turbulence event based on the preset initial neighborhood radius and minimum sample number parameters, to obtain a plurality of clustering clusters; S13, calculate turbulence event density of each clustering cluster, and identify a geographical range represented by a clustering cluster with the highest turbulence event density as a target training area; S14, extract data corresponding to the spatiotemporal range of the target training area from a historical meteorological database and a time-series flight state data set, and use the data as target historical meteorological data and target time-series flight state data.

3. The method of claim 1, wherein, In S100, the process of constructing the first-scale meteorological risk assessment model includes: S101, spatiotemporally match the target historical meteorological data and the target time-series flight state data, calculate a plurality of turbulence prediction indexes based on the matched data, and select turbulence prediction indexes with a prediction performance representation value reaching a preset threshold to form a basic feature set; S102, pre-process the basic feature set to obtain a pre-processed feature time series; the pre-processing includes outlier processing, data standardization, and data transformation; S103, perform signal decomposition and screening on the pre-processed feature time series to extract key dynamic modes related to turbulence generation and dissipation mechanisms; S104, train an integrated classification model based on the screened key dynamic modes to form the first-scale meteorological risk assessment model.

4. The method of claim 1, wherein, In S100, the process of constructing the second-scale turbulence disturbance model includes: S110, identify turbulence events from the target time-series flight state data, and divide process data of the turbulence events into pre-event stage process data, event stage process data, and post-event stage process data based on the peak time of vertical acceleration; S120, extract time-domain statistical features and frequency-domain statistical features of each stage process data, and average the same features of each stage process data to construct average feature templates representing typical states of each stage; S130, determine at least one key indicator from the time-domain statistical features and the frequency-domain statistical features by comparing the pre-event average feature template and the event average feature template; S140, use the pre-event average feature template, the event average feature template, the post-event average feature template, and the determined key indicator to jointly form the second-scale turbulence disturbance model.

5. The method of claim 1, wherein, Fusing the background trend term and the nonlinear correction term through the ensemble Kalman filtering algorithm includes a state propagation step, which is realized through the following state propagation equation: PropagatedState = x t • a persist +(Trend LS + Adjust SS ) • (1 - a persist ); wherein PropagatedState represents the prior estimate of the vertical wind speed state at future time instant after fusing the background trend term and the non-linear adjustment term, x t denotes the set-membership state at current time instant t; a persist is a pre-defined persistence weight factor; Trend LS is the background trend term; Adjust SS is the non-linear adjustment term.

6. The method of claim 5, wherein, When generating the nonlinear correction term, a cross-scale coupling factor is introduced to modulate the correction strength; The cross-scale coupling factor is calculated from the first-scale meteorological features and the second-scale turbulence features.

7. An electronic device, comprising: It includes a processor and a memory; The processor is configured to execute the steps of the method according to any one of claims 1 to 6 by invoking the program or instructions stored in the memory.

8. A computer-readable storage medium, characterized in that, The computer readable storage medium is configured to store the program or instructions, which cause the computer to execute the steps of the method according to any one of claims 1 to 6.

Images