A wind speed adaptive eddy flux frequency response correction method

Abstract

The application discloses a wind speed adaptive vortex flux frequency response correction method and relates to the technical field of electric digital data processing. The wind speed adaptive vortex flux frequency response correction method comprises the following steps: obtaining high-frequency observation data containing three-dimensional wind speed and gas concentration, calculating initial flux after pretreatment and coordinate rotation attack angle correction; dynamically updating three groups of transfer functions of path averaging, sensor separation and instrument response with real-time average wind speed, and integrating frequency response correction factors according to the matching of the atmospheric stability and the co-spectrum model; simultaneously, adaptively selecting block averaging, linear detrending or exponential filtering according to the turbulence state parameters, and completing high-frequency loss compensation by combining the correction factors, and finally outputting the final flux and mass level after WPL correction. The wind speed driving transfer function updating and filter adaptive switching are adopted, so that the fixed parameter correction deviation problem is solved, and the automatic connection of the whole process and the synchronous output of the flux and mass are realized.

Description

Technical Field

[0001] This invention relates to the field of electronic digital data processing technology, specifically to a wind speed adaptive eddy flux frequency response correction method. Background Technology

[0002] Eddy covariance is a commonly used data processing method for calculating fluxes of mass and energy exchange between the Earth's surface and the atmosphere. In practical flux data processing, physical effects such as path averaging, sensor separation, dynamic response delay, and high-frequency attenuation in ultrasonic anemometers and gas analyzers can cause high-frequency turbulence signal loss, affecting the accuracy of calculated results for carbon dioxide flux, water vapor flux, sensible heat flux, and latent heat flux. Therefore, frequency response correction is usually required during flux calculation to compensate for high-frequency signal loss and improve the accuracy of the flux calculation results.

[0003] Currently, common eddy flux calculation algorithms typically use fixed instrument parameters and transfer functions during frequency response correction, without dynamically adjusting the relevant correction parameters based on real-time wind speed. This makes it difficult to guarantee the correction effect under low and high wind speed conditions. High-pass filtering methods usually execute according to a preset fixed mode without adaptive matching based on the actual turbulence state, which can easily lead to insufficient or excessive high-frequency loss compensation.

[0004] The limitations of existing technologies include at least the following problems: Eddy flux frequency response correction methods typically use fixed instrument parameters or fixed transfer functions, failing to dynamically adjust parameters related to path averaging, sensor spatial separation, instrument dynamic response delay, and high-frequency attenuation based on real-time wind speed. This results in excessively large high-frequency loss correction errors under low or high wind speed conditions. High-pass filtering methods typically use fixed modes, failing to adaptively select from block averaging, linear detrending, and exponential filtering based on the on-site turbulence conditions. This easily leads to insufficient or excessive high-frequency loss compensation. Under low wind speed, complex terrain, and high humidity observation conditions, these errors are further amplified, causing deviations between flux calculation results and actual flux values, affecting the accuracy, stability, and comparability of flux data. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a wind speed-adaptive eddy flux frequency response correction method. This method solves the problems of fixed parameters and filtering modes in existing eddy flux frequency response correction techniques, which lack wind speed and filtering adaptation capabilities, leading to large correction errors and reduced accuracy and stability of flux calculation under different wind speeds and complex observation conditions.

[0006] To achieve the above objectives, the present invention provides the following technical solution: a wind speed adaptive eddy flux frequency response correction method, comprising the following steps: acquiring raw high-frequency eddy flux observation data containing three-dimensional wind speed, carbon dioxide concentration, water vapor concentration, and atmospheric temperature; obtaining preprocessed observation data through outlier removal and anomaly repair; sequentially performing coordinate rotation and angle of attack correction on the preprocessed observation data to obtain corrected observation data, and calculating the initial carbon dioxide flux, initial water vapor flux, initial sensible heat flux, and initial latent heat flux; extracting the real-time average wind speed from the corrected observation data, and updating the path average effect transfer function, sensor spatial separation effect transfer function, and instrument dynamic response based on the real-time average wind speed. The frequency response correction factor is obtained by delaying the transfer function, selecting a cospectral model based on atmospheric stability parameters, and integrating the results. The target filtering method is determined from block average filtering, linear detrending filtering, and exponential filtering based on the field turbulence state parameters. The corrected observation data is processed using the target filtering method, and the frequency response correction factor is used to compensate for high-frequency losses in the initial fluxes, resulting in frequency-response corrected fluxes. WPL correction is applied to the frequency-response corrected carbon dioxide and water vapor fluxes, and combined with the frequency-response corrected sensible heat flux and the latent heat flux determined by the WPL corrected water vapor flux, the corrected fluxes are formed. A quality evaluation is completed based on stability and turbulence characteristic indices, and the corrected fluxes and their quality levels are output.

[0007] Furthermore, the specific steps for outlier removal and outlier repair of the original high-frequency eddy flux observation data are as follows: calculate the mean and standard deviation of each variable within the sliding window; mark the observation values ​​that exceed the mean plus or minus a preset multiple of the standard deviation as outliers, and mark the outlier positions as missing values; perform interpolation repair on the missing positions formed by the outlier markings to generate preprocessed observation data.

[0008] Furthermore, the specific steps for performing coordinate rotation and angle of attack correction on the preprocessed observation data are as follows: calculate the prevailing turbulent wind direction and perform coordinate rotation to bring the average lateral wind speed and average vertical wind speed to zero; extract the airflow injection angle deviation of the data after coordinate rotation; based on the pre-calibrated correspondence between the angle of attack and the vertical wind speed component deviation, compensate and correct the vertical wind speed component to generate the corrected observation data.

[0009] Furthermore, the specific steps for updating the path average effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function are as follows: obtain the instrument optical path length, the sensor horizontal separation distance, and the instrument dynamic response time constant; update the path average effect transfer function based on the real-time average wind speed and the instrument optical path length; update the sensor spatial separation effect transfer function based on the real-time average wind speed and the sensor horizontal separation distance; update the instrument dynamic response delay transfer function based on the real-time average wind speed and the instrument dynamic response time constant.

[0010] Furthermore, the specific steps for constructing the path average effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function are as follows: the path average effect transfer function is constructed based on the correlation between the instrument optical path length, real-time average wind speed, and frequency; the sensor spatial separation effect transfer function is constructed based on the correlation between the sensor horizontal separation distance, real-time average wind speed, and frequency; and the instrument dynamic response delay transfer function is constructed based on the correlation between the instrument dynamic response time constant, real-time average wind speed, and frequency.

[0011] Furthermore, the specific steps for selecting a cospectral model and calculating the frequency response correction factor by integrating atmospheric stability parameters are as follows: Obtain atmospheric stability parameters including the ratio of observation height to Moning-Obukhov length; Select a neutral cospectral model, a stable cospectral model, or an unstable cospectral model based on the comparison results of atmospheric stability parameters and preset thresholds; Multiply the updated path-averaged effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function to obtain the total transfer function; Perform a frequency-domain weighted integration of the total transfer function and the selected cospectral model to determine the frequency response correction factor.

[0012] Furthermore, the expression for the total transfer function is as follows: In the formula, For the updated path-average effect transfer function, For the updated sensor spatial separation effect transfer function, This is the updated instrument dynamic response delay transfer function.

[0013] Furthermore, the specific steps for determining the target filtering method based on the on-site turbulence state parameters are as follows: Turbulence intensity and frictional wind speed are collected as on-site turbulence state parameters; the turbulence intensity is compared with a preset turbulence intensity threshold, and the frictional wind speed is compared with a preset frictional wind speed threshold; when the on-site turbulence state meets the steady turbulence condition, block average filtering is determined as the target filtering method; when the on-site turbulence state meets the low-frequency trend condition, linear detrending filtering is determined as the target filtering method; when the on-site turbulence state meets the low-turbulence or weak turbulence condition, exponential filtering is determined as the target filtering method; when the on-site turbulence state simultaneously meets multiple filtering conditions, the target filtering method is determined according to a preset priority.

[0014] Furthermore, the specific steps for performing WPL correction on the frequency response corrected carbon dioxide flux and water vapor flux are as follows: extract temperature fluctuation data, water vapor fluctuation data, and average gas concentration data; calculate density perturbation calibration parameters based on the temperature fluctuation data, water vapor fluctuation data, and average gas concentration data; and use the density perturbation calibration parameters to perform density effect correction on the frequency response corrected carbon dioxide flux and water vapor flux.

[0015] Furthermore, the specific steps for determining the flux quality level based on the stationarity index and the turbulence characteristic index are as follows: divide the observation period into multiple sub-periods and calculate the average value of each flux after correction in each sub-period; determine the stationarity index based on the deviation between the average value of each sub-period and the overall average value of the observation period; calculate the turbulence characteristic index based on the turbulence statistical parameters in the observation period; and determine the flux quality level based on the joint judgment result of the stationarity index and the turbulence characteristic index.

[0016] The present invention has the following beneficial effects:

[0017] (1) The wind speed adaptive eddy flux frequency response correction method updates the path average effect transfer function, sensor spatial separation effect transfer function and instrument dynamic response delay transfer function by using real-time average wind speed as the driving variable, replacing the existing fixed instrument parameters or fixed transfer function processing method. When the observed wind speed changes, the path average effect, sensor spatial separation effect and instrument dynamic response delay can be adjusted according to the current wind conditions, so that the frequency response correction process is matched with the actual wind speed conditions. This can reduce the problem of insufficient compensation for high frequency loss under low wind speed conditions and compensation deviation from reality under high wind speed conditions, so that the calculation results of carbon dioxide flux, water vapor flux, sensible heat flux and latent heat flux are closer to the actual observation state, and improve the usability of flux correction results under different wind conditions.

[0018] (2) The wind speed adaptive eddy flux frequency response correction method introduces the field turbulence state parameters and determines the target filtering mode among block average filtering, linear detrending filtering and exponential filtering. This avoids the long-term use of a single fixed mode for high-pass filtering. When the turbulence is relatively stable during the observation period, block average filtering can be used to retain the main turbulence information. When there are low-frequency trend terms in the observation data, linear detrending filtering can be used to weaken drift interference. When the observation period is in a low-turbulence or weak-turbulence state, exponential filtering can be used to reduce the impact of abnormal fluctuations. This method can make the filtering mode adapt to the field turbulence state, reduce the error caused by insufficient or excessive high-frequency loss compensation, and make the corrected flux result closer to the changes in the field observation signal.

[0019] (3) The wind speed adaptive eddy flux frequency response correction method introduces atmospheric stability parameters into the cospectral model selection process, determines the current atmospheric stratification state based on the ratio of observation height to Moning-Obukhov length, selects between corresponding cospectral models, and then performs frequency domain weighted integration by combining the updated path average effect transfer function, sensor spatial separation effect transfer function and instrument dynamic response delay transfer function to obtain the frequency response correction factor. This method can make the cospectral model match the actual turbulence spectrum shape, reduce the correction deviation caused by spectrum mismatch under complex terrain, high humidity and wind speed fluctuation conditions, and make the frequency response correction result more consistent with the field observation environment.

[0020] Of course, any product implementing this invention does not necessarily need to achieve all of the advantages described above at the same time. Attached Figure Description

[0021] Figure 1 This is a flowchart of a wind speed adaptive eddy flux frequency response correction method according to the present invention.

[0022] Figure 2 This is a flowchart illustrating the specific steps of performing coordinate rotation and angle of attack correction sequentially on preprocessed observation data in a wind speed adaptive eddy flux frequency response correction method of the present invention. Detailed Implementation

[0023] Please see Figure 1 This invention provides a technical solution: a wind speed adaptive eddy flux frequency response correction method, comprising the following steps: acquiring raw high-frequency eddy flux observation data containing three-dimensional wind speed, carbon dioxide concentration, water vapor concentration, and atmospheric temperature; obtaining preprocessed observation data through outlier removal and anomaly repair; sequentially performing coordinate rotation and angle of attack correction on the preprocessed observation data to obtain corrected observation data, and calculating the initial carbon dioxide flux, initial water vapor flux, initial sensible heat flux, and initial latent heat flux; extracting the real-time average wind speed from the corrected observation data, and updating the path average effect transfer function, sensor spatial separation effect transfer function, and instrument dynamic response delay transfer function based on the real-time average wind speed. The frequency response correction factor is obtained by integrating a recursive function and selecting a cospectral model based on atmospheric stability parameters. The target filtering method is determined from block average filtering, linear detrending filtering, and exponential filtering based on the field turbulence state parameters. The corrected observation data is processed using the target filtering method, and the frequency response correction factor is used to compensate for high-frequency losses in the initial fluxes, resulting in frequency-response corrected fluxes. WPL correction is applied to the frequency-response corrected carbon dioxide and water vapor fluxes, and the corrected fluxes are formed by combining the frequency-response corrected sensible heat flux and the latent heat flux determined by the WPL corrected water vapor flux. A quality evaluation is completed based on stability and turbulence characteristic indices, and the corrected fluxes and their quality levels are output.

[0024] The sampling frequency of the original high-frequency eddy flux observation data is configured to be 10Hz to 100Hz, which is determined according to the performance parameters of the observation instrument and the requirements of the observation scenario.

[0025] The size of the sliding window is adjusted according to the sampling frequency to match the processing window with the sampling frequency;

[0026] The initial carbon dioxide flux, initial water vapor flux, initial sensible heat flux, and initial latent heat flux were all calculated based on the three-dimensional wind speed, carbon dioxide concentration, water vapor concentration, and atmospheric temperature in the corrected observation data. The calculations were performed using the basic principle of the eddy covariance method, which is to obtain the initial flux value by averaging the product of the pulsating components of each variable over time.

[0027] Specifically, the steps for outlier removal and outlier repair of the original high-frequency eddy flux observation data are as follows:

[0028] The mean and standard deviation of each variable within the sliding window are calculated as follows:

[0029] The length of the sliding window can be set to 100-1000 data points, and can be adjusted according to the sampling frequency;

[0030] In one implementation, when the sampling frequency is 10Hz, the sliding window length is set to 300 data points, corresponding to 30 seconds of observation data;

[0031] For the three-dimensional wind speed components, carbon dioxide concentration, water vapor concentration and atmospheric temperature in the original high-frequency eddy flux observation data, the mean of all data points in each window is calculated by the arithmetic mean method within their respective sliding windows, and the standard deviation of the data in the window is obtained by the calculation logic of the sample standard deviation.

[0032] Observations exceeding the range of the mean plus or minus a preset multiple of the standard deviation are marked as outliers and processed accordingly.

[0033] The preset multiplier is set to 2 to 3 times, which can be adjusted according to the dispersion of the observed data.

[0034] In one implementation, the preset multiplier is set to 2.5 times, that is, when an observation value exceeds the range of mean minus 2.5 times the standard deviation to mean plus 2.5 times the standard deviation, the observation value is marked as an outlier;

[0035] The processing method is to set the data positions marked as outliers to empty or invalid values, while retaining the corresponding time index to avoid disrupting the time synchronization relationship of high-frequency multivariate data, and at the same time record the position and value of the outliers;

[0036] Interpolation is performed to repair the gaps caused by outlier markers, generating preprocessed observation data, specifically as follows:

[0037] Interpolation repair uses either linear interpolation or cubic spline interpolation, with cubic spline interpolation being preferred.

[0038] In one implementation, for a single isolated missing position, linear interpolation is performed using two adjacent valid data points. That is, the repair value of the missing position is calculated by linear proportional relationship based on the values ​​of valid data points before and after the missing position and the corresponding observation time.

[0039] For multiple consecutive missing positions, cubic spline interpolation is used to repair them, so that the repaired data can be connected with the trend of change of the previous and subsequent valid data.

[0040] After all missing locations are repaired, all valid data (including original valid data and repaired data) are integrated to form preprocessed observation data.

[0041] In this implementation scheme, abrupt changes in each observed variable are identified by a sliding window statistical method. Outlier locations are marked as missing values ​​and then interpolated for repair. This reduces the interference of abnormal pulses on subsequent flux calculations and preserves the time index of the original high-frequency observation data. It also avoids disrupting the synchronization relationship between three-dimensional wind speed, carbon dioxide concentration, water vapor concentration and atmospheric temperature, thus providing continuous and aligned input data for subsequent coordinate rotation, angle of attack correction and frequency response correction.

[0042] Specifically, such as Figure 2 As shown, the specific steps for performing coordinate rotation and angle of attack correction on the preprocessed observation data are as follows:

[0043] The prevailing turbulent wind direction is calculated using a quadratic coordinate rotation method, and the wind speed coordinate system is adjusted to bring the average lateral wind speed and average vertical wind speed to zero. Specifically:

[0044] First, the three-dimensional wind speeds (u, v, w, where u is the longitudinal wind speed, v is the lateral wind speed, and w is the vertical wind speed) in the preprocessed observation data are averaged over time to obtain the average longitudinal wind speed, average lateral wind speed, and average vertical wind speed.

[0045] The angle between the prevailing turbulent wind direction and the initial longitudinal coordinate axis is calculated by using the ratio of the average longitudinal wind speed to the average lateral wind speed, combined with the arctangent function.

[0046] Based on this included angle, a second coordinate rotation is performed. The first rotation rotates the lateral wind speed v and the longitudinal wind speed u to the direction of the prevailing turbulent wind, so that the average lateral wind speed after rotation is zero. The second rotation adjusts the coordinate axis of the vertical wind speed w, so that the average vertical wind speed after rotation is zero, thus completing the coordinate rotation operation.

[0047] In one implementation, the average longitudinal wind speed during the observation period is 2.3 m / s, the average lateral wind speed is 0.2 m / s, and the calculated turbulent prevailing wind angle is 4.8°. After rotation, both the average lateral wind speed and the average vertical wind speed approach 0.

[0048] The airflow injection angle deviation of the data after coordinate rotation is extracted as follows:

[0049] After the coordinate rotation is completed, the three-dimensional wind speed data (u', v', w') after rotation are obtained, where u' is the longitudinal wind speed along the main turbulent wind direction after rotation, v' is the lateral wind speed after rotation, and w' is the vertical wind speed after rotation.

[0050] The airflow injection angle deviation, also known as the angle of attack deviation, is calculated by combining the ratio of the vertical wind speed w' to the longitudinal wind speed u' after rotation with the arctangent function. It is used to characterize the deviation angle between the airflow injection direction and the horizontal plane of the turbulent prevailing wind direction.

[0051] In one implementation, after rotation, at a certain moment, u'=2.2m / s and w'=0.15m / s, and the calculated angle of attack deviation is 3.9°;

[0052] Based on the pre-calibrated correspondence between the angle of attack and the vertical wind speed component deviation, the vertical wind speed component is compensated and corrected to generate corrected observation data, specifically as follows:

[0053] Beforehand, indoor calibration experiments were conducted to obtain the vertical wind speed component deviations corresponding to different angles of attack, and a correspondence table between the angle of attack and the vertical wind speed deviation was established (as shown in Table 1). The data in the table can be adjusted according to the model of the observation instrument.

[0054] In the actual calibration process, based on the extracted angle of attack deviation, the corresponding vertical wind speed component deviation is determined by looking up a table or linear interpolation, and the rotated vertical wind speed w' is compensated and corrected. The correction logic is to subtract the corresponding vertical wind speed component deviation from the rotated vertical wind speed to obtain the corrected vertical wind speed.

[0055] The corrected vertical wind speed is combined with the rotated u' and v', while retaining the carbon dioxide concentration, water vapor concentration and atmospheric temperature data in the rotated observation data. The combined data are then used to generate the corrected observation data.

[0056] Table 1. Example of establishing the correspondence between angle of attack and vertical wind speed deviation.

[0057] Angle of attack deviation α (°) Vertical wind speed component deviation Δw (m / s) 0 0.00 1 0.02 2 0.04 3 0.07 4 0.09

[0058] In this implementation scheme, the wind speed coordinate system is adjusted to the prevailing turbulent wind direction by a quadratic coordinate rotation, and the vertical wind speed component is compensated and corrected by combining the angle of attack calibration relationship. This can reduce the impact of sensor installation angle, airflow injection angle and terrain disturbance on the vertical wind speed calculation. Since eddy flux mainly depends on the covariance between vertical wind speed fluctuation and scalar fluctuation, this processing can improve the reliability of the wind speed component used in the initial flux calculation and reduce the input error in the subsequent frequency response correction.

[0059] Specifically, the update process for the path average effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function involves the following steps:

[0060] The instrument's optical path length, sensor horizontal separation distance, and instrument dynamic response time constant are obtained, specifically as follows:

[0061] The optical path length of the instrument is the observation optical path length of the ultrasonic anemometer and the gas analyzer, which can be obtained by consulting the instrument manual or by field measurement.

[0062] The optical path length of the instrument and the horizontal separation distance of the sensor can be determined according to the specific instrument model and installation method.

[0063] In one embodiment, the instrument's optical path is 0.1m, and the sensor's horizontal separation distance is 0.15m;

[0064] The dynamic response time constant of an instrument is the delay time of the instrument's response to a signal. It is obtained through instrument calibration experiments. The dynamic response time constant of an ultrasonic anemometer is usually 0.01 to 0.1 s, and the dynamic response time constant of a gas analyzer is usually 0.1 to 1.0 s.

[0065] In one embodiment, the dynamic response time constant of the ultrasonic anemometer is set to 0.05 s, and the dynamic response time constant of the gas analyzer is set to 0.3 s.

[0066] Based on the real-time average wind speed and the instrument optical path update path average effect transfer function, the specific details are as follows:

[0067] The real-time average wind speed is extracted from the corrected observation data and is obtained by averaging the corrected longitudinal wind speed u' over time. The calculation period is consistent with the length of the sliding window.

[0068] The update logic for the path average effect transfer function is as follows:

[0069] The higher the real-time average wind speed, the smaller the high-frequency signal loss caused by the path averaging effect, and the closer the amplitude of the transfer function is to 1.

[0070] Based on real-time average wind speed With respect to the instrument's optical path S, adjust the parameters of the path-average effect transfer function so that the transfer function can adaptively match the path-average effect under the current wind speed conditions.

[0071] In one implementation, when the real-time average wind speed When the optical path length of the instrument is S=0.3m, adjust the attenuation coefficient of the transfer function to 0.92;

[0072] When the real-time average wind speed When the optical path length S of the instrument is 0.3m, the attenuation coefficient is adjusted to 0.98 to achieve dynamic updating of the transfer function;

[0073] The sensor spatial separation effect transfer function is updated based on the real-time average wind speed and the horizontal separation distance of the sensor, specifically as follows:

[0074] The signal delay caused by the spatial separation effect of the sensor is inversely proportional to the real-time average wind speed and directly proportional to the horizontal separation distance of the sensor;

[0075] Based on real-time average wind speed The horizontal separation distance d from the sensor is adjusted by adjusting the time delay parameter of the sensor's spatial separation effect transfer function. The larger the real-time average wind speed, the smaller the time delay parameter and the smaller the phase shift of the transfer function.

[0076] In one implementation, the sensor's horizontal separation distance d = 0.5m, when the real-time average wind speed... At that time, the time delay parameter is set to 0.25s;

[0077] When the real-time average wind speed At that time, the time delay parameter was set to 0.10s so that the transfer function could characterize the sensor spatial separation effect under the current wind speed conditions;

[0078] The instrument dynamic response delay transfer function is updated based on the real-time average wind speed and the instrument dynamic response time constant, specifically as follows:

[0079] The cutoff frequency of the instrument's dynamic response delay transfer function is related to the real-time average wind speed. The higher the real-time average wind speed, the higher the frequency of the high-frequency turbulence signal. Therefore, the cutoff frequency of the transfer function needs to be increased accordingly to reduce the loss of high-frequency signals caused by the instrument's dynamic response delay.

[0080] Based on the real-time average wind speed and the instrument's dynamic response time constant, the cutoff frequency of the transfer function is adjusted. The adjustment logic is that the cutoff frequency is proportional to the reciprocal of the instrument's dynamic response time constant, and is multiplied by an adjustment coefficient related to the real-time average wind speed. The larger the real-time average wind speed, the larger the adjustment coefficient.

[0081] In one implementation, the instrument's dynamic response time constant is 0.3s, and when the real-time average wind speed is 1.0m / s, the adjustment coefficient is 0.8 and the cutoff frequency is 0.42Hz.

[0082] When the real-time average wind speed is 3.0 m / s, the adjustment coefficient is 1.2 and the cutoff frequency is 0.64 Hz, so that the wind speed is adaptively updated in the instrument's dynamic response delay transfer function.

[0083] The specific steps for constructing the path average effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function are as follows:

[0084] The path averaging effect transfer function is constructed based on the correlation between instrument optical path length, real-time average wind speed, and frequency. Specifically:

[0085] The path-average effect transfer function is used to characterize the high-frequency signal attenuation caused by the instrument's optical path length. Its construction requires consideration of the instrument's optical path length S and the real-time average wind speed. The relationship between the signal frequency f and the three is as follows:

[0086] The higher the signal frequency, the longer the optical path of the instrument, and the smaller the real-time average wind speed, the more obvious the signal attenuation caused by the path averaging effect.

[0087] The transfer function can be constructed based on the instrument type, observation optical path, sensor spacing, real-time average wind speed and frequency. The specific function form can be determined according to the instrument calibration results or the preset frequency response model.

[0088] In one implementation, when the frequency f = 1 Hz, the instrument optical path S = 0.3 m, and the real-time average wind speed... At that time, the transfer function magnitude is 0.95;

[0089] When the frequency f=10Hz, the instrument optical path S=0.3m, and the real-time average wind speed... At that time, the transfer function amplitude was 0.88, consistent with the attenuation trend of high-frequency signals;

[0090] Based on the correlation between sensor horizontal separation distance, real-time average wind speed, and frequency, a sensor spatial separation effect transfer function is constructed, which is as follows:

[0091] The sensor spatial separation effect transfer function is used to characterize the signal phase shift and amplitude attenuation caused by the spatial distance between two sensors. The correlation is as follows:

[0092] The higher the signal frequency, the greater the horizontal separation distance of the sensor, and the smaller the real-time average wind speed, the more obvious the phase shift and amplitude attenuation.

[0093] The transfer function is used to characterize the relationship between the sensor's horizontal separation distance, real-time average wind speed, and frequency.

[0094] In one implementation, the sensor's horizontal separation distance d = 0.5m and the real-time average wind speed... At a frequency f=5Hz, the transfer function amplitude is 0.93 and the phase shift is 0.12rad;

[0095] When the frequency f=15Hz, the transfer function amplitude is 0.82 and the phase shift is 0.36rad, which is consistent with the trend of the spatial separation effect.

[0096] The instrument dynamic response delay transfer function is constructed based on the relationship between the instrument's dynamic response time constant, real-time average wind speed, and frequency. Specifically:

[0097] The instrument dynamic response delay transfer function is used to characterize the signal delay and attenuation caused by the instrument's own response speed. The correlation is as follows:

[0098] The higher the signal frequency, the larger the instrument's dynamic response time constant, and the smaller the real-time average wind speed, the more obvious the signal delay and attenuation.

[0099] The transfer function is constructed to characterize the interaction among the three elements;

[0100] In one implementation, the instrument's dynamic response time constant Real-time average wind speed At a frequency f=2Hz, the transfer function amplitude is 0.96;

[0101] When the frequency f=8Hz, the transfer function amplitude is 0.78, which characterizes the effect of the instrument's dynamic response delay on signals of different frequencies.

[0102] In this implementation scheme, by incorporating real-time average wind speed into the transfer function construction process of path averaging effect, sensor spatial separation effect, and instrument dynamic response delay, the three types of transfer functions can be updated with changes in current wind conditions. Compared with fixed transfer functions, this method can better fit the attenuation of high-frequency turbulence signals under different wind speed conditions, reducing the problems of insufficient compensation at low wind speeds or compensation deviating from reality at high wind speeds.

[0103] Specifically, the steps for selecting the cospectral model based on atmospheric stability parameters and calculating the frequency response correction factor by integration are as follows:

[0104] The atmospheric stability parameters, including the ratio of observation altitude to the Moning-Obukhov length, are obtained as follows:

[0105] Atmospheric stability parameters are mainly expressed as the ratio of the observation height z to the Moning-Obukhov length L. Wherein, the observation height z is the height of the observation instrument above the ground, which is obtained through on-site measurement and has a value range of 2 to 10 m;

[0106] In one implementation, the observation height z = 5m;

[0107] The Moning-Obukhov length L was calculated using corrected observation data. The calculation process is as follows:

[0108] First calculate the friction speed sensible heat flux air density Specific heat at constant pressure Atmospheric temperature Then through the formula Calculated;

[0109] In the formula, k is the von Kármán constant (value is 0.41), and g is the gravitational acceleration (value is 9.8 m / s²).

[0110] The corresponding atmospheric stability parameters are obtained through the above calculation process;

[0111] Based on the comparison results between atmospheric stability parameters and preset thresholds, a neutral cospectral model, a stable cospectral model, or an unstable cospectral model are selected, specifically as follows:

[0112] Three thresholds are preset: a neutral threshold, a stable threshold, and an unstable threshold. The neutral threshold is set to [-0.05, 0.05], the stable threshold is set to (0.05, +∞), and the unstable threshold is set to (-∞, -0.05).

[0113] The calculated atmospheric stability parameters Compared with three thresholds, when When the atmospheric condition is determined to be neutral, a neutral cospectral model is selected.

[0114] when When the atmospheric state is stable, a stable cospectral model is selected.

[0115] when When the atmospheric condition is determined to be unstable, an unstable cospectral model is selected.

[0116] In one implementation, If the value is within the neutral threshold range, a neutral cospectral model is selected.

[0117] In another implementation, If the value is within the stable threshold range, a stable cospectral model is selected.

[0118] The updated path-averaged effect transfer function, sensor spatial separation effect transfer function, and instrument dynamic response delay transfer function are multiplied together to obtain the total transfer function, which is as follows:

[0119] The total transfer function G(f) is the product of the three transfer functions, that is:

[0120] ;

[0121] In the formula, For the updated path-average effect transfer function, For the updated sensor spatial separation effect transfer function, The updated instrument dynamic response delay transfer function;

[0122] The product operation involves multiplying the amplitudes of three transfer functions at the same frequency f and adding their phases to obtain the amplitude and phase of the total transfer function, which is used to characterize the total high-frequency signal loss caused by all physical effects.

[0123] In one implementation, when the frequency f = 5Hz, , Then the magnitude of the total transfer function ;

[0124] The frequency response correction factor is determined by performing a frequency-domain weighted integral of the total transfer function and the selected cospectral model, specifically as follows:

[0125] First, determine the integration frequency range, which is the effective frequency range of the observed data, i.e., from 0.01Hz to 1 / 2 of the sampling frequency (Nyquist frequency).

[0126] In one embodiment, the sampling frequency is 10Hz, and the integration frequency range is 0.01Hz to 5Hz.

[0127] The cospectral model is a selected neutral, stable, or unstable cospectral model used to characterize the energy distribution of the turbulent signal;

[0128] The formula for calculating the weighted integral is:

[0129] ;

[0130] In the formula, This is the frequency response correction factor. For the selected cospectral model, The frequency is the lower limit of integration. The upper limit frequency for integration;

[0131] The frequency response correction factor is obtained through this integral calculation. This is used to compensate for high-frequency signal loss.

[0132] In one implementation, the integral calculation yields... This indicates that the high-frequency signal loss is approximately 19%, which needs to be compensated for by this correction factor.

[0133] In this implementation scheme, by incorporating atmospheric stability parameters into the cospectral model selection process and performing frequency-domain weighted integration of the selected cospectral model and the total transfer function, the calculation of the frequency response correction factor can simultaneously consider instrument response, sensor deployment, and atmospheric stratification. This avoids applying a single cospectral model under different stability conditions, reduces frequency response correction deviations caused by turbulence spectrum mismatch, and is particularly suitable for flux data processing under observation conditions such as low wind speed, complex terrain, and high humidity.

[0134] Specifically, the steps for determining the target filtering method based on the on-site turbulence state parameters are as follows:

[0135] Turbulence intensity and frictional wind speed were collected as on-site turbulence state parameters, specifically:

[0136] The field turbulence state parameters used for filtering determination include turbulence intensity I and friction wind speed, wherein turbulence intensity I is calculated as the ratio of the standard deviation of the corrected longitudinal wind speed fluctuation to the real-time average wind speed.

[0137] Frictional wind speed is calculated from corrected observation data and is used to characterize the intensity of turbulence. Its value range is usually 0.05 to 1.0 m / s.

[0138] In one implementation, the calculated turbulence intensity I = 0.12 and the frictional wind speed is 0.35 m / s, indicating that the on-site turbulence is relatively stable and of moderate intensity.

[0139] The turbulence intensity is compared with a preset turbulence intensity threshold, and the friction wind speed is compared with a preset friction wind speed threshold, respectively.

[0140] The preset stability threshold is set to 0.2. When the turbulence intensity I ≤ 0.2, the turbulence state is judged to be stable.

[0141] When I > 0.2, the turbulent state is determined to be unstable;

[0142] The preset frictional wind speed threshold is set to 0.2 m / s. At that time, it was determined that the turbulence intensity was relatively strong;

[0143] when When the turbulence intensity is low (low turbulence or weak turbulence), it is determined to be relatively weak.

[0144] In one implementation, the turbulence intensity I = 0.12 ≤ 0.2 is determined to be turbulent and stable;

[0145] Frictional wind speed The turbulence intensity was determined to be relatively strong.

[0146] In another implementation, the turbulence intensity I = 0.25 > 0.2, which is considered as turbulent instability;

[0147] Frictional wind speed This is determined to be a low-turbulence state;

[0148] When the on-site turbulence condition meets the steady turbulence condition, block averaging filtering is determined as the target filtering method, specifically as follows:

[0149] The conditions for steady turbulence are turbulence intensity I ≤ 0.2 and frictional wind speed. At this point, the turbulence signal has no obvious low-frequency trend term and the data fluctuation is relatively stable. Using block averaging filtering can reduce the influence of high-frequency noise while retaining the effective information of the turbulence signal.

[0150] The block length for block averaging filtering is set to 10–30 data points, which can be adjusted according to the sampling frequency.

[0151] In one implementation, the sampling frequency is 10Hz, the block length is set to 20 data points, corresponding to a block averaging time of 2 seconds, and filtering is achieved by performing an arithmetic average on the data within each block.

[0152] When the on-site turbulence condition meets the low-frequency trend condition, linear detrending filtering is determined as the target filtering method, specifically as follows:

[0153] The low-frequency trend condition is that the turbulence intensity I > 0.2 and the frictional wind speed is... At this time, the turbulence signal has obvious low-frequency trend terms (such as the trend caused by instrument drift and slow environmental changes). Linear detrending filtering can reduce the influence of low-frequency trend terms while retaining high-frequency turbulence signals.

[0154] The process of linear detrending filtering is as follows:

[0155] A linear trend line is obtained by performing linear fitting on a continuous range of observation data. The corresponding value of the linear trend line is then subtracted from the original data to obtain the detrended data.

[0156] In one implementation, 100 data points are selected as a segment, and a trend line is obtained by linear fitting. The detrended data can effectively eliminate low-frequency trend interference.

[0157] When the on-site turbulence conditions meet the criteria for low or weak turbulence, exponential filtering is determined as the target filtering method, specifically as follows:

[0158] Low or weak turbulence conditions are defined as frictional wind speeds less than 0.2 m / s (regardless of turbulence intensity). At this time, the turbulence signal is weak, and high-frequency signals are easily masked by noise. Exponential filtering can reduce the impact of abnormal fluctuations on flux calculation under weak turbulence conditions and retain the main change information in the weak turbulence signal.

[0159] The filter coefficients for the exponential filter are set to 0.1 to 0.3.

[0160] In one implementation, the filter coefficient is set to 0.2, and the filtering logic is that the current filtered data value is equal to the product of the filter coefficient and the current original data value, plus the product of (1 minus the filter coefficient) and the previous filtered data value.

[0161] When the on-site turbulence condition simultaneously meets multiple filtering conditions, the target filtering method is determined according to a preset priority, specifically as follows:

[0162] The preset priorities, from highest to lowest, are as follows:

[0163] Low-turbulence or weak-turbulence conditions > Low-frequency trend conditions > Steady-flow turbulence conditions;

[0164] That is, when the on-site turbulence conditions simultaneously meet the low turbulence condition and other conditions, exponential filtering should be preferred;

[0165] When both the low-frequency trend condition and the steady turbulence condition are met, linear detrending filter should be preferred.

[0166] In one implementation, the turbulence intensity I = 0.18 ≤ 0.2 (satisfying the stationary condition), and the frictional wind speed... (Satisfying the low turbulence condition), at this point both conditions are met simultaneously. Based on priority, exponential filtering is determined as the target filtering method.

[0167] In another implementation, the turbulence intensity I = 0.22 > 0.2, and the frictional wind speed... The on-site turbulence state was determined to meet the low-frequency trend condition, and linear detrending filtering was selected as the target filtering method.

[0168] The priority settings are shown in Table 2:

[0169] Table 2 Priority Setting Examples

[0170] Priority Corresponding filtering conditions Target filtering method 1 (highest) <0.2m / s (Low turbulence / weak turbulence) Exponential filtering 2 I>0.2 and ≥0.2m / s (low frequency trend) Linear detrending filter 3 (lowest) I≤0.2 and ≥0.2m / s (stationary turbulent flow) Block average filtering

[0171] In this implementation scheme, the on-site turbulence state is determined by turbulence intensity and frictional wind speed. The target filtering method is selected from block average filtering, linear detrending filtering and exponential filtering according to preset priority. This can avoid the processing deviation caused by fixed filtering mode. For different data states such as steady turbulence, low-frequency trend and weak turbulence, matching filtering methods are used respectively. This helps to reduce the impact of trend term, weak turbulence noise or abnormal fluctuations on flux correction results while retaining effective turbulence signal.

[0172] Specifically, the steps for applying WPL correction to the frequency response-corrected carbon dioxide flux and water vapor flux are as follows:

[0173] Extracting temperature fluctuation data, water vapor fluctuation data, and average gas concentration data, specifically as follows:

[0174] Temperature fluctuation data is the pulsating component of atmospheric temperature in the corrected observation data, which is obtained by subtracting its time average from the corrected atmospheric temperature.

[0175] Water vapor fluctuation data is the pulsating component of water vapor concentration in the corrected observation data, that is, the corrected water vapor concentration minus its time average.

[0176] The average gas concentration data includes average carbon dioxide concentration and average water vapor concentration, both of which are obtained by time averaging of the corresponding concentration data in the corrected observation data.

[0177] In one implementation, the average atmospheric temperature is 293K, and the atmospheric temperature at a certain moment is 293.2K, resulting in a temperature fluctuation of 0.2K.

[0178] The average water vapor concentration is 0.015 kg / kg, and the water vapor concentration at a certain moment is 0.016 kg / kg, resulting in a water vapor fluctuation of 0.001 kg / kg.

[0179] The average carbon dioxide concentration is 380 ppm. When converting to mass concentration or mass mixing ratio, it can be calculated based on the molar mass of carbon dioxide, the molar mass of air, and the air density.

[0180] The density perturbation calibration parameters are calculated based on temperature fluctuation data, water vapor fluctuation data, and average gas concentration data, specifically as follows:

[0181] The density perturbation calibration parameter is used to characterize the impact of temperature fluctuations and water vapor fluctuations on the density perturbation results of gas flux calculations. The calculation process is as follows:

[0182] First, calculate the air density based on atmospheric pressure, dry air gas constant, and virtual temperature. Atmospheric pressure is obtained through on-site observation, dry air gas constant is taken as a preset constant, and virtual temperature is determined based on atmospheric temperature and water vapor content.

[0183] Based on air density, average carbon dioxide concentration, average water vapor concentration, average atmospheric temperature, temperature fluctuation data, and water vapor fluctuation data, the density perturbation calibration parameters used for WPL correction are calculated.

[0184] The density perturbation calibration parameters include temperature perturbation calibration parameters and water vapor perturbation calibration parameters. The temperature perturbation calibration parameters are used to correct the impact of air density changes caused by temperature fluctuations on carbon dioxide flux and water vapor flux, while the water vapor perturbation calibration parameters are used to correct the impact of the dilution effect caused by water vapor fluctuations on carbon dioxide flux.

[0185] Density effect correction was applied to the frequency-response corrected carbon dioxide and water vapor fluxes using density perturbation calibration parameters, specifically as follows:

[0186] The core of WPL correction is to correct the impact of density disturbances caused by temperature and water vapor fluctuations on flux. The correction formula is:

[0187] ;

[0188] In the formula, This represents the WPL-corrected carbon dioxide flux. This represents the frequency-response corrected carbon dioxide flux. It is the time average of the product of vertical wind speed fluctuations and temperature fluctuations. This is the time average of the product of vertical wind speed fluctuations and water vapor fluctuations, where... and The density perturbation coefficient is determined by the average temperature, average water vapor concentration, average carbon dioxide concentration, and air density.

[0189] The WPL correction for water vapor flux can be made based on the density disturbance term corresponding to temperature fluctuations. The correction formula is as follows:

[0190] ;

[0191] In the formula, This is the water vapor flux corrected for WPL. This is the frequency response corrected water vapor flux;

[0192] In one implementation, the frequency response corrected carbon dioxide flux Calculations yielded The WPL-corrected carbon dioxide flux ;

[0193] The initial latent heat flux is used in the frequency response correction process;

[0194] After completing the WPL correction of the water vapor flux, the final output latent heat flux is obtained by converting the WPL-corrected water vapor flux with the latent heat of vaporization.

[0195] In this implementation scheme, by further performing WPL correction on the frequency response-corrected carbon dioxide flux and water vapor flux, it is possible to compensate for air density disturbances caused by temperature fluctuations and water vapor fluctuations. This ensures that the gas flux results no longer rely solely on the original concentration fluctuation calculation values. By combining the frequency response-corrected sensible heat flux and the latent heat flux determined by the WPL-corrected water vapor flux, a complete corrected flux result can be formed, making the generation processes of carbon dioxide flux, water vapor flux, sensible heat flux, and latent heat flux correspond before and after.

[0196] Specifically, the steps for determining the flux quality level based on stability and turbulence characteristic indices are as follows:

[0197] The observation period is divided into multiple sub-periods, and the average value of each corrected flux within each sub-period is calculated, specifically as follows:

[0198] The observation period is usually set at 30 minutes, which is divided into multiple equal-length sub-periods, with the length of each sub-period set to 1 to 5 minutes.

[0199] In one implementation, the 30-minute observation period is divided into 10 sub-periods, each sub-period being 3 minutes long;

[0200] For the corrected carbon dioxide flux, water vapor flux, sensible heat flux and latent heat flux, calculate the arithmetic mean for each sub-period to obtain the flux mean for each sub-period.

[0201] For example, the corrected average values ​​of carbon dioxide flux for each sub-period within a certain observation period are 0.000022, 0.000023, 0.000021, and 0.000024 kg / (m²·s), respectively.

[0202] The stationarity index is determined based on the degree of deviation between the average value of each sub-period and the overall average value of the observation period. Specifically:

[0203] The stationarity indicator uses the stationarity index IST, and its calculation logic is as follows:

[0204] First, calculate the overall average value of each flux after correction within the observation period. Then, calculate the relative deviation between the average value of each sub-period and the overall average value. Squaring all relative deviations and taking the average value, then taking the square root, yields the IST.

[0205] The smaller the IST value, the higher the throughput quality level.

[0206] In one implementation, the overall average value of the corrected carbon dioxide flux is 0.0000225 kg / (m²·s), and the root mean square of the relative deviations of the 10 sub-periods is calculated to be IST=0.08, indicating that the data has good stationarity.

[0207] Turbulence characteristic indices are calculated based on turbulence statistical parameters during the observation period, specifically as follows:

[0208] The turbulence characteristic index (ITC) is calculated based on one or more of the following: frictional wind speed, turbulence intensity, wind speed fluctuation statistics, integral time scale, and observation height.

[0209] In one implementation, three turbulence statistical parameters—turbulence intensity I, frictional wind speed, and turbulence energy dissipation rate—are selected and calculated by weighted summation. The weighting logic involves multiplying each of the three parameters by its corresponding weighting coefficient and then summing the results. The sum of the three weighting coefficients is 1, where the three weighting coefficients are 0.3, 0.4, and 0.3, respectively.

[0210] Turbulent energy dissipation rate can be estimated based on corrected wind speed fluctuation data, wind speed spectrum, or structure function method;

[0211] The higher the ITC value, the higher the availability of the data;

[0212] In one implementation, the turbulence intensity I = 0.12, the frictional wind speed is 0.35 m / s, and the turbulence energy dissipation rate is 0.002 m² / s³, so ITC = 0.3 × 0.12 + 0.4 × 0.35 + 0.3 × 0.002 ≈ 0.176 is calculated.

[0213] The flux quality level is determined based on the joint assessment results of stability and turbulence characteristic indices, specifically as follows:

[0214] The quality levels are divided into three levels by setting the classification thresholds for the stability index IST and the turbulence characteristic index ITC:

[0215] The joint judgment rules for Level 1 (Excellent), Level 2 (Qualified), and Level 3 (Unqualified) are shown in Table 3.

[0216] Based on the calculated IST and ITC values, determine the corresponding throughput quality level by referring to the table;

[0217] In one implementation, IST=0.08, ITC=0.176, and according to the table, it is judged as Level 1 (Excellent).

[0218] In another implementation, IST=0.22, ITC=0.10, and it is judged as Level 3 (unqualified).

[0219] When the levels corresponding to the stability index and the turbulence characteristic index are inconsistent, the flux quality level shall be determined according to the lower level.

[0220] Table 3 Examples of Joint Judgment Rules

[0221] Quality grade Stationarity index IST Turbulence Characteristic Index (ITC) Level 1 (High Quality) IST≤0.10 ITC≥0.15 Level 2 (Qualified) 0.10<IST≤0.20 0.10≤ITC<0.15 Level 3 (Unqualified) IST>0.20 ITC<0.10

[0222] In this implementation scheme, by combining the stationarity index with the turbulence characteristic index for flux quality level determination, the corresponding data level can be given at the same time as the output flux calculation results. Users can not only obtain the corrected carbon dioxide flux, water vapor flux, sensible heat flux and latent heat flux, but also judge whether the data of the observation period is suitable for direct use, needs to be reviewed or should be rejected according to the quality level, thereby improving the convenience of subsequent data screening, comparative analysis and long-term observation application.

[0223] In one application verification example, several sets of high-frequency eddy flux observation data are selected as input data, including three-dimensional wind speed, carbon dioxide concentration, water vapor concentration and atmospheric temperature;

[0224] The present invention is used to preprocess the input data, perform coordinate rotation, angle of attack correction, wind speed adaptive frequency response correction, adaptive selection of filtering method, and WPL correction to obtain the correction results of water vapor flux and carbon dioxide flux, and compare them with the calculation results of EddyPro software.

[0225] The comparison results are shown in Tables 4 and 5.

[0226] Table 4 Comparison of water vapor flux results between the present invention and EddyPro

[0227] Serial Number <![CDATA[EddyPro water vapor flux before correction / mmol·m -2 ·s -1 > <![CDATA[EddyPro corrected water vapor flux / mmol·m -2 ·s -1 > <![CDATA[The water vapor flux before the correction of the present invention / mmol·m -2 ·s -1 > <![CDATA[The water vapor flux after correction of the present invention / mmol·m -2 ·s -1 > 1 0.107263 0.111644 0.105847 0.111091 2 0.0215706 0.0224067 0.022631 0.020998 3 -0.0465513 -0.04886 -0.047121 -0.049676 4 -0.0546232 -0.0569981 -0.05194 -0.050869 5 0.202858 0.211957 0.205953 0.21982

[0228] Table 5 Comparison of CO2 flux results between the method of this invention and EddyPro

[0229] Serial Number <![CDATA[EddyPro corrected carbon dioxide flux / μmol·m -2 ·s -1 > <![CDATA[EddyPro corrected carbon dioxide flux / μmol·m -2 ·s -1 > <![CDATA[Before the correction of the present invention, the carbon dioxide flux / μmol·m -2 ·s -1 > <![CDATA[The corrected carbon dioxide flux of the present invention / μmol·m -2 ·s -1 > 1 0.217026 0.225034 0.212536 0.229425 2 -0.956069 -0.992941 -1.188335 -0.996097 3 -0.300327 -0.312483 -0.175958 -0.314902 4 0.956882 1.07038 0.925196 1.070502 5 -0.0503556 -0.053302 2.194712 -0.058312

[0230] As can be seen from Table 4, the water vapor flux modified by the present invention and the water vapor flux modified by EddyPro have the same overall trend, and the positive and negative directions and the magnitude of the numerical changes are basically corresponding.

[0231] As can be seen from Table 5, the carbon dioxide flux corrected by the present invention is generally close to the carbon dioxide flux corrected by EddyPro. For some samples with large deviations before correction, after frequency response correction and WPL correction, the corrected results are close to the comparison results.

[0232] The above results demonstrate that the present invention can correct high-frequency losses and density disturbances while maintaining the flux change trend, thus meeting the data processing requirements of eddy flux observation.

[0233] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0234] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A wind speed adaptive eddy flux frequency response correction method, characterized in that, Includes the following steps: Raw high-frequency eddy flux observation data containing three-dimensional wind speed, carbon dioxide concentration, water vapor concentration and atmospheric temperature were obtained, and preprocessed observation data were obtained after outlier removal and outlier repair. The preprocessed observation data were subjected to coordinate rotation and angle of attack correction in sequence to obtain the corrected observation data, and the initial carbon dioxide flux, initial water vapor flux, initial sensible heat flux and initial latent heat flux were calculated. The real-time average wind speed of the corrected observation data is extracted, and the path average effect transfer function, sensor spatial separation effect transfer function and instrument dynamic response delay transfer function are updated based on the real-time average wind speed. The cospectral model is selected and integrated to obtain the frequency response correction factor by combining atmospheric stability parameters. Based on the on-site turbulence state parameters, the target filtering method is determined from block average filtering, linear detrending filtering and exponential filtering. The target filtering method is used to process the corrected observation data, and the frequency response correction factor is used to compensate for the high-frequency loss of each initial flux to obtain the frequency response corrected flux. The frequency response-corrected carbon dioxide flux and water vapor flux are corrected by WPL, and the corrected fluxes are formed by combining the frequency response-corrected sensible heat flux and the latent heat flux determined by the WPL-corrected water vapor flux. The quality evaluation is completed based on the stability index and the turbulence characteristic index, and the corrected fluxes and flux quality levels are output.

2. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for outlier removal and outlier repair of the original high-frequency eddy flux observation data are as follows: Calculate the mean and standard deviation of each variable within the sliding window; Observations that exceed the range of the mean plus or minus a preset multiple of the standard deviation are marked as outliers, and the outlier locations are marked as missing values; Interpolation is used to repair the missing locations caused by outlier markers, generating preprocessed observation data.

3. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for performing coordinate rotation and angle of attack correction on the preprocessed observation data are as follows: Calculate the prevailing turbulent wind direction and perform coordinate rotation to bring the mean lateral wind speed and mean vertical wind speed to zero; Extract the airflow injection angle deviation from the data after coordinate rotation; Based on the pre-calibrated correspondence between the angle of attack and the vertical wind speed component deviation, the vertical wind speed component is compensated and corrected to generate corrected observation data.

4. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for updating the path average effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function are as follows: Obtain the instrument's optical path length, sensor horizontal separation distance, and instrument dynamic response time constant; The path average effect transfer function is updated based on the real-time average wind speed and the instrument optical path length. The sensor spatial separation effect transfer function is updated based on the real-time average wind speed and the horizontal separation distance of the sensor. The instrument dynamic response delay transfer function is updated based on the real-time average wind speed and the instrument dynamic response time constant.

5. The wind speed adaptive eddy flux frequency response correction method of claim 4, wherein, The specific steps for constructing the path average effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function are as follows: A path-average effect transfer function is constructed based on the correlation between instrument optical path length, real-time average wind speed and frequency; A sensor spatial separation effect transfer function is constructed based on the correlation between sensor horizontal separation distance, real-time average wind speed and frequency. The instrument dynamic response delay transfer function is constructed based on the correlation between the instrument dynamic response time constant, real-time average wind speed and frequency.

6. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for selecting the cospectral model based on atmospheric stability parameters and calculating the frequency response correction factor by integration are as follows: Obtain atmospheric stability parameters, including the ratio of observation altitude to Moning-Obukhov length; Based on the comparison results between atmospheric stability parameters and preset thresholds, a neutral cospectral model, a stable cospectral model, or an unstable cospectral model is selected. The total transfer function is obtained by multiplying the updated path-averaged effect transfer function, the sensor spatial separation effect transfer function, and the instrument dynamic response delay transfer function. The frequency response correction factor is determined by performing a frequency-domain weighted integral of the total transfer function and the selected cospectral model.

7. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The expression for the total transfer function is as follows: ； wherein is the updated path average effect transfer function, is the updated sensor spatial separation effect transfer function, is the updated instrument dynamic response delay transfer function.

8. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for determining the target filtering method based on the on-site turbulence state parameters are as follows: Turbulence intensity and frictional wind speed were collected as on-site turbulence state parameters. Compare the turbulence intensity with the preset turbulence intensity threshold and the friction wind speed with the preset friction wind speed threshold respectively; When the on-site turbulence condition meets the steady turbulence condition, block average filtering is determined as the target filtering method. When the on-site turbulence condition meets the low-frequency trend condition, linear detrending filtering is determined as the target filtering method. When the on-site turbulence condition meets the low turbulence or weak turbulence conditions, exponential filtering is determined as the target filtering method. When the on-site turbulence condition simultaneously meets multiple filtering conditions, the target filtering method is determined according to the preset priority.

9. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for applying WPL correction to the frequency response-corrected carbon dioxide and water vapor fluxes are as follows: Extract temperature fluctuation data, water vapor fluctuation data, and average gas concentration data; Density perturbation calibration parameters are calculated based on temperature fluctuation data, water vapor fluctuation data, and average gas concentration data. Density effect correction was performed on the frequency response-corrected carbon dioxide flux and water vapor flux using density perturbation calibration parameters.

10. The wind speed adaptive eddy flux frequency response correction method of claim 1, wherein, The specific steps for determining flux quality level based on stationarity and turbulence characteristic indices are as follows: The observation period was divided into multiple sub-periods, and the average value of each flux after correction was calculated in each sub-period. The stability index is determined based on the degree of deviation between the average value of each sub-period and the overall average value of the observation period; Calculate turbulence characteristic indices based on turbulence statistical parameters during the observation period; Flux quality level is determined based on the combined judgment results of stability index and turbulence characteristic index.

Images