An inversion method, system, medium and program product for atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar

Abstract

The present invention discloses a method, system, medium and program product for retrieving atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar, which relates to the technical field of atmospheric remote sensing and boundary layer detection. By controlling the lidar to emit 6 beams at a set elevation angle and azimuth angle, collecting the multi-directional radial wind speed time series, and calculating the radial velocity variance of each beam on the basis of data quality control. Further, establish the relationship between it and the velocity variances and covariances in three directions, construct an observation equation set and retrieve the turbulence statistics, and combine coordinate rotation to obtain the turbulence statistics in the along-wind, cross-wind and vertical directions. Finally, calculate the turbulent kinetic energy and output the vertical profiles of the turbulence statistics at each height layer. The present invention can realize continuous monitoring and high-precision characterization of the atmospheric boundary layer turbulence statistics, and meet the requirements of the business system for the boundary layer turbulence inversion products in terms of stability, accuracy and deployability.

Description

Technical Field

[0001] The present invention belongs to the technical field of atmospheric remote sensing and boundary layer detection, and relates to lidar wind field detection and turbulence parameter inversion technology. Specifically, it is an atmospheric boundary layer turbulence statistic inversion method, system, medium and program product based on the DBS scanning mode of a 6-beam ground-based optical quantum wind measurement lidar. Background Art

[0002] The atmospheric boundary layer is a key area for the exchange of matter and energy between the earth's surface and the free atmosphere. Its turbulence structure and evolution have important impacts on meteorological forecasting, air quality assessment, wind energy resource development, and aviation safety. Therefore, accurately and continuously obtaining the vertical distribution characteristics of turbulence statistics (such as velocity variance, momentum flux, turbulent kinetic energy, etc.) in the atmospheric boundary layer is of great theoretical and practical significance for understanding boundary layer dynamics, developing boundary layer parameterization schemes, and establishing high-precision numerical models.

[0003] Traditional atmospheric boundary layer observation methods mainly include sounding, boundary layer towers, radar wind profilers, aircraft platforms, and satellite remote sensing. Among them, sounding and tower observations can obtain the vertical distribution information of basic meteorological elements such as wind speed, temperature, and humidity. However, limited by low time resolution and limited spatial coverage, it is difficult to meet the need for continuous inversion of turbulence statistics. Although aircraft observations have high accuracy, due to problems such as high cost and low frequency, they are only applicable to short-term and local area research. Satellite remote sensing plays a role in large-scale meteorological monitoring, but its spatial resolution and penetration ability for the low-altitude boundary layer are still limited.

[0004] In recent years, with the development of lidar technology, especially the wide application of optical quantum Doppler wind measurement lidar, it has become possible to obtain a three-dimensional wind field with high spatio-temporal resolution. This technology emits laser beams and receives the scattered signals of the laser by atmospheric aerosols, and uses the Doppler frequency shift to invert the wind speed information. In actual observations, Doppler lidar can adopt different scanning modes to achieve wind field inversion and turbulence parameter estimation. Common modes include Velocity Azimuth Display (VAD), Range Height Indicator (RHI), and Doppler Beam Swinging (DBS), etc. The VAD and RHI methods can be used to invert the mean wind field and vertical velocity structure, but they have limited accuracy in turbulence statistic inversion. The DBS method obtains the radial velocities in multiple directions by setting multiple laser beams with fixed azimuths and elevation angles, and thus inverses parameters such as velocity variance, momentum flux, and turbulent kinetic energy.

[0005] However, the traditional DBS method mostly uses a 4-beam configuration, and its inversion results are easily affected by the angle distribution between beams, wind field non-uniformity and atmospheric disturbances, resulting in insufficient calculation accuracy of turbulence statistics. In the prior art, Chinese patent CN109814131B proposes a laser radar turbulence parameter inversion method based on discontinuous conical scanning, which reduces the speed measurement error of radial wind speed by pausing the scan, and uses the Kolmogorov model combined with the velocity structure function to eliminate the noise term, but its scanning mode is still mainly based on a single fixed zenith angle, and it is difficult to obtain the three-dimensional turbulence structure of the entire boundary layer. Its discontinuous scanning strategy causes the laser beam to have a window area in time, which limits the ability to fully capture high-frequency turbulent disturbances and is not conducive to reconstructing the continuous vertical profile of turbulence statistics in the entire boundary layer. In addition, some studies have attempted to increase the number of beams to 5 or 6 to enhance the inversion capability, but there is still a lack of standardized guidance on scanning strategies and parameter settings (such as beam azimuth, elevation angle, time interval, and spatial resolution). Especially in the deployment of operational observations, there is no mature solution for how to balance scanning efficiency and inversion accuracy. In addition, existing studies have focused on products such as average wind speed and boundary layer height. High-precision inversion products for turbulence statistics have not been systematically established, and data quality control and result reliability assessment are still insufficient.

[0006] In summary, the existing atmospheric boundary layer turbulence observation methods have obvious deficiencies in terms of coverage height, spatiotemporal resolution, inversion accuracy and product systematization. Therefore, how to build a turbulence statistics inversion method and system with high accuracy, strong stability and suitable for business deployment based on a ground-based photon wind lidar system with high spatiotemporal resolution is a technical problem that needs to be solved urgently. Summary of the invention

[0007] 1. Purpose of the invention

[0008] In view of the above-mentioned defects and shortcomings of the prior art, the present invention aims to provide a method, system, medium and program product for inverting turbulence statistics in the atmospheric boundary layer based on the DBS scanning mode of a 6-beam ground-based photon wind measurement lidar. The method adopts 6-beam laser synchronous measurement and DBS scanning mode to obtain the spatiotemporal distribution of the atmospheric three-dimensional wind field in the boundary layer with high spatiotemporal resolution, and constructs a radial velocity variance matrix containing multi-beam observation directions to invert the turbulence statistics of the atmospheric boundary layer such as velocity variance, momentum flux, turbulent kinetic energy, etc., and introduces coordinate rotation to unify the results, thereby outputting the turbulence vertical profile with high spatiotemporal resolution, realizing continuous monitoring and high-precision characterization of the turbulence structure in the atmospheric boundary layer, and meeting the business system's requirements for the stability, accuracy and deployability of the boundary layer turbulence inversion products.

[0009] (II) Technical solution

[0010] To achieve the object of the present invention and solve its technical problems, the present invention adopts the following technical solutions:

[0011] The first object of the present invention is to provide a method for retrieving turbulence statistics in the atmospheric boundary layer based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar, for obtaining the vertical profiles of turbulence statistics such as velocity variance, momentum flux, and turbulent kinetic energy in the atmospheric boundary layer. When the retrieval method is implemented, it at least includes the following steps:

[0012] S100. Configuration of the 6-beam DBS scanning mode:

[0013] Control the ground-based optical quantum wind lidar to observe in the Doppler Beam Swinging (DBS) scanning mode, and set 6 detection beams, including 1 zenith beam vertically upward and 5 inclined beams distributed at a set elevation angle φ and equally spaced azimuth angles along the conical surface θ ;

[0014] S200. Acquisition and processing of radial wind speed data:

[0015] Synchronously obtain the time series data of the radial wind speed in the direction of each detection beam, record the corresponding azimuth angle θ and elevation angle φ parameters, and perform quality control on the acquired radial wind speed data, screening out data points with a signal-to-noise ratio lower than the preset threshold, and removing abnormal data points that are significantly deviated from the average value or show mutations;

[0016] S300. Calculate the radial velocity variance of each beam:

[0017] Perform statistical processing on the time series data of the radial wind speed in the direction of each beam collected v ri ( t ) to calculate the radial velocity variance corresponding to each beam , and construct an observation vector of the radial velocity variance accordingly , is the radial velocity variance of the 1st to 6th beams, i = 1, 2, …, 6;

[0018] S400. Establish a relationship between the radial velocity variance and turbulence statistics:

[0019] According to the elevation angle φ i and azimuth angle θ i corresponding to each beam, respectively establish the radial velocity variance of each beam and 6 turbulence statistics Relationship between:

[0020]

[0021] wherein, are the pulsating wind speed components in the horizontal east - west direction, horizontal north - south direction, and vertical direction respectively are the velocity variances in three directions respectively, are the covariances of the velocity components in three directions, namely the momentum flux terms;

[0022] S500. Construct the radial velocity variance observation equations and parameter matrix:

[0023] Based on the relationship between the radial velocity variances in 6 beam directions and the turbulence statistics, construct the linear observation equations of the radial velocity variance in matrix form , where S is the radial velocity variance observation vector, M is the 6×6 coefficient matrix calculated according to the elevation angles φ i and azimuth angles θ i of each beam, Σ is the vector of turbulence statistics to be solved and ;

[0024] S600. Inversion solution of turbulence statistics and coordinate rotation:

[0025] Based on the radial velocity variance observation equations , use the inverse matrix method or the least - squares fitting method to inversely solve the vector Σ of turbulence statistics, and perform coordinate rotation on the inversion result based on the mean wind direction to ensure that the velocity perturbation variance and covariance are expressed in the coordinate system aligned with the mean wind direction;

[0026] S700. Turbulent kinetic energy TKE calculation and result output:

[0027] Based on the inverted velocity variance, calculate the turbulent kinetic energy , and finally output the vertical profiles of turbulence statistics including velocity variance, momentum flux, and turbulent kinetic energy at each height layer to characterize the vertical structure of the atmospheric boundary layer turbulence.

[0028] The second object of the present invention is to provide an atmospheric boundary layer turbulence statistics inversion system for implementing the above - mentioned atmospheric boundary layer turbulence statistics inversion method based on the 6 - beam ground - based optical quantum wind lidar DBS scanning mode of the present invention. The system includes at least the following modules:

[0029] Lidar observation module, which controls the ground - based optical quantum wind lidar to emit 6 detection beams in the DBS scanning mode, including 1 vertical zenith beam and 5 inclined beams evenly distributed along the conical surface;

[0030] A data acquisition and processing module, which is used to synchronously acquire the radial wind speed time series data in the directions of each detection beam, record the corresponding azimuth and elevation angle parameters, and perform quality control on the radial wind speed data;

[0031] A statistic calculation and matrix construction module, which is used to calculate the radial velocity variance of each beam at each height layer, and establish a linear equation set between it and the turbulence statistics based on the beam geometric angle;

[0032] An inversion solution and coordinate rotation module, which is used to solve the turbulence statistic vector by using the inverse matrix method or the least square fitting method, and perform coordinate transformation based on the mean wind direction;

[0033] A result calculation and output module, which is used to calculate the turbulent kinetic energy based on the inversely obtained velocity variance, and output the vertical profiles of the turbulence statistics at each height layer.

[0034] The third invention object of the present invention is to provide a computer program product, including computer instructions, which are used to execute the above-mentioned method for inverting the turbulence statistics of the atmospheric boundary layer based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar of the present invention.

[0035] The fourth invention object of the present invention is to provide a computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the above-mentioned method for inverting the turbulence statistics of the atmospheric boundary layer based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar of the present invention is realized.

[0036] (III) Technical effects

[0037] Compared with the prior art, the method, system, medium and program product for inverting the turbulence statistics of the atmospheric boundary layer based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar of the present invention have the following beneficial and remarkable technical effects:

[0038] (1) By introducing the 6-beam DBS scanning mode (1 zenith beam and 5 inclined beams), the present invention can synchronously acquire the radial wind speed data of multiple height layers in the atmospheric boundary layer, realize the high spatio-temporal resolution observation of the turbulence statistics, and at the same time, by accurately constructing the mathematical relationship between the radial velocity variance and the turbulence statistics and optimizing the solution of the matrix equation, the systematic error caused by the underdetermined equation of the traditional 4-5 beams is avoided, and the inversion accuracy of the turbulence statistics such as the velocity variance and the momentum flux is significantly improved.

[0039] (2) The present invention constructs an analytical relationship model between radial velocity variance and six turbulence statistics, and uniformly solves them in matrix form, avoiding the uncertainty of empirical estimation item by item, and realizing the systematic inversion of velocity variance and momentum flux. In addition, the present invention adopts a strategy combining quality control and coordinate rotation to establish a stable mapping process from the original observation to the turbulence component in the main wind direction coordinate system, which improves the physical interpretability of the inversion results and is particularly suitable for boundary layer environments with drastic wind direction changes.

[0040] (3) The present invention is based on ground-based quantum wind laser radar technology, which has the advantages of small size, light weight, and easy deployment. It can be widely used in atmospheric boundary layer turbulence observation in different scenarios such as urban environments, mountainous terrains, and ocean areas. At the same time, the present invention can realize continuous vertical profile observation of atmospheric boundary layer turbulence statistics, provide key parameters for atmospheric diffusion models, numerical weather forecasts, and climate models, and has significant scientific value and application prospects. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 It is a flow chart of the atmospheric boundary layer turbulence statistics inversion method based on the 6-beam ground-based photon wind laser radar DBS scanning mode of the present invention;

[0042] Figure 2 It is a schematic diagram of the atmospheric boundary layer turbulence statistics inversion system of the present invention;

[0043] Figure 3 The time-height profiles of various turbulence statistics obtained by inverting the six-beam detection mode in Hefei area from 11-16 00:00 to 11-21 23:30 Beijing time in 2023, where: (a) u Directional velocity variance , (b) is v Directional velocity variance , (c) is the vertical velocity variance , (d) is the turbulent kinetic energy , pink and cyan are the boundary layer heights obtained from wind lidar and ERA5 reanalysis, respectively.

[0044] Description of reference numerals:

[0045] LiDAR observation module 100, data acquisition and processing module 200, statistics calculation and matrix construction module 300, inversion solution and coordinate rotation module 400, result calculation and output module 500. DETAILED DESCRIPTION

[0046] The present invention aims to provide a method, system, medium and program product for inverting the atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar, so as to obtain the vertical profiles of turbulence statistics such as velocity variance, momentum flux and turbulent kinetic energy in the atmospheric boundary layer. To make the purpose, technical solution and advantages of the present invention clearer, the technical solutions in the embodiments of the present invention will be described in more detail below with reference to the accompanying drawings in the embodiments of the present invention. The described embodiments are part of the embodiments of the present invention, rather than all of the embodiments, and the described embodiments are exemplary and are intended to explain the present invention, and should not be construed as a limitation to the present invention.

[0047] Embodiment 1: Inversion method

[0048] As a specific example, the method for inverting the atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar provided by the embodiment of the present invention is as Figure 1 shown, and when implemented, it mainly includes the following steps:

[0049] S100. Configuration of 6-beam DBS scanning mode:

[0050] Control the ground-based optical quantum wind lidar to perform observations according to the Doppler Beam Swinging (DBS) scanning mode, and set 6 detection beams, including 1 vertical upward zenith beam and 5 inclined beams with a set elevation angle φ and evenly spaced azimuth angles along the conical surface θ distributed.

[0051] In the embodiment of the present invention, the 5 inclined beams are evenly distributed along the azimuth angle with the zenith beam as the center on the conical surface formed with the vertical beam, and the azimuth angles of the inclined beams θ i are arranged at intervals of 72° in sequence starting from the due north direction to ensure the symmetry of the laser beam on the horizontal plane, that is, the azimuth angles are set as θ 1 = 0°, θ 2 = 72°, θ 3 = 144°, θ 4 = 216°, θ 5 = 288°. The elevation angles of the inclined beams are the same and are set in the range of 45° to 75° to ensure that the detection height covers the main range of the atmospheric boundary layer, and the 6 beams are emitted in sequence at fixed time intervals to achieve omnidirectional coverage of the typical scale turbulent perturbation structure.

[0052] S200. Collection and processing of radial wind speed data:

[0053] Synchronously obtain the time series data of the radial wind speed in the direction of each detection beam, and record the corresponding azimuth angleθ with the elevation angle φ parameters, and perform quality control on the acquired radial wind speed data, screening out data points with a signal-to-noise ratio lower than a preset threshold, and eliminating abnormal data points that significantly deviate from the average value or show mutations.

[0054] Preferably, the sampling period for traversing the radial wind speed time series data of 6 beams is 20 - 30 seconds, and the spatial resolution is set to 15 - 30 meters to ensure the ability to capture mesoscale and small-scale turbulent perturbation processes in the atmospheric boundary layer.

[0055] In addition, in step S200, the quality control of the radial wind speed data preferably includes the following sub-steps:

[0056] S201. Preprocess the original radial wind speed data, eliminating invalid or incorrect data caused by too low detection echo intensity, incorrect time stamps, echo oversaturation, or abnormal system return codes;

[0057] S202. Calculate the signal-to-noise ratio CNR of each detection point based on the characteristic parameters of the laser pulse echo signal, and screen out data points with a signal intensity insufficient to support reliable velocity estimation whose signal-to-noise ratio is lower than a preset threshold (preferably between -18 dB and -20 dB);

[0058] S203. Apply the sliding window statistical analysis method to the retained data for local mean and standard deviation calculation, and combine the Z-score method or the empirical upper and lower limit method to identify drastic jump points, extreme deviation points, or data spikes, and eliminate data points that do not conform to the characteristics of stationary perturbations;

[0059] S204. Perform linear interpolation or time series reconstruction on the data with local missing measurements caused by screening to ensure that all 6 beams have a consistent sampling time reference within the same observation period.

[0060] S300. Calculate the radial velocity variance of each beam:

[0061] For the radial wind speed time series data in the directions of each beam collected v ri ( t ) perform statistical processing to calculate the radial velocity variance corresponding to each beam , and construct a radial velocity variance observation vector accordingly is the radial velocity variance of the 1st - 6th beams, i = 1, 2, …, 6.

[0062] In the embodiments of the present invention, the calculation of the radial velocity variance of each beam includes the following sub-steps:

[0063] S301. For each beam direction iThe height layer corresponding to each distance k , for the radial wind speed time series data calculate its mean value within the time window T , where i = 1, 2, …, 6, k = 1, 2, …, N , N is the number of effective height layers that can be covered in the vertical direction;

[0064] S302. Based on the calculated mean value solve the sample variance of the radial velocity within the time window T , where n is the number of effective samples, t j is the sampling time point;

[0065] S303. For each beam direction and each height layer, organize and form a radial velocity variance matrix of multiple heights and multiple directions, as the basic input data for subsequent inversion of turbulence statistics. When necessary, linearly interpolate and complete the height layers with data gaps to ensure the structural continuity and calculation stability of the vertical profile.

[0066] S400. Establish a relationship between the radial velocity variance and the turbulence statistics:

[0067] According to the elevation angle φ i and azimuth angle θ i corresponding to each beam, establish the relationship between the radial velocity variance of each beam and six turbulence statistics respectively:

[0068]

[0069] where are the pulsating wind speed components in the horizontal east-west, horizontal north-south, and vertical directions respectively are the velocity variances in the three directions respectively, are the covariances of the velocity components in the three directions, that is, the momentum flux terms.

[0070] It should be noted that when establishing the relationship between the radial velocity variance and the turbulence statistics, based on the assumptions of the homogeneity and isotropy of atmospheric turbulence (that is, assuming that within the atmospheric boundary layer, the statistical characteristics of turbulence are uniformly distributed in space and the same in all directions), the radial velocity variance is expressed as a linear combination of six turbulence statistics (including the velocity variances in the three directions and the momentum fluxes in the three directions); the coefficients of this linear combination are determined by the elevation angle of each beamφ i and azimuth angle θ i is determined and reflects the contribution of the wind speed components in different directions to the radial velocity variance.

[0071] S500. Construct the observation equations and parameter matrix of the radial velocity variance:

[0072] Based on the relationship between the radial velocity variance and the turbulence statistics in 6 beam directions, construct a linear observation equation of the radial velocity variance in matrix form , where S is the observation vector of the radial velocity variance, M is the 6×6 coefficient matrix calculated according to the elevation angle φ i and azimuth angle θ i of each beam, Σ is the vector of turbulence statistics to be solved, and .

[0073] It should be noted that when determining the coefficient matrix M of the observation equations, according to the elevation angle φ i and azimuth angle θ i of each beam, calculate the projection coefficients of each beam in the three coordinate axes directions, and arrange these projection coefficients in a certain order to form the coefficient matrix M, whose dimension is 6×6, each row corresponds to a beam, and each column corresponds to a turbulence statistic. For example, the 6 matrix elements corresponding to the i row are respectively

[0074] . This coefficient matrix reflects the sensitivity of each beam to different turbulence statistics, and its numerical value directly affects the accuracy of the inversion of turbulence statistics; in actual calculations, it is necessary to ensure the measurement accuracy of the elevation angle φ i and azimuth angle θ i of each beam to avoid the influence of the calculation error of the coefficient matrix M on the inversion result.

[0075] S600. Inversion solution of turbulence statistics and coordinate rotation:

[0076] Based on the observation equations of the radial velocity variance , use the inverse matrix method or the least squares fitting method to inversely solve the vector Σ of turbulence statistics, and perform coordinate rotation on the inversion result based on the mean wind direction to ensure that the velocity perturbation variance and covariance are expressed in the coordinate system aligned with the mean wind direction.

[0077] Preferably, the coordinate rotation can adopt the double coordinate rotation method, including the following sub-steps:

[0078] S621. Calculate the horizontal average wind speed direction at each altitude layer based on the average component of the wind speed obtained by inversion , where and are the east-west and north-south components of the average wind speed respectively, and based on this, construct a first rotation matrix that rotates from the radar local coordinate system to the average wind direction coordinate system to eliminate the influence of the velocity distribution projection in the non-principal axis direction;

[0079] S622. Based on the covariance term of the vertical velocity component and the horizontal wind speed component (such as ),

[0080] construct a second rotation matrix to align the vertical direction completely with the ground normal direction, project the turbulent perturbation tensor onto the principal axis system in the sense of atmospheric dynamics, realize the physical consistent expression of the turbulent flux term, and at the same time eliminate the cross-covariance components of non-physical solutions.

[0081] In addition, when using the inverse matrix method to inversely solve the turbulent statistic vector Σ, it includes the following sub-steps:

[0082] S611. Judge whether the coefficient matrix M of the observation equation system is invertible. If M is not invertible, use the singular value decomposition (SVD) method for processing;

[0083] S612. Calculate the inverse matrix M -1 of the coefficient matrix M, and then multiply the radial velocity variance observation vector S by the inverse matrix M -1 to obtain the solution of the turbulent statistic vector Σ.

[0084] It should be noted that when using the least squares fitting method to inversely solve the turbulent statistic vector Σ, first, a target function needs to be constructed, and this target function represents the sum of the squared residuals between the observed value of the radial velocity variance and the calculated value of the model; then, use an optimization algorithm (such as the gradient descent method or the Levenberg-Marquardt algorithm) to minimize the target function to obtain the optimal solution of the turbulent statistic vector Σ.

[0085] Furthermore, in step S600, the ill-conditioned matrix analysis and regularization processing method is used to solve the possible ill-conditioned problem of the coefficient matrix M, including:

[0086] Calculate the condition number of the coefficient matrix M to judge the ill-conditioned degree of the observation equation system;

[0087] When the condition number of the coefficient matrix M exceeds the preset threshold, use the regularization method to solve the corrected equation system , where λ is the regularization parameter and I is the identity matrix;

[0088] Determine the optimal regularization parameter by the L-curve method or the generalized cross-validation method λ ;

[0089] Solve for the turbulent statistics Σ based on the corrected equations, and perform uncertainty analysis on the inversion results through the Monte Carlo simulation method, and output the confidence interval of each turbulent statistic.

[0090] S700. Turbulent kinetic energy TKE calculation and result output:

[0091] Calculate the turbulent kinetic energy based on the velocity variance obtained by inversion , and finally output the vertical profiles of turbulent statistics including velocity variance, momentum flux, and turbulent kinetic energy at each height layer, which are used to characterize the vertical structure of turbulent perturbations in the atmospheric boundary layer.

[0092] S800. Test the physical validity of the turbulence parameters, including:

[0093] Verify the physical rationality of the inversion results based on the constraint conditions of atmospheric turbulence theory, including testing the non-negativity condition of the velocity variance ; verify that the covariance satisfies the Cauchy-Schwarz inequality , , ; calculate the turbulent anisotropy index and verify whether its variation range conforms to the characteristics of atmospheric boundary layer turbulence; test whether the sign and vertical variation trend of the momentum flux under different atmospheric stability conditions conform to physical laws, and mark or correct the results that do not meet the constraint conditions.

[0094] S900. Evaluate the errors of the inverted turbulent statistics, including:

[0095] S901. Multi-source comparison error evaluation: Compare and verify the turbulent statistics inverted at each height layer with other existing observational data (such as the temperature and wind field profiles provided by the radiosonde system, the turbulent kinetic energy and momentum flux measured by the boundary layer tower mast ultrasonic anemometer, and the high-resolution wind speed perturbation obtained by the airborne platform), calculate the deviation, root mean square error, and / or correlation coefficient of the inversion results, and evaluate the physical consistency and observational credibility of the inversion results under the actual atmospheric background;

[0096] S902. Error source identification: Identify the main sources of inversion errors, including lidar system accuracy limitations, echo signal SNR fluctuations, insufficient observation time window, statistical biases caused by the non-stationarity of boundary layer turbulence, and / or beam angle projection errors, and label the potential error dominant factors for different height layers and different beam directions respectively during the identification process, providing a basis for subsequent propagation modeling and weighting processing;

[0097] S903. Propagation Error Modeling and Quantitative Estimation: Matrix Structure Based on Inversion Algorithm Conduct quantitative analysis on the inversion error of turbulence statistics, and calculate the covariance matrix of turbulence statistics , obtain the estimated error variances and covariances of each inversion component, as well as the cross-correlation uncertainty between turbulence statistics, C S is the covariance matrix of the uncertainty of the observed radial velocity variance;

[0098] S904. Evaluation Quantity Output and Stratified Diagnosis: Output the relative error percentages and confidence interval boundaries of each turbulence statistic component in each height layer, identify the noise-dominated regions and high-confidence sections, and visualize the inversion results and error envelopes for subsequent quality control and data release reference.

[0099] In summary, Example 1 elaborates in detail the complete technical process of the method for inverting the atmospheric boundary layer turbulence statistics based on the 6-beam ground-based optical quantum wind lidar DBS scanning mode of the present invention. From the observation configuration, data processing to parameter inversion and error evaluation, the vertical profile of turbulence statistics with high spatio-temporal resolution is finally obtained, reflecting the comprehensive advantages of the present invention in terms of accuracy, stability and practicality.

[0100] Example 2: Inversion System

[0101] Based on the method for inverting the atmospheric boundary layer turbulence statistics based on the 6-beam ground-based optical quantum wind lidar DBS scanning mode proposed in Example 1, this Example 2 further provides an inversion system structure that can be actually deployed and run. This system combines the six-beam DBS scanning observation ability of the ground-based optical quantum wind lidar and the multi-module collaborative calculation architecture, and can efficiently complete the extraction and output of turbulence statistics, as specifically Figure 2 shown. Through modular design, this system decomposes the complex inversion process into several independent modules with clear functions, which is convenient for the development, maintenance and upgrade of the system. At the same time, standardized interfaces are used for data exchange between modules, ensuring the flexibility and scalability of the system. This system mainly includes the following modules:

[0102] The lidar observation module 100 is used to control the ground-based optical quantum wind lidar to emit 6 detection beams in the DBS scanning mode, including 1 vertical zenith beam and 5 inclined beams evenly distributed along the conical surface. In terms of hardware implementation, this module includes a laser emission sub-module, a beam control sub-module, and a timing synchronization sub-module. The laser emission sub-module adopts optical quantum detection technology in the 1.5μm band to generate laser pulses that meet the requirements of atmospheric detection through a fiber amplifier; the beam control sub-module realizes the rapid switching of six beams through a high-precision galvanometer system, where the elevation angle of the inclined beam can be dynamically set within the range of 45°-75°, and the azimuth angle is evenly distributed at intervals of 72°, forming a complete ring-cone structure; the timing synchronization sub-module ensures that the emission and reception of each beam strictly follow the preset time sequence and maintain time reference synchronization with the subsequent data processing module. This module realizes the sequential emission of 6 beams at fixed time intervals through a preset scanning program and the lidar control system, and can adjust parameters such as the scanning period and beam elevation angle according to different meteorological conditions or application requirements to meet the observation requirements of multi-scale turbulent perturbations in the atmospheric boundary layer, providing stable, continuous, and well-covered original wind field information for subsequent data processing and parameter inversion.

[0103] The data acquisition and processing module 200 is used to synchronously obtain the time series data of the radial wind speed in the direction of each detection beam, record the corresponding azimuth angle and elevation angle parameters, and perform quality control on the radial wind speed data. In terms of hardware implementation, this module includes a multi-channel data acquisition card, a signal processing sub-module, and a quality control sub-module. The data acquisition card synchronously records the echo signals in the directions of six beams at a sampling rate not lower than 1MHz. The signal processing sub-module extracts the radial wind speed values of each range gate (with a resolution of 15-30m) through fast Fourier transform FFT and calculates quality indicators such as the carrier-to-noise ratio CNR. The quality control sub-module first performs system preprocessing on the original data, including abnormal code identification, time alignment verification, and range gate channel screening; then calculates the signal-to-noise ratio CNR based on the laser echo signal and eliminates unreliable sample points below the preset threshold (–18 dB to –20 dB); further uses the sliding window method and the Z-score method to identify local wind speed mutation points and deviation outliers to ensure that only data points that meet stable statistical characteristics are retained. In addition, to solve the data discontinuity problem caused by quality control, this module also supports linear interpolation or time series reconstruction of the time series to ensure that the sampling time reference of all beams is consistent.

[0104] The statistic calculation and matrix construction module 300 is used to calculate the radial velocity variance of each beam at each altitude layer and establish a linear equation system between it and the turbulent statistics based on the beam geometric angles. This module corresponds to the calculation process of steps S300-S500 in Embodiment 1 and consists of three sub-modules working together:

[0105] The variance calculation sub-module performs statistical analysis on the wind speed time series (30 - 60 s window) of each beam and each altitude layer to calculate the radial velocity variance; the matrix construction sub-module, based on the preset beam geometric parameters ( φ i , θ i ), automatically generates a 6×6 coefficient matrix M according to the trigonometric function relation given in Embodiment 1; the equation system assembly sub-module combines the observation vector S and the coefficient matrix M into a complete linear equation system

[0106] S = M·Σ. To improve the calculation efficiency, this module adopts a parallel computing architecture, processes data of multiple altitude layers simultaneously, and realizes fast access to large amounts of data through the memory mapping technology. There is also a cache mechanism inside the module. When some beam data is missing, it can be intelligently filled based on historical data to ensure the continuity of the equation system solution. In addition, this module also has the function of interpolating and complementing the missing altitude layers to ensure the continuity of the vertical profile, and supports batch construction of multi-temporal equation system structures to realize automatic preprocessing of a large amount of data under a sliding window.

[0107] The inversion solution and coordinate rotation module 400 is used to solve the turbulent statistic vector by using the inverse matrix method or the least squares fitting method, and perform coordinate transformation based on the mean wind direction. This module is the core calculation unit of the system. The solution operator sub-module in it provides two algorithm options: directly using the LU decomposition method to find the inverse for a well-conditioned equation system; for ill-conditioned cases (condition number > 100), the Tikhonov regularization algorithm is enabled, and the optimal regularization parameter is automatically determined by the L-curve method. The coordinate rotation sub-module uses the quaternion method to realize three-dimensional coordinate system transformation: first, calculate the first rotation angle based on the mean horizontal wind to eliminate the wind direction influence; then calculate the second rotation angle according to the vertical flux to align the coordinate system with the main axis of the mean wind field. This module also includes a result verification sub-module to check whether the inverted turbulent statistics satisfy and other physical constraint conditions, and automatically trigger the recalculation process for abnormal results.

[0108] The result calculation and output module 500 calculates the turbulent kinetic energy based on the inverted velocity variance and outputs the vertical profiles of the turbulent statistics of each altitude layer. The turbulent kinetic energy calculation sub-module is responsible for calculating the turbulent kinetic energy values of each altitude layer; the boundary layer height detection sub-module automatically identifies the position of the boundary layer top through the vertical velocity variance threshold method (w′² < 0.3 m² / s²); the data output sub-module generates standard data products including velocity variance, momentum flux, turbulent kinetic energy, and boundary layer height, supporting multiple scientific data formats such as NetCDF and HDF5. In addition, it also includes a quality control sub-module to check the spatial consistency and time continuity of the inversion results, mark suspicious data points, and generate metadata information including error estimates.

[0109] Through the collaborative operation of the above modules, the system can completely and efficiently implement all the key links of the inversion method described in Embodiment 1, and has good adaptability, scalability, and engineering deployment capabilities.

[0110] Embodiment 3: Example Verification

[0111] To verify the practicability and reliability of the atmospheric boundary layer turbulence statistic inversion method based on the 6-beam ground-based optical quantum wind lidar DBS scanning mode proposed in the present invention, a six-day radar field measurement and observation experiment was continuously carried out from 00:00 on November 16, 2023 to 23:30 on November 21, 2023 (Beijing time) at the experimental site in Hefei, Anhui. The observation equipment used a ground-based optical quantum wind lidar with 6-beam DBS scanning ability. The beam configuration was the same as that described in Embodiment 1. The data sampling period was 1 second, the sliding time window was 60 seconds, the vertical resolution was set to 30 meters, and the coverage height was up to 2.7 km at most. After data collection, steps such as radial wind speed quality control, velocity variance calculation, turbulence statistic matrix inversion, coordinate rotation, and turbulent kinetic energy solution were completed according to the technical process of Embodiment 1.

[0112] Figure 3 The time-height profile distributions of typical turbulence statistics inverted by the method of the present invention during the above observation period are shown. In the figure, (a)-(d) respectively represent the east-west velocity variance , the north-south velocity variance , the vertical velocity variance , and the turbulent kinetic energy . The pink and cyan curves respectively represent the boundary layer heights Z - z i determined based on the wind lidar and ERA5 reanalysis data. It can be seen from the figure that the four turbulence statistics show good spatio-temporal structure characteristics in different diurnal cycles, and show a highly coupled correlation with the change of the boundary layer height. Especially during the daytime convective development stage, and both increase significantly in the height range of 500 - 1500 m, and can clearly depict the strong turbulent activities inside the typical convective boundary layer. When the nocturnal stable layer appears, the three-component velocity variances quickly weaken and show a characteristic of concentrating near the ground, reflecting the sensitivity of this method to capturing the boundary layer turbulence suppression process. In addition, in the figure Z - z i and ERA5- z iThe diurnal variation trends are generally consistent. Especially on days with severe convection (such as November 16th, 18th, and 21st), the peak heights of the two are close, verifying the effectiveness of the present invention in boundary layer height identification. It should be noted that under weak wind and stable conditions (such as the night of November 19th), ERA5- z i shows an obvious lag phenomenon, while Z - z i can still identify the local turbulence enhancement process, indicating that the inversion mechanism of the present invention has a higher timeliness response ability and vertical structure resolution ability. In terms of boundary layer height estimation, the dual-criterion method adopted by the present invention, that is, the vertical velocity variance is less than 0.3 m² / s² for the first time or the signal-to-noise ratio CNR is lower than 0.4, is used as the boundary layer top criterion. From Figure 3 the continuity and boundary clarity, this method has strong stability and discrimination ability in practical applications and can be used for the automatic extraction of boundary layer height in business operations and the real-time release of turbulence products.

[0113] In summary, the measured analysis results of Example 3 fully show that the inversion method proposed by the present invention can not only accurately reflect the three-dimensional structure evolution process of turbulent disturbances, but also has good boundary layer structure identification ability and consistency in comparison with reanalysis data, and has strong engineering adaptability and scientific research promotion value.

[0114] Through the above embodiments, the purpose of the present invention is completely and effectively achieved. Those skilled in the art can understand that the present invention includes but is not limited to the content described in the drawings and the above specific embodiments. Although the present invention has been described with respect to the currently considered most practical and preferred embodiments, it should be understood that the present invention is not limited to the disclosed embodiments, and any modification that does not deviate from the functional and structural principles of the present invention will be included in the scope of the claims.

Claims

1. An inversion method for atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar, characterized in that Including: S100. Set the ground-based quantum wind lidar to the DBS scanning mode, and configure 1 zenith beam vertically upward and 5 tilted beams distributed at equal azimuth angles with a set elevation angle; S200. Synchronously obtain the radial wind speed time series data in each beam direction, and record the corresponding azimuth angle θ and elevation angle φ parameters, and perform quality control on the wind speed data; S300. Statistically analyze the radial wind speed data of each beam, and calculate the variance of each radial velocity respectively , i = 1, 2, …, 6, and construct an observation vector of the radial velocity variance accordingly , where , , …, are the variances of the radial velocities of the 1st to 6th beams; S400. According to the elevation angle corresponding to each beam φ i and the azimuth angle θ i , establish the relational expressions between the radial velocity variances and the six turbulence statistics respectively: wherein, are the pulsating wind speed components in the east-west, north-south, and vertical directions respectively, and are the velocity variances and covariances in the three directions respectively; S500. Construct the Observation Equation System for Radial Velocity Variance , where S is the observation vector of radial velocity variance, and M is a 6×6 coefficient matrix calculated based on the elevation angle φ i and azimuth angle θ i of each beam, Σ is the vector of turbulence statistics to be solved, and ; S600. Observation equation set based on radial velocity variance , the inverse matrix method or the least squares fitting method is used to inversely solve the turbulence statistic vector Σ, and coordinate rotation is performed on the inversion result based on the mean wind direction to ensure that the velocity perturbation variance and covariance are expressed in the coordinate system aligned with the mean wind direction; S700. Calculate the turbulent kinetic energy based on the velocity variance obtained by inversion Finally, output the vertical profiles of velocity variance, momentum flux, and turbulent kinetic energy at each height layer, which are used to characterize the vertical structure of the turbulence in the atmospheric boundary layer.

2. The method for retrieving the atmospheric boundary layer turbulence statistic based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar according to claim 1, characterized in that In step S100, five tilted beams are evenly distributed along the azimuth angle centered on the zenith beam on the conical surface formed with the vertical beam, and the azimuth angles of the respective tilted beams θ i are arranged at intervals of 72° successively starting from the due north direction to ensure the symmetry of the laser beam on the horizontal plane. The elevation angles of the respective tilted beams are the same and are set within the range of 45° to 75° to ensure that the detection height covers the main range of the atmospheric boundary layer. Moreover, the six beams are emitted successively at fixed time intervals to achieve omnidirectional coverage of the typical-scale turbulent perturbation structure.

3. The method for retrieving the atmospheric boundary layer turbulence statistic based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar according to claim 1, wherein In step S200, the sampling period for traversing the radial wind speeds of the 6 beams is 20 - 30 seconds, and the spatial resolution is set to 15 - 30 meters.

4. The method for retrieving the atmospheric boundary layer turbulence statistic based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar according to claim 1 or 3, characterized in that, In step S200, the quality control of the radial wind speed data includes the following sub-steps: S201. Preprocess the original radial wind speed data, and eliminate invalid data or error data caused by too low detection echo intensity, disordered timestamps, echo oversaturation, or abnormal system return codes; S202. Calculate the signal-to-noise ratio CNR of each detection point based on the characteristic parameters of the laser pulse echo signal, and screen out the data points with signal strength lower than the preset threshold and insufficient to support reliable speed estimation; S203. Apply the sliding window statistical analysis method to the retained data for local mean and standard deviation calculation, and combine the Z-score method or the empirical upper and lower limit method to identify sharp jump points, extreme deviation points, or data spikes, and eliminate the data points that do not conform to the characteristics of stationary disturbances; S204. Perform linear interpolation or time series reconstruction on the data with local missing measurements caused by screening to ensure that all 6 beams have a consistent sampling time reference within the same observation period.

5. The method for retrieving the atmospheric boundary layer turbulence statistic based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar according to claim 1, characterized in that In step S300, the calculation of the radial velocity variance of each beam includes the following sub-steps: S301. For each beam direction i and each altitude layer corresponding to each distance k , calculate the mean value of the radial wind speed time series data within the time window T, where i i = 1, 2, …, 6, k j = 1, 2, …, N , N and N is the number of effective altitude layers that can be covered in the vertical direction; S302. Based on the calculated mean value Solve for the sample variance of the radial velocity within the time window T , where n is the number of valid samples, t j is the sampling time point; S303. For each beam direction and each altitude layer, organize and form a radial velocity variance matrix of multiple altitudes and multiple directions as the basic input data for subsequent inversion of turbulence statistics.

6. The method for inverting the atmospheric boundary layer turbulence statistic according to the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar as claimed in claim 1, wherein In step S600, the coordinate rotation adopts the double coordinate rotation method, including the following sub-steps: Based on the average component of the wind speed obtained by inversion, calculate the horizontal average wind speed direction of each altitude layer , where and are the east-west and north-south components of the average wind speed respectively, and based on this, construct a first rotation matrix that rotates from the radar local coordinate system to the average wind direction coordinate system to eliminate the influence of the velocity distribution projection in the non-principal axis direction; S622. Construct a second rotation matrix based on the covariance term between the vertical velocity component and the horizontal wind speed component to completely align the vertical direction with the ground normal direction, project the turbulence perturbation tensor onto the main axis system in the sense of atmospheric dynamics, realize the physical consistency expression of the turbulence flux term, and eliminate the cross-covariance components of non-physical solutions at the same time.

7. The method for retrieving the atmospheric boundary layer turbulence statistic according to the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar as claimed in claim 1, wherein In step S600, when using the inverse matrix method to inversely solve the turbulence statistic vector Σ, it includes the following sub-steps: S611. Judge whether the coefficient matrix M of the observation equation system is invertible. If M is not invertible, the singular value decomposition method is used for processing; S612. Calculate the inverse matrix M of the coefficient matrix M -1 , and then multiply the radial velocity variance observation vector S by the inverse matrix M -1 to obtain the solution of the turbulence statistic vector Σ.

8. The method for retrieving the atmospheric boundary layer turbulence statistic quantity based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar according to claim 1, characterized in that In step S600, the ill-conditioned matrix analysis and regularization processing method are used to solve the possible ill-conditioned problem of the coefficient matrix M: Calculate the condition number of the coefficient matrix M to judge the ill-conditioned degree of the observation equation system; When the condition number of the coefficient matrix M exceeds a preset threshold, a regularization method is used to solve the modified system of equations , where λ is the regularization parameter and I is the identity matrix; Determine the optimal regularization parameter by the L-curve method or the generalized cross-validation method λ ; Solve the turbulence statistic Σ based on the modified equation system, and perform uncertainty analysis on the inversion result through the Monte Carlo simulation method, and output the confidence interval of each turbulence statistic.

9. The method for retrieving the atmospheric boundary layer turbulence statistic according to the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar as claimed in claim 1, wherein It also includes step S800 to test the physical validity of the turbulence parameters, including the following sub-steps: S801. Verify the physical rationality of the inversion result based on the constraint conditions of atmospheric turbulence theory, including checking the non-negativity condition of the velocity variance ; S802. Verify that the covariance satisfies the Cauchy - Schwarz inequality ; S803. Calculate the turbulence anisotropy index and verify whether the variation range conforms to the characteristics of the atmospheric boundary layer turbulence; S804. Test whether the sign and vertical variation trend of the momentum flux under different atmospheric stability conditions conform to physical laws, and mark or correct the results that do not meet the constraint conditions.

10. The method for retrieving the atmospheric boundary layer turbulence statistic based on the DBS scanning mode of the 6-beam ground-based optical quantum wind lidar according to claim 9, characterized in that, It also includes step S900 for error evaluation of the retrieved turbulence statistics, and the error evaluation includes the following sub-steps: S901. Compare and verify the turbulence statistics retrieved at each altitude layer with other existing observational data, calculate the deviation, root mean square error, and / or correlation coefficient of the retrieval results, and evaluate the physical consistency and observational credibility of the retrieval results under the actual atmospheric background; S902. Identify the main sources of retrieval errors, including lidar system accuracy limitations, signal-to-noise ratio fluctuations of echo signals, insufficient observation time window, statistical biases caused by the non-steadiness of boundary layer turbulence, and / or beam angle projection errors; Matrix Structure Based on Inversion Algorithm , quantitatively analyze the inversion error of turbulent statistics, and calculate the covariance matrix of turbulent statistics , obtain the estimated error variances and covariances of each inversion component and the cross-correlation uncertainty between turbulent statistics, where C S is the uncertainty covariance matrix of radial velocity variance observations; S904. Output the relative error percentages and confidence interval boundaries of each turbulence statistical component in each altitude layer, identify the noise-dominated regions and high-confidence sections, and visualize the retrieval results and error envelopes.

11. An atmospheric boundary layer turbulence statistic inversion system for implementing the atmospheric boundary layer turbulence statistic inversion method based on the DBS scanning mode of a 6-beam ground-based optical quantum wind lidar according to any one of claims 1 to 10, characterized in that, It includes: A lidar observation module that controls a ground-based optical quantum wind measurement lidar to emit 6 detection beams in the DBS scanning mode, including 1 vertical zenith beam and 5 inclined beams evenly distributed along the conical surface; A data acquisition and processing module for synchronously acquiring the radial wind speed time series data in the directions of each detection beam, recording the corresponding azimuth and elevation angle parameters, and performing quality control on the radial wind speed data; A statistic calculation and matrix construction module for calculating the radial velocity variance of each beam at each altitude layer and establishing a linear equation set between it and the turbulence statistics based on the beam geometric angles; An inversion solution and coordinate rotation module for solving the turbulence statistic vector by using the inverse matrix method or the least square fitting method and performing coordinate transformation based on the mean wind direction; A result calculation and output module for calculating the turbulent kinetic energy based on the retrieved velocity variance and outputting the vertical profiles of the turbulence statistics at each altitude layer.

12. A computer program product comprising computer instructions, characterized in that, It is used to execute the method for retrieving the atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind measurement lidar according to any one of claims 1 to 10.

13. A computer-readable storage medium, on which a computer program is stored, characterized in that, When the computer program is executed by a processor, it realizes the method for retrieving the atmospheric boundary layer turbulence statistics based on the DBS scanning mode of a 6-beam ground-based optical quantum wind measurement lidar according to any one of claims 1 to 10.

Images