A water vapor profile correction method and system based on multi-source data

Abstract

The present application relates to a kind of water vapor profile correction method and system based on multi-source data, belong to atmospheric remote sensing and meteorological data processing technical field.The method includes: obtaining multi-source water vapor profile and unified to microwave radiometer time spatial scale;With sounding water vapor profile combination history statistical profile generates dynamic reference profile;When cloudy, cloud radar parameter is used to carry out cloud layer saturation correction;The deviation time-varying rate of dynamic reference profile and inversion profile in sliding time window and vertical gradient anomaly index are combined to determine suspected deviation layer;Variational fusion correction is carried out to suspected deviation layer, and incremental relaxation correction is carried out to non-deviation layer;Disturbance is applied to input profile and error variance to obtain correction set, calculate uncertainty and compare with climatological standard deviation, filter credible layer output set mean value, the rest layer is back to dynamic reference profile and is marked low confidence.The present application improves correction accuracy and reliability.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric remote sensing and meteorological data processing technology. More specifically, this invention relates to a method and system for correcting water vapor profiles based on multi-source data. Background Technology

[0002] The water vapor profile retrieved by microwave radiometers is an important observational constraint for numerical weather prediction assimilation systems. However, under strong advection transport weather conditions, the retrieved profile often exhibits spurious extreme values ​​or systematic shifts at certain altitudes, requiring correction by combining information such as dynamic baseline profiles.

[0003] Existing correction methods typically set a fixed deviation threshold to identify abnormal height layers, and then apply a uniform correction strategy to all height layers, such as using variational assimilation or weighted fusion of dynamic baseline profile errors and observation errors.

[0004] This approach is effective when the deviation evolves slowly, but when the water vapor profile deviation changes rapidly, the sign reverses, or the vertical gradient deviates significantly from the historical climatological norm, the fixed threshold is difficult to balance detection sensitivity and specificity, and is prone to causing missed detection or misjudgment of the anomalous layer.

[0005] Meanwhile, the uniform correction method across all layers does not distinguish between deviated and undeviated layers, nor does it consider the differences in water vapor variability between the boundary layer and the free troposphere, as well as the different causes of deviation in each layer. This results in insufficient correction for some deviated layers, while undeviated layers may be overcorrected, introducing additional anthropogenic disturbances.

[0006] Furthermore, existing correction processes generally only output deterministic correction solutions, lacking quantitative analysis of the propagation of input profile errors and statistical parameter errors during the correction process, and thus cannot assess the reliability of the correction results. If low-confidence correction layers are directly introduced into the assimilation system, it will lead to error accumulation and affect forecast accuracy.

[0007] Therefore, how to dynamically identify abnormal deviation layers under strong advection weather, how to implement differentiated corrections for altitude layers with different characteristics, and how to quantify correction uncertainties and control low-confidence outputs have become urgent technical problems to be solved. Summary of the Invention

[0008] To address the aforementioned technical problems, the present invention provides solutions in the following aspects.

[0009] In a first aspect, the present invention provides a water vapor profile correction method based on multi-source data, employing the following technical solution: A water vapor profile correction method based on multi-source data, comprising: Microwave radiometer data and multi-source inverted water vapor profiles were acquired, and then unified to the microwave radiometer spatial scale through unit conversion, time alignment, and spatial interpolation to serve as inverted profiles; historical statistical water vapor profiles were also acquired. Using the total water vapor volume of the entire layer obtained by balloon sounding as a constraint, a dynamic baseline profile is generated by combining the water vapor obtained by balloon sounding with historical statistical profiles; when there are clouds, cloud radar parameters are used to perform cloud saturation correction on the inversion profile to obtain an optimized inversion profile. Based on the time-varying rate of the deviation between the dynamic baseline profile and the optimized inversion profile within the sliding time window and the vertical gradient anomaly index, suspected deviation layers are jointly identified; variational fusion correction is performed on the suspected deviation layers using the dynamic baseline error variance and the observation error variance; and incremental relaxation correction is performed on the non-deviation layers by scaling the deviation amount using the relaxation coefficient. Gaussian perturbations of the error standard deviation are applied to the inversion profile and the dynamic baseline profile, and amplitude perturbations are applied to the error variance to obtain the correction set. The set standard deviation is calculated as the correction uncertainty and compared with the climatological standard deviation. The set mean is output for the confidence layer, and the remaining layers are backed to the dynamic baseline profile and marked with low confidence.

[0010] Furthermore, by time matching and linear interpolation, the inversion profile, historical statistical profile, and dynamic baseline profile are unified to a preset vertical height grid, so that water vapor at each altitude layer has the same height resolution.

[0011] Furthermore, the multi-source data includes: water vapor profiles retrieved by microwave radiometers, historical statistical water vapor profiles, water vapor profiles from balloon sounding, and cloud height and cloud water content parameters detected by cloud radar.

[0012] Furthermore, the calculation of the deviation time-varying rate includes: defining the difference between the dynamic baseline profile and the optimized inversion profile at corresponding times as the deviation amount; obtaining the deviation amount at the current time and the deviation amount at the start time of the sliding time window; scaling the difference between the two according to the duration of the sliding time window to obtain the deviation time-varying rate.

[0013] Furthermore, the calculation of the vertical gradient anomaly index includes: obtaining the climatological mean and standard deviation of the vertical gradient at each altitude layer, which are obtained in advance based on historical radiosonde samples; calculating the real-time vertical gradient according to the optimized inversion profile; and dividing the absolute value of the difference between the real-time gradient and the climatological mean by the climatological standard deviation to obtain the vertical gradient anomaly index.

[0014] Furthermore, the labeling of the vertical height layer as a suspected deviation layer includes: if the absolute value of the deviation time-varying rate exceeds a preset time-varying rate threshold and the vertical gradient anomaly index exceeds a preset gradient anomaly threshold, it is labeled as a suspected deviation layer; if the deviation amount at the current time is opposite in sign to the deviation amount at the start of the sliding time window, it is labeled as a suspected deviation layer; otherwise, it is labeled as a non-deviation layer.

[0015] Furthermore, solving for the optimal correction increment includes: constructing a cost function with the correction increment as the variable, which includes a dynamic baseline profile constraint term based on the dynamic baseline error variance and an observation constraint term based on the observation error variance; obtaining the optimal correction increment by minimizing the cost function, wherein the optimal correction increment is determined by weighted fusion based on the dynamic baseline error variance, the observation error variance, and the deviation at the current time.

[0016] Furthermore, the incremental relaxation correction for the non-biased layer includes: multiplying the deviation at the current moment by a preset relaxation coefficient to obtain the relaxation correction increment; adding the value of the optimized inversion profile at the corresponding height layer to the relaxation correction increment to obtain the corrected water vapor.

[0017] Furthermore, obtaining a set consisting of multiple correction profiles includes: setting a preset number of members; the mean of the Gaussian perturbation is zero, and the standard deviations are the error standard deviations of the inverted profile and the dynamic baseline profile, respectively; applying a perturbation of a preset magnitude to the error variance includes multiplying the dynamic baseline error variance and the observation error variance by random coefficients uniformly distributed within a preset percentage range; the preset multiplier is a preset amplification factor; the remaining layers reverting to the dynamic baseline profile and marking low confidence includes setting the water vapor of the unconfidential layer as the value of the dynamic baseline profile in the corresponding layer, and attaching a low confidence identifier; the low confidence identifier is a first quality code, a second quality code is a reliable layer, and a third quality code is a data missing layer.

[0018] Secondly, the present invention provides a water vapor profile correction system based on multi-source data, which adopts the following technical solution: A water vapor profile correction system based on multi-source data includes a processor and a memory, wherein the memory stores computer program instructions, and when the computer program instructions are executed by the processor, the above-mentioned water vapor profile correction method based on multi-source data is implemented.

[0019] The present invention has the following beneficial effects: 1. By using the time-varying rate of the deviation between the dynamic baseline profile and the inverted profile within a sliding time window, combined with the vertical gradient anomaly index based on historical sounding statistics, a dual-condition judgment of time-varying abrupt changes and morphological distortion is adopted to identify suspected deviation layers layer by layer. This can sensitively capture anomalies such as rapid changes in deviation and sign reversal under strong advection weather, solve the defects of fixed threshold methods that are prone to missed detection or misjudgment, and improve the dynamic adaptability and accuracy of anomaly layer identification.

[0020] 2. Correction based on the difference in the bias layer identifier: For the bias layer, a cost function is constructed using the error covariance to achieve variational fusion, effectively integrating dynamic benchmarks and observational information to eliminate systematic bias; For the non-bias layer, incremental relaxation is used to gradually approximate with a controlled step size to avoid overcorrection and the introduction of false structures; Combining the difference in error variance between the boundary layer and the free troposphere, the correction intensity is adaptively matched with the atmospheric vertical variability, resulting in a smoother and more reasonable correction profile.

[0021] 3. By applying statistically consistent random perturbations to the inversion profile, dynamic baseline profile, and error variance, a correction set is obtained. The set dispersion is used as the layer-by-layer correction uncertainty. The credible layer is screened by comparing with historical climatological variability and the set mean is output. The uncredible layer is backed to the dynamic baseline profile and a low confidence level is added. This achieves explicit quantification and propagation blocking of correction errors, improving the reliability of water vapor profile applications.

[0022] 4. This invention integrates multiple observation data, among which microwave radiometer and GNSS / MET have the ability to work continuously in all weather conditions, lidar and cloud radar provide supplementary constraints under their respective applicable conditions (night or when clouds are present), and combined with balloon radiosonde statistical constraints, it extends to periods without radiosonde through sliding time windows and adaptive weights, so that high-quality water vapor profile correction results can be output under all weather conditions (day and night, sunny and cloudy), which is superior to the time-limited limitations of single sensors. Attached Figure Description

[0023] Figure 1 The schematic diagram illustrates the steps of a water vapor profile correction method and system based on multi-source data according to the present invention. Detailed Implementation

[0024] The following will refer to the appendices in the embodiments of the present invention. Figure 1 The technical solutions in the embodiments of the present invention are clearly and completely described herein. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0025] S1: Acquire water vapor profiles from microwave radiometer inversion, historical statistical water vapor profiles, and other multi-source data.

[0026] First, acquire multi-source observation data, including: (1) Water vapor profile retrieved by microwave radiometer The brightness temperature is obtained by receiving multi-channel microwave radiation and processing it through an inversion algorithm. The time resolution is determined by the radiometer observation period, which is usually several minutes. The vertical coverage extends from near the ground to a height of 10 km. (2) Water vapor profile or backscatter ratio profile retrieved by lidar, optional input, skip the cloud-assisted identification function when missing; (3) GNSS / MET inversion of whole-layer precipitation Optional input; skips global total auxiliary constraint function when missing. (4) Cloud boundary height and cloud liquid water content profile retrieved by cloud radar. Optional input. Skip cloud interior correction function if missing. (5) Water vapor profile measured by balloon sounding and the total amount of water vapor in the entire layer Optional input; if missing, total constraints in the dynamic baseline profile cannot be generated. (6) Historical statistical water vapor profiles obtained based on historical sounding data from the same period This data reflects the long-term climatological average vertical water vapor distribution and is essential.

[0027] All the above profiles are unified to a preset vertical height grid through linear interpolation, ranging from 0m to 10km above ground. The interlayer spacing is not exactly the same at each height layer, ranging from 100-250m, so that water vapor at each height layer has the same height resolution. Among them, 10km and interlayer spacing are empirical values ​​and can be adjusted by the implementer according to the specific implementation scenario.

[0028] When a certain type of optional data is missing, subsequent steps that rely on that data are automatically skipped or downgraded, without affecting the core correction process.

[0029] In addition, historical statistical water vapor profiles are obtained by calculating the arithmetic mean of water vapor at each altitude using radiosonde observation data within the same season, month, or number of days over many consecutive years. They reflect the long-term climatological average vertical water vapor distribution, and their height resolution and number of vertical cover layers depend on the statistical sample density of the historical radiosonde data used.

[0030] After obtaining the above profiles, the inverted profiles and historical statistical profiles are unified to the same preset vertical height grid through linear interpolation.

[0031] The preset vertical height grid starts from 0m above ground and ends at 10km, with interlayer spacing... -250m, where the maximum height and interlayer spacing are preset parameters that can be adjusted by the implementer according to the specific implementation scenario.

[0032] This preprocessing ensures that each height layer has the same interlayer spacing and height index, and subsequent deviation time-varying rate, vertical gradient anomaly index, variational fusion correction, and uncertainty screening are all performed on this unified grid.

[0033] S2: Based on the time-varying rate of the deviation between the dynamic baseline profile and the inverted profile within the sliding time window, and the vertical gradient anomaly index obtained from the historical statistical profile, the suspected deviation layer identifier is obtained by joint determination.

[0034] Under strong advection transport weather conditions, the water vapor profile retrieved by microwave radiometers exhibits rapid changes in deviation, sign reversal, and vertical gradient deviation from historical climatic norms at certain altitudes.

[0035] If these phenomena are detected solely by fixed thresholds, it is easy for them to be missed or misjudged due to improper threshold settings.

[0036] To dynamically capture the aforementioned evolutionary characteristics, a time-varying deviation rate within a sliding time window and a vertical gradient anomaly index obtained based on historical statistical profiles are introduced to perform joint determination for each height layer.

[0037] Dynamic baseline profile The dynamic baseline profile is generated by combining the total water vapor volume of the entire layer obtained from balloon sounding with historical statistical profiles. The calculation method is as follows: calculate the total amount of the entire layer of the historical statistical profile. Let the scaling factor ,but .

[0038] Based on a grid that has been unified to the same height and inversion profile The deviation at each height level at each moment is defined as The unit is grams per kilogram (g / kg).

[0039] Then set a sliding time window, the duration of which is... Take an empirical value of 1 hour, in hours (h). These are preset parameters, but can also be adjusted by the implementer based on the actual data update frequency.

[0040] For the current moment and height Get the start time of the sliding time window deviation If the starting time of the sliding time window deviation If data is missing and cannot be obtained, the bias time-varying rate condition is abandoned, and the determination is made solely based on the vertical gradient anomaly index and sign reversal condition. A data missing marker and quality code are then added to the final output data for that moment. , This is represented as the third quality code.

[0041] and the deviation at the current moment And calculate the time-varying rate of the deviation accordingly. The calculation formula is: In the formula, For height Water vapor deviation time variability, in grams per kilogram per hour (g / kg·h) ).

[0042] This represents the change in deviation within the sliding time window. The ratio of the two values ​​represents the duration of the sliding time window and reflects the average rate of change of the deviation at that altitude over time.

[0043] The larger the value, the more drastic the deviation between the dynamic baseline profile and the inverted profile at that altitude level.

[0044] Meanwhile, a vertical gradient anomaly index is introduced to quantify the degree of distortion in the vertical shape of the current inverted profile.

[0045] Obtain the climatological mean of the vertical gradient at each altitude level. and climatological standard deviation Both of these statistics are measured in grams per kilogram per cubic meter (g / kg·m ... ).

[0046] They utilize historical statistical water vapor profiles The set of historical sounding samples from the same source was obtained during the initialization phase using the following method: First, the water vapor profile of each historical sounding sample is linearly interpolated to the same preset vertical height grid as in step S1, with the interlayer spacing... -250m, then for each height level Calculate the difference in water vapor between adjacent layers and divide by The vertical gradient of the sample at that layer is obtained, and finally the arithmetic mean of the vertical gradients of all samples at the same height layer is calculated. and sample standard deviation .

[0047] If the number of historical samples for a certain height layer is less than 30, then the samples from the adjacent height layer will be used. and Perform linear interpolation supplementation for the highest layer ( ) and the bottom layer ( The gradient values ​​are directly taken from the values ​​of adjacent layers. The value of 30 is an empirical value that can be adjusted by the implementer based on the specific implementation scenario. Subsequently, the real-time vertical gradient of each height layer is calculated based on the inversion profile. .

[0048] Specifically, for adjacent height layers, the difference in water vapor concentration between the two layers is divided by the interlayer spacing (100-250m) of the uniform height grid to obtain the result in g / kg·m ... Real-time vertical gradient.

[0049] For the highest level Take it from the lower layer The gradient values ​​are the same for the bottom layer; Take it from the upper layer The gradient values ​​are the same.

[0050] The vertical boundary layers (the highest and lowest layers) do not participate in the anomaly determination of the vertical gradient anomaly index; that is, they are directly marked as unbiased layers.

[0051] Next, the vertical gradient anomaly index is calculated. : In the formula, It is a dimensionless quantity that describes the current inversion profile at a certain height. The relative degree to which the vertical gradient deviates from the historical climatological average gradient.

[0052] like Then define This indicates that the historical gradient of this layer has not changed, and the current gradient is not considered an anomaly.

[0053] The absolute deviation between the real-time vertical gradient and the climatological mean is normalized by dividing by the standard deviation. A larger value indicates significant distortion of the vertical gradient of that layer.

[0054] Obtaining the time-varying rate of the deviation and gradient anomaly index Then, cloud layer identification is introduced: if both lidar and cloud radar data are available, the lidar depolarization ratio and cloud radar echo intensity are used to jointly determine each altitude layer. Is it located inside the clouds?

[0055] Specifically, if the following conditions are met simultaneously: (a) the depolarization ratio of the lidar (a) Dimensionless; (b) Measured echo intensity from cloud radar dBZ, then this layer is labeled as a cloud layer, and recorded as... ;otherwise .

[0056] If one type of data is missing, then A value of 0 indicates that no special cloud processing is performed.

[0057] for For the cloud interior height layer, the following adjustments are made: the threshold for determining the vertical gradient anomaly index remains unchanged (it will still be used in subsequent joint determinations). However, this layer needs to be recorded as a cloud layer so that its observation error variance can be included in the subsequent variational fusion correction step. Multiply by cloud magnification factor (This reflects an increase in the inversion error of the radiometer within the cloud).

[0058] Then, for each height layer Perform joint judgment to obtain layer-by-layer suspected deviation layer identifiers. The value of this identifier is 1 or 0, which corresponds to the suspected deviation layer and the non-deviation layer, respectively.

[0059] The decision rule includes two conditions; if either condition is met, the decision is made... Set to 1 otherwise set to 0.

[0060] Condition 1 is Exceeding the preset time-varying rate threshold and Exceeding the preset gradient anomaly threshold Condition two is and The opposite sign means the product is negative.

[0061] If there is at least one intermediate moment within the sliding time window Make and The signs are opposite, and the absolute value of the deviation at this intermediate moment exceeds the preset time-varying rate threshold multiplied by [a certain value]. This is also considered to satisfy condition two. This enhanced rule is used to capture cases where signs are reversed within the window but have the same sign at the beginning and end.

[0062] in, Take an empirical value of 0.5 g / kg. , The empirical value of 3.0 is used. Both are preset parameters and can be adjusted by the implementer according to the specific implementation scenario.

[0063] Global constraint steps (based on total water vapor content in the entire layer from balloon sounding): If the current observation time Compared to the most recent balloon sounding release satisfy The time frame is 1 hour, which is a preset threshold that can be adjusted by the implementer. The total amount of water vapor in the entire atmosphere is then calculated using the balloon's sounding capability. Microwave radiometer inversion of total layer And calculate the relative deviation. .

[0064] like (Preset threshold), then the above threshold will be... Temporarily reduced by 40%, that is This is used for determining the suspected deviation layer at this moment. Otherwise, keep... .

[0065] If balloon sounding data is missing or the time difference is greater than 1 hour, no adjustment will be made. Keep the default value of 0.5.

[0066] Condition 1 captures height layers with drastic time-varying deviations and significant vertical structural distortions, while condition 2 captures height layers with reversed deviation directions. Both conditions indicate the possibility of anomalies in the water vapor information of these layers.

[0067] After performing the above judgment layer by layer, a suspected deviation layer identification profile is obtained that corresponds one-to-one with the height layer. .

[0068] S3: Based on the identification, variational fusion correction is performed on the suspected deviation layer using the dynamic baseline error variance and the observation error variance to solve for the optimal correction increment.

[0069] For suspected deviation layer identification Taking the height layer with an empirical value of 1, due to strong advection weather, the inverted water vapor may show false extreme values ​​or systematic shifts. Relying solely on the inverted profile or dynamic baseline profile cannot reliably describe the true water vapor state of this layer.

[0070] Therefore, during correction, variational fusion of the two information sources—dynamic baseline error variance and observation error variance—is necessary to obtain the optimal correction increment. The dynamic baseline error variance is pre-statistically calculated. and observation error variance The units are all (g / kg)². These two variance profiles can be obtained by statistically analyzing the differences between historical radiosonde data from the same period and the dynamic baseline profile and radiometer inversion values.

[0071] The specific statistical method is as follows: Collect a sounding sample set that originates from the same source as the historical statistical profile, with no fewer than 100 samples. 100 is an empirical value and can be adjusted by the implementer according to the specific implementation scenario. For each sample time... The dynamic baseline profile of water vapor was obtained at the same time. and radiometer inversion of water vapor and sounding observations of water vapor The variance of the dynamic benchmark error is then: The variance of the observation error is If based on the height of the top of the atmospheric boundary layer m-layered statistics, then calculate the boundary layer ( ) and outside the boundary layer ( The variance value is used for all heights within the same layer.

[0072] For cloud identification The height level, and the observation error variance used in subsequent calculations. It needs to be multiplied by the cloud magnification factor. That is, the variance of the actual observation error involved in the calculation is For non-cloudy areas, the value remains unchanged.

[0073] After obtaining the error variance of each layer, for the current time... and height Take the deviation of this layer The unit is g / kg, and a system for constructing correction increments is used. Cost function This is used to express the weighted deviation between the dynamic benchmark and the observed information under different increment values, and its expression is: In the formula, The cost function is dimensionless. The correction increment to be solved is expressed in g / kg. For dynamic baseline profile constraint terms, The dynamic baseline error variance is used to penalize the difference between the correction increment and the current deviation. The smaller the dynamic baseline error variance, the greater the penalty weight.

[0074] make right The partial derivatives are zero, which gives the condition for minimizing the cost function. From this, the optimal correction increment for this height layer can be derived. The calculation formula is as follows: In the formula, This represents the optimal correction increment for this layer, expressed in g / kg.

[0075] The optimal correction increment for this layer Add to inversion profile value Above, water vapor after correction of suspected bias layer was obtained. The unit is g / kg, which is the output of this layer in this step.

[0076] S4: Incremental relaxation correction is achieved by scaling the deviation amount using a preset relaxation coefficient on the unbiased layer.

[0077] for The non-biased layer, its deviation amount No sign reversal occurred within the sliding time window and and None of them exceeded their respective thresholds at the same time, indicating that there are no abnormal features such as systematic shifts or false extrema in the inversion values ​​of this layer.

[0078] Since the deviation magnitude of such height layers is small and they do not exhibit anomalous evolution, using the same variational fusion correction as for suspected deviation layers would introduce a dynamic baseline error variance and observation error variance weighting that does not match the actual error characteristics of the layer, and in When the value is relatively large, the observation noise will be further amplified. Therefore, incremental relaxation correction is used for these layers to limit the unnecessary correction intensity.

[0079] Incremental relaxation correction introduces a preset relaxation coefficient. This is achieved by scaling the deviation. The empirical value of 0.3 is used. It is dimensionless and is a preset parameter. It can also be adjusted by the implementer according to the correction step size requirements of the non-biased layer.

[0080] The deviation of this layer at the current moment and Multiplying them together yields the relaxation correction increment. The unit is grams per kilogram (g / kg). This product operation allows the correction increment to take only a portion of the original deviation, thereby avoiding the additional jumps that might be introduced by completely eliminating the deviation while approaching the dynamic baseline profile value.

[0081] Increment of relaxation correction Add to the value of the inverted profile at the corresponding height Above, the water vapor after correction of the unbiased layer was obtained. The unit is grams per kilogram (g / kg). This output value, together with the correction results of suspected bias layers, constitutes a correction profile that includes all height layers.

[0082] In cases where suspected and unbiased layers coexist within a continuous height layer, at the boundary between the two types of layers (i.e., The layer is a suspected deviation layer. (If the boundary layer is an unbiased layer, or vice versa), to avoid discontinuous jumps in the correction results in the vertical direction, the correction results of each boundary layer are weighted and smoothed. Specifically: let the boundary layer be... (Suspected deviation layer) and (Non-biased layer), then the final output is in For variational fusion correction values, These are relaxation correction values. If there are multiple consecutive boundary layers (more than 2 layers), linear interpolation is used for full-segment smoothing. 0.75 and 0.25 are empirical values ​​that can be adjusted by the implementer according to the specific implementation scenario, but their sum must be 1.

[0083] The error variance adaptive adjustment during the non-sonarization period is used to adjust the basic error variance used in subsequent ensemble perturbation processes, without affecting the already completed S2~S4 correction results.

[0084] Define the current time Compared with the most recent effective balloon sounding time Time difference (Hour).

[0085] like Hours: No additional adjustments will be made.

[0086] like Hour: The variance of the dynamic benchmark error Multiply by a factor of 1.2 to moderately increase confidence in the dynamic baseline profile.

[0087] like Missing hourly or sounding data: Multiply by 1.5, and simultaneously include the observation error variance. Multiplying by 1.2 reduces the confidence in inversion observations and increases reliance on dynamic baseline profiles.

[0088] The time thresholds (1 hour, 3 hours) and scaling factors (1.2, 1.5, 1.2) mentioned above are all preset parameters.

[0089] When all balloon sounding data is missing, all All judgments are considered to be greater than 3 hours, and the last case will be executed directly.

[0090] S5: Apply Gaussian perturbations corresponding to their respective error standard deviations to the inversion profile and the dynamic baseline profile, and apply a perturbation of a preset amplitude to the error variance to obtain a set of multiple correction profiles.

[0091] After completing the variational fusion correction of the suspected biased layer and the incremental relaxation correction of the unbiased layer, a preliminary correction profile was obtained.

[0092] However, the inversion profile, dynamic baseline profile, and two error variance profiles used in this correction process all have uncertainties caused by radiometer observation noise, model forecast bias, and the limited statistical nature of historical samples. A single deterministic correction cannot reflect the propagation and cumulative effects of these uncertainties in the correction calculation.

[0093] Therefore, it is necessary to construct a set of disturbance profiles and disturbance parameters that can reasonably simulate the above-mentioned input and parameter uncertainties, and to re-execute the correction process for each set of disturbance inputs to obtain a set consisting of multiple correction profiles.

[0094] First, set the number of members in the collection. , The empirical value of 20 is used as the preset parameter, but it can also be adjusted by the implementer according to the stability of the required uncertainty estimate.

[0095] For each set member ,in Two independent Gaussian random perturbation sequences were obtained and applied to the inversion profile and the dynamic baseline profile, respectively.

[0096] Inversion profile standard deviation of error The values ​​were obtained in advance from the differences between historical radiosonde and radiometer inversion values ​​from the same period, and the unit is grams per kilogram (g / kg).

[0097] For each height layer The generated data follows a mean of 0 and a standard deviation of 0. Gaussian distribution of random numbers Its unit and The random number is then superimposed onto the inversion profile to obtain the inversion profile after the member perturbation. .

[0098] Correspondingly, dynamic baseline profile standard deviation of error It was also obtained in advance from the statistical differences between historical radiosonde data and dynamic baseline profiles, with the unit being grams per kilogram (g / kg).

[0099] For each height layer The result is that the sample follows a mean of 0 and a standard deviation of 0. Gaussian distribution of random numbers Its unit and Consistent with the data, and superimposed onto the dynamic baseline profile, we obtain the dynamic baseline profile after the member disturbance. .

[0100] The application of the two sets of Gaussian perturbations causes the inversion profile and dynamic baseline profile of each member to generate random fluctuations around the original profile that are consistent with the error statistics of the profile itself, thereby simulating the uncertainty of the input profile.

[0101] Based on this, due to the dynamic baseline error variance used in variational fusion correction and observation error variance Since it is derived from finite sample statistics, its values ​​also have estimation biases, so it is necessary to further apply amplitude perturbations to the two variance profiles.

[0102] For each set member Two independent uniformly distributed random numbers are obtained. and ,in, and All within the range The value is internal, dimensionless, and can be adjusted by the implementer according to the specific implementation scenario.

[0103] Order No. The variance of the dynamic reference error of the disturbance of each member is The variance of the disturbance observation error is All units are (g / kg)².

[0104] This uniform disturbance reflects the engineering experience that there is uncertainty in the estimation of error variance in practical applications. Moreover, the sample size is smaller and the uncertainty is greater at higher altitudes, so the disturbance amplitude increases with increasing height.

[0105] For example, for height m (boundary layer), uniformly distributed interval is For height m, the interval is For height m, the interval is Thus, each set member independently obtains random coefficients for each height level. and .

[0106] The inversion profile after the above disturbance Dynamic baseline profile after disturbance The variance of the dynamic reference error after disturbance and the variance of observation errors after perturbation Substituting into the consistent correction process of S3 and S4, variational fusion correction or incremental relaxation correction is implemented layer by layer to obtain the correction profile corresponding to that member. .

[0107] right After each member performs the above operations, a total of A complete set of correction profiles, which constitute the output set. .

[0108] S6: Calculate the cascaded standard deviations as the correction uncertainty.

[0109] In obtaining by The set of corrected profiles consists of a set of members. Subsequently, the corrected water vapor values ​​of different members at the same altitude layer exhibit numerical dispersion corresponding to the disturbance. This dispersion originates from random disturbances in the input profile and error variance, directly reflecting the uncertainty of the correction process at this layer. Therefore, it is necessary to calculate the set standard deviation for each layer to quantify this uncertainty.

[0110] For each height level First, calculate the set mean of this layer. The calculation formula is: In the formula, The corrected ensemble mean of water vapor in this layer, expressed in grams per kilogram (g / kg). The number of members in the set is taken as an empirical value of 20. For the first The correction values ​​of each set member at this level.

[0111] Then calculate the ensemble standard deviation of this layer based on the ensemble mean. Its expression is: In the formula, Units and The same applies, expressed in grams per kilogram (g / kg). The summation term is the sum of squares of the deviations of each member's correction from the mean of that stratum, and the denominator uses... The unbiased standard deviation is obtained, and this standard deviation is used as the value of the correction uncertainty of this layer.

[0112] While calculating the correction uncertainty, retrieve the data based on historical statistical profiles from the initialization phase. Climatic water vapor standard deviation obtained from pre-statistical analysis of sounding samples from the same source The unit is grams per kilogram (g / kg). This statistic reflects the degree of fluctuation of water vapor at various altitudes under the historical natural weather background.

[0113] Climatic water vapor standard deviation The statistical method is based on historical statistical profiles. A set of historical sounding samples from the same source (consistent with the samples used for error variance statistics) is used to linearly interpolate the water vapor profile of each sample to a preset vertical height grid. Then, for each height layer... Calculate the sample standard deviation of water vapor in this layer for all samples, which is... .

[0114] If the number of samples at a certain height layer is less than 30, then the samples from adjacent height layers will be used. Perform linear interpolation supplementation.

[0115] Introduce preset multiples The empirical value of 1.5 is taken as a preset parameter, which is dimensionless and can also be adjusted by the implementer. The confidence test threshold is calculated layer by layer. In the formula, The unit is grams per kilogram (g / kg), which represents the upper limit of the allowable correction uncertainty after being amplified by the climatological standard deviation for that layer.

[0116] The empirical value of 1.5 is used. It is dimensionless and is a preset parameter. Implementers can adjust it according to the specific implementation scenario.

[0117] For each height level ,like If the correction uncertainty of this layer does not exceed the allowable range determined by the natural variability of climatology, then this layer is determined to be a reliable layer, and the ensemble mean is directly used. This represents the final water vapor value output by this layer.

[0118] like If the correction uncertainty of that layer is too large, the correction result is considered unreliable. In this case, the layer is judged as an untrustworthy layer, the correction result of that layer is discarded, and the output value is rolled back to the value of the dynamic baseline profile at the corresponding height. The unit is grams per kilogram (g / kg), with a low confidence level marker.

[0119] Quality Label The first quality code is defined as follows: This indicates a low confidence level, an untrusted layer, and outputs a dynamic baseline profile value; the second quality code is: This indicates high confidence, a trusted layer, and the mean of the output set; the third quality code is: This indicates that data is missing or out of time, making correction impossible.

[0120] The label is output along with the water vapor profile. The layer is assigned an observation weight of 0.

[0121] After completing the above determination layer by layer, the output values ​​of all height layers are determined according to the corresponding rules to form a final water vapor profile. The profile uses the ensemble mean in the credible layers and the dynamic baseline profile value in the uncredible layers. Each layer carries a corresponding confidence level label to distinguish the credibility of the correction results.

[0122] This invention also discloses a water vapor profile correction system based on multi-source data, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement the water vapor profile correction method based on multi-source data of this invention. The system also includes other components well-known to those skilled in the art, such as a communication bus and a communication interface; their configuration and functions are known in the art and will not be described further here.

[0123] The above are merely preferred embodiments of the present invention. It should be noted that those skilled in the art can make several improvements and substitutions without departing from the technical principles of the present invention, and these improvements and substitutions should also be considered within the scope of protection of the present invention.

Claims

1. A method for correcting water vapor profiles based on multi-source data, characterized in that, include: The water vapor profile obtained from microwave radiometer and multi-source data inversion is unified to the microwave radiometer spatial scale through unit conversion, time alignment and spatial interpolation, and used as the inversion profile; historical statistical water vapor profile is also obtained. Using the total water vapor volume of the entire layer obtained by balloon sounding as a constraint, a dynamic baseline profile is generated by combining the water vapor profile obtained by balloon sounding with historical statistical profiles; when there are clouds, cloud radar parameters are used to perform cloud saturation correction on the inverted profile to obtain an optimized inverted profile. Based on the time-varying rate of the deviation between the dynamic baseline profile and the optimized inversion profile within the sliding time window and the vertical gradient anomaly index, suspected deviation layers are jointly identified; variational fusion correction is performed on the suspected deviation layers using the dynamic baseline error variance and the observation error variance; and incremental relaxation correction is performed on the non-deviation layers by scaling the deviation amount using the relaxation coefficient. Gaussian perturbations of the error standard deviation are applied to the inversion profile and the dynamic baseline profile, and amplitude perturbations are applied to the error variance to obtain the correction set. The set standard deviation is calculated as the correction uncertainty and compared with the climatological standard deviation. The set mean is output for the confidence layer, and the remaining layers are backed to the dynamic baseline profile and marked with low confidence.

2. The water vapor profile correction method based on multi-source data according to claim 1, characterized in that, By using time matching and linear interpolation, the inversion profile, historical statistical profile, and dynamic baseline profile are unified to a preset vertical height grid, so that water vapor at each altitude layer has the same height resolution.

3. The water vapor profile correction method based on multi-source data according to claim 1, characterized in that, The multi-source data includes: water vapor profiles retrieved from microwave radiometers, historical statistical water vapor profiles, water vapor profiles from balloon sounding, and cloud height and cloud water content parameters detected by cloud radar.

4. The method according to claim 1, characterized in that, The calculation of the time-varying rate of the deviation includes: The difference between the dynamic baseline profile and the optimized inversion profile at the corresponding time is defined as the deviation. Obtain the deviation at the current moment and the deviation at the start of the sliding time window; The time-varying rate of the deviation is obtained by scaling the difference between the two values ​​according to the duration of the sliding time window.

5. The method according to claim 1, characterized in that, The calculation of the vertical gradient anomaly index includes: Obtain the mean and standard deviation of the vertical gradient climatology at each altitude layer, which are pre-calculated based on historical sounding samples. Calculate the real-time vertical gradient based on the optimized inversion profile; The vertical gradient anomaly index is obtained by dividing the absolute value of the difference between the real-time gradient and the climatological mean by the climatological standard deviation.

6. The method according to claim 1, characterized in that, Marking suspected deviation layers for vertical height layers includes: If the absolute value of the deviation time-varying rate exceeds the preset time-varying rate threshold and the vertical gradient anomaly index exceeds the preset gradient anomaly threshold, it is marked as a suspected deviation layer. If the deviation at the current time is opposite in sign to the deviation at the start of the sliding time window, it is marked as a suspected deviation layer; otherwise, it is marked as a non-deviation layer.

7. The method according to claim 1, characterized in that, Solving for the optimal correction increment includes: Construct a cost function with the correction increment as the variable, which includes a dynamic baseline profile constraint term based on the dynamic baseline error variance and an observation constraint term based on the observation error variance; The optimal correction increment is obtained by minimizing the cost function. The optimal correction increment is determined by weighted fusion based on the dynamic baseline error variance, the observation error variance, and the current time deviation.

8. The method according to claim 1, characterized in that, Incremental relaxation correction for unbiased layers includes: Multiply the deviation at the current moment by the preset relaxation coefficient to obtain the relaxation correction increment; The value of the optimized inversion profile at the corresponding height layer is added to the relaxation correction increment to obtain the corrected water vapor.

9. The method according to claim 1, characterized in that, The set of multiple corrected profiles includes: Set the preset number of members; The mean of the Gaussian perturbation is 0, and the standard deviations are the error standard deviations of the inverted profile and the dynamic baseline profile, respectively. Applying a perturbation of a preset magnitude to the error variance includes multiplying the dynamic baseline error variance and the observation error variance by random coefficients that are uniformly distributed within a preset percentage range, respectively. The preset magnification factor is a preset magnification factor; the regression of the remaining layers to the dynamic reference profile and marking low confidence includes setting the water vapor of the unconfidential layer to the value of the dynamic reference profile in the corresponding layer and adding a low confidence mark; The low confidence level identifier is the first quality code, the second quality code is the trusted layer, and the third quality code is the data missing layer.

10. A water vapor profile correction system based on multi-source data, characterized in that, include: A processor and a memory, the memory storing computer program instructions that, when executed by the processor, implement a water vapor profile correction method based on multi-source data according to any one of claims 1-9.

Images