A high-resolution modeling method and system for low-altitude turbulence intensity fields

Abstract

This invention relates to the field of low-altitude airspace meteorological monitoring and support technology, and provides a high-resolution modeling method for low-altitude turbulence intensity fields, including: S1: data preprocessing and spatiotemporal alignment, performing time synchronization, coordinate unification, data cleaning, missing value imputation, and outlier removal on multi-source heterogeneous detection data; S2: three-level fusion of multi-source data based on turbulence dissipation rate, obtaining a priori background field based on the data preprocessed in step S1; S3: constructing a variational assimilation objective function with turbulent kinetic energy (TKE) as the core state variable based on various key parameters inverted in step S2 and the obtained prior background field, and solving the optimal three-dimensional turbulence intensity analysis field using the L-BFGS-B optimization algorithm; S4: result output and visualization. This method can effectively fuse meteorological observation data from different sources and at different resolutions, and reconstruct a high-resolution, high-integrity three-dimensional turbulence intensity field that meets the safety requirements of low-altitude flight using variational assimilation theory.

Description

Technical Field

[0001] This invention relates to the technical field of low-altitude airspace meteorological monitoring and support, and particularly to a high-resolution modeling method and system for low-altitude turbulence intensity fields. Specifically, it involves fusing multi-source heterogeneous meteorological observation data, based on three-level fusion and variational assimilation, and using variational assimilation theory to reconstruct a high-resolution three-dimensional low-altitude turbulence intensity field. Background Technology

[0002] Low-altitude airspace (typically below 1000 meters above sea level) serves as the primary operating space for emerging aircraft such as unmanned aerial vehicles (UAVs) and electric vertical takeoff and landing (eVTOL) aircraft, and its atmospheric environment is highly complex and unique. Unlike the relatively stable and homogeneous atmospheric conditions of mid- and high-altitude regions, the low-altitude atmosphere is directly and strongly influenced by surface friction, uneven heating (such as urban heat islands and temperature differences between farmland and roads), complex terrain (such as hills and canyons), and man-made structures (such as building complexes and towers). These factors work together to generate turbulent vortices ranging in size from tens to hundreds of meters, roughly equivalent to the wingspan or fuselage size of a small aircraft. This scale resonance effect makes aircraft highly susceptible to being engulfed by turbulence, resulting in severe turbulence and even loss of control, seriously threatening flight safety and operational stability. This has become a key meteorological bottleneck restricting the large-scale and safe development of low-altitude airspace.

[0003] Existing turbulence monitoring technologies rely on an observation network composed of various sensors, primarily including wind profiler radar, Doppler weather radar, and automatic weather stations. Wind profiler radar provides high-resolution wind field information in single-point vertical profiles, but its horizontal coverage is limited and stations are sparse. Doppler weather radar (such as S-band and X-band) has broad horizontal coverage and can detect radial wind speed and spectral width, but its data exhibits a conical distribution, its detection capability decreases under clear skies, and its vertical resolution is insufficient. Automatic weather station networks provide high-precision near-surface meteorological data, but completely lack vertical detection capabilities in the air. These devices exhibit significant differences in spatiotemporal resolution (from 1 minute to 10 minutes, from 60 meters to several kilometers), detected elements, and coverage, resulting in a severe "data heterogeneity" problem and creating information silos.

[0004] The aforementioned data heterogeneity means that no single data source can independently meet the stringent requirements of low-altitude flight turbulence monitoring, demanding "high resolution (100 meters horizontally, 50 meters vertically), high precision, and three-dimensional continuity." Therefore, effectively fusing these spatiotemporally mismatched and physically diverse observational data to generate a physically plausible, accurate three-dimensional continuous field that characterizes small-scale turbulence structures has become a core scientific problem and technical challenge. Traditional simple interpolation or single-point extrapolation methods generate numerous non-physical assumptions in sparse data regions, failing to accurately reflect the intermittency and spatial structure of turbulence, and thus hindering refined flight path planning and real-time safety warnings. This necessitates the development of an advanced algorithmic model capable of deeply fusing multi-source heterogeneous data while conforming to the physical laws of atmospheric turbulence. Summary of the Invention

[0005] To address the aforementioned problems, the present invention aims to provide a high-resolution modeling method and system for low-altitude turbulence intensity fields. This method effectively integrates meteorological observation data from different sources and at different resolutions, and utilizes variational assimilation theory to reconstruct a high-resolution, high-integrity three-dimensional turbulence intensity field that meets the safety requirements for low-altitude flight.

[0006] The above-mentioned objective of this invention is achieved through the following technical solutions: A high-resolution modeling method for low-altitude turbulence intensity fields includes the following steps: S1: Perform data preprocessing and spatiotemporal alignment, acquire detection data from ground automatic weather stations, wind profiler radar, S-band and X-band Doppler weather radar, and perform time synchronization, coordinate unification, data cleaning, missing value imputation and outlier removal on multi-source heterogeneous detection data to achieve spatiotemporal alignment; S2: Three-level fusion of multi-source data based on turbulent dissipation rate. Based on the data preprocessed in step S1, the turbulent dissipation rate ε is inverted using the structure function method and the vertical profile of turbulent kinetic energy (TKE) is calculated using wind profiler radar data. The spectral width component caused by pure turbulence is extracted using Doppler weather radar spectral width data and a statistical relationship model between it and TKE is established. The near-surface gradient Richardson number is calculated using data from ground automatic weather stations and a statistical calibration relationship between it and near-surface TKE is constructed. The prior background field is obtained through three-level fusion. S3: Variational assimilation optimization modeling and solution. Based on the various key parameters inverted in step S2 and the obtained prior background field, a variational assimilation objective function with turbulent kinetic energy TKE as the core state variable is constructed. By introducing background field constraints, observation operators, and various physical constraints, the optimal three-dimensional turbulence intensity analysis field is solved using the L-BFGS-B optimization algorithm. The three-dimensional analysis field meets the grid specifications of 100 meters horizontal resolution, 50 meters vertical resolution, and 0-2 kilometers vertical range. S4: Results Output and Visualization. Outputs a three-dimensional turbulence intensity analysis field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters. This three-dimensional turbulence intensity field is overlaid with geographic information data to generate isosurface maps, horizontal slice maps, and vertical profile maps, realizing a three-dimensional spatial visualization of low-altitude turbulence intensity.

[0007] Furthermore, step S1 specifically includes: Step S11: Doppler radar data processing: Valid echoes are screened based on the signal-to-noise ratio threshold, and the missing values ​​of radial velocity and velocity spectrum width are filled by a two-dimensional interpolation method based on range and azimuth angle. Abnormal radial velocity values ​​that exceed the physical range are identified and corrected using the 3σ criterion. Step S12: Data processing of ground automatic weather stations: unify near-surface temperature, wind speed, air pressure and other data into WGS84 coordinate system and UTC time format, use inverse distance weighted interpolation to spatially interpolate missing values, and use the 3σ criterion to identify and correct outliers. Step S13: Wind profiler radar data processing: Wavelet denoising method is used to smooth the wind speed and wind direction time series. Based on the assumption of vertical continuity of atmospheric motion, linear interpolation and nonlinear interpolation methods are used to fill in the missing data at different height levels caused by weak echoes. Step S14: Spatiotemporal unification: Align all observation data in time to a 6-10 minute update frequency, and project and interpolate the radar polar coordinate data to the target Cartesian grid system in space.

[0008] Furthermore, step S2 specifically includes: Step S21: Vertical baseline construction: Using the high-resolution vertical wind profile provided by the wind profiler radar, the turbulent dissipation rate ε at each height layer is calculated using the second-order structure function method; combined with the shear generation term P calculated from the vertical wind speed gradient, and the buoyancy term B determined based on the imaginary potential temperature gradient and the gradient Richardson number, the vertical turbulent kinetic energy TKE profile is obtained by inversion through the turbulent energy balance equation: Where z is the target height, This is the near-ground reference height. Data provided by ground stations; the vertical resolution of this TKE profile is 50m, but it can only reflect single-point features in the horizontal direction, and needs to be extended to a two-dimensional plane through a second-level fusion. The shear generation term P represents the shear of the mean wind field, that is, the transfer of energy from wind speed or direction to turbulence through dynamic instability as wind speed or direction changes with height. This is the primary source of TKE (Total Kinematic Energy). For a given target height z, its calculation formula is: in, It is the vertical gradient of the horizontal wind speed *u* at height *z*, i.e., the vertical wind shear. It is the momentum turbulence exchange coefficient; the buoyancy term B characterizes the promoting or inhibiting effect of thermal instability or stable stratification on turbulence development. When the atmosphere is unstable, B > 0, and buoyancy provides energy to the turbulence; when the atmosphere is stable, B < 0, and buoyancy consumes turbulence energy. It is determined by the vertical thermal structure of the atmosphere, i.e., the stratification stability, and its calculation formula is: Where g is the acceleration due to gravity. This is the potential temperature after taking into account air humidity, also known as the virtual potential temperature, which more accurately reflects the buoyancy effect of the atmosphere. It is the vertical gradient of the virtual potential temperature at height z. It is the turbulent heat exchange coefficient; Step S22: Horizontal spatial constraint correction: Perform multi-factor calibration on the original Doppler radar spectral width, eliminate instrument noise, wind shear contribution and precipitation interference, and extract the spectral width component caused by pure turbulence. Based on the power function statistical relationship, the spectral width is converted into a TKE value: Where k is the scaling factor, with a value of 0.85, and m is the exponent, with a value of 1.2; Step S23: Adaptive Weight Fusion: The TKE results retrieved from radar are matched to the vertical reference grid using Kriging interpolation. Fusion is achieved through an adaptive weighting function. The fusion formula is as follows: Where x is the horizontal coordinate of the target Cartesian grid, and z is the height; For the turbulent kinetic energy after fusion, For the vertical reference profile based on wind profiler radar, Turbulent kinetic energy obtained by Doppler radar inversion and interpolation; weighting coefficients Wind profiler radar data quality factor Doppler radar range attenuation factor Joint decision; Step S24: Near-Earth Boundary Calibration: Calculate the gradient Richardson number using data from ground-based automatic weather stations. The calculation formula is as follows: in, denoted by potential temperature, g is the acceleration due to gravity, z is the altitude, and u and v are the components of the horizontal wind in the east-west and north-south directions, respectively. Construct the statistical calibration formula: in, For the gradient Richardson number, The average wind speed at a height of 10m is given. The calculation result of the above formula is corrected for deviation from the near-surface TKE obtained by the second-level fusion, i.e., <200m, with a correction amount of . The 20m height layer is the near-ground height layer. The correction amount is transferred to the 200m height layer through vertical linear interpolation, and finally the three-dimensional turbulence intensity field is obtained.

[0009] Furthermore, in step S23, the formula for calculating the adaptive weight function is: in, This refers to the quality factor of wind profiler radar data. Doppler radar range attenuation factor, wind profiler radar data quality factor The signal-to-noise ratio (SNR) of the radar echo at altitude z is the primary factor determining the data quality. A higher SNR indicates more reliable data quality, and its weight should be increased accordingly. The calculation formula is as follows: in, Here is the measured signal-to-noise ratio at height z. and The effective threshold for signal-to-noise ratio is set based on radar performance. Doppler radar range attenuation factor This characterizes the attenuation of radar detection capability with distance, and is related to the spatial distance R(x, z) from the target grid point to the radar station. The greater the distance, the lower the radar detection accuracy and resolution, and the lower the data weight should be. The calculation formula is as follows: Where R(x,z) is the straight-line distance from the grid point to the radar station. γ is the effective maximum detection range of the radar, and γ is the attenuation exponent, which is usually taken as 1.5 to 2.5 and is used to control the rate at which the weight decays with distance.

[0010] Furthermore, step S3 specifically includes: Step S31: Cost function construction: Construct a cost function that includes background and observation terms to constrain the deviation between the analysis field and the prior background field, and to measure the fit between the analysis field and the observation data; Step S32: Definition of Observation Operators: Define observation operators for Doppler radar, wind profiler radar, and ground observation stations respectively, to establish the physical relationship between the state variable TKE and the original observation values ​​of various sensors, and map the gridded analysis field to the observation space of the k-th type of sensor: Doppler radar observation operator: in, To analyze the turbulent kinetic energy in the field, a and b are model parameters obtained based on statistical analysis of historical data. This power function form can effectively characterize the nonlinear relationship between turbulent energy and radar spectral width observations. Wind profiler radar observation operator: Where c is the proportionality coefficient. Let z be the turbulent mixing length at height z. This operator reflects the physical nature that the greater the turbulent kinetic energy, the stronger the dissipation, and takes into account the variation of the mixing length with height. Ground observation operator: This operator directly extracts the analysis field. The turbulent kinetic energy values ​​at the bottom layer of the grid, corresponding to the near-ground height, are compared with the measured turbulence estimates from the ground station. Step S33: Physical Constraint Settings: Set Non-negativity Constraints Smoothness constraints, boundary constraints, and resolution constraints are used to ensure that the solution conforms to physical laws; Nonnegativity constraint: Turbulent kinetic energy is the sum of squares of fluctuating velocities, and its value must be greater than or equal to zero; Smoothness constraint: The real turbulent field is usually continuously changing in space and will not have an infinitely large gradient. This is achieved through the background field error covariance matrix. Boundary constraints: Based on the fundamental laws of turbulent motion, constraints are imposed on the numerical range and variation trend of ε and TKE; Resolution constraint: Ensure that the model output meets the required spatiotemporal resolution; Step S34: Optimization Solution: The cost function is minimized using the L-BFGS-B algorithm. The analysis field I is iteratively updated using gradient descent. The process terminates when the gradient norm is less than the set tolerance or when the maximum number of iterations is reached. Specific steps are as follows: S341: Initialization: Given an initial guess I0, usually the background field I b ; S342: Iterative Update: In the nth iteration, calculate the gradient of the cost function. : The L-BFGS-B algorithm is used to generate the search direction d. n ; S343: Line search: along d n Perform a one-dimensional search to determine the direction. Ensure that the updated solution satisfies I. i≥0; S344: Convergence judgment: When the gradient norm is less than the set tolerance or the maximum number of iterations is reached, the iteration terminates and the analysis field I is output; Step S35: Model parameter setting and initialization settings: Grid system: horizontal resolution of 100 meters, vertically divided into 40 layers from ground to 2000 meters at 50-meter intervals, horizontal range covers the smallest bounding rectangle of all ground stations, i.e. longitude 118.4°E–119.0°E, latitude 31.8°N–32.6°N; Background field initialization: The climate-averaged turbulent kinetic energy or preliminary interpolation results are used as the initial estimate; Background error covariance matrix B: Based on the Gaussian model, the horizontal correlation length is set to 2 kilometers and the vertical correlation length to 200 meters. The standard deviation is determined based on climate statistics and observation data. Observation error covariance matrix R: Assuming that each observation is independent, it is a diagonal matrix, and the diagonal elements are set according to the sensor accuracy. Numerical Implementation: Based on Python 3.8 or later, relying on the core functions of NumPy and SciPy libraries, the L-BFGS-B algorithm is used for constrained optimization. To address the large grid size, recursive filters or spectral methods are used to quickly calculate the inverse of the background error covariance. Efficient spherical to Cartesian coordinate transformation and interpolation model validity verification are performed on radar data.

[0011] Furthermore, in step S31, the specific form of the cost function based on variational assimilation theory is as follows: in, As a background term, the project's constraint analysis field Do not deviate too much from a priori background field. The background field can be obtained based on climate mean, low-resolution numerical model forecasts or simple interpolation results. It introduces prior knowledge about the spatial structure of the turbulent field. B is the background field error covariance matrix, which adopts a Gaussian correlation function. Its diagonal elements represent the error variance of the background TKE value, while the off-diagonal elements characterize the correlation of background errors at different spatial points, reflecting the spatial continuity and smoothness scale of the TKE field. The observation term is used to measure the relationship between the analysis field and all K types of observation data. The degree of fit; It is an observation operator; The three-dimensional turbulence intensity analysis field to be solved is the state variable, and the variables on its mesh are turbulent kinetic energy or turbulent dissipation rate. B represents the prior background field obtained through three-level fusion, where B is the background error covariance matrix, constructed using a Gaussian correlation function, and its diagonal elements represent the background field. The error variance at different spatial locations, and the off-diagonal elements quantify the correlation between background errors at different spatial points, are represented by a Gaussian correlation function model. The function's form is as follows: Where C is the covariance between two points in space. It is the background error variance. and These are the vertical and horizontal distances between two points, respectively. and These are the correlation scales in the vertical and horizontal directions, respectively, used to control the smoothness of the analysis field; : The kth type of observation data; The observation operator for the k-th type of observation data is a physical model responsible for mapping the state variable I from the model grid space to the observation space, thus allowing it to be compared with the actual observation values. Compare; The observation error covariance matrix of the k-th type of observation data is usually assumed to be a diagonal matrix. The elements on its diagonal represent the uncertainty of various observation instruments, i.e., the error variance, which determines the relative weight of observation data from different sources and with different precision in the cost function. In step S33, the boundary constraints specifically include: Numerical range constraint: Based on the statistical measurements of low-altitude turbulence, the reasonable range for ε is 10. −6 ≤ε≤10 −2 m 2 / s 3 The reasonable range for TKE is 0.1 ≤ TKE ≤ 10 m 2 / s 2 If the calculation exceeds this range, it will be automatically truncated to the boundary value. Vertical variation constraint: The attenuation rate of near-ground TKE with altitude must not exceed 0.002 m. 2 / ( s 2 The absolute value of the vertical gradient of the TKE at mid-to-high altitudes must not exceed 0.005 m. 2 / ( s 2 m), to avoid non-physical vertical jumps; Energy balance constraint: For any grid node, P+B≥0.5ε must be satisfied, and the generation of turbulent energy must not be much less than the dissipation, otherwise the energy conservation law will be violated; In step S33, the resolution constraint specifically includes: Horizontal resolution constraint: The difference between ε and TKE between any two adjacent horizontal grid nodes shall not exceed 20% of the mean of the region. If it exceeds this, Gaussian smoothing shall be used for correction with a smoothing radius of 100m. Vertical resolution constraint: The ε and TKE values ​​of each 50m height layer need to be verified by linear interpolation based on the data of the upper and lower layers. If the deviation exceeds 10%, the coupling weight is recalculated to ensure the continuity of the vertical layering.

[0012] Furthermore, step S4 specifically includes: Step S41: Output of 3D turbulence intensity analysis field: The mesh system corresponding to the output 3D turbulence intensity analysis field satisfies: Horizontal resolution is The value is 100 meters; Vertical direction within the height range Internal settings There are 1 vertical layers, each with a thickness of [missing information]. ; Step S42: Result visualization processing: Extract horizontal slice data from the three-dimensional turbulence intensity field calculated in step S3 at 50-meter vertical intervals; generate contour lines in the horizontal direction at a 100-meter grid resolution; overlay the above gridded data with the digital elevation model and administrative division vector data to generate a three-dimensional turbulence intensity spatial distribution map with geographic reference, a horizontal distribution map of a specified height layer, and a typical profile vertical structure map; Step S43: Results and Verification: The turbulence intensity calculated by the model exhibits a vertically stratified distribution, with strong turbulence near the ground and weak turbulence at higher altitudes. The near-surface layer is affected by both ground friction and thermal convection, resulting in significant airflow fluctuations and ample turbulent energy supply. The upper atmosphere, far from the ground, relies mainly on wind shear to maintain weaker turbulence. The model results accurately capture this physical process, verifying the calculation logic based on the vertical energy balance equation and proving that the model effectively simulates the dissipation of turbulent energy with altitude. To more clearly quantify the stratified variation of turbulence intensity with altitude calculated by the model, and to verify its consistency with atmospheric boundary layer theory and model construction logic, statistical results of turbulence intensity at different altitudes are summarized.

[0013] A high-resolution modeling system for low-altitude turbulence intensity fields, used to perform the high-resolution modeling method for low-altitude turbulence intensity fields as described above, includes: The data preprocessing and alignment module is used to perform data preprocessing and spatiotemporal alignment, acquire detection data from ground automatic weather stations, wind profiler radar, and S-band and X-band Doppler weather radar, and perform time synchronization, coordinate unification, data cleaning, missing value imputation and outlier removal on multi-source heterogeneous detection data to achieve spatiotemporal alignment; The multi-source data three-level fusion module is used for the three-level fusion of multi-source data based on the turbulence dissipation rate. Based on the data preprocessed in step S1, the turbulence dissipation rate ε is inverted using the structure function method and the vertical profile of turbulence kinetic energy (TKE) is calculated using wind profile radar data. The spectral width component caused by pure turbulence is extracted using Doppler weather radar spectral width data and a statistical relationship model between it and TKE is established. The near-surface gradient Richardson number is calculated using ground automatic weather station data and a statistical calibration relationship between it and near-surface TKE is constructed. The prior background field is obtained through three-level fusion. The variational assimilation modeling and solution module is used for variational assimilation optimization modeling and solution. Based on the various key parameters inverted in step S2 and the obtained prior background field, it constructs a variational assimilation objective function with turbulent kinetic energy TKE as the core state variable. By introducing background field constraints, observation operators, and various physical constraints, the optimal three-dimensional turbulence intensity analysis field is solved using the L-BFGS-B optimization algorithm. The three-dimensional analysis field meets the grid specifications of 100 meters horizontal resolution, 50 meters vertical resolution, and 0-2 kilometers vertical range. The results output visualization module is used for results output and visualization. It outputs a three-dimensional turbulence intensity analysis field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters. The three-dimensional turbulence intensity field is overlaid with geographic information data to generate isosurface maps, horizontal slice maps and vertical profile maps, so as to realize the three-dimensional spatial visualization of low-altitude turbulence intensity.

[0014] A computer device, characterized in that it includes a memory and one or more processors, wherein the memory stores computer code, and when the computer code is executed by the one or more processors, causes the one or more processors to perform the method as described above.

[0015] A computer-readable storage medium, characterized in that the computer-readable storage medium stores computer code, which, when executed, is performed as described above.

[0016] Compared with the prior art, the beneficial effects of the present invention are: By combining a core architecture of "three-level fusion" and "variable assimilation," this method successfully solves the challenge of heterogeneous fusion and high-precision reconstruction of multi-source data in low-altitude turbulence monitoring. Its main advantages are: First, it innovatively proposes a three-level fusion strategy of "vertical reference - horizontal constraint - near-ground calibration," achieving effective complementarity and deep fusion of the advantages of heterogeneous data from wind profiler radar, Doppler radar, and ground stations. Second, based on optimal modeling using variational assimilation theory, it ensures that the final generated three-dimensional turbulence field is mathematically optimal and physically continuous and reasonable. Its high-resolution output of 100 meters horizontally and 50 meters vertically can accurately capture turbulent vortices that threaten low-altitude flight safety. Finally, while ensuring high accuracy, this method also possesses excellent engineering practicality, high computational efficiency, and good data integrity, providing real-time and reliable meteorological environment products for route planning and flight safety of UAVs and urban air traffic. Attached Figure Description

[0017] Figure 1 This is an overall flowchart of a high-resolution modeling method for low-altitude turbulence intensity fields according to the present invention; Figure 2 This is a spatial distribution map of the multi-source detection sites of this invention; Figure 3 This is a flowchart illustrating the principle of the multi-source data three-level fusion algorithm of the present invention; Figure 4 This is a three-dimensional spatial distribution map of turbulence intensity calculated by the present invention; Figure 5 This is a horizontal distribution diagram of turbulence intensity at different height layers calculated by the present invention; Figure 6 This is a structural diagram of the high-resolution modeling system for low-altitude turbulence intensity field of the present invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0019] Those skilled in the art will understand that, unless specifically stated otherwise, the singular forms “a,” “an,” “the,” and “the” used herein may also include the plural forms. It should be further understood that the term “comprising” as used in this specification means the presence of the stated features, integers, steps, operations, elements, and / or components, but does not exclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and / or groups thereof.

[0020] This invention discloses a method for modeling high-resolution three-dimensional turbulence intensity fields at low altitudes, aiming to solve the problem of constructing high-precision, high-resolution three-dimensional turbulence fields from multi-source meteorological observation data due to strong heterogeneity and spatiotemporal resolution mismatch. The method includes: preprocessing and spatiotemporally aligning multi-source data such as Doppler radar, wind profiler radar, and automatic weather stations; constructing a vertical benchmark based on turbulence dissipation rate, inverting turbulent kinetic energy using radar spectral width, and achieving three-level fusion of multi-source data through adaptive weighted fusion to form a priori background field; constructing an objective function using variational assimilation theory, combining various physical constraints, and solving for the optimal three-dimensional turbulence intensity analysis field using the L-BFGS-B algorithm; finally outputting a three-dimensional turbulence field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters, and visualizing it using geographic information. This invention integrates the advantages of multi-source data, possessing high precision, high integrity, and strong practicality, and can provide reliable meteorological environment product support for low-altitude flight safety.

[0021] The following is an illustration through specific examples: First Embodiment like Figure 1 As shown, this embodiment provides a high-resolution modeling method for low-altitude turbulence intensity fields, characterized by the following steps: S1: Perform data preprocessing and spatiotemporal alignment, acquiring detection data from ground-based automatic weather stations, wind profiler radar, and S-band and X-band Doppler weather radar. Perform time synchronization, coordinate unification, data cleaning, missing value imputation, and outlier removal on the multi-source heterogeneous detection data to achieve spatiotemporal alignment. For example... Figure 2 As shown, taking the detection stations in the vicinity of Nanjing as an example, the distribution of detection source stations is as follows: Figure 2 As shown.

[0022] In this embodiment, step S1 specifically includes: Step S11: Doppler radar data processing: Valid echoes are screened based on the signal-to-noise ratio threshold, and the missing values ​​of radial velocity and velocity spectrum width are filled by a two-dimensional interpolation method based on range and azimuth angle. Abnormal radial velocity values ​​that exceed the physical range are identified and corrected using the 3σ criterion. Step S12: Data processing of ground automatic weather stations: unify near-surface temperature, wind speed, air pressure and other data into WGS84 coordinate system and UTC time format, use inverse distance weighted interpolation to spatially interpolate missing values, and use the 3σ criterion to identify and correct outliers. Step S13: Wind profiler radar data processing: Wavelet denoising method is used to smooth the wind speed and wind direction time series. Based on the assumption of vertical continuity of atmospheric motion, linear interpolation and nonlinear interpolation methods are used to fill in the missing data at different height levels caused by weak echoes. Step S14: Spatiotemporal unification: Align all observation data in time to a 6-10 minute update frequency, and project and interpolate the radar polar coordinate data to the target Cartesian grid system in space.

[0023] S2: Three-level fusion of multi-source data based on turbulent dissipation rate. Based on the data preprocessed in step S1, the turbulent dissipation rate ε is inverted using the structure function method and the vertical profile of turbulent kinetic energy (TKE) is calculated using wind profiler radar data. The spectral width component caused by pure turbulence is extracted using Doppler weather radar spectral width data and a statistical relationship model between it and TKE is established. The near-surface gradient Richardson number is calculated using data from ground automatic weather stations and a statistical calibration relationship between it and near-surface TKE is constructed. The prior background field is obtained through three-level fusion.

[0024] In this embodiment, step S2 specifically includes: Step S21: Vertical baseline construction: Using the high-resolution vertical wind profile provided by the wind profiler radar, the turbulent dissipation rate ε at each height layer is calculated using the second-order structure function method; combined with the shear generation term P calculated from the vertical wind speed gradient, and the buoyancy term B determined based on the imaginary potential temperature gradient and the gradient Richardson number, the vertical turbulent kinetic energy (TKE) profile is obtained by inversion through the turbulent energy balance equation: Where z is the target height, This is the near-ground reference height. Data provided by ground stations; the vertical resolution of this TKE profile is 50m, but it can only reflect single-point features in the horizontal direction, and needs to be extended to a two-dimensional plane through a second-level fusion. The shear generation term P represents the shear of the mean wind field, that is, the transfer of energy from wind speed or direction to turbulence through dynamic instability as wind speed or direction changes with height. This is the primary source of TKE (Total Kinematic Energy). For a given target height z, its calculation formula is: in, It is the vertical gradient of the horizontal wind speed *u* at height *z*, i.e., the vertical wind shear. It is the momentum turbulence exchange coefficient; the buoyancy term B characterizes the promoting or inhibiting effect of thermal instability or stable stratification on turbulence development. When the atmosphere is unstable, B > 0, and buoyancy provides energy to the turbulence; when the atmosphere is stable, B < 0, and buoyancy consumes turbulence energy. It is determined by the vertical thermal structure of the atmosphere, i.e., the stratification stability, and its calculation formula is: Where g is the acceleration due to gravity. This is the potential temperature after taking into account air humidity, also known as the virtual potential temperature, which more accurately reflects the buoyancy effect of the atmosphere. It is the vertical gradient of the virtual potential temperature at height z. It is the turbulent heat exchange coefficient; Step S22: Horizontal spatial constraint correction: Perform multi-factor calibration on the original Doppler radar spectral width, eliminate instrument noise, wind shear contribution and precipitation interference, and extract the spectral width component caused by pure turbulence. Based on the power function statistical relationship, the spectral width is converted into a TKE value: Where k is the scaling factor, with a value of 0.85, and m is the exponent, with a value of 1.2; Step S23: Adaptive Weight Fusion: The TKE results retrieved from radar are matched to the vertical reference grid using Kriging interpolation. Fusion is achieved through an adaptive weighting function. The fusion formula is as follows: Where x is the horizontal coordinate of the target Cartesian grid, and z is the height; For the turbulent kinetic energy after fusion, For the vertical reference profile based on wind profiler radar, Turbulent kinetic energy obtained by Doppler radar inversion and interpolation; weighting coefficients Wind profiler radar data quality factor Doppler radar range attenuation factor Joint decision; Step S24: Near-Earth Boundary Calibration: Calculate the gradient Richardson number using data from ground-based automatic weather stations. The calculation formula is as follows: in, denoted by potential temperature, g is the acceleration due to gravity, z is the altitude, and u and v are the components of the horizontal wind in the east-west and north-south directions, respectively. Construct the statistical calibration formula: in, For the gradient Richardson number, The average wind speed at a height of 10m is given. The calculation result of the above formula is corrected for deviation from the near-surface TKE obtained by the second-level fusion, i.e., <200m, with a correction amount of . The 20m height layer is the near-ground height layer. The correction amount is transferred to the 200m height layer through vertical linear interpolation, and finally the three-dimensional turbulence intensity field is obtained.

[0025] Furthermore, in step S23, the formula for calculating the adaptive weight function is: in, This refers to the quality factor of wind profiler radar data. Doppler radar range attenuation factor, wind profiler radar data quality factor The signal-to-noise ratio (SNR) of the radar echo at altitude z is the primary factor determining the data quality. A higher SNR indicates more reliable data quality, and its weight should be increased accordingly. The calculation formula is as follows: in, Here is the measured signal-to-noise ratio at height z. and The effective threshold for signal-to-noise ratio is set based on radar performance. Doppler radar range attenuation factor This characterizes the attenuation of radar detection capability with distance, and is related to the spatial distance R(x, z) from the target grid point to the radar station. The greater the distance, the lower the radar detection accuracy and resolution, and the lower the data weight should be. The calculation formula is as follows: Where R(x,z) is the straight-line distance from the grid point to the radar station. γ is the effective maximum detection range of the radar, and γ is the attenuation exponent, which is usually taken as 1.5 to 2.5 and is used to control the rate at which the weight decays with distance.

[0026] S3: Variational Assimilation Optimization Modeling and Solution. Based on the key parameters retrieved in step S2 and the obtained prior background field, a variational assimilation objective function with turbulent kinetic energy (TKE) as the core state variable is constructed. By introducing background field constraints, observation operators, and various physical constraints, the optimal three-dimensional turbulence intensity analysis field is solved using the L-BFGS-B optimization algorithm. The principle diagram of its fusion modeling is shown below. Figure 3 As shown, the three-dimensional analysis field meets the grid specifications of 100 meters horizontal resolution, 50 meters vertical resolution, and 0-2 kilometers vertical range.

[0027] In this embodiment, step S3 specifically includes: Step S31: Cost function construction: Construct a cost function that includes background and observation terms to constrain the deviation between the analysis field and the prior background field, and to measure the fit between the analysis field and the observation data; Step S32: Definition of Observation Operators: Define observation operators for Doppler radar, wind profiler radar, and ground observation stations respectively, to establish the physical relationship between the state variable TKE and the original observation values ​​of various sensors, and map the gridded analysis field to the observation space of the k-th type of sensor: Doppler radar observation operator: in, To analyze the turbulent kinetic energy in the field, a and b are model parameters obtained based on statistical analysis of historical data. This power function form can effectively characterize the nonlinear relationship between turbulent energy and radar spectral width observations. Wind profiler radar observation operator: Where c is the proportionality coefficient. Let z be the turbulent mixing length at height z. This operator reflects the physical nature that the greater the turbulent kinetic energy, the stronger the dissipation, and takes into account the variation of the mixing length with height. Ground observation operator: This operator directly extracts the analysis field. The turbulent kinetic energy values ​​at the bottom layer of the grid, corresponding to the near-ground height, are compared with the measured turbulence estimates from the ground station. Step S33: Physical Constraint Settings: Set Non-negativity Constraints Smoothness constraints, boundary constraints, and resolution constraints are used to ensure that the solution conforms to physical laws; Nonnegativity constraint: Turbulent kinetic energy is the sum of squares of fluctuating velocities, and its value must be greater than or equal to zero; Smoothness constraint: The real turbulent field is usually continuously changing in space and will not have an infinitely large gradient. This is achieved through the background field error covariance matrix. Boundary constraints: Based on the fundamental laws of turbulent motion, constraints are imposed on the numerical range and variation trend of ε and TKE; Resolution constraint: Ensure that the model output meets the required spatiotemporal resolution; Step S34: Optimization Solution: The cost function is minimized using the L-BFGS-B algorithm. The analysis field I is iteratively updated using gradient descent. The process terminates when the gradient norm is less than the set tolerance or when the maximum number of iterations is reached. Specific steps are as follows: S341: Initialization: Given an initial guess I0, usually the background field Ib ; S342: Iterative Update: In the nth iteration, calculate the gradient of the cost function. : The L-BFGS-B algorithm is used to generate the search direction d. n ; S343: Line search: along d n Perform a one-dimensional search in the direction to determine the step size. Ensure that the updated solution satisfies I. i ≥0; S344: Convergence judgment: When the gradient norm is less than the set tolerance or the maximum number of iterations is reached, the iteration terminates and the analysis field I is output; Step S35: Model parameter setting and initialization settings: Grid system: horizontal resolution of 100 meters, vertically divided into 40 layers from ground to 2000 meters at 50-meter intervals, horizontal range covers the smallest bounding rectangle of all ground stations, i.e. longitude 118.4°E–119.0°E, latitude 31.8°N–32.6°N; Background field initialization: The climate-averaged turbulent kinetic energy or preliminary interpolation results are used as the initial estimate; Background error covariance matrix B: Based on the Gaussian model, the horizontal correlation length is set to 2 kilometers and the vertical correlation length to 200 meters. The standard deviation is determined based on climate statistics and observation data. The observation error covariance matrix R: Assuming that each observation is independent, it is a diagonal matrix, and the diagonal elements are set according to the sensor accuracy; for example, the ground station wind speed error variance is 0.25 square meters per square second, and the radar spectral width error variance is negatively correlated with the signal-to-noise ratio.

[0028] Numerical Implementation: Based on Python 3.8 or later, relying on the core functions of NumPy and SciPy libraries, the L-BFGS-B algorithm is used for constrained optimization. To address the large grid size, recursive filters or spectral methods are used to quickly calculate the inverse of the background error covariance. Efficient spherical to Cartesian coordinate transformation and interpolation model validity verification are performed on radar data.

[0029] Furthermore, in step S31, the specific form of the cost function constructed based on variational assimilation theory is as follows: in, As a background term, the project's constraint analysis field Do not deviate too much from a priori background field. The background field can be obtained based on climate mean, low-resolution numerical model forecasts or simple interpolation results. It introduces prior knowledge about the spatial structure of the turbulent field. B is the background field error covariance matrix, which adopts a Gaussian correlation function. Its diagonal elements represent the error variance of the background TKE value, while the off-diagonal elements characterize the correlation of background errors at different spatial points, reflecting the spatial continuity and smoothness scale of the TKE field. The observation term is used to measure the relationship between the analysis field and all K types of observation data. The degree of fit; It is an observation operator; The three-dimensional turbulence intensity analysis field to be solved is the state variable, and the variables on its mesh are turbulent kinetic energy or turbulent dissipation rate. B represents the prior background field obtained through three-level fusion, where B is the background error covariance matrix, constructed using a Gaussian correlation function, and its diagonal elements represent the background field. The error variance at different spatial locations, and the off-diagonal elements quantify the correlation between background errors at different spatial points, are represented by a Gaussian correlation function model. The function's form is as follows: Where C is the covariance between two points in space. It is the background error variance. and These are the vertical and horizontal distances between two points, respectively. and These are the correlation scales in the vertical and horizontal directions, respectively, used to control the smoothness of the analysis field; : The kth type of observation data; The observation operator for the k-th type of observation data is a physical model responsible for mapping the state variable I from the model grid space to the observation space, thus allowing it to be compared with the actual observation values. Compare; The observation error covariance matrix of the k-th type of observation data is usually assumed to be a diagonal matrix. The elements on its diagonal represent the uncertainty of various observation instruments, i.e., the error variance, which determines the relative weight of observation data from different sources and with different precision in the cost function. In step S33, the boundary constraints specifically include: Numerical range constraint: Based on the statistical measurements of low-altitude turbulence, the reasonable range for ε is 10. −6 ≤ε≤10 −2 m 2 / s 3The reasonable range for TKE is 0.1 ≤ TKE ≤ 10 m 2 / s 2 If the calculation exceeds this range, it will be automatically truncated to the boundary value. Vertical variation constraint: The attenuation rate of near-ground TKE with altitude must not exceed 0.002 m. 2 / ( s 2 The absolute value of the vertical gradient of the TKE at mid-to-high altitudes must not exceed 0.005 m. 2 / ( s 2 m), to avoid non-physical vertical jumps; Energy balance constraint: For any grid node, P+B≥0.5ε must be satisfied, and the generation of turbulent energy must not be much less than the dissipation, otherwise the energy conservation law will be violated; In step S33, the resolution constraint specifically includes: Horizontal resolution constraint: The difference between ε and TKE between any two adjacent horizontal grid nodes shall not exceed 20% of the mean of the region. If it exceeds this, Gaussian smoothing shall be used for correction with a smoothing radius of 100m. Vertical resolution constraint: The ε and TKE values ​​of each 50m height layer need to be verified by linear interpolation based on the data of the upper and lower layers. If the deviation exceeds 10%, the coupling weight is recalculated to ensure the continuity of the vertical layering.

[0030] S4: Results Output and Visualization. Outputs a three-dimensional turbulence intensity analysis field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters. This three-dimensional turbulence intensity field is overlaid with geographic information data to generate isosurface maps, horizontal slice maps, and vertical profile maps, realizing a three-dimensional spatial visualization of low-altitude turbulence intensity.

[0031] In this embodiment, step S4 specifically includes: Step S41: Output of 3D turbulence intensity analysis field: The mesh system corresponding to the output 3D turbulence intensity analysis field satisfies: Horizontal resolution is The value is 100 meters; Vertical direction within the height range Internal settings There are 1 vertical layers, each with a thickness of [missing information]. ; Step S42: Result visualization processing: Extract horizontal slice data from the three-dimensional turbulence intensity field calculated in step S3 at 50-meter vertical intervals; generate contour lines in the horizontal direction at a 100-meter grid resolution; overlay the above gridded data with the digital elevation model and administrative division vector data to generate a three-dimensional turbulence intensity spatial distribution map with geographic reference, a horizontal distribution map of a specified height layer, and a typical profile vertical structure map; Step S43: Result Characterization and Validation: The vertical stratification of turbulence intensity calculated by the model is shown in... Figure 4 The turbulence intensity calculated by the model is shown in a three-dimensional spatial distribution map. The turbulence intensity is strong near the ground and weak at high altitudes, exhibiting a vertical stratified distribution. The near-surface layer is affected by the combined effects of ground friction and thermal convection, resulting in significant airflow pulsation and ample turbulent energy supply. The upper atmosphere, far from the ground, relies mainly on wind shear to maintain weaker turbulence. The model results accurately capture this physical process, verifying the calculation logic based on the vertical energy balance equation and proving that the model effectively simulates the dissipation of turbulent energy with altitude. To more clearly quantify the stratified variation of turbulence intensity with altitude calculated by the model, and to verify its consistency with atmospheric boundary layer theory and model construction logic, the statistical results and correlation analysis of turbulence intensity at different altitudes are summarized in Table 1, and in... Figure 5 The figures show the horizontal distribution of turbulence intensity at heights of 100m, 500m, 1000m, and 2000m.

[0032] Table 1 Distribution of turbulence intensity in the upper layer Second Embodiment like Figure 6 As shown, this embodiment provides a high-resolution modeling system for low-altitude turbulence intensity fields for performing the high-resolution modeling method for low-altitude turbulence intensity fields as described in the first embodiment, comprising: The data preprocessing and alignment module is used to perform data preprocessing and spatiotemporal alignment, acquire detection data from ground automatic weather stations, wind profiler radar, and S-band and X-band Doppler weather radar, and perform time synchronization, coordinate unification, data cleaning, missing value imputation and outlier removal on multi-source heterogeneous detection data to achieve spatiotemporal alignment; The multi-source data three-level fusion module is used for the three-level fusion of multi-source data based on the turbulence dissipation rate. Based on the data preprocessed in step S1, the turbulence dissipation rate ε is inverted using the structure function method and the vertical profile of turbulence kinetic energy (TKE) is calculated using wind profile radar data. The spectral width component caused by pure turbulence is extracted using Doppler weather radar spectral width data and a statistical relationship model between it and TKE is established. The near-surface gradient Richardson number is calculated using ground automatic weather station data and a statistical calibration relationship between it and near-surface TKE is constructed. The prior background field is obtained through three-level fusion. The variational assimilation modeling and solution module is used for variational assimilation optimization modeling and solution. Based on the various key parameters inverted in step S2 and the obtained prior background field, it constructs a variational assimilation objective function with turbulent kinetic energy TKE as the core state variable. By introducing background field constraints, observation operators, and various physical constraints, the optimal three-dimensional turbulence intensity analysis field is solved using the L-BFGS-B optimization algorithm. The three-dimensional analysis field meets the grid specifications of 100 meters horizontal resolution, 50 meters vertical resolution, and 0-2 kilometers vertical range. The results output visualization module is used for results output and visualization. It outputs a three-dimensional turbulence intensity analysis field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters. The three-dimensional turbulence intensity field is overlaid with geographic information data to generate isosurface maps, horizontal slice maps and vertical profile maps, so as to realize the three-dimensional spatial visualization of low-altitude turbulence intensity.

[0033] A computer-readable storage medium stores computer code that, when executed, performs the methods described above. Those skilled in the art will understand that all or part of the steps in the various methods of the above embodiments can be implemented by a program instructing related hardware. This program can be stored in a computer-readable storage medium, which may include: read-only memory (ROM), random access memory (RAM), a magnetic disk, or an optical disk, etc.

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

[0035] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0036] It should be noted that the above embodiments can be freely combined as needed. The above description is only a preferred embodiment of the present invention. It should be pointed out that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A high-resolution modeling method for low-altitude turbulence intensity fields, characterized in that, Includes the following steps: S1: Perform data preprocessing and spatiotemporal alignment, acquire detection data from ground automatic weather stations, wind profiler radar, S-band and X-band Doppler weather radar, and perform time synchronization, coordinate unification, data cleaning, missing value imputation and outlier removal on multi-source heterogeneous detection data to achieve spatiotemporal alignment; S2: Three-level fusion of multi-source data based on turbulent dissipation rate. Based on the data preprocessed in step S1, the turbulent dissipation rate ε is inverted using the structure function method and the vertical profile of turbulent kinetic energy (TKE) is calculated using wind profiler radar data. The spectral width component caused by pure turbulence is extracted using Doppler weather radar spectral width data and a statistical relationship model between it and TKE is established. The near-surface gradient Richardson number is calculated using data from ground automatic weather stations and a statistical calibration relationship between it and near-surface TKE is constructed. The prior background field is obtained through three-level fusion. S3: Variational assimilation optimization modeling and solution. Based on the various key parameters inverted in step S2 and the obtained prior background field, a variational assimilation objective function with turbulent kinetic energy TKE as the core state variable is constructed. By introducing background field constraints, observation operators, and various physical constraints, the optimal three-dimensional turbulence intensity analysis field is solved using the L-BFGS-B optimization algorithm. The three-dimensional analysis field meets the grid specifications of 100 meters horizontal resolution, 50 meters vertical resolution, and 0-2 kilometers vertical range. S4: Results Output and Visualization. Outputs a three-dimensional turbulence intensity analysis field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters. This three-dimensional turbulence intensity field is overlaid with geographic information data to generate isosurface maps, horizontal slice maps, and vertical profile maps, realizing a three-dimensional spatial visualization of low-altitude turbulence intensity.

2. The high-resolution modeling method for low-altitude turbulence intensity field according to claim 1, characterized in that, Step S1 is as follows: Step S11: Doppler radar data processing: Valid echoes are screened based on the signal-to-noise ratio threshold, and the missing values ​​of radial velocity and velocity spectrum width are filled by a two-dimensional interpolation method based on range and azimuth angle. Abnormal radial velocity values ​​that exceed the physical range are identified and corrected using the 3σ criterion. Step S12: Data processing of ground automatic weather stations: unify near-surface temperature, wind speed, air pressure and other data into WGS84 coordinate system and UTC time format, use inverse distance weighted interpolation to spatially interpolate missing values, and use the 3σ criterion to identify and correct outliers. Step S13: Wind profiler radar data processing: Wavelet denoising method is used to smooth the wind speed and wind direction time series. Based on the assumption of vertical continuity of atmospheric motion, linear interpolation and nonlinear interpolation methods are used to fill in the missing data at different height levels caused by weak echoes. Step S14: Spatiotemporal unification: Align all observation data in time to a 6-10 minute update frequency, and project and interpolate the radar polar coordinate data to the target Cartesian grid system in space.

3. The high-resolution modeling method for low-altitude turbulence intensity field according to claim 1, characterized in that, Step S2 is as follows: Step S21: Vertical baseline construction: Using the high-resolution vertical wind profile provided by the wind profiler radar, the turbulent dissipation rate ε at each height layer is calculated using the second-order structure function method; combined with the shear generation term P calculated from the vertical wind speed gradient, and the buoyancy term B determined based on the imaginary potential temperature gradient and the gradient Richardson number, the vertical turbulent kinetic energy TKE profile is obtained by inversion through the turbulent energy balance equation: Where z is the target height, This is the near-ground reference height. Data provided by ground stations; the vertical resolution of this TKE profile is 50m, but it can only reflect single-point features in the horizontal direction, and needs to be extended to a two-dimensional plane through a second-level fusion. The shear generation term P represents the shear of the mean wind field, that is, the transfer of energy from wind speed or direction to turbulence through dynamic instability as wind speed or direction changes with height. This is the primary source of TKE (Total Kinematic Energy). For a given target height z, its calculation formula is: in, It is the vertical gradient of the horizontal wind speed *u* at height *z*, i.e., the vertical wind shear. It is the momentum turbulence exchange coefficient; the buoyancy term B characterizes the promoting or inhibiting effect of thermal instability or stable stratification on turbulence development. When the atmosphere is unstable, B > 0, and buoyancy provides energy to the turbulence; when the atmosphere is stable, B < 0, and buoyancy consumes turbulence energy. It is determined by the vertical thermal structure of the atmosphere, i.e., the stratification stability, and its calculation formula is: Where g is the acceleration due to gravity. This is the potential temperature after taking into account air humidity, also known as the virtual potential temperature, which more accurately reflects the buoyancy effect of the atmosphere. It is the vertical gradient of the virtual potential temperature at height z. It is the turbulent heat exchange coefficient; Step S22: Horizontal spatial constraint correction: Perform multi-factor calibration on the original Doppler radar spectral width, eliminate instrument noise, wind shear contribution and precipitation interference, and extract the spectral width component caused by pure turbulence. Based on the power function statistical relationship, the spectral width is converted into a TKE value: Where k is the scaling factor, with a value of 0.85, and m is the exponent, with a value of 1.2; Step S23: Adaptive Weight Fusion: The TKE results retrieved from radar are matched to the vertical reference grid using Kriging interpolation. Fusion is achieved through an adaptive weighting function. The fusion formula is as follows: Where x is the horizontal coordinate of the target Cartesian grid, and z is the height; For the turbulent kinetic energy after fusion, For the vertical reference profile based on wind profiler radar, Turbulent kinetic energy obtained by Doppler radar inversion and interpolation; weighting coefficients Wind profiler radar data quality factor Doppler radar range attenuation factor Joint decision; Step S24: Near-Earth Boundary Calibration: Calculate the gradient Richardson number using data from ground-based automatic weather stations. The calculation formula is as follows: in, denoted by potential temperature, g is the acceleration due to gravity, z is the altitude, and u and v are the components of the horizontal wind in the east-west and north-south directions, respectively. Construct the statistical calibration formula: in, For the gradient Richardson number, The average wind speed at a height of 10m is given. The calculation result of the above formula is corrected for deviation from the near-surface TKE obtained by the second-level fusion, i.e., <200m, with a correction amount of . The 20m height layer is the near-ground height layer. The correction amount is transferred to the 200m height layer through vertical linear interpolation, and finally the three-dimensional turbulence intensity field is obtained.

4. The high-resolution modeling method for low-altitude turbulence intensity field according to claim 3, characterized in that, In step S23, the formula for calculating the adaptive weight function is: in, This refers to the quality factor of wind profiler radar data. Doppler radar range attenuation factor, wind profiler radar data quality factor The signal-to-noise ratio (SNR) of the radar echo at altitude z is the primary factor determining the data quality. A higher SNR indicates more reliable data quality, and its weight should be increased accordingly. The calculation formula is as follows: in, Here is the measured signal-to-noise ratio at height z. and The effective threshold for signal-to-noise ratio is set based on radar performance. Doppler radar range attenuation factor This characterizes the attenuation of radar detection capability with distance, and is related to the spatial distance R(x, z) from the target grid point to the radar station. The greater the distance, the lower the radar detection accuracy and resolution, and the lower the data weight should be. The calculation formula is as follows: Where R(x,z) is the straight-line distance from the grid point to the radar station. γ is the effective maximum detection range of the radar, and γ is the attenuation exponent, which is usually taken as 1.5 to 2.5 and is used to control the rate at which the weight decays with distance.

5. The high-resolution modeling method for low-altitude turbulence intensity field according to claim 1, characterized in that, Step S3 is as follows: Step S31: Cost function construction: Construct a cost function that includes background and observation terms to constrain the deviation between the analysis field and the prior background field, and to measure the fit between the analysis field and the observation data; Step S32: Definition of Observation Operators: Define observation operators for Doppler radar, wind profiler radar, and ground observation stations respectively, to establish the physical relationship between the state variable TKE and the original observation values ​​of various sensors, and map the gridded analysis field to the observation space of the k-th type of sensor: Doppler radar observation operator: in, To analyze the turbulent kinetic energy in the field, a and b are model parameters obtained based on statistical analysis of historical data. This power function form can effectively characterize the nonlinear relationship between turbulent energy and radar spectral width observations. Wind profiler radar observation operator: Where c is the proportionality coefficient. Let z be the turbulent mixing length at height z. This operator reflects the physical nature that the greater the turbulent kinetic energy, the stronger the dissipation, and takes into account the variation of the mixing length with height. Ground observation operator: This operator directly extracts the analysis field. The turbulent kinetic energy values ​​at the bottom layer of the grid, corresponding to the near-ground height, are compared with the measured turbulence estimates from the ground station. Step S33: Physical Constraint Settings: Set Non-negativity Constraints Smoothness constraints, boundary constraints, and resolution constraints are used to ensure that the solution conforms to physical laws; Nonnegativity constraint: Turbulent kinetic energy is the sum of squares of fluctuating velocities, and its value must be greater than or equal to zero; Smoothness constraint: The real turbulent field is usually continuously changing in space and will not have an infinitely large gradient. This is achieved through the background field error covariance matrix. Boundary constraints: Based on the fundamental laws of turbulent motion, constraints are imposed on the numerical range and variation trend of ε and TKE; Resolution constraint: Ensure that the model output meets the required spatiotemporal resolution; Step S34: Optimization Solution: The cost function is minimized using the L-BFGS-B algorithm. The analysis field I is iteratively updated using gradient descent. The process terminates when the gradient norm is less than the set tolerance or when the maximum number of iterations is reached. Specific steps are as follows: S341: Initialization: Given an initial guess I0, usually the background field I b ; S342: Iterative Update: In the nth iteration, calculate the gradient of the cost function. : The L-BFGS-B algorithm is used to generate the search direction d. n ; S343: Line search: along d n Perform a one-dimensional search to determine the direction. Ensure that the updated solution satisfies I. i ≥0; S344: Convergence judgment: When the gradient norm is less than the set tolerance or the maximum number of iterations is reached, the iteration terminates and the analysis field I is output; Step S35: Model parameter setting and initialization settings: Grid system: horizontal resolution of 100 meters, vertically divided into 40 layers from ground to 2000 meters at 50-meter intervals, horizontal range covers the smallest bounding rectangle of all ground stations, i.e. longitude 118.4°E–119.0°E, latitude 31.8°N–32.6°N; Background field initialization: The climate-averaged turbulent kinetic energy or preliminary interpolation results are used as the initial estimate; Background error covariance matrix B: Based on the Gaussian model, the horizontal correlation length is set to 2 kilometers and the vertical correlation length to 200 meters. The standard deviation is determined based on climate statistics and observation data. Observation error covariance matrix R: Assuming that each observation is independent, it is a diagonal matrix, and the diagonal elements are set according to the sensor accuracy. Numerical Implementation: Based on Python 3.8 or later, relying on the core functions of NumPy and SciPy libraries, the L-BFGS-B algorithm is used for constrained optimization. To address the large grid size, recursive filters or spectral methods are used to quickly calculate the inverse of the background error covariance. Efficient spherical to Cartesian coordinate transformation and interpolation model validity verification are performed on radar data.

6. The high-resolution modeling method for low-altitude turbulence intensity field according to claim 5, characterized in that, In step S31, the specific form of the cost function based on variational assimilation theory is as follows: in, As a background term, the project's constraint analysis field Do not deviate too much from a priori background field. The background field can be obtained based on climate mean, low-resolution numerical model forecasts or simple interpolation results. It introduces prior knowledge about the spatial structure of the turbulent field. B is the background field error covariance matrix, which adopts a Gaussian correlation function. Its diagonal elements represent the error variance of the background TKE value, while the off-diagonal elements characterize the correlation of background errors at different spatial points, reflecting the spatial continuity and smoothness scale of the TKE field. The observation term is used to measure the relationship between the analysis field and all K types of observation data. The degree of fit; It is an observation operator; The three-dimensional turbulence intensity analysis field to be solved is the state variable, and the variables on its mesh are turbulent kinetic energy or turbulent dissipation rate. B represents the prior background field obtained through three-level fusion, where B is the background error covariance matrix, constructed using a Gaussian correlation function, and its diagonal elements represent the background field. The error variance at different spatial locations, and the off-diagonal elements quantify the correlation between background errors at different spatial points, are represented by a Gaussian correlation function model. The function's form is as follows: Where C is the covariance between two points in space. It is the background error variance. and These are the vertical and horizontal distances between two points, respectively. and These are the correlation scales in the vertical and horizontal directions, respectively, used to control the smoothness of the analysis field; : The kth type of observation data; The observation operator for the k-th type of observation data is a physical model responsible for mapping the state variable I from the model grid space to the observation space, thus allowing it to be compared with the actual observation values. Compare; The observation error covariance matrix of the k-th type of observation data is usually assumed to be a diagonal matrix. The elements on its diagonal represent the uncertainty of various observation instruments, i.e., the error variance, which determines the relative weight of observation data from different sources and with different precision in the cost function. In step S33, the boundary constraints specifically include: Numerical range constraint: Based on the statistical measurements of low-altitude turbulence, the reasonable range for ε is 10. −6 ≤ε≤10 −2 m 2 / s 3 The reasonable range for TKE is 0.1 ≤ TKE ≤ 10 m 2 / s 2 If the calculation exceeds this range, it will be automatically truncated to the boundary value. Vertical variation constraint: The attenuation rate of near-ground TKE with altitude must not exceed 0.002 m. 2 / ( s 2 The absolute value of the vertical gradient of the TKE at mid-to-high altitudes must not exceed 0.005 m. 2 / ( s 2 m), to avoid non-physical vertical jumps; Energy balance constraint: For any grid node, P+B≥0.5ε must be satisfied, and the generation of turbulent energy must not be much less than the dissipation, otherwise the energy conservation law will be violated; In step S33, the resolution constraint specifically includes: Horizontal resolution constraint: The difference between ε and TKE between any two adjacent horizontal grid nodes shall not exceed 20% of the mean of the region. If it exceeds this, Gaussian smoothing shall be used for correction with a smoothing radius of 100m. Vertical resolution constraint: The ε and TKE values ​​of each 50m height layer need to be verified by linear interpolation based on the data of the upper and lower layers. If the deviation exceeds 10%, the coupling weight is recalculated to ensure the continuity of the vertical layering.

7. The high-resolution modeling method for low-altitude turbulence intensity field according to claim 1, characterized in that, Step S4 is as follows: Step S41: Output of 3D turbulence intensity analysis field: The mesh system corresponding to the output 3D turbulence intensity analysis field satisfies: Horizontal resolution is The value is 100 meters; Vertical direction within the height range Internal settings There are 1 vertical layers, each with a thickness of [missing information]. ; Step S42: Result visualization processing: Extract horizontal slice data from the three-dimensional turbulence intensity field calculated in step S3 at 50-meter vertical intervals; In the horizontal direction, contour lines are generated at a grid resolution of 100 meters; the above gridded data are overlaid with digital elevation model and administrative division vector data to generate a three-dimensional spatial distribution map of turbulence intensity with geographic reference, a horizontal distribution map of a specified height layer, and a vertical structure map of a typical profile. Step S43: Results and Verification: The turbulence intensity calculated by the model exhibits a vertically stratified distribution, with strong turbulence near the ground and weak turbulence at higher altitudes. The near-surface layer is affected by both ground friction and thermal convection, resulting in significant airflow fluctuations and ample turbulent energy supply. The upper atmosphere, far from the ground, relies mainly on wind shear to maintain weaker turbulence. The model results accurately capture this physical process, verifying the calculation logic based on the vertical energy balance equation and proving that the model effectively simulates the dissipation of turbulent energy with altitude. To more clearly quantify the stratified variation of turbulence intensity with altitude calculated by the model, and to verify its consistency with atmospheric boundary layer theory and model construction logic, statistical results of turbulence intensity at different altitudes are summarized.

8. A high-resolution modeling system for low-altitude turbulence intensity fields for performing the high-resolution modeling method for low-altitude turbulence intensity fields as described in any one of claims 1-7, characterized in that, include: The data preprocessing and alignment module is used to perform data preprocessing and spatiotemporal alignment, acquire detection data from ground automatic weather stations, wind profiler radar, and S-band and X-band Doppler weather radar, and perform time synchronization, coordinate unification, data cleaning, missing value imputation and outlier removal on multi-source heterogeneous detection data to achieve spatiotemporal alignment; The multi-source data three-level fusion module is used for the three-level fusion of multi-source data based on the turbulence dissipation rate. Based on the data preprocessed in step S1, the turbulence dissipation rate ε is inverted using the structure function method and the vertical profile of turbulence kinetic energy (TKE) is calculated using wind profile radar data. The spectral width component caused by pure turbulence is extracted using Doppler weather radar spectral width data and a statistical relationship model between it and TKE is established. The near-surface gradient Richardson number is calculated using ground automatic weather station data and a statistical calibration relationship between it and near-surface TKE is constructed. The prior background field is obtained through three-level fusion. The variational assimilation modeling and solution module is used for variational assimilation optimization modeling and solution. Based on the various key parameters inverted in step S2 and the obtained prior background field, it constructs a variational assimilation objective function with turbulent kinetic energy TKE as the core state variable. By introducing background field constraints, observation operators, and various physical constraints, the optimal three-dimensional turbulence intensity analysis field is solved using the L-BFGS-B optimization algorithm. The three-dimensional analysis field meets the grid specifications of 100 meters horizontal resolution, 50 meters vertical resolution, and 0-2 kilometers vertical range. The results output visualization module is used for results output and visualization. It outputs a three-dimensional turbulence intensity analysis field with a horizontal resolution of 100 meters and a vertical resolution of 50 meters. The three-dimensional turbulence intensity field is overlaid with geographic information data to generate isosurface maps, horizontal slice maps and vertical profile maps, so as to realize the three-dimensional spatial visualization of low-altitude turbulence intensity.

9. A computer device, characterized in that, The device includes a memory and one or more processors, wherein the memory stores computer code that, when executed by the one or more processors, causes the one or more processors to perform the method as described in any one of claims 1 to 7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer code, and when the computer code is executed, the method as described in any one of claims 1 to 7 is performed.

Images