Kilometer-level hour-by-hour multi-element forecasting method and system for meteorological high-altitude field generation

Abstract

The invention provides a kilometer-level hour-by-hour multi-element forecasting method and system for meteorological high-altitude field generation, and relates to the technical field of meteorological forecasting, and the method comprises the steps: obtaining multi-source high-altitude meteorological observation data and numerical mode forecasting output data, and carrying out the time-space alignment processing; constructing a multi-dimensional feature space based on a meteorological element physical constraint relationship to recognize uncertainty distribution features; giving a dynamic weight by using a probability forecast generation model and carrying out probability density estimation; and reliability evaluation and uncertainty quantification are carried out in combination with a historical similar weather process. According to the invention, the precision and temporal-spatial resolution of high-altitude meteorological field forecasting are improved, and the reliability of forecasting is enhanced.

Description

Technical Field

[0001] This invention relates to meteorological forecasting technology, and more particularly to a method and system for generating kilometer-level hourly multi-element forecasts of upper-air meteorological fields. Background Technology

[0002] Upper-air meteorological field forecasting has significant applications in weather forecasting, aviation transportation, and national defense. Traditional upper-air meteorological field forecasting mainly relies on numerical weather prediction models (NWPs), which simulate the future state of the atmosphere by solving atmospheric motion equations. With the improvement of computing power and the development of observation technology, multi-source observational data such as satellite remote sensing, weather radar, and upper-air sounding provide rich data sources for upper-air field forecasting. At the same time, various numerical models, such as global models and regional models, can also output forecast products with different spatiotemporal resolutions.

[0003] The spatiotemporal resolution of multi-source observational data and numerical model forecast data is inconsistent. Differences in acquisition frequency, coverage, and vertical hierarchy among data sources make data fusion difficult, hindering the effective integration of the advantages of multi-source data for high-precision forecasts. Traditional deterministic forecasting methods struggle to accurately characterize forecast uncertainties, especially in complex weather systems and extreme weather events, where forecast results lack reliable probabilistic information, failing to provide users with risk assessment and decision-making references. Existing forecasting methods generally lack sufficient consideration of the physical constraints between meteorological elements, leading to inconsistencies in forecast results across multiple elements, affecting the physical rationality and practicality of the forecasts.

[0004] As meteorological services develop towards greater precision and intelligence, users are placing higher demands on the accuracy, spatiotemporal resolution, and reliability of upper-air meteorological field forecasts. This necessitates the development of new technologies and methods to provide kilometer-level hourly multi-element upper-air field forecast products with reasonable uncertainty quantification capabilities, thereby providing more accurate meteorological information services for various application scenarios. Summary of the Invention

[0005] This invention provides a method and system for generating kilometer-level hourly multi-element forecasts of meteorological upper-air fields, which can solve the problems in the prior art.

[0006] A first aspect of the present invention provides a method for kilometer-level hourly multi-element forecasting of meteorological upper-air fields, comprising: Acquire multi-source upper-air meteorological observation data and numerical model forecast output data, perform spatiotemporal alignment processing on the multi-source upper-air meteorological observation data and the numerical model forecast output data, construct a multi-dimensional feature space based on the physical constraint relationship of meteorological elements through the spatiotemporal aligned data, and identify the uncertainty distribution characteristics of each set member in the multi-dimensional feature space through error covariance analysis between set members, and generate a set of set samples with confidence level labels; The set of samples is input into the probability forecast generation model. The probability forecast generation model assigns corresponding dynamic weights to each set member based on the confidence index. At the same time, it estimates the probability density of the set sample with dynamic weights in combination with the spatial continuity constraint of the meteorological field, and obtains a probability forecast result containing multiple probability quantiles. Based on the spatiotemporal evolution characteristics of each meteorological element in the probabilistic forecast results, and combined with similar weather process cases in the historical database, the reliability of the probabilistic forecast results is assessed and the uncertainty is quantified, generating a probability distribution product containing the forecast value and its probability of occurrence, providing decision support information for forecast services.

[0007] The multi-source upper-air meteorological observation data and the numerical model forecast output data are spatiotemporally aligned. A multidimensional feature space is constructed based on the physical constraints of meteorological elements using the spatiotemporally aligned data, including: The time dimension of the multi-source upper-air meteorological observation data and the numerical model forecast output data is standardized and mapped to a unified time reference. The spatial dimension is then resampled in a gridded manner, and the data with different spatial resolutions are projected onto a unified spatial grid to obtain spatiotemporally aligned data. Based on the spatiotemporally aligned data, a cross-element correlation feature matrix is ​​constructed according to the thermodynamic balance and dynamic conservation relationships between meteorological elements. The coupling constraint relationship between pressure gradient and wind field is used to characterize the physical consistency features between vertical levels. A multi-dimensional feature space is constructed by fusing the cross-element correlation feature matrix and the physical consistency features.

[0008] In the multidimensional feature space, the uncertainty distribution characteristics of each set member are identified through error covariance analysis among set members, generating a set of set samples with confidence level labels, including: Based on the feature vectors of each set member in the multidimensional feature space, an error covariance matrix is ​​constructed between set members in each feature dimension. The error covariance matrix is ​​then optimized by hierarchical feature decomposition. The dominant feature vectors and corresponding feature value spectra obtained by the decomposition are used to construct an error subspace with a hierarchical structure. In the error subspace, the projection components of the error vector of each set member on each dominant eigenvector are calculated based on a nonlinear projection algorithm. The prediction deviation of the set members in the multidimensional feature space is analyzed in combination with the energy distribution characteristics of the eigenvalue spectrum, thus characterizing the uncertainty distribution characteristics of the set members. Based on the aforementioned uncertainty distribution characteristics, the nonlinear projection algorithm is used to diagnose the Mahalanobis distance between each set member and the center of the overall error distribution of the set. The condition number of the error covariance matrix is ​​then fused to reconstruct the Mahalanobis distance. The reconstructed distance values ​​are mapped to confidence indicators that characterize the reliability of the forecast, and finally, a set of set samples with confidence indicators is generated.

[0009] The projection components of the error vector of each set member onto each dominant eigenvector are calculated based on a nonlinear projection algorithm. The prediction deviation of the set members in the multidimensional feature space is analyzed by combining the energy distribution characteristics of the eigenvalue spectrum. Obtain the error vector of each set member in the multidimensional feature space, and calculate the orthogonal projection components of the error vector on each dominant feature vector through a nonlinear projection algorithm. The nonlinear projection algorithm calculates the projection basis vector by introducing a random perturbation term that follows a Gaussian distribution, optimizes the projection basis vector to determine the optimal projection direction, and obtains the spatiotemporal evolution characteristics of the orthogonal projection components and the corresponding projection weight coefficients based on the optimal projection direction. Based on the spatiotemporal evolution characteristics of the orthogonal projection components and the projection weight coefficients, a weight calculation method incorporating variational entropy is constructed. The weight calculation method is then convolved with the spatiotemporal evolution characteristics, and the convolution result is adjusted using the projection weight coefficients to generate a distribution curve reflecting system stability. The spatiotemporal scale of the prediction deviation is identified by analyzing the local feature points of the distribution curve using the nonlinear projection algorithm, and the degree of prediction deviation of the set members in the multidimensional feature space is quantitatively evaluated.

[0010] The ensemble sample set is input into the probabilistic forecast generation model. The model assigns corresponding dynamic weights to each ensemble member based on the confidence level identifier. Simultaneously, it estimates the probability density of the dynamically weighted ensemble sample set in conjunction with the spatial continuity constraint of the meteorological field, yielding a probabilistic forecast result containing multiple probability quantiles, including: The set of samples is input into the probability prediction generation model, which calculates the dynamic weight of each set member in the probability density estimation based on the confidence index of each set member. A weighted kernel density function is constructed based on the dynamic weights. The predicted values ​​of each set member are then fused using the dynamic weights according to the weighted kernel density function to form a preliminary probability density distribution. Based on the spatial continuity constraint of the meteorological field, the preliminary probability density distribution is spatially smoothed. By introducing a spatial correlation regularization term to constrain the difference in probability density distribution between adjacent grid points, the preliminary probability density distribution satisfies the physical continuity of meteorological field evolution in the spatial dimension, resulting in a spatially consistent probability density distribution. The probability density distribution is solved by quantile calculation, and the quantile values ​​corresponding to multiple preset probability levels are calculated. Based on the multiple probability quantile values, a probability prediction result containing different confidence intervals is constructed.

[0011] Based on the spatial continuity constraint of the meteorological field, the preliminary probability density distribution is spatially smoothed. By introducing a spatial correlation regularization term to constrain the difference in probability density distribution between adjacent grid points, the preliminary probability density distribution satisfies the physical continuity of meteorological field evolution in the spatial dimension, resulting in a spatially consistent probability density distribution, including: Based on the spatial continuity constraint of the meteorological field, a spatial correlation regularization term is constructed. The spatial correlation regularization term calculates the probability density distribution gradient between the target grid point and its neighboring grid points, and sets the allowable range of gradient change based on the spatial diffusion characteristics of meteorological elements, thus obtaining a spatial constraint index with physical meaning. The spatial correlation regularization term and its corresponding spatial constraint index are introduced into the optimization objective function of the probability density distribution. By adding the penalty term corresponding to the spatial correlation regularization term to the optimization objective function, the difference in probability density distribution that violates the spatial continuity constraint is suppressed, thus forming an optimization objective that takes into account both physical constraints and numerical stability. The initial probability density distribution is iteratively adjusted based on the optimization objective. In each iteration, the deviation of the current probability density distribution from the spatial continuity constraint is calculated according to the optimization objective function. The probability density distribution values ​​of each grid point are adjusted by gradient descent so that the difference in probability density distribution between adjacent grid points gradually converges to the allowable range, thus obtaining a spatially consistent probability density distribution.

[0012] Based on the spatiotemporal evolution characteristics of each meteorological element in the probabilistic forecast results, and combined with similar weather process cases in the historical database, the reliability of the probabilistic forecast results is assessed and the uncertainty is quantified, generating a probability distribution product containing the forecast value and its probability of occurrence, including: The evolution trend and spatial distribution of each meteorological element in the time series are extracted from the probability forecast results to construct a spatiotemporal evolution feature vector representing the current forecast weather process. The spatiotemporal evolution feature vector is matched with the feature vector of historical weather processes stored in the historical database. The distance between the feature vectors is calculated to filter out historical weather process cases that are similar to the current forecast weather process. Based on the deviation statistics between the actual observation results and the historical forecast results in the historical weather process cases, the reliability index of the forecast value corresponding to each probability quantile in the current probability forecast result is calculated, and the uncertainty level of each probability quantile is quantified by analyzing the distribution characteristics of the deviation statistics at different probability quantiles. The multiple probability quantiles, the forecast values ​​corresponding to each probability quantile, the reliability index, and the uncertainty level in the current probability forecast result are correlated and organized to generate a probability distribution product containing the forecast values ​​and their probabilities of occurrence.

[0013] A second aspect of the present invention provides a kilometer-level hourly multi-element forecasting system for generating meteorological upper-air fields, comprising: The first unit is used to acquire multi-source upper-air meteorological observation data and numerical model forecast output data, perform spatiotemporal alignment processing on the multi-source upper-air meteorological observation data and the numerical model forecast output data, construct a multi-dimensional feature space based on the physical constraint relationship of meteorological elements through the spatiotemporal aligned data, and identify the uncertainty distribution characteristics of each set member through error covariance analysis among set members in the multi-dimensional feature space, and generate a set of set samples with confidence level labels. The second unit is used to input the set of aggregate samples into the probability forecast generation model. The probability forecast generation model assigns corresponding dynamic weights to each aggregate member based on the confidence index, and at the same time, combines the spatial continuity constraint of the meteorological field to estimate the probability density of the aggregate samples with dynamic weights to obtain a probability forecast result containing multiple probability quantile values. The third unit is used to assess the reliability and quantify the uncertainty of the probability forecast results based on the spatiotemporal evolution characteristics of each meteorological element in the probability forecast results and in combination with similar weather process cases in the historical database, and to generate a probability distribution product containing the forecast value and its probability of occurrence, so as to provide decision support information for forecast services.

[0014] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0015] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0016] The beneficial effects of this application are as follows: By performing spatiotemporal alignment processing on multi-source upper-air meteorological observation data and numerical model forecast data, and constructing a multi-dimensional feature space based on the physical constraints of meteorological elements, the effective fusion of heterogeneous data was achieved, improving the integrity and consistency of the data.

[0017] Innovatively, the uncertainty distribution characteristics of each ensemble member are identified through error covariance analysis among ensemble members, and an ensemble sample set with confidence level labels is generated, which effectively quantifies the uncertainty of the forecast and improves the reliability of the forecast.

[0018] The probability forecast generation model of this invention assigns dynamic weights to each ensemble member based on confidence level labels, which breaks through the limitations of traditional equal-weighted ensemble forecasts and more accurately reflects the differences in forecasting capabilities among different ensemble members. At the same time, combined with the spatial continuity constraint of the meteorological field, it ensures the physical rationality of the forecast results.

[0019] By combining similar weather event cases from historical databases to assess the reliability and quantify the uncertainty of probabilistic forecast results, a probability distribution product containing forecast values ​​and their probabilities of occurrence is generated, providing users with more comprehensive decision support information and significantly improving the accuracy and practicality of hourly multi-element forecasts at the kilometer level for upper-air meteorological fields. Attached Figure Description

[0020] Figure 1 This is a flowchart illustrating the method for generating kilometer-level hourly multi-element forecasts of meteorological upper-air fields according to an embodiment of the present invention. Figure 2 This is a flowchart of the probability density estimation spatial distribution consistency processing in an embodiment of the present invention; Figure 3 This is a schematic diagram of the structure of the deep learning weather forecasting model in an embodiment of the present invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. 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.

[0022] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0023] Figure 1 This is a flowchart illustrating the method for generating kilometer-level hourly multi-element forecasts of the upper-air meteorological field according to an embodiment of the present invention. Figure 1 As shown, the method includes: Acquire multi-source upper-air meteorological observation data and numerical model forecast output data, perform spatiotemporal alignment processing on the multi-source upper-air meteorological observation data and the numerical model forecast output data, construct a multi-dimensional feature space based on the physical constraint relationship of meteorological elements through the spatiotemporal aligned data, and identify the uncertainty distribution characteristics of each set member in the multi-dimensional feature space through error covariance analysis between set members, and generate a set of set samples with confidence level labels; The set of samples is input into the probability forecast generation model. The probability forecast generation model assigns corresponding dynamic weights to each set member based on the confidence index. At the same time, it estimates the probability density of the set sample with dynamic weights in combination with the spatial continuity constraint of the meteorological field, and obtains a probability forecast result containing multiple probability quantiles. Based on the spatiotemporal evolution characteristics of each meteorological element in the probabilistic forecast results, and combined with similar weather process cases in the historical database, the reliability of the probabilistic forecast results is assessed and the uncertainty is quantified, generating a probability distribution product containing the forecast value and its probability of occurrence, providing decision support information for forecast services.

[0024] In one optional implementation, the multi-source upper-air meteorological observation data and the numerical model forecast output data undergo spatiotemporal alignment processing. A multidimensional feature space is constructed based on the physical constraints of meteorological elements using the spatiotemporally aligned data, including: The time dimension of the multi-source upper-air meteorological observation data and the numerical model forecast output data is standardized and mapped to a unified time reference. The spatial dimension is then resampled in a gridded manner, and the data with different spatial resolutions are projected onto a unified spatial grid to obtain spatiotemporally aligned data. Based on the spatiotemporally aligned data, a cross-element correlation feature matrix is ​​constructed according to the thermodynamic balance and dynamic conservation relationships between meteorological elements. The coupling constraint relationship between pressure gradient and wind field is used to characterize the physical consistency features between vertical levels. A multi-dimensional feature space is constructed by fusing the cross-element correlation feature matrix and the physical consistency features.

[0025] Standardized mapping was performed on the time dimension of multi-source upper-air meteorological observation data and numerical model forecast output data. Because meteorological observation data comes from diverse sources, including radiosonde, satellite remote sensing, and weather radar, these devices have different observation intervals; for example, radiosonde data is observed twice daily, satellite remote sensing data once per hour, and numerical model forecast output data once every 6 hours. To achieve time dimension unification, a linear interpolation method was used to map data with different time resolutions to a unified time reference. Specifically, the data source with the highest time resolution was selected as the reference, and interpolation was performed on other lower time resolution data. For example, if satellite remote sensing data is observed once per hour, the radiosonde and numerical model forecast output data are mapped to the hourly time points through linear interpolation, thus achieving time dimension alignment.

[0026] The spatial dimension of the data is resampled using a grid. Due to differences in spatial coverage and resolution among different data sources—such as uneven distribution of meteorological sounding stations, spatial resolutions of satellite remote sensing data ranging from several kilometers to tens of kilometers, and numerical model forecast grids of tens of kilometers—it is necessary to project these data with different spatial resolutions onto a unified spatial grid. A bilinear interpolation method is used to resample the data from each data source onto a predefined unified spatial grid, typically an equally spaced latitude and longitude grid, such as a spatial resolution of 0.1 degrees by 0.1 degrees. For point data (such as meteorological sounding station data), Kriging interpolation is used to perform spatial interpolation, generating a continuous spatial distribution field, which is then resampled onto the unified grid using bilinear interpolation. In this way, data from different sources are aligned in the spatial dimension.

[0027] After completing spatiotemporal alignment, a cross-element correlation feature matrix is ​​constructed based on the aligned data, according to the thermodynamic equilibrium and kinetic conservation relationships among meteorological elements. Meteorological elements exhibit close physical correlations, such as the thermodynamic equilibrium relationships between temperature, humidity, and air pressure, and the geostrophic equilibrium relationship between wind and pressure fields. First, based on the ideal gas law, the characteristic values ​​of the relationships between temperature, air pressure, and specific volume at each grid point are calculated, forming a thermodynamic correlation matrix. Next, based on the atmospheric hydrostatic equilibrium equation, the characteristic values ​​of the relationship between air pressure change and density at each level are calculated, forming vertical hydrostatic features. Then, based on the geostrophic wind relationship, the degree of coordination between the horizontal pressure gradient and the wind field is calculated, forming geostrophic equilibrium features. Furthermore, the saturation relationship characteristics between humidity and temperature are calculated, forming humidity thermodynamic features. Finally, these features are combined into a cross-element correlation feature matrix, where each row represents a grid point and each column represents a type of inter-element correlation feature.

[0028] This study utilizes the coupling constraint relationship between pressure gradient and wind field to characterize the physical consistency characteristics between vertical layers. The pressure gradient force and wind field in the atmosphere are closely coupled, exhibiting continuity and consistency in the vertical direction. First, the horizontal pressure gradient vector field for each layer is calculated. Second, based on the geostrophic equilibrium relationship between pressure gradient and wind field, the geostrophic wind field for each layer is calculated. Then, the difference between the actual wind field and the theoretical geostrophic wind field is calculated, forming the geostrophic wind deviation field. Next, the variation of the geostrophic wind deviation field between adjacent layers is analyzed, and the vertical coherence index is calculated. Finally, combining the relationship between the vertical rate of change of pressure and the vertical shear of the wind field, a physical consistency characteristic vector between vertical layers is constructed, which describes the physical constraint relationship of the atmospheric vertical structure.

[0029] A multi-dimensional feature space is constructed by integrating cross-factor correlation feature matrices and physical consistency features. An enhanced feature matrix is ​​formed by merging the cross-factor correlation feature matrix and the physical consistency feature vector using a feature concatenation method. Each row of this matrix represents a spatiotemporal grid point, and the columns contain all correlation and consistency features. To handle scale differences between features, the feature matrix is ​​normalized to ensure that all feature values ​​are distributed within the same numerical range. Principal component analysis is used to reduce the dimensionality of the feature matrix, retaining key feature information and reducing redundancy, resulting in a compact multi-dimensional feature space. This feature space preserves the physical constraints between meteorological elements and reflects the consistency of the atmospheric vertical structure, providing a physical basis for subsequent meteorological forecast bias correction.

[0030] In one optional implementation, identifying the uncertainty distribution characteristics of each set member in the multidimensional feature space through error covariance analysis among set members, and generating a set of set samples with confidence level labels includes: Based on the feature vectors of each set member in the multidimensional feature space, an error covariance matrix is ​​constructed between set members in each feature dimension. The error covariance matrix is ​​then optimized by hierarchical feature decomposition. The dominant feature vectors and corresponding feature value spectra obtained by the decomposition are used to construct an error subspace with a hierarchical structure. In the error subspace, the projection components of the error vector of each set member on each dominant eigenvector are calculated based on a nonlinear projection algorithm. The prediction deviation of the set members in the multidimensional feature space is analyzed in combination with the energy distribution characteristics of the eigenvalue spectrum, thus characterizing the uncertainty distribution characteristics of the set members. Based on the aforementioned uncertainty distribution characteristics, the nonlinear projection algorithm is used to diagnose the Mahalanobis distance between each set member and the center of the overall error distribution of the set. The condition number of the error covariance matrix is ​​then fused to reconstruct the Mahalanobis distance. The reconstructed distance values ​​are mapped to confidence indicators that characterize the reliability of the forecast, and finally, a set of set samples with confidence indicators is generated.

[0031] Based on the eigenvectors of each set member in the multidimensional feature space, the error covariance matrix between set members in each feature dimension is constructed by calculating the numerical deviation covariance between any two set members in the same feature dimension. The eigenvectors contain the numerical distribution of meteorological elements such as temperature, humidity, wind speed, and air pressure at 3D spatial grid points. Each set member corresponds to a high-dimensional vector structure containing the meteorological element values ​​of all grid points. The value of the element in the i-th row and j-th column of the error covariance matrix is ​​obtained by calculating the expected value of the product of the deviations of set member i and set member j across all feature dimensions. The deviation is defined as the eigenvalue of a single member minus the arithmetic mean of all members in that feature dimension. The covariance matrix is ​​a symmetric positive definite matrix, with dimensions equal to the total dimension of the feature space, typically in the range of 10,000 to 100,000. The numerical range of the matrix elements is determined according to the physical dimensions and variability of the meteorological elements; temperature-related elements are typically in the range of 0.1 to 10.0, and wind speed-related elements are in the range of 0.01 to 1.0.

[0032] Hierarchical eigenvalue decomposition optimizes the error covariance matrix using a singular value decomposition algorithm, decomposing the n-dimensional covariance matrix into a product of orthogonal eigenvector matrices, diagonal eigenvalue matrices, and transpose matrices. Eigenvalues ​​are arranged in descending order of value, and the corresponding eigenvectors form a standard orthogonal basis system. The optimization process achieves dimensionality reduction by retaining the top k dominant eigenvectors with a cumulative contribution rate of over 95%, while eliminating minor components whose eigenvalues ​​are less than 0.1% of the total eigenvalue sum. The number of dominant eigenvectors, k, is typically controlled between 10% and 20% of the original dimension n, preserving the main variation information while significantly reducing computational complexity. The eigenvalue spectrum reflects the intensity distribution of error variation along each dominant direction, and the eigenvector corresponding to the m-th eigenvalue forms a hierarchical error subspace basis. The decomposition process uses the Lanzos iterative algorithm to handle large-scale sparse matrices, with an iteration convergence accuracy set to 1e-8 and a maximum iteration count limited to 1000.

[0033] In the error subspace, a nonlinear projection algorithm is used to calculate the projection components of the error vector of each set member onto each dominant eigenvector. Linear projection coefficients are obtained by performing an inner product operation between the original high-dimensional error vector and each dominant eigenvector. The error vector is defined as the eigenvector of the i-th set member minus the average eigenvector of all members. The nonlinear projection algorithm employs a radial basis function kernel transformation method, using a Gaussian kernel function to nonlinearly map the linear projection results. The kernel parameter is set to the square root of the standard deviation of the projection components. The nonlinear projection component of the i-th member on the j-th dominant eigenvector is obtained by transforming the original linear projection coefficients using the kernel function. The absolute value of the transformed projection component reflects the degree of deviation of the set member in the corresponding dominant direction, and the sign indicates the positive or negative direction of the deviation. The numerical stability of the projection algorithm is guaranteed by a singular value truncation method, and regularization is automatically enabled when the condition number of the eigenvector exceeds 1e6.

[0034] The predicted deviation of set members in the multidimensional feature space is quantitatively assessed by weighted summation of the contributions of each projected component, based on the energy distribution characteristics of the eigenvalue spectrum. The weighting coefficient is set as the proportion of the j-th eigenvalue to the sum of the first k eigenvalues, ensuring that deviations in the dominant direction receive higher weight. The predicted deviation of the i-th set member is calculated as the square root of the sum of squares of the weighted projected components; a larger value indicates a more significant difference between the set member and the overall set distribution. The uncertainty distribution characteristics are obtained by statistically analyzing the probability distribution of the deviations of all set members, including fourth-order statistical moment characteristics such as the mean, standard deviation, skewness, and kurtosis. The distribution shape is determined using the Kolmogorov-Smirnov test to determine whether it conforms to the normal distribution assumption, with a significance level set at 0.05. When the deviation distribution significantly deviates from the normal distribution, the Box-Cox transform is used for distribution normalization, and the transform parameters are determined using the maximum likelihood estimation method.

[0035] Based on the characteristics of uncertainty distribution, a nonlinear projection algorithm is used to diagnose the Mahalanobis distance between each set member and the center of the overall error distribution of the set. The distance is measured by calculating the quadratic form of the difference between the error vector of each set member and the vector of the error distribution center under the inverse transformation of the error covariance matrix. The center of the error distribution is defined as the weighted average of the error vectors of all set members, with the weights determined by sorting the reciprocals of the root mean square errors of each member's historical 24-hour forecasts. The Mahalanobis distance of the i-th set member is calculated as the square root of the quadratic form of the product of the error vector difference and the inverse of the covariance matrix. This distance metric considers the covariance relationship between feature dimensions and the differences in the degree of variation of each dimension. The Mahalanobis distance is calculated by solving the inverse of the covariance matrix using the Choleski decomposition method, avoiding the numerical instability problem of direct matrix inversion. When the covariance matrix is ​​close to singular, the Moore-Penrose pseudo-inverse is used instead of the conventional matrix inversion operation, and the truncation threshold for the pseudo-inverse calculation is set to 1e-10.

[0036] The condition number-based reconstruction of the Mahalanobis distance from the fusion error covariance matrix achieves distance correction by introducing a condition number-related correction factor. The condition number is defined as the ratio of the largest to the smallest eigenvalue of the covariance matrix, reflecting the matrix's numerical stability and ill-conditioned nature. When the condition number exceeds 1e4, ridge regression regularization is used to correct the covariance matrix by adding a regularization term equal to 1% of the smallest eigenvalue to the diagonal elements. The corrected matrix is ​​the original matrix plus the regularized identity matrix. The condition number normalization factor is set to 1 divided by 1 plus the logarithm of the condition number, ensuring a stronger correction as the condition number increases. The reconstructed Mahalanobis distance is equal to the original Mahalanobis distance multiplied by the normalization factor. The reconstruction process also considers the local density distribution in the feature space, estimating the local density around each set member using the k-nearest neighbor method, where k is the square root of the total number of set members rounded down.

[0037] The reconstructed distance values ​​are mapped to confidence levels, characterizing forecast reliability, using a three-segment piecewise linear mapping function for probabilistic transformation. Distance values ​​within the range of 0 to the mean distance are linearly mapped to a confidence level of 90% to 75%. Distance values ​​within the range of the mean distance to the mean distance plus two standard deviations are linearly mapped to a confidence level of 75% to 60%. Distance values ​​exceeding the mean distance plus two standard deviations correspond to confidence levels below 60%. Confidence levels are expressed as percentages, with precision retained to one decimal place. The segmentation points and slope parameters of the mapping function are determined based on statistical results from 1000 historical validation data periods. Reliability chart analysis ensures that the confidence level maintains a monotonically increasing relationship with the actual forecast accuracy. An adaptive adjustment mechanism is also incorporated into the mapping function: the distance distribution parameters are re-statistically analyzed every 100 periods, updating the segmentation thresholds of the mapping function. Confidence calculation also considers seasonal variations, using independent sets of mapping parameters for different seasons.

[0038] A sample set with confidence level labels is generated to augment data by attaching a corresponding confidence level label and uncertainty quantification information to each original set member. The sample set retains all original meteorological element values, adding confidence level, Mahalanobis distance, and projection component fields as metadata attributes. The data structure uses a nested dictionary format, with the outer key being the member number and the inner values ​​containing four subfields: feature vector, confidence level, Mahalanobis distance, and projection component. The storage format of the sample set supports efficient random access and batch read operations, employing a hierarchical data format to support large-scale data compression and fast input / output operations. The confidence level label is used for dynamic weight allocation in the subsequent probabilistic forecast generation process, with higher confidence members receiving larger weight coefficients. The sample set also includes global statistical information, recording metadata such as the uncertainty distribution parameters, dominant feature vectors, and eigenvalue spectra of the entire set.

[0039] In one optional implementation, the projection components of the error vector of each set member onto each dominant eigenvector are calculated based on a nonlinear projection algorithm. The prediction deviation of the set members in the multidimensional feature space is then analyzed by combining the energy distribution characteristics of the eigenvalue spectrum. Obtain the error vector of each set member in the multidimensional feature space, and calculate the orthogonal projection components of the error vector on each dominant feature vector through a nonlinear projection algorithm. The nonlinear projection algorithm calculates the projection basis vector by introducing a random perturbation term that follows a Gaussian distribution, optimizes the projection basis vector to determine the optimal projection direction, and obtains the spatiotemporal evolution characteristics of the orthogonal projection components and the corresponding projection weight coefficients based on the optimal projection direction. Based on the spatiotemporal evolution characteristics of the orthogonal projection components and the projection weight coefficients, a weight calculation method incorporating variational entropy is constructed. The weight calculation method is then convolved with the spatiotemporal evolution characteristics, and the convolution result is adjusted using the projection weight coefficients to generate a distribution curve reflecting system stability. The spatiotemporal scale of the prediction deviation is identified by analyzing the local feature points of the distribution curve using the nonlinear projection algorithm, and the degree of prediction deviation of the set members in the multidimensional feature space is quantitatively evaluated.

[0040] Obtain the error vectors of each set member in the multidimensional feature space. In practical applications, these error vectors can represent the difference between the output of a weather forecast model and the observed values, the deviation between a financial market prediction model and the actual market trend, etc. Assuming there are n set members, and the error vector of each member has a dimension of m, an error matrix E can be constructed, where each row represents the error vector of a set member.

[0041] The orthogonal projection components of the error vector onto the dominant eigenvectors are calculated using a nonlinear projection algorithm. This algorithm first generates a random perturbation term ξ following a Gaussian distribution, with a standard deviation set to a small value between 0.05 and 0.1 to ensure moderate perturbation. This random perturbation term is then introduced into the initial projection basis vector: v_initial = v_0 + ξ, where v_0 is the initially estimated basis vector.

[0042] The projection basis vectors are optimized to determine the optimal projection direction. Gradient descent is used to iteratively optimize the projection direction. The objective function is designed to maximize the variance after projection and minimize the projection error. In each iteration, the gradient of the objective function with respect to the current projection direction is calculated, and the projection direction is updated: v_t+1 = v_t - η·∇f(v_t), where η is the learning rate, initially set to 0.01 and gradually decreasing with the number of iterations. When the change in projection direction between two consecutive iterations is less than a preset threshold (e.g., 10), the optimal projection direction is determined. -5 When the algorithm converges, or when the maximum number of iterations (e.g., 500) is reached, the optimal projection direction v_opt is obtained.

[0043] Calculate the orthogonal projection components of the error vector based on the optimal projection direction. For each error vector e_i, calculate its projection onto the optimal projection direction v_opt: p_i = (e_i·v_opt)·v_opt / ||v_opt|| 2 By analyzing the temporal variation patterns and spatial distribution characteristics of these projection components, the spatiotemporal evolution characteristics of the projection components can be obtained. Simultaneously, based on the relationship between the projection direction and the original feature space, the weighting coefficient w_j for each projection direction is calculated, reflecting the importance of different projection directions to the error interpretation.

[0044] Based on the spatiotemporal evolution characteristics of orthogonal projection components and projection weight coefficients, a weight calculation method incorporating variational entropy is constructed. Variational entropy H(X) is used to measure the uncertainty of the system and is calculated by analyzing the distribution characteristics of the projection components at different spatiotemporal scales. The weight function is constructed as: W(t, s) = α·exp(-H(t, s) / β), where t represents the time dimension, s represents the spatial dimension, and α and β are adjustment parameters that control the overall magnitude of the weights and their sensitivity to entropy changes, respectively.

[0045] The weight calculation method is convolved with spatiotemporal evolution features. An appropriate convolution kernel K(t, s) is selected, and the convolution operation is performed: F(t, s) = ∫∫W(τ, σ)·E(t-τ, s-σ)·K(τ, σ)dτdσ, where E(t, s) represents the spatiotemporal evolution features. Convolution operations can capture feature changes at local spatiotemporal scales, enhancing pattern recognition capabilities.

[0046] The convolution result is adjusted by incorporating projection weight coefficients: F_adj(t, s) = F(t, s)·∑w_j·φ_j(t, s), where φ_j(t, s) is the spatiotemporal feature function corresponding to the j-th projection direction. This step considers the relative importance of different projection directions, making the final result more accurately reflect the dynamic characteristics of the system.

[0047] A distribution curve reflecting system stability is generated. The adjusted result F_adj(t, s) is integrated along the time and spatial dimensions to obtain a one-dimensional distribution curve C(x), where x can represent variables such as forecast lead time or spatial location. The peak value, trough value, and inflection point of this curve reflect the stability characteristics of the system under different states.

[0048] By analyzing local feature points on the distribution curve using a nonlinear projection algorithm, key points on the curve are identified, including extreme points, inflection points, and abrupt change points. These points often correspond to critical moments or locations of forecast deviations. Then, the density distribution and clustering characteristics of these feature points are analyzed to identify the main spatiotemporal scales of forecast deviations. For example, a dense appearance of feature points within a certain time period indicates poor forecast stability during that period.

[0049] The degree of forecast deviation of ensemble members in the multidimensional feature space is quantitatively evaluated. Based on the performance of each ensemble member at key feature points, its deviation index is calculated: D_i = ∑w_k·|F_i(k) - F_ref(k)| / σ_k, where F_i(k) represents the value of the i-th ensemble member at the k-th feature point, F_ref(k) represents the reference value (which can be the ensemble average or observed value), σ_k is the standard deviation at that point, and w_k is the weighting coefficient. A larger deviation index indicates a higher degree of forecast deviation and lower reliability for that ensemble member.

[0050] In practical applications, ensemble members can be sorted according to the deviation index. Members with smaller deviations can be selected to build a more reliable forecast ensemble, or members with larger deviations can be specifically improved to enhance the overall forecast quality. This method can effectively identify and evaluate forecast deviations in complex systems, providing a scientific basis for subsequent forecast corrections and system optimization.

[0051] In one optional implementation, the ensemble sample set is input into a probabilistic forecast generation model. The probabilistic forecast generation model assigns corresponding dynamic weights to each ensemble member based on the confidence level identifier. Simultaneously, it performs probability density estimation on the ensemble sample with dynamic weights, incorporating the spatial continuity constraint of the meteorological field, to obtain a probabilistic forecast result containing multiple probability quantiles, including: The set of samples is input into the probability prediction generation model, which calculates the dynamic weight of each set member in the probability density estimation based on the confidence index of each set member. A weighted kernel density function is constructed based on the dynamic weights. The predicted values ​​of each set member are then fused using the dynamic weights according to the weighted kernel density function to form a preliminary probability density distribution. Based on the spatial continuity constraint of the meteorological field, the preliminary probability density distribution is spatially smoothed. By introducing a spatial correlation regularization term to constrain the difference in probability density distribution between adjacent grid points, the preliminary probability density distribution satisfies the physical continuity of meteorological field evolution in the spatial dimension, resulting in a spatially consistent probability density distribution. The probability density distribution is solved by quantile calculation, and the quantile values ​​corresponding to multiple preset probability levels are calculated. Based on the multiple probability quantile values, a probability prediction result containing different confidence intervals is constructed.

[0052] The ensemble sample set is input into the probabilistic forecast generation model, which calculates the dynamic weights of each ensemble member in the probability density estimation based on their confidence level. Specifically, assume the ensemble forecast system contains N ensemble members, each corresponding to a set of forecast values ​​and a corresponding confidence level. The forecast value of ensemble member i is denoted as Xi, and its confidence level is denoted as Ci. During the dynamic weight calculation, the original confidence levels are first normalized to ensure the total weights are equal to 1. The original confidence level label Ci is nonlinearly mapped using an exponential transformation to enhance the weight contribution of high-confidence members. For example, for ensemble member i, its dynamic weight Wi can be calculated as a normalized exponential transformation value, ensuring that the sum of weights satisfies the normalization constraint. This approach effectively enhances the influence of ensemble members with better forecasting skills in probability density estimation while suppressing the negative impact of members with poorer skills.

[0053] A weighted kernel density function is constructed based on dynamic weights. The predicted values ​​of each ensemble member are fused using dynamic weights to form a preliminary probability density distribution. In practice, a Gaussian kernel function is selected as the basic kernel function to expand the probability density of the predicted values ​​of each ensemble member. For any grid point v within the prediction target area, its preliminary probability density distribution p(v) is obtained through weighted kernel density estimation. That is, the predicted value Xi of each ensemble member is expanded into a local probability density using the kernel function and then weighted and averaged according to the corresponding dynamic weights Wi.

[0054] The choice of bandwidth parameter for the kernel function has a significant impact on the probability density estimation results. In practice, an adaptive bandwidth determination method can be adopted to dynamically adjust the bandwidth parameter according to the dispersion of the ensemble samples. For example, for regions with low dispersion, a smaller bandwidth value can be used to preserve detailed features; while for regions with high dispersion, a larger bandwidth is used to obtain a smoother probability distribution.

[0055] Based on the spatial continuity constraint of the meteorological field, the preliminary probability density distribution is spatially smoothed. Meteorological element fields typically exhibit continuous spatial variation, and the weather conditions of adjacent grid points should maintain a certain degree of consistency. To satisfy this physical constraint, a spatial correlation regularization term is introduced to constrain the differences in probability density distributions between adjacent grid points.

[0056] For a grid point v and its neighboring grid point set N(v), a spatial smoothing objective function is defined, which includes an initial probability density fidelity term and a spatial regularization term. By minimizing this objective function, an optimized probability density distribution that respects the original set of prediction information and satisfies the spatial continuity constraint is obtained. During the optimization process, an iterative solution method is used to gradually adjust the probability density distribution of each grid point until the convergence condition is met.

[0057] For temperature field forecasting, the gradient of the temperature probability distribution between adjacent grid points can be set to not exceed a reasonable threshold. This threshold can be dynamically adjusted based on factors such as terrain features and seasonal changes. In plain areas, spatial continuity constraints can be stronger; while in complex terrain areas, the constraint strength can be appropriately relaxed.

[0058] The probability density distribution is solved by quantile calculation, and the quantile values ​​corresponding to multiple preset probability levels are calculated to form a probability forecast result containing different confidence intervals. Specifically, for the preset probability level set P={p1, p2, ..., pm}, the corresponding quantile values ​​Q={q1, q2, ..., qm} are solved, where qi represents the forecast variable value when the cumulative probability distribution function reaches level pi.

[0059] In practical applications, commonly used probability levels include 10%, 25%, 50%, 75%, and 90%, each corresponding to the boundaries of different confidence intervals. For example, the 25% and 75% quantiles constitute the central 50% confidence interval, indicating that the forecast variable has a 50% probability of falling within this interval. These quantiles can be directly used to generate weather forecast products, providing users with comprehensive uncertainty information.

[0060] Through case studies, a 24-hour precipitation forecast for a certain region was tested. The probability forecast results generated by the above method improved the Brier score by 15% compared with the traditional equal-weighted method. The reliability map showed that the calibration effect was better, especially in the forecast of extreme weather events.

[0061] The above method achieves dynamic weighting of ensemble members by introducing confidence level indicators, effectively utilizing the differences in forecasting capabilities among different members; at the same time, combined with the spatial continuity constraint of the meteorological field, it ensures the physical consistency of probabilistic forecast results and improves the accuracy and reliability of probabilistic forecasts.

[0062] In one optional implementation, the preliminary probability density distribution is spatially smoothed according to the spatial continuity constraint of the meteorological field. By introducing a spatial correlation regularization term to constrain the difference in probability density distributions between adjacent grid points, the preliminary probability density distribution satisfies the physical continuity of meteorological field evolution in the spatial dimension, resulting in a spatially consistent probability density distribution, including: Based on the spatial continuity constraint of the meteorological field, a spatial correlation regularization term is constructed. The spatial correlation regularization term calculates the probability density distribution gradient between the target grid point and its neighboring grid points, and sets the allowable range of gradient change based on the spatial diffusion characteristics of meteorological elements, thus obtaining a spatial constraint index with physical meaning. The spatial correlation regularization term and its corresponding spatial constraint index are introduced into the optimization objective function of the probability density distribution. By adding the penalty term corresponding to the spatial correlation regularization term to the optimization objective function, the difference in probability density distribution that violates the spatial continuity constraint is suppressed, thus forming an optimization objective that takes into account both physical constraints and numerical stability. The initial probability density distribution is iteratively adjusted based on the optimization objective. In each iteration, the deviation of the current probability density distribution from the spatial continuity constraint is calculated according to the optimization objective function. The probability density distribution values ​​of each grid point are adjusted by gradient descent so that the difference in probability density distribution between adjacent grid points gradually converges to the allowable range, thus obtaining a spatially consistent probability density distribution.

[0063] The spatial correlation regularization term is constructed based on the spatial continuity constraint of the meteorological field and is implemented by analyzing the physical propagation characteristics of meteorological elements between adjacent grid points. The spatial correlation regularization term is calculated based on the probability density distribution gradient between the target grid point and its neighboring grid points. The neighborhood is defined as a 3×3 or 5×5 grid area centered on the target grid point. The probability density distribution gradient is obtained by calculating the difference between the probability density value of the target grid point and the probability density values ​​of each neighboring grid point. The difference is calculated using a Euclidean distance weighted method, with the weighting coefficient being the reciprocal of the actual physical distance between grid points. The gradient calculation considers the spatial derivatives in eight directions, including the four main directions (east, west, south, and north) and the four diagonal directions (northeast, southeast, northwest, and southwest). The gradients in each direction are numerically approximated using a central difference scheme.

[0064] The permissible range for gradient changes based on the spatial diffusion characteristics of meteorological elements is determined by analyzing the physical diffusion rate and influence range of different meteorological elements. The permissible gradient range for temperature is set at 0.5 to 2.0 degrees Celsius per kilometer, for humidity at 5% to 15% relative humidity per kilometer, for wind speed at 0.1 to 1.0 meters per second per kilometer, and for air pressure at 0.01 to 0.1 hPa per kilometer. The upper and lower limits of the permissible range are determined based on statistical analysis of historical observation data, obtained by calculating the 95% confidence interval of the gradient distribution of each meteorological element between adjacent grid points. The permissible gradient range also considers the influence of terrain complexity; in areas with drastic terrain changes, such as mountainous areas or land-sea junctions, the permissible range is expanded by 1.5 to 2.0 times. The spatial constraint index is defined as the ratio of the actual gradient value to the median value of the permissible range. A ratio greater than 1.0 indicates a violation of spatial continuity constraints; the larger the ratio, the more severe the violation.

[0065] The spatial correlation regularization term is mathematically quantified by calculating the weighted sum of squares of gradient violations at all grid points. The violation degree is calculated using the numerical value of the portion exceeding the allowable range. When the actual gradient is within the allowable range, the violation degree is 0; when it exceeds the allowable range, the violation degree is equal to the absolute value of the excess portion. The weighting coefficients are determined based on the geographical location of the grid points and the importance of meteorological elements: 0.5 for ocean areas, 1.0 for plains, 1.5 for mountains, and 2.0 for densely populated urban areas. The weighting coefficients for temperature and air pressure are 1.0, for humidity 0.8, and for wind speed 1.2. The regularization term also includes a smoothness penalty component, which measures the spatial curvature change of the probability density distribution by calculating the sum of squares of the second-order gradient.

[0066] The spatial correlation regularization term is introduced into the objective function of the probability density distribution optimization by constructing a multi-objective optimization framework. The original objective function includes a data fitting term and a model complexity control term. The data fitting term measures the consistency between the probability density distribution and the observed data, while the model complexity control term prevents overfitting. The spatial correlation regularization term is added as a third objective component to the optimization function, with the weight coefficient λ set to be an adjustable parameter between 0.1 and 0.5. The weight coefficient is determined using cross-validation, dividing the historical data into training and validation sets, and selecting the λ value that minimizes the validation set error. The objective function adopts a weighted summation form, with the sum of the weight coefficients of the three components constrained to 1.0.

[0067] The penalty term is designed as a piecewise linear function. When the gradient violation is less than a threshold, the penalty term is zero; when it exceeds the threshold, the penalty term increases proportionally to the violation severity. The threshold is set to 10% of the allowable range width to ensure that slight gradient deviations do not trigger excessive penalties. The penalty term also incorporates an adaptive mechanism, dynamically adjusting the penalty strength based on the convergence progress of the current iteration. When convergence is slow, the penalty coefficient is increased to strengthen spatial constraints; when convergence is too fast, the penalty coefficient is decreased to avoid over-smoothing. The calculation of the penalty term is implemented in parallel, dividing the grid region into blocks to improve computational efficiency, and overlapping the boundaries between blocks to ensure spatial continuity.

[0068] The initial probability density distribution is iteratively adjusted based on the optimization objective using a gradient descent algorithm for numerical solution. In each iteration, the deviation of the current probability density distribution from the spatial continuity constraint is calculated, and this deviation is measured by the gradient vector of the objective function. Automatic differentiation is employed for gradient calculation, supporting efficient gradient solving for high-dimensional probability distributions. An adaptive adjustment strategy is used for the learning rate. The initial learning rate is set to 0.01 and dynamically adjusted according to the descent trend of the objective function. The learning rate is halved when the objective function decreases by less than 1e-6 for five consecutive iterations, and increased by 1.2 times when the decrease exceeds 1e-3 for five consecutive iterations.

[0069] The gradient descent process employs momentum acceleration to improve convergence speed, with a momentum coefficient set to 0.9. The iterative update formula considers historical gradient information, and the current update direction is a weighted average of the current gradient and historical momentum. To prevent gradient explosion, a gradient clipping mechanism is introduced, performing normalization when the gradient magnitude exceeds a threshold of 10.0. The iterative process also employs an early stopping strategy, automatically terminating the iteration when the objective function improvement after 20 consecutive iterations is less than 1e-7. The probability density distribution is normalized every 10 iterations to ensure that the constraint that the probability integral equals 1.0 is always satisfied.

[0070] The convergence of probability density distribution differences between adjacent grid points is determined by calculating the statistical characteristics of the spatial gradient. The convergence criterion is set at 95% or more of the grid point gradient values ​​being within the allowable range, and the remaining 5% exceeding the allowable range by no more than 20% of its width. A multi-scale strategy is employed in the convergence process, iteratively optimizing from a coarse grid resolution, gradually refining to the target resolution. The optimization result of the coarse grid serves as the initial value for the fine grid. The multi-scale strategy comprises three levels: a coarse grid resolution four times the target resolution, a medium grid resolution twice the target resolution, and finally, refinement to the target resolution. Bicubic spline interpolation is used between each level to ensure a smooth transition in the probability density distribution.

[0071] The generation of spatially consistent probability density distributions ensures the physical plausibility of the results through post-processing steps. Post-processing includes normalization checks, non-negativity constraints, and boundary condition handling. The normalization check ensures that the integral of the probability density function at each grid point equals 1.0, with an error tolerance set to 1e-6. Non-negativity constraints are implemented using a projection method, projecting negative values ​​to 0 and then re-normalizing. Boundary condition handling employs either periodic or natural boundary conditions, selecting the appropriate boundary type based on the geographical characteristics of the forecast domain. Periodic conditions are used for oceanic boundaries, while natural boundary conditions are used for land boundaries, meaning the gradient at the boundary is 0.

[0072] The data structure design employs a hierarchical storage approach, with probability density distribution data stored as a 4-dimensional array. The dimensions are longitude, latitude, meteorological elements, and probability quantiles. The data access interface supports multiple query modes by grid point, by element, and by quantile. Memory management utilizes a block loading strategy, dividing large-scale grid data into 128×128 sub-blocks and loading them into memory for processing as needed. The computation module is implemented using multi-threaded parallelism, with the number of threads set to 80% of the CPU cores to avoid excessive contention. Lock-free data structures are used between threads, and atomic operations ensure data consistency.

[0073] In one optional implementation, based on the spatiotemporal evolution characteristics of each meteorological element in the probabilistic forecast results, and combined with similar weather process cases in historical databases, the reliability of the probabilistic forecast results is assessed and the uncertainty is quantified to generate a probability distribution product containing forecast values ​​and their probabilities of occurrence, including: The evolution trend and spatial distribution of each meteorological element in the time series are extracted from the probability forecast results to construct a spatiotemporal evolution feature vector representing the current forecast weather process. The spatiotemporal evolution feature vector is matched with the feature vector of historical weather processes stored in the historical database. The distance between the feature vectors is calculated to filter out historical weather process cases that are similar to the current forecast weather process. Based on the deviation statistics between the actual observation results and the historical forecast results in the historical weather process cases, the reliability index of the forecast value corresponding to each probability quantile in the current probability forecast result is calculated, and the uncertainty level of each probability quantile is quantified by analyzing the distribution characteristics of the deviation statistics at different probability quantiles. The multiple probability quantiles, the forecast values ​​corresponding to each probability quantile, the reliability index, and the uncertainty level in the current probability forecast result are correlated and organized to generate a probability distribution product containing the forecast values ​​and their probabilities of occurrence.

[0074] The evolution trends of various meteorological elements over time are extracted from probabilistic forecast results using a sliding window analysis method. The time series evolution trends include descriptors with four dimensions: monotonicity, periodicity, volatility, and extreme value characteristics. The monotonicity descriptor calculates the increasing or decreasing trend of values ​​over a continuous time period, using the slope of linear regression as a quantitative indicator. An absolute slope value greater than 0.1 indicates a significant trend; a positive value indicates an upward trend, and a negative value indicates a downward trend. The periodicity descriptor identifies the main frequency components using Fast Fourier Transform, extracting the amplitude and phase information of the three main periods: 24 hours, 12 hours, and 6 hours. The volatility descriptor calculates the coefficient of variation using the ratio of the standard deviation to the mean; a coefficient of variation greater than 0.3 indicates high volatility, and less than 0.1 indicates low volatility. Extreme value characteristics include four sub-features: the time of occurrence of the maximum value, the time of occurrence of the minimum value, the duration of the extreme value, and the intensity of the extreme value.

[0075] Spatial distribution morphological feature extraction employs a combination of spatial gradient analysis and morphological pattern recognition to construct descriptive vectors. Spatial gradient analysis calculates the first and second derivatives of each grid point in the two main east-west and north-south directions. The first derivative reflects the spatial rate of change, while the second derivative reflects the spatial curvature characteristics. The main features of the gradient field include four statistical quantities: the location of the gradient maximum, the distribution of gradient zeros, gradient directional divergence, and gradient curl. Morphological pattern recognition identifies spatial structural patterns such as high-value regions, low-value regions, saddle points, and ridges through connected component analysis. The area proportion, centroid, aspect ratio, and orientation angle of each pattern constitute the core elements of the morphological descriptive vector. Spatial correlation features are calculated using a variogram; the variogram values ​​at different lag distances reflect the spatial autocorrelation structure.

[0076] The spatiotemporal evolution feature vectors are constructed using a hierarchical, cascaded approach to organize features across dimensions. The temperature feature vector contains 48 temporal dimensions and 32 spatial dimensions; the humidity feature vector contains 36 temporal dimensions and 28 spatial dimensions; the wind speed feature vector contains 42 temporal dimensions and 35 spatial dimensions; and the air pressure feature vector contains 40 temporal dimensions and 30 spatial dimensions. Principal component analysis (PCA) is used to reduce the dimensionality of each feature vector, retaining principal components with a cumulative variance contribution rate of 95%. The resulting feature vectors have a dimensionality between 20 and 30. Z-score normalization is employed to ensure that each dimension has the same numerical scale. Multi-factor feature vectors are weighted and concatenated to form a comprehensive feature vector. The weighting coefficients are set according to the forecast importance of each element: temperature and air pressure have a weight of 0.3, while humidity and wind speed have a weight of 0.2.

[0077] The historical weather event feature vectors in the historical database are stored using a categorized index structure to improve query efficiency. The database contains nearly 20 years of historical weather event cases, each case including a 72-hour forecast period and corresponding actual observation results. Feature vectors are stored according to three levels of classification: weather type, season, and geographical region. Weather types include six main categories: sunny, cloudy, rain, snow, thunderstorm, and typhoon. Seasons are divided into spring, summer, autumn, and winter. Geographical regions are divided into four climate zones: tropical, subtropical, temperate, and frigid. The index structure uses a multi-dimensional B+ tree, supporting fast retrieval based on feature vector dimensions. The storage format of historical cases includes six fields: case identifier, timestamp, feature vector, forecast result, observation result, and error statistics.

[0078] Similarity matching calculation employs a hybrid metric combining weighted Euclidean distance and cosine similarity. During distance calculation, each feature dimension is assigned a different weight based on its importance in describing the weather process: temporal evolution features are weighted at 0.6, and spatial distribution features at 0.4. Euclidean distance is calculated as the square root of the weighted sum of squares of the differences between corresponding dimensions, while cosine similarity is calculated as the cosine of the angle between two feature vectors. The hybrid similarity index is defined as a weighted average of the reciprocal of the Euclidean distance and the cosine similarity, with a weight ratio of 3:7. A similarity threshold of 0.85 is set; historical cases exceeding this threshold are considered similar weather processes. The number of similar cases is controlled between 5 and 20; if the number is insufficient, the similarity threshold is lowered to 0.8, and if the number is excessive, the threshold is increased to 0.9.

[0079] The selected historical weather event cases were ranked according to their similarity scores, with the top 10 cases receiving the highest similarity serving as the primary reference for reliability assessment. Case selection also considered temporal proximity constraints, prioritizing cases that were closer in time, with the time weight using an exponential decay function and a decay coefficient of 0.95. Geographical similarity was also considered, with cases within the same climate zone given higher priority, and cross-climate zone similarity scores multiplied by a penalty factor of 0.8. Seasonal similarity was weighted at 0.9, with cases within the same season given higher priority than cross-seasonal cases. The final selection results underwent multiple validations to ensure quality, including statistical significance tests, time series stationarity tests, and spatial distribution consistency tests.

[0080] Reliability indices based on historical weather case studies are quantitatively assessed by analyzing the statistical deviations between historical forecasts and actual observations. These deviation statistics include three dimensions: systematic bias, random bias, and absolute bias. Systematic bias is calculated as the average difference between historical forecasts and observed values; a positive value indicates an overestimation, a negative value indicates an underestimation, and an absolute value greater than 10% of the observed standard deviation indicates significant systematic bias. Random bias is measured using the standard deviation of the bias sequence, reflecting the random component of forecast uncertainty. Absolute bias is calculated using the mean absolute error, providing an unbiased estimate of the bias magnitude. Reliability indices for each probability quantile are calculated using conditional bias statistics, grouping historical cases according to their forecast probability quantiles and calculating the deviation statistics for each group.

[0081] Uncertainty level quantification achieves accurate estimation by analyzing the distribution characteristics of deviation statistics at different probability quantiles. The deviation distribution characteristics are modeled using the kernel density estimation method, with a Gaussian kernel chosen as the kernel function, and the bandwidth parameter adaptively determined through cross-validation. The uncertainty level of each probability quantile is defined as the standard deviation of the corresponding deviation distribution; a larger standard deviation indicates a higher level of uncertainty. Quantile-specific uncertainty is identified by comparing the differences in deviation distributions at different quantiles, and the Kolmogorov-Smirnov test is used to determine the significance of differences between distributions. The uncertainty level of extreme quantiles is usually higher than that of middle quantiles; the uncertainty levels at the 5th and 95th percentiles are generally 1.5 to 2.0 times that of the 50th percentile.

[0082] The generation of probability distribution products employs a structured data organization method to achieve the associated storage of multi-dimensional information. The data structure includes eight dimensions of information fields: forecast time, spatial location, meteorological elements, probability quantiles, forecast value, probability of occurrence, reliability index, and uncertainty level. Forecast time is stored in ISO 8601 format, spatial location is represented using latitude and longitude coordinates, meteorological elements are encoded using standard codes, and probability quantiles are set at 19 quantiles at 5% intervals from the 5th to the 95th percentile. Forecast values ​​are retained to two decimal places, the probability of occurrence is expressed as a percentage with one decimal place, the reliability index uses dimensionless values ​​between 0 and 1, and the uncertainty level is expressed as an absolute value in the same unit as the forecast value.

[0083] The hierarchical information management system employs a nested structure. The top-level structure is organized according to forecast time and spatial location, the middle-level structure is categorized by meteorological elements, and the bottom-level structure contains detailed information for each probability quantile. Each quantile record includes five subfields: forecast value, probability of occurrence, confidence interval, reliability level, and uncertainty indicator. The confidence interval is calculated based on the uncertainty level, using the normal distribution assumption to determine the upper and lower limits of the 95% confidence interval. The reliability level is divided into three levels: high, medium, and low. A reliability index greater than 0.8 is high, 0.6 to 0.8 is medium, and less than 0.6 is low. The uncertainty indicator uses color coding: green indicates low uncertainty, yellow indicates medium uncertainty, and red indicates high uncertainty.

[0084] The data format uses JSON for cross-platform compatibility, supporting network transmission and data exchange. The probability distribution product also includes metadata information, recording auxiliary information such as generation time, data version, processing parameters, and quality identifiers. The quality identifiers include evaluation results across three quality dimensions: data integrity, timeliness, and accuracy level. The access interface supports querying and filtering based on time range, spatial region, meteorological elements, probability thresholds, and other conditions. The caching strategy uses the LRU algorithm to manage hot data, with a cache capacity of 1GB and a cache hit rate maintained above 80%.

[0085] like Figure 3 As shown, the method further includes: Large-scale upper-air meteorological field input data for the forecast start time is automatically acquired through a data interface module. This module employs a RESTful API design based on the HTTP protocol, supporting integration with the data publishing service of the Global Numerical Weather Prediction Center. The module configures connection attributes such as the data source address, authentication key, and timeout parameters, with a connection timeout set to 30 seconds, a read timeout set to 300 seconds, and a maximum of three retries. Data requests are asynchronous to avoid blocking the main thread. The input data has a 6-hour time interval and a horizontal resolution of 0.25 degrees, with the time interval precisely controlled at 21,600 seconds and the horizontal resolution corresponding to a grid spacing of approximately 28 kilometers. The data covers the latitude and longitude range of the forecast target area, with boundary extensions of at least one grid point to provide boundary condition support.

[0086] The data structure is defined using a four-dimensional array storage format, with dimensions of time, vertical hierarchy, latitude, and longitude. The time dimension includes the starting time and the first three time intervals, providing a 24-hour historical evolution background. The vertical hierarchy covers isobaric surfaces from 1000 hPa to 100 hPa, selecting eight standard isobaric surfaces: 1000, 925, 850, 700, 500, 300, 200, and 100 hPa. The latitude and longitude dimensions are determined based on the target area, and the number of grid points is controlled within 200 x 200 to balance computational efficiency and accuracy requirements.

[0087] Meteorological field variables comprising multiple vertical levels adopt standard meteorological variable coding specifications. The wind speed field includes three components: u, v, and w. The u component represents east-west wind speed, the v component represents north-south wind speed, and the w component represents vertical wind speed. Geometric height z is expressed in geopotential meters, representing the height of the isobaric surface relative to sea level. Temperature t is expressed in Kelvin, and specific humidity q is expressed in kilograms of water vapor per kilogram of dry air. Numerical range constraints are set as follows: u and v components between -100 and +100 m / s; w component between -10 and +10 Pa / s; geopotential height between 0 and 30,000 geopotential meters; temperature between 180 and 330 Kelvin; and specific humidity between 0 and 0.03.

[0088] Data quality control ensures the integrity of input data through outlier detection and missing value handling mechanisms. Outlier detection uses a 3-standard-deviation rule, marking values ​​outside a reasonable range as outliers and replacing them with interpolated values. Missing value handling prioritizes temporal interpolation, using linear interpolation to complete the missing values ​​from adjacent time intervals. Spatial interpolation is used as an alternative, employing bilinear interpolation to calculate missing value values ​​based on the values ​​of the surrounding four grid points. Data preprocessing includes normalization, standardizing each variable according to its physical dimensions, setting the mean to 0 and the standard deviation to 1.

[0089] The deep learning prediction model employs an encoder-decoder architecture to implement an end-to-end neural network structure. The encoder is responsible for extracting high-dimensional feature representations of the input data and consists of four 3D convolutional layers and two temporal attention layers. The kernel size of the 3D convolutional layers is set to 3x3x3 with a stride of 1 and uniform padding. The first layer has 64 output channels, the second layer has 128, the third layer has 256, and the fourth layer has 512. Each convolutional layer is followed by batch normalization and a ReLU activation function. The momentum parameter of the batch normalization is set to 0.9, and the epsilon parameter is set to 1e-5.

[0090] The temporal attention layer employs a multi-head attention mechanism to process time series features, with 8 attention heads and a hidden layer dimension of 512. Attention calculation involves linear transformations of the query matrix, key matrix, and value matrix, and weight initialization uses the Xavier initialization method. The dropout rate is set to 0.1 to prevent overfitting. Positional encoding uses sine and cosine encoding to generate a unique encoding vector for each position in the time dimension.

[0091] The decoder maps the encoded features to the target high-resolution prediction results and consists of four transposed convolutional layers and a multi-task output header. The transposed convolutional layers upsample the feature maps. Layer 1 has 512 input channels and 256 output channels, with an upsampling factor of 2. Layer 2 has 256 input channels and 128 output channels, with an upsampling factor of 2. Layer 3 has 128 input channels and 64 output channels, with an upsampling factor of 4. Layer 4 has 64 input channels and 32 output channels, with an upsampling factor of 2. Overall, the upsampling factor reaches 64 times, improving the resolution from 0.25 degrees to approximately 1 kilometer.

[0092] The multi-task output head sets up independent output branches for different meteorological elements. The precipitation output branch uses a 1x1 convolutional layer with one output channel and a ReLU activation function to ensure non-negative output. The temperature output branch has one output channel and a linear activation function. The humidity output branch has one output channel and a Sigmoid activation function, limiting the output range to 0 to 1. The wind speed output branch has two output channels, corresponding to the u and v components respectively, and uses a hyperbolic tangent function. The radiation output branch has one output channel and a ReLU activation function to ensure non-negative output.

[0093] The model training employs supervised learning, with a weighted multi-task loss function. Precipitation loss uses a mean squared error loss with a weighting coefficient of 2.0 to enhance precipitation forecasting capabilities. Temperature loss uses a mean squared error loss with a weighting coefficient of 1.0. Humidity loss uses a binary cross-entropy loss with a weighting coefficient of 1.5. Wind speed loss uses a vector loss function, considering both wind speed magnitude and direction errors, with a weighting coefficient of 1.0. Radiation loss uses a mean squared error loss with a weighting coefficient of 0.8. The total loss is the weighted sum of the losses for each component.

[0094] The Adam optimizer was chosen, with an initial learning rate of 0.001, beta1 parameter set to 0.9, beta2 parameter set to 0.999, and epsilon parameter set to 1e-8. Cosine annealing was used for learning rate scheduling, reducing the rate to 0.5 times its current value every 50 epochs, with a minimum learning rate of 1e-6. The batch size was set to 8, the number of training epochs was 200, and an early stopping mechanism was triggered when the validation loss showed no improvement for 20 consecutive epochs.

[0095] The inference computation process is accelerated by GPU for efficient forecasting. Model loading utilizes pre-trained weight files in PyTorch's pth format, approximately 500MB in size. Approximately 8GB of GPU memory is required for model parameter storage and intermediate result caching. The input data batch size is set to 1, and a single inference time is approximately 10 seconds. The output contains 72 time steps, each corresponding to a 1-hour interval, with a spatial resolution of 1024 x 1024 grid points, corresponding to a grid spacing of approximately 1 kilometer.

[0096] Continuous forecasting of hourly multi-meteorological element forecasts for the next 72 hours is achieved through a temporal unfolding mechanism. The temporal unfolding employs an autoregressive approach, forecasting the state one hour later each time, and then using the forecast result as input for the next time step. The initial state is derived from feature representations extracted by the encoder, containing complete information about the atmospheric state. The hidden state dimension is 512, and temporal information is transmitted through LSTM units. The parameters of the LSTM's forget gate, input gate, and output gate are learned independently.

[0097] The horizontal resolution of the forecast field is precisely controlled to a spacing of 1 kilometer using bilinear interpolation. The target grid adopts an equidistant projected coordinate system, with the origin set at the center of the target area, the x-axis pointing east, and the y-axis pointing north. The grid point spacing is precisely 1000 meters, with an error tolerance of 10 meters. Mirror boundary conditions are used for boundary processing to avoid the influence of boundary effects on the forecast results.

[0098] Output element physical units and range constraints ensure the reasonableness of the forecast results. Precipitation is measured in millimeters per hour, with a range limited to 0 to 100; values ​​outside this range are truncated. Temperature is measured in degrees Celsius, with a range limited to -50 to 50. Humidity is measured as a percentage of relative humidity, with a range limited to 0 to 100. Wind speed is measured in meters per second, with the u and v components limited to a range of -50 to 50. Radiation is measured in watts per square meter, with a range limited to 0 to 1500.

[0099] A second aspect of the present invention provides a kilometer-level hourly multi-element forecasting system for generating meteorological upper-air fields, comprising: The first unit is used to acquire multi-source upper-air meteorological observation data and numerical model forecast output data, perform spatiotemporal alignment processing on the multi-source upper-air meteorological observation data and the numerical model forecast output data, construct a multi-dimensional feature space based on the physical constraint relationship of meteorological elements through the spatiotemporal aligned data, and identify the uncertainty distribution characteristics of each set member through error covariance analysis among set members in the multi-dimensional feature space, and generate a set of set samples with confidence level labels. The second unit is used to input the set of aggregate samples into the probability forecast generation model. The probability forecast generation model assigns corresponding dynamic weights to each aggregate member based on the confidence index, and at the same time, combines the spatial continuity constraint of the meteorological field to estimate the probability density of the aggregate samples with dynamic weights to obtain a probability forecast result containing multiple probability quantile values. The third unit is used to assess the reliability and quantify the uncertainty of the probability forecast results based on the spatiotemporal evolution characteristics of each meteorological element in the probability forecast results and in combination with similar weather process cases in the historical database, and to generate a probability distribution product containing the forecast value and its probability of occurrence, so as to provide decision support information for forecast services.

[0100] A third aspect of the present invention provides an electronic device, comprising: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the aforementioned method.

[0101] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0102] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0103] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for kilometer-level hourly multi-element forecasting of upper-air meteorological fields, characterized in that, include: Acquire multi-source upper-air meteorological observation data and numerical model forecast output data, perform spatiotemporal alignment processing on the multi-source upper-air meteorological observation data and the numerical model forecast output data, construct a multi-dimensional feature space based on the physical constraint relationship of meteorological elements through the spatiotemporal aligned data, and identify the uncertainty distribution characteristics of each set member in the multi-dimensional feature space through error covariance analysis between set members, and generate a set of set samples with confidence level labels; The set of samples is input into the probability forecast generation model. The probability forecast generation model assigns corresponding dynamic weights to each set member based on the confidence index. At the same time, it estimates the probability density of the set sample with dynamic weights in combination with the spatial continuity constraint of the meteorological field, and obtains a probability forecast result containing multiple probability quantiles. Based on the spatiotemporal evolution characteristics of each meteorological element in the probabilistic forecast results, and combined with similar weather process cases in the historical database, the reliability of the probabilistic forecast results is assessed and the uncertainty is quantified, generating a probability distribution product containing the forecast value and its probability of occurrence, providing decision support information for forecast services.

2. The method according to claim 1, characterized in that, The multi-source upper-air meteorological observation data and the numerical model forecast output data are spatiotemporally aligned. A multidimensional feature space is constructed based on the physical constraints of meteorological elements using the spatiotemporally aligned data, including: The time dimension of the multi-source upper-air meteorological observation data and the numerical model forecast output data is standardized and mapped to a unified time reference. The spatial dimension is then resampled in a gridded manner, and the data with different spatial resolutions are projected onto a unified spatial grid to obtain spatiotemporally aligned data. Based on the spatiotemporally aligned data, a cross-element correlation feature matrix is ​​constructed according to the thermodynamic balance and dynamic conservation relationships between meteorological elements. The coupling constraint relationship between pressure gradient and wind field is used to characterize the physical consistency features between vertical levels. A multi-dimensional feature space is constructed by fusing the cross-element correlation feature matrix and the physical consistency features.

3. The method according to claim 1, characterized in that, In the multidimensional feature space, the uncertainty distribution characteristics of each set member are identified through error covariance analysis among set members, generating a set of set samples with confidence level labels, including: Based on the feature vectors of each set member in the multidimensional feature space, an error covariance matrix is ​​constructed between set members in each feature dimension. The error covariance matrix is ​​then optimized by hierarchical feature decomposition. The dominant feature vectors and corresponding feature value spectra obtained by the decomposition are used to construct an error subspace with a hierarchical structure. In the error subspace, the projection components of the error vector of each set member on each dominant eigenvector are calculated based on a nonlinear projection algorithm. The prediction deviation of the set members in the multidimensional feature space is analyzed in combination with the energy distribution characteristics of the eigenvalue spectrum, thus characterizing the uncertainty distribution characteristics of the set members. Based on the aforementioned uncertainty distribution characteristics, the nonlinear projection algorithm is used to diagnose the Mahalanobis distance between each set member and the center of the overall error distribution of the set. The condition number of the error covariance matrix is ​​then fused to reconstruct the Mahalanobis distance. The reconstructed distance values ​​are mapped to confidence indicators that characterize the reliability of the forecast, and finally, a set of set samples with confidence indicators is generated.

4. The method according to claim 3, characterized in that, The projection components of the error vector of each set member onto each dominant eigenvector are calculated based on a nonlinear projection algorithm. The prediction deviation of the set members in the multidimensional feature space is analyzed by combining the energy distribution characteristics of the eigenvalue spectrum. Obtain the error vector of each set member in the multidimensional feature space, and calculate the orthogonal projection components of the error vector on each dominant feature vector through a nonlinear projection algorithm. The nonlinear projection algorithm calculates the projection basis vector by introducing a random perturbation term that follows a Gaussian distribution, optimizes the projection basis vector to determine the optimal projection direction, and obtains the spatiotemporal evolution characteristics of the orthogonal projection components and the corresponding projection weight coefficients based on the optimal projection direction. Based on the spatiotemporal evolution characteristics of the orthogonal projection components and the projection weight coefficients, a weight calculation method incorporating variational entropy is constructed. The weight calculation method is then convolved with the spatiotemporal evolution characteristics, and the convolution result is adjusted using the projection weight coefficients to generate a distribution curve reflecting system stability. The spatiotemporal scale of the prediction deviation is identified by analyzing the local feature points of the distribution curve using the nonlinear projection algorithm, and the degree of prediction deviation of the set members in the multidimensional feature space is quantitatively evaluated.

5. The method according to claim 1, characterized in that, The ensemble sample set is input into the probabilistic forecast generation model. The model assigns corresponding dynamic weights to each ensemble member based on the confidence level identifier. Simultaneously, it estimates the probability density of the dynamically weighted ensemble sample set in conjunction with the spatial continuity constraint of the meteorological field, yielding a probabilistic forecast result containing multiple probability quantiles, including: The set of samples is input into the probability prediction generation model, which calculates the dynamic weight of each set member in the probability density estimation based on the confidence index of each set member. A weighted kernel density function is constructed based on the dynamic weights. The predicted values ​​of each set member are then fused using the dynamic weights according to the weighted kernel density function to form a preliminary probability density distribution. Based on the spatial continuity constraint of the meteorological field, the preliminary probability density distribution is spatially smoothed. By introducing a spatial correlation regularization term to constrain the difference in probability density distribution between adjacent grid points, the preliminary probability density distribution satisfies the physical continuity of meteorological field evolution in the spatial dimension, resulting in a spatially consistent probability density distribution. The probability density distribution is solved by quantile calculation, and the quantile values ​​corresponding to multiple preset probability levels are calculated. Based on the multiple probability quantile values, a probability prediction result containing different confidence intervals is constructed.

6. The method according to claim 5, characterized in that, Based on the spatial continuity constraint of the meteorological field, the preliminary probability density distribution is spatially smoothed. By introducing a spatial correlation regularization term to constrain the difference in probability density distribution between adjacent grid points, the preliminary probability density distribution satisfies the physical continuity of meteorological field evolution in the spatial dimension, resulting in a spatially consistent probability density distribution, including: Based on the spatial continuity constraint of the meteorological field, a spatial correlation regularization term is constructed. The spatial correlation regularization term calculates the probability density distribution gradient between the target grid point and its neighboring grid points, and sets the allowable range of gradient change based on the spatial diffusion characteristics of meteorological elements, thus obtaining a spatial constraint index with physical meaning. The spatial correlation regularization term and its corresponding spatial constraint index are introduced into the optimization objective function of the probability density distribution. By adding the penalty term corresponding to the spatial correlation regularization term to the optimization objective function, the difference in probability density distribution that violates the spatial continuity constraint is suppressed, thus forming an optimization objective that takes into account both physical constraints and numerical stability. The initial probability density distribution is iteratively adjusted based on the optimization objective. In each iteration, the deviation of the current probability density distribution from the spatial continuity constraint is calculated according to the optimization objective function. The probability density distribution values ​​of each grid point are adjusted by gradient descent so that the difference in probability density distribution between adjacent grid points gradually converges to the allowable range, thus obtaining a spatially consistent probability density distribution.

7. The method according to claim 1, characterized in that, Based on the spatiotemporal evolution characteristics of each meteorological element in the probabilistic forecast results, and combined with similar weather process cases in the historical database, the reliability of the probabilistic forecast results is assessed and the uncertainty is quantified, generating a probability distribution product containing the forecast value and its probability of occurrence, including: The evolution trend and spatial distribution of each meteorological element in the time series are extracted from the probability forecast results to construct a spatiotemporal evolution feature vector representing the current forecast weather process. The spatiotemporal evolution feature vector is matched with the feature vector of historical weather processes stored in the historical database. The distance between the feature vectors is calculated to filter out historical weather process cases that are similar to the current forecast weather process. Based on the deviation statistics between the actual observation results and the historical forecast results in the historical weather process cases, the reliability index of the forecast value corresponding to each probability quantile in the current probability forecast result is calculated, and the uncertainty level of each probability quantile is quantified by analyzing the distribution characteristics of the deviation statistics at different probability quantiles. The multiple probability quantiles, the forecast values ​​corresponding to each probability quantile, the reliability index, and the uncertainty level in the current probability forecast result are correlated and organized to generate a probability distribution product containing the forecast values ​​and their probabilities of occurrence.

8. A kilometer-level hourly multi-element forecasting system for generating upper-air meteorological fields, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to acquire multi-source upper-air meteorological observation data and numerical model forecast output data, perform spatiotemporal alignment processing on the multi-source upper-air meteorological observation data and the numerical model forecast output data, construct a multi-dimensional feature space based on the physical constraint relationship of meteorological elements through the spatiotemporal aligned data, and identify the uncertainty distribution characteristics of each set member through error covariance analysis among set members in the multi-dimensional feature space, and generate a set of set samples with confidence level labels. The second unit is used to input the set of aggregate samples into the probability forecast generation model. The probability forecast generation model assigns corresponding dynamic weights to each aggregate member based on the confidence index, and at the same time, combines the spatial continuity constraint of the meteorological field to estimate the probability density of the aggregate samples with dynamic weights to obtain a probability forecast result containing multiple probability quantile values. The third unit is used to assess the reliability and quantify the uncertainty of the probability forecast results based on the spatiotemporal evolution characteristics of each meteorological element in the probability forecast results and in combination with similar weather process cases in the historical database, and to generate a probability distribution product containing the forecast value and its probability of occurrence, so as to provide decision support information for forecast services.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images