Meteorological large model prediction method based on data correction model

Abstract

The invention relates to the technical field of numerical weather forecast, in particular to a meteorological large model prediction method based on a data correction model, which comprises the following steps of: firstly acquiring multi-source atmospheric observation, screening observation by topology-optimal transmission quantum annealing, and constructing a weighted error covariance; applying mass, energy and earth rotation gradient, and generating a conservation assimilation field through diffusion implicit sampling; calculating a mutual information mask and coupling a cloud top optical flow fine tuning phase; cloud motion consistent field pulse codes are sent to the symplectic decomposition pulse neural network for neural form hardware reasoning, a pulse threshold is adjusted in a closed loop to control energy drift, and an uncertainty field is output through parallel disturbance reasoning. The method has the advantages of high resolution, low power consumption and probability prediction capability, and the extreme weather path and intensity prediction precision is obviously improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of numerical weather prediction, and in particular to a meteorological large model prediction method based on a data correction model. BACKGROUND

[0002] Meteorological large model prediction is shifting from traditional equation set integration to a "physics-data dual driving" framework. High-frequency occultation, radar, lightning, and stationary satellite brightness temperature observations provide rich constraints for rapidly evolving weather, but existing assimilation methods mainly rely on three-dimensional variation or ensemble Kalman filtering: observation quality control relies on semi-empirical thresholds, making it difficult to automatically identify spatial skeletons, and redundant or isolated observations often introduce noise; the cost function is usually quadratic, which only gives a local optimum for highly nonlinear systems; after assimilation, the initial value is inferred by a floating-point GPU network, which is power-hungry and accumulates energy drift, leading to instability in long-time integration; uncertainty is obtained through post-processing Monte Carlo sampling, which is not real-time.

[0003] These defects result in larger path and intensity errors in extreme convection and typhoon rapid intensification scenarios, and limit the timeliness of probabilistic forecast products. SUMMARY

[0004] To address the many problems existing in the prior art, the present application provides a meteorological large model prediction method based on a data correction model, which outputs a weighted covariance matrix based on topological-optimal transport quantum annealing observation screening, then obtains a conserved assimilated field through diffusion reverse integration driven by a physical-observation composite gradient, and then uses a symplectic decomposition pulse network to infer the master and perturbed initial values in parallel on a neuromorphic chip, real-time monitors energy and self-adapts thresholds, achieving high-precision, low-power, and probabilistic uncertainty prediction.

[0005] A meteorological large model prediction method based on a data correction model, comprising the following steps: Collecting multi-source atmospheric observation data and completing time alignment to form a time-aligned observation set; converting the initial field of the previous cycle into a coarse background field based on a graph structure neural weather model; Performing topological analysis on the time-aligned observation set, using an optimal transport-topology preserving binary optimization model and obtaining an observation selection vector through quantum annealing, constructing a weighted observation error covariance matrix, and obtaining an assimilated initial field through diffusion implicit sampling under the joint action of a mass conservation gradient, an energy conservation gradient, a geostrophic balance gradient, and a weighted observation gradient; Calculating the mutual information between the assimilated initial field and the geopotential height of the previous cycle to generate a mutual information mask, estimating a cloud top optical flow field, and injecting the gradient of the difference between the mutual information mask weighted geopotential height gradient and the cloud top optical flow field into the diffusion fine-tuning to obtain an initial field consistent with cloud motion; The initial field pulse consistent with cloud motion is encoded and input into the symplectic decomposition pulse neural network deployed on neuromorphic hardware for inference to generate a forecast field; the total energy of the system is monitored at preset intervals, and when the energy drift exceeds the threshold, the pulse firing threshold is adjusted and the topology preservation weight is simultaneously increased; the initial field is perturbed to generate an uncertainty field, and the forecast field, uncertainty field and energy drift index are output.

[0006] Preferably, before time alignment, the cloud pixel removal processing, the neighborhood consistency verification processing based on spatial first-order difference and the fixed time interval resampling processing are sequentially performed on the multi-source atmospheric observation data.

[0007] Preferably, the topology analysis adopts a sparse kernel persistent homology algorithm to calculate the connected persistence of each observation on a zero-dimensional topological scale and the ring-shaped persistence on a one-dimensional topological scale, and the two are weighted and summed to obtain the observation topological importance.

[0008] Preferably, the objective function of the optimal transmission-topology preservation binary optimization model is composed of a second-order Wasserstein distance item between the observation empirical distribution and the coarse background field distribution and an observation topological importance weighted penalty item, and the quantum annealing solves the objective function according to a decreasing simulated temperature sequence and outputs an observation selection vector.

[0009] Preferably, the diffusion implicit sampling adopts a numerical solution method of fixed step backward integration, each step simultaneously applies mass conservation gradient, energy conservation gradient and geostrophic balance gradient, and the observation gradient determined by the weighted observation error covariance matrix is used to correct the state variable.

[0010] Preferably, the mutual information mask calculates the mutual information between the potential height of the initial field after assimilation and the potential height of the previous cycle by using a kernel density estimation method, and selects the grids with mutual information greater than a preset threshold to form the mask.

[0011] Preferably, the cloud top optical flow field extracts multi-scale features of satellite infrared brightness temperature sequences by a space-time transformer model, and obtains the multi-scale features by simultaneous regression with the time-aligned lightning cluster flow observation.

[0012] Preferably, the convolution weights of the symplectic decomposition pulse neural network are first decomposed into regular terms and potential terms according to the symplectic mapping rule, then fixed-point quantization is performed, and then the pulse frequency and pulse polarity coding method are used to map to the pulse core of the neuromorphic hardware.

[0013] Preferably, during the pulse neural network inference process, when the total energy of the system drifts beyond the preset threshold, the pulse firing threshold is automatically increased, and the observation topological importance weighted penalty item weight in the objective function is simultaneously increased to maintain energy conservation.

[0014] Preferably, the uncertainty field is calculated by parallel inference of a fixed number of perturbed initial fields and then calculating the grid variance at each time.

[0015] Compared with the prior art, the application has the advantages and beneficial effects that: By topological sparse kernel analysis combined with optimal transmission quantum annealing screening observation, the effect of retaining spatial skeleton and suppressing redundant noise is realized. By physical conservation-observation gradient coupled diffusion implicit sampling, the effect of fast convergence of nonlinear state and energy conservation is realized. By symplectic mapping fixed-point quantization pulse convolution network in neuromorphic hardware inference, the effects of low power consumption and energy closed-loop control of high-resolution prediction are realized. By parallel perturbation inference, the effect of real-time output of grid-level uncertainty is realized. BRIEF DESCRIPTION OF DRAWINGS

[0016] Figure 1 is a flowchart of the method of the application; Figure 2 is a schematic diagram of energy closed-loop adaptive adjustment in the application. DETAILED DESCRIPTION

[0017] Hereinafter, embodiments of the present disclosure will be described with reference to the accompanying drawings. However, it should be understood that these descriptions are merely exemplary and are not intended to limit the scope of the present disclosure. In the following detailed description, numerous specific details are set forth in order to provide a thorough understanding of the embodiments of the present disclosure. However, it will be apparent to one skilled in the art that one or more embodiments can be practiced without these specific details. In addition, in the following description, descriptions of well-known structures and techniques have been omitted to avoid unnecessarily obscuring the concept of the present disclosure.

[0018] The terms used herein are merely used to describe specific embodiments and are not intended to limit the present disclosure. The terms "include", "comprise" and the like used herein indicate the presence of the described features, steps, operations and / or components, but do not exclude the presence or addition of one or more other features, steps, operations or components.

[0019] All terms used herein (including technical and scientific terms) have meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein should be interpreted to have meanings consistent with the context of the present specification, and should not be interpreted in an idealized or overly formal manner.

[0020] As Figure 1 shown, a meteorological large model prediction method based on a data correction model includes the following steps: Collect multi-source atmospheric observation data and complete time alignment to form a time-aligned observation set; convert the initial field of the previous cycle into a coarse background field based on a graph-structured neural weather model. In the overall process of this invention, the first step, "collecting multi-source atmospheric observation data and completing time alignment to form a time-aligned observation set; converting the initial field of the previous cycle into a coarse background field based on a graph-structured neural weather model," undertakes the dual responsibility of unifying the observation time and constructing the subsequent assimilation reference field.

[0021] In principle, multi-source observation data refers to a combination of ground-based, airborne, and satellite observations, each with its own spatiotemporal resolution and systematic error. To ensure that data from different sources can jointly constrain the assimilation field during subsequent calculations, time alignment must first be performed. This invention employs a linear interpolation method based on a reference time axis to map all observations to a unified time. Suppose that a certain observation record is at the original time. The corresponding observation error variance is After interpolation, the observations at the same time are obtained. and error variance .in Calculated through error propagation:

[0022] Indicates the interpolation time difference. This represents the first-order gradient of the observation in the time dimension. This correction ensures that the error is reasonably expanded over the time mapping process, eliminates the time asynchrony between observations, and provides a time-consistent "time-aligned observation set" for subsequent assimilation.

[0023] After time alignment is completed, a reliable and computationally efficient prior field needs to be constructed as the background for subsequent diffusion assimilation. This invention employs a graph-structured neural weather model to generate the coarse background field. The graph-structured neural weather model uses geographic grid points as nodes on a spherical grid, with edge weights between nodes set according to a distance function, and achieves spatial information exchange through a graph attention mechanism. Its core hidden layer uses a multi-head attention operator.

[0024] in , , These are query tensors, key tensors, and value tensors, respectively. The dimension of each attention head. Explicit distance weighting is included between nodes to maintain rotational consistency in the irregular spherical mesh. The model input is the initial 3D atmospheric field tensor from the previous loop, and the output is the forward-progressed time step. The coarse-resolution prediction field is obtained. To simultaneously cover information in both the horizontal and vertical directions, vertical layer position signals are embedded in the node features, and the output is reconstructed to a unified vertical coordinate through bilinear interpolation.

[0025] In terms of implementation, the multi-source observation includes five observation channels: occultation refractive index, radar reflectivity, infrared brightness temperature, lightning cluster flow, and airborne temperature, humidity, and wind. The system first performs time window filtering according to the channel's native time stamp, and then performs linear interpolation using the center time of the window as the target interpolation time. The interpolated records are stored in a unified format: {observation value, observation error variance, observation operator}. The observation operator describes the correspondence between the observation value and the assimilation grid; for example, radar reflectivity is mapped to the grid using three-dimensional bilinear coefficients, improving the efficiency of subsequent residual calculations.

[0026] In the coarse background field generation stage, the system reads the initial three-dimensional atmospheric field tensor stored in the previous loop as input, and outputs five variables—temperature, specific humidity, east-west wind, north-south wind, and air pressure—in the horizontal direction through one forward inference. and vertical direction The predicted values ​​on the layer are calculated in less than 20 seconds. This coarse background field does not require strict conservation at the energy closure level, but the potential height and wind field phase must be within an acceptable error range to provide good initial values ​​for subsequent diffusion assimilation.

[0027] After time alignment, the variance of the observation residuals converges to a uniform order of magnitude, reducing the weighting imbalance of multi-source observations and thus mitigating the ill-conditioned nature of the subsequent weighted observation error covariance matrix. Compared to the static reanalysis field, the coarse background field generated by the graph-structured neural weather model can reduce the mean square error of the initial geopotential height by approximately 10%, significantly accelerating the convergence speed of diffusion implicit sampling. The entire process ensures the spatiotemporal consistency between the observation information and the model's initial state, laying a reliable foundation for subsequent topological analysis, quantum annealing observation selection, and diffusion assimilation.

[0028] Preferably, before time alignment, the multi-source atmospheric observation data are sequentially processed by cloud pixel removal, neighborhood consistency check based on spatial first-order difference, and resampling at fixed time intervals.

[0029] This invention performs cloud pixel removal, neighborhood consistency checks, and fixed-time resampling on multi-source atmospheric observation data before assimilation. The purpose is to remove invalid pixels, balance spatial errors, and unify temporal resolution, so that subsequent topological analysis and diffusion assimilation are based on a quality-balanced input set.

[0030] The cloud pixel removal process relies on a dual judgment method: cloud detection flags and radiation thresholds built into the channel. First, cloud flag bits in the observation file are read, directly excluding pixels flagged as clouds. Then, for unmarked data, an infrared brightness temperature threshold method is used; if the brightness temperature is below a preset threshold, it is considered a cloud pixel and removed. The infrared threshold is determined by historical clearing field probability statistics, with different thresholds set for nighttime and daytime to account for the influence of sunlight. This step effectively removes pixels with strong attenuation from occultation detection cloud penetration observations and radar, significantly reducing artifact residuals.

[0031] The neighborhood consistency check is designed to address the spatial discretization error of the remaining pixels. Let grid points be used. The observed value is Its eight-neighbor mean is Calculate the first-order difference error:

[0032] in Representing grid points The absolute difference from its neighbors. If If a pixel exceeds the percentile threshold of its historical distribution within the same observation range, it is considered to have local inconsistency and is recorded as a suspicious pixel. Suspicious pixels are not directly deleted; instead, the neighborhood mean is used to determine their classification. Replacement is used to maintain spatial continuity. This avoids gradient spikes caused by point anomalies and also avoids information loss due to over-masking.

[0033] Fixed-interval resampling maps all channels to a unified time axis using linear interpolation. ,in The preset observation time resolution is used. Interpolation employs linear weights for both left and right frames, and time propagation correction is applied to the observation error variance: if the original observation is in and At that moment, insert The variance of the error at that time is:

[0034] in and As time weight, and This represents the original error variance. After resampling, all observations have consistent times and corresponding error estimates, and the subsequent covariance matrix can be directly assembled according to the observation order.

[0035] Example 1: The above process was applied to radar and occultation observations during a typhoon landfall. Approximately 12% of low-brightness-temperature pixels were removed during the cloud pixel removal stage; 3% of isolated outliers were corrected using a neighborhood consistency test; and resampling mapped observations from different times to a unified time axis every ten minutes. Results showed that the variance of the observation residuals before assimilation was reduced by approximately 20% compared to the original data, and the mean square error of the geopotential height of the coarse background field and the initial field after assimilation was reduced by approximately 8%, demonstrating that preprocessing has a significant effect on improving the prediction accuracy of the data correction model.

[0036] Topological analysis was performed on the time-aligned observation set. An optimal transport-topology-preserving binary optimization model was adopted and the observation selection vector was obtained through quantum annealing. A weighted observation error covariance matrix was constructed. Under the combined effect of the mass conservation gradient, energy conservation gradient, geostrophic balance gradient and weighted observation gradient, the assimilated initial field was obtained through diffusion implicit sampling. Given a time-aligned observation set, this invention completes the assimilation correction of the initial field through a three-stage cascaded process: topology analysis and observation selection, diffusion implicit sampling with consistent physical-observational multi-gradients, and adaptive update of error covariance. This process ensures the representativeness of the selected observations in the spatial topology and uses a diffusion generation mechanism to couple physical conservation conditions and observational constraints to the same optimization trajectory, thereby establishing a low-bias, high-smoothness, and convergently stable assimilated initial field.

[0037] The first stage involves topological analysis and observation selection. A sparse kernel persistent cohomology algorithm is used to calculate the persistent strip length of each observation at both the zero-dimensional and one-dimensional topological scales. The zero-dimensional strip reflects the persistence scale of isolated connected components, while the one-dimensional strip reflects the persistence scale of ring structures. Let the strip length of a certain observation at the zero-dimensional scale be... The strip length on a one-dimensional scale is Both are weighted by coefficients and Linear composition yields topological importance To reduce the need for adjusting human wave parameters, this invention will and Set as a constant related to the observation type, determined by historical cohomology statistics of occultation, radar, and infrared data.

[0038] Subsequently, an optimal transport-topology-preserving binary optimization model is constructed. The optimization objective consists of two terms: first, the second-order Wasserstein distance between the empirical observation distribution and the coarse background field distribution, used to measure whether the selected observations can cover the dominant mode of background error; second, a topology-preserving penalty term, used to maintain the previously calculated topological importance. The objective function is written as:

[0039] in Choose a vector for binary observations. Represents the second-order Wasserstein distance. Here are the topological constraint weight constants. This model is solved using quantum annealing hardware, leveraging quantum randomness to quickly escape local optima. The quantum annealing returns... The decision on whether an observation is adopted is made, and the removal of observations with high topological importance will result in a significant penalty, thereby driving the solution to maintain the representativeness of the spatial skeleton.

[0040] The second stage involves constructing the weighted observation error covariance matrix. Based on and original error variance ,definition:

[0041] in It is a very small positive number to prevent matrix singularities. This matrix directly participates in subsequent observation gradient calculations, reflecting the observation selection results.

[0042] The third stage involves multi-gradient consistent diffusion implicit sampling. First, noise is added to the coarse background field to obtain the diffusion starting point, then the diffusion stochastic differential equation is integrated along the reverse time direction. In each inverse step, four types of gradients are simultaneously accumulated: mass conservation gradient, energy conservation gradient, geostrophic equilibrium gradient, and the observation gradient derived from the weighted covariance matrix. Taking the energy conservation gradient as an example, let the total system energy be:

[0043] in air density, For constant volume specific heat capacity, For temperature, For the horizontal wind component, It is the acceleration due to gravity. This represents the potential height. The energy conservation gradient is achieved through... The partial derivative with respect to the state vector is obtained and combined with other physical gradients to generate the overall physical constraint term. The inverse integral uses an implicit Heun scheme, which can reduce the number of steps while maintaining numerical stability.

[0044] Example 2: Assimilation of a rapidly developing tropical cyclone. In the topology analysis stage, approximately 70% of the original observations were selected, reducing the second-order Wasserstein distance by about 15%. Implicit diffusion sampling converged in 50 inverse steps, reducing the iteration count by about one-third compared to unselected observations. After final assimilation, the initial field showed a reduction of approximately 13% in the mean square error of the 500 hPa geopotential height compared to the coarse background field, and subsequently reduced the cyclone path deviation by approximately 20 km in the 72-hour forecast, demonstrating the assimilation gain of this invention in extreme weather scenarios.

[0045] Preferably, the topology analysis uses a sparse kernel persistent cohomology algorithm to calculate the connectivity persistence of each observation at the zero-dimensional topology scale and the ring persistence at the one-dimensional topology scale. The topological importance of the observation is obtained by weighted summation of the two.

[0046] One of the core innovations of this invention is the "topology analysis-quantum annealing observation screening-multigradient diffusion implicit sampling" link, which is used to transform the time-aligned observation set into an assimilated initial field. This link balances the representativeness of the observation space, the consistency of physical conservation, and the stability of numerical convergence.

[0047] First, the topology analysis step uses a sparse kernel persistent cohomology algorithm to measure the stability of observations at low-dimensional cohomology scales. Specifically, a sparse kernel persistent cohomology algorithm is applied to the observation points in the geographic coordinate space with a radius of... A Gaussian kernel is used to record the generation and demise of topological structures as the kernel radius gradually increases from zero. Zero-dimensional homology focuses on connected components, reflecting whether an observation is isolated; one-dimensional homology focuses on ring structures, capturing large-scale morphologies such as typhoon eyes and frontal closures. For a single observation... The length of the zero-dimensional strip is denoted as The length of a one-dimensional strip is denoted as Therefore, topological importance is defined as follows:

[0048] in and where are constant weights, representing the relative contributions of zero-dimensional and one-dimensional topological scales to the assimilation effect, respectively. The weights are determined through offline sensitivity experiments: when Take 0.7, When the value is 0.3, it can balance the isolated observation rejection rate and the large-scale ring feature retention rate in global data assimilation system examples, thus obtaining stable benefits under most weather patterns.

[0049] After obtaining the topological importance, this invention incorporates it into the optimal transport-topology-preserving binary optimization model. An observational empirical distribution is then constructed. With coarse background field distribution Both are discretized to a unified histogram interval. The observation-selected vector... Expressing whether or not to adopt the observation The model objective function is written as:

[0050] in The second-order Wasserstein distance, weights The degree of topology preservation is controlled. When an observation with high importance is excluded, the penalty term is significantly increased, thereby driving the optimization result to retain key structures. This binary optimization problem can have tens of thousands of dimensions, and traditional simulated annealing is too time-consuming. To address this, this invention utilizes quantum annealing hardware to quickly search for feasible solutions within the energy decay sequence, outputting the optimal selection vector in an average time of approximately 2.5 milliseconds. .

[0051] according to Construct a weighted observation error covariance matrix. For the selected observations, the diagonal elements of the covariance are set as the original error variance. The rejected observations are approximated as having no information using the maximum diagonal element, mathematically expressed as:

[0052] in The value is a very small positive number, avoiding division by zero. The observed gradient can be directly calculated from this matrix.

[0053] The assimilation core uses diffuse implicit sampling. First, Gaussian noise is injected into the coarse background field to obtain the diffused terminal field. Then, during the reverse time integration, four gradients are accumulated at each step: the mass conservation gradient... Energy conservation gradient Geostrophic gradient and observation gradient The observation gradient is determined by the aforementioned covariance matrix:

[0054] For the observation operator, Let this be the current state vector. Let be the observation vector. The inverse integral uses a second-order implicit Heun scheme. This ensures stability while reducing the number of iterations. The entire sampling process typically converges in 50 steps.

[0055] Example 3: In a case study of a Meiyu front weather event, the original number of occultation refractive index observations was 18,000. After topological analysis, 13,000 observations were retained. Observations with zero maintenance duration significantly higher than one dimension were identified as isolated noise and removed. The quantum annealing solution took 2.3 milliseconds. After 50 steps of diffuse implicit sampling, the mean square error of the assimilated initial field at a 500 hPa geopotential height was reduced by 12% compared to the coarse background field. Simultaneously, the error at the front inflection point was reduced by approximately 15 km in the 24-hour path forecast, confirming the suppressive effect of the observation-diffusion coupling scheme on phase errors.

[0056] The combination of topological analysis and quantum annealing ensures that key observational structures are not filtered out; multi-gradient diffusion implicit sampling utilizes physical conservation and observational residuals to avoid the problems of over-smoothing or iterative divergence in traditional three-dimensional variational methods. The overall process reduces the initial mean square error by an average of about 10% in multiple playback experiments and provides significant improvements in dynamic-statistical joint evaluation indicators (energy closure error, path deviation, and precipitation peak location), which is the core supporting component of this invention.

[0057] Preferably, the objective function of the optimal transport-topology-preserving binary optimization model consists of a second-order Wasserstein distance term between the observation empirical distribution and the coarse background field distribution and an observation topology importance-weighted penalty term. Quantum annealing solves the objective function according to a decreasing simulation temperature sequence and outputs an observation selection vector.

[0058] In the assimilation link of this invention, the optimal transport-topology-preserving binary optimization model undertakes the key task of "filtering redundancy and preserving the skeleton from observations". The model principle is based on the mass conservation distribution matching idea and the sparse constraint of topological stability, coupling the two types of information into a single objective function, and then using a quantum annealing solver to output the observation selection vector in an extremely short time, thereby providing a high signal-to-noise ratio observation gradient for diffusion implicit sampling.

[0059] First, we need to construct the observation experience distribution. The empirical distribution is obtained by dividing the time-aligned observation set into histogram intervals according to the variable and normalizing the counts. The coarse background field is obtained by the same histogram division. The statistical difference between the two is measured by the second-order Wasserstein moment, whose discrete approximation is written as:

[0060] in and The first The center value of each histogram interval This represents the weight of the interval. To incorporate topological information, this invention utilizes a sparse kernel persistent cohomology algorithm to calculate the length of the connected strip at the zero-dimensional scale for each observation. Length of one-dimensional annular strip Then, according to experience weighting and Linear composition as topological importance When an observation forms a stable bridge or loop in a spatial network, it will have a high topological importance. The final objective function is denoted as:

[0061] vector For binary variables, elements Indicates adoption of observation Penalty coefficient The tradeoff between statistical matching and topology preservation is addressed. The objective function is encoded as a quantum Ising Hamiltonian and fed into a quantum annealing chip, where a decreasing temperature sequence is simulated internally. Gradually shrink the energy landscape and quickly search for the near-optimal. The annealing time is approximately 2.5 milliseconds, which meets the requirements for rapid cyclic assimilation.

[0062] When discretizing the objective function, the number of histogram intervals is adaptively determined based on the variable type; for example, 32 intervals are used for temperature and 24 intervals for wind speed, balancing resolution and computational cost. The weights for topological importance are determined through offline cross-validation. The quantum annealing process uses a 50-step temperature table, with the initial temperature set at 2.5°C and the final temperature at 0.1°C. The decay at each step follows an exponential law.

[0063] In obtaining Then, the covariance diagonal elements are arranged according to... Perform adaptive scaling. This represents the variance of the original observation error. for Small constants of the order of magnitude ensure that the matrix is ​​nonsingular. The updated covariance matrix is ​​directly used for observation gradient calculation, and the gradient formula is written as:

[0064] in To observe the operator matrix, Let this be the current state vector. This is the observation vector.

[0065] The experiments were conducted on five typhoons, three severe convective weather events, and two Meiyu fronts. Compared to random observation elimination without topological constraints, this invention reduces the second-order Wasserstein distance by an average of approximately 10%, while preserving key loop structures such as typhoon eyewall radar echo closures and frontal deflection lines. After quantum annealing screening, the assimilation convergence steps were reduced from 70 to 50; the computation time was shortened by two orders of magnitude compared to traditional simulated annealing. Ultimately, continuous improvements were observed in indicators such as 72-hour geopotential height forecast error, typhoon track deviation, and precipitation peak location error.

[0066] Example 4: During a typhoon in the Northwest Pacific, 18,000 original observations were collected, and the model retained 13,000. The retention rate of observations with zero maintenance greater than 0.8 and one maintenance greater than 0.4 exceeded 95%. After assimilation, the root mean square error of the 500 hPa geopotential height decreased from 18.2 meters to 15.9 meters, and the three-day cumulative path deviation of the typhoon center decreased from 78 kilometers to 63 kilometers. The results show that the optimal transport-topology preservation model effectively avoids the loss of spatial skeleton information while ensuring statistical consistency, providing higher-quality observation gradients for diffuse implicit sampling.

[0067] Preferably, the diffusion implicit sampling adopts a numerical solution method with fixed step-size inverse integration. At each step, the mass conservation gradient, energy conservation gradient and geostrophic balance gradient are applied simultaneously, and the state variables are corrected by the observation gradient determined by the weighted observation error covariance matrix.

[0068] Implicit diffusion sampling is the core component of the "data correction-physical conservation" closed loop in this invention. The idea is to treat the coarse background field as the final state of the diffusion process, gradually eliminating random noise through inverse integration, and simultaneously injecting the observation information and three types of physical conservation constraints into the same gradient field, thereby obtaining a self-consistent assimilated initial field. The implementation method can be divided into three parts: step size setting, composite gradient construction, and implicit update.

[0069] The step size is set using a fixed reverse time interval. Let the total number of steps in the diffusion process be... The time step for inverse integration is Define discrete time series ,in Corresponding to the maximum state of diffuse noise, This corresponds to the initial state after assimilation. A fixed step size facilitates batch vectorization in a hardware-accelerated environment and avoids convergence oscillations caused by adaptive step sizes in high-dimensional state spaces.

[0070] The composite gradient consists of four parts. First, the state vector is given. It includes variables such as temperature, specific humidity, horizontal wind speed, and air pressure. The physical gradient part includes the mass conservation gradient. Energy conservation gradient With the geostrophic gradient Taking the energy conservation gradient as an example, the total energy is defined as:

[0071] in air density, For constant volume specific heat capacity, For temperature, For the horizontal wind component, It is the acceleration due to gravity. This refers to the potential height. (Regarding...) about We can obtain the result by taking the partial derivative. The other two conserved gradients are similar. The observed gradient is determined by the weighted covariance matrix. With observation operator Together, we determine that the expression is:

[0072] in These are the observation vectors after observation and filtering. The diagonal element is ,in The variance of the observation error. To ensure numerical stability, the smallest positive number is used. The total gradient is formed by a linear combination of the four types of gradients. Weight , , , Cross-validation was used to determine that different gradients are comparable in terms of scale.

[0073] Implicit updates employ a second-order Heun back-integral scheme. Let the current state be... First, calculate the predicted gradient. To obtain the predicted state Subsequently Recalculate the gradient at the point And corrected by the average of the two. Iterate in reverse until The initial field after assimilation can then be obtained. The advantage of implicit schemes lies in their large stability region, which can be maintained within a fixed range. Relatively large for internal use To avoid numerical explosion.

[0074] Example 5: Testing a rapid squall line development process, setting... , After incorporating the physical gradient, the width of the water vapor convergence zone in the 6-hour forecast was reduced by approximately 12% compared to the control experiment without the conserved gradient, demonstrating smooth and consistent conservation characteristics. Incorporating the observation gradient reduced the peak position error of the simulated radar reflectivity by 11 kilometers. Overall, the diffusion implicit sampling method reduced the convergence steps by 35% compared to the three-dimensional variational method, while also lowering the energy closure error by approximately 0.4%, demonstrating the significant advantages of this invention in fusing observation and physical constraints.

[0075] The mutual information between the assimilated initial field and the potential height of the previous cycle is calculated to generate a mutual information mask. The cloud top optical flow field is estimated. The gradient of the difference between the weighted potential height gradient and the cloud top optical flow field is injected into the mutual information mask for diffusion fine-tuning to obtain an initial field consistent with cloud motion. The purpose of the synchronous phase correction step is to further reduce the positional error of the cloud system, ensuring that the geopotential height is consistent with the observed cloud motion direction, provided that the initial field after assimilation already satisfies the observation and physical conservation constraints. This step consists of three parts: mutual information mask generation, cloud top optical flow estimation, and diffusion fine-tuning. These three parts are tightly coupled and directly affect the accuracy and energy closure performance of subsequent symplectic form predictions.

[0076] Mutual information mask generation is used to identify areas where the geopotential height field generally exhibits a strong correlation with historical geopotential height fields, but rapid phase shifts often occur in areas of strong local convection or tropical cyclones. This invention quantifies the similarity between the "new" and "old" fields across grid cells using mutual information, selecting reliable regions as subsequent gradient weighting factors. Mutual information is defined as:

[0077] in Let be the joint probability density of the assimilated potential height and the potential height of the previous cycle. and This corresponds to the edge density. Density estimation uses the Gaussian kernel method, and the bandwidth is automatically determined by the silver standard. The mutual information threshold is set to 0.8 bits; grid cells exceeding the threshold are marked as valid, and their binary masks are denoted as [insert values ​​here]. This reflects the "potential structure stability". Using a mask instead of point-by-point numerical difference can suppress the influence of outliers on gradient amplification.

[0078] Cloud top optical flow estimation reflects the translational velocity of clouds in the upper troposphere and is directionally coupled with the geopotential height gradient under the hydrostatic equilibrium approximation. This invention utilizes a spacetime transformer model to characterize the infrared brightness temperature sequence of a geostationary satellite and lightning cluster flow, outputting a horizontal optical flow vector field. The model includes four attention phases and can capture multi-scale motion at a resolution of 4 kilometers. Lightning cluster flow is projected onto the brightness temperature grid through time alignment as an additional channel, improving the motion capture capability of rapidly evolving cloud towers at low brightness temperatures.

[0079] Mutual information mask and optical flow coupling gradient, potential height gradient With Cloudtop Light Flow Theoretically, the geostrophic approximation direction should be consistent. This is because the greater the geopotential height, the higher the isobaric surface, and the balance between the horizontal pressure gradient force and the Coriolis force causes the airflow to advance along the isopotential lines. If the two directions are inconsistent, it indicates a phase shift in geopotential height. This invention constructs a differential gradient:

[0080] in For element-wise multiplication, As a proportionality constant, the optical flow velocity dimension is linearly mapped to the height gradient dimension. The difference gradient is then subjected to a first partial derivative. Then, it is added to the existing physics-observation composite gradient to form a new gradient. .only The grid cells will be weighted to avoid over-correction of low mutual information regions.

[0081] Diffusion fine-tuning uses the assimilated initial field as the final state of the diffusion process, and completes the fine-tuning in 10 steps of back-integration. Each time step has a size of 0.1. Predicted state. With correction status The calculation method is the same as the aforementioned implicit Heun format, but the gradient is replaced with... Because of the addition of cloud top optical flow constraints, the fine-tuning is focused on the high cloud top motion region, and the average change is less than 2% of the overall field, which will not violate the previous conservation conditions.

[0082] Example 6: In a subtropical front, cloud top optical flow showed that the cold-season convection belt advanced rapidly from west to east, while the geopotential gradient of the initial field after coarse assimilation was biased southward. The mutual information mask marked the frontal region as 1 and the remaining background region as 0. After 10 fine-tuning steps, the geopotential line shifted northward by approximately 30 km, coinciding with the subsequent 3-hour evolution of satellite cloud top temperature, reducing the predicted frontal line position error by 12 km. The experiment shows that the optical flow-potential coupling gradient guided by the mutual information mask effectively corrects the phase error of the fast system without interfering with the statically stable backfield. Statistical analysis of all experimental cases showed that after adding this step, the 24-hour mean square error of the phase at 500 hPa altitude decreased by approximately 8% compared to the control experiment without fine-tuning, and the 72-hour precipitation peak position error decreased by approximately 10 km. Since the gradient update only takes effect in the high mutual information region, the energy conservation residual only increases by less than 0.05%, and its impact on the overall energy balance is negligible. In summary, this step significantly improves cloud-scale phase consistency, provides a more accurate initial cloud field for symplectic spiking neural networks, and thus enhances the forecasting performance of short-term severe convective paths and integrated precipitation.

[0083] Preferably, the mutual information mask uses the kernel density estimation method to calculate the mutual information between the initial field potential height after assimilation and the potential height of the previous cycle, and selects grids with mutual information greater than a preset threshold to form the mask.

[0084] Mutual information masks are used to screen for a "phase confidence region" between the assimilated initial field and the previous cycle's potential height field. The design principle is that only grids with statistically high correlation between the two time-series potential height fields are suitable as the region for cloud motion phase correction. Forcibly applying the optical flow-potential difference gradient in low mutual information regions can easily lead to overcorrection or even numerical oscillations. Therefore, this invention limits the gradient effect to a high signal-to-noise ratio subdomain by calculating mutual information at the grid scale and setting a threshold, thereby simultaneously ensuring local phase consistency and overall conservation balance.

[0085] Mutual information measures the difference between the joint and independent distributions of two random variables, capturing both linear and nonlinear correlations. Let random variables... This represents the grid value of the assimilated potential height field at the target time. This represents the corresponding grid value of the potential height field in the previous cycle. Mutual information is defined as... ,in For joint probability density, and The edge density is denoted as . If the two fields are locally translated or scaled, the mutual information will be significantly greater than zero; conversely, if the two fields have low correlation in the grid, the mutual information will approach zero.

[0086] Density estimation is performed using the Gaussian kernel method. For a single grid point... ,collect One historical assimilation-prediction sample pair Gaussian kernel density expression:

[0087] bandwidth According to Silver Standard Automatically determined, among which The standard deviation is the joint sample value. Since the geopotential height field is typically smoothed on a global structural scale, but local gradients are large in convective regions, to avoid excessive bandwidth that could smooth out details, the actual bandwidth is capped at 500 meters above the automatic value. The edge density is obtained by integrating the joint density along the direction of another variable.

[0088] The mutual information integral is performed using a two-dimensional Gauss-Legend de Gauss numerical integral with 64 integration nodes. The result is a mutual information lattice matrix. This invention sets a threshold of 0.8 bits, sets grid cells above the threshold to 1, and those below the threshold to 0, thus forming a binary mask. The threshold was derived from cross-validation of multiple cases: in 20 cases of tropical cyclones and severe convection, when the threshold was set to 0.8 bits, more than 35% of the low-correlation regions could be eliminated while maintaining 90% of the potential structure.

[0089] The estimation of the cloud top optical flow field relies on a spatiotemporal transformer model. The model takes as input a combination of six frames of geostationary satellite infrared brightness temperature and lightning cluster flow projections from the past 30 minutes, and outputs an optical flow vector. The resolution is 4 kilometers. The model employs four self-attention encoding layers, each containing eight attention heads, each with a dimension of 64. The output optical flow is upscaled to the same mesh as the potential height mask through bicubic interpolation.

[0090] In a high mutual information grid, the potential height gradient is... With mapped optical flow Difference calculation of gradient:

[0091] The constant scaling factor. The difference gradient is converted into a potential height correction term through a single spatial gradient operation. And add it to the original diffusion gradient, written as:

[0092] These are weighting coefficients. This correction only applies to... The grid plays a role in ensuring that the gradient in the low mutual information region maintains the original physical-observation balance.

[0093] The introduction of mutual information masks significantly reduced unnecessary optical flow driving. In a strong convective burst experiment, in the unmasked control scheme, the optical flow gradient mistakenly stretched parts of the clear sky, leading to local discontinuities in the geopotential height. The masked scheme, however, corrected this only in the high mutual information region, reducing the mean square angle between the geopotential height and the optical flow direction from 37 degrees to 21 degrees. The resulting downstream benefits are reflected in symplectic forecasts: the positional deviation of the 24-hour cloud top brightness temperature forecast decreased from 32 km to 25 km, while the increase in energy drift was controlled within 0.05%. This demonstrates that mutual information masks effectively achieve the goal of "local phase alignment and overall conservation and stability."

[0094] Example 7 uses a squall line event in North China that occurred at 12:00 on June 12, 2024. The mutual information between the assimilated initial geopotential height field and the previous cycle field was calculated, with the mask coverage area accounting for 58% of the total area. After five steps of fine-tuning using the mask superimposed with the optical flow difference gradient, the peak radar reflectivity simulation showed a 10-kilometer reduction in the 3-hour forecast peak area, and an improvement of the observation-simulation correlation coefficient of 0.06. This example verifies the effectiveness and universality of mutual information masks in rapidly developing weather conditions.

[0095] Preferably, the cloud top optical flow field is obtained by extracting multi-scale features of the satellite infrared brightness temperature sequence from the spatiotemporal converter model, and then combining the multi-scale features with time-aligned lightning cluster flow observations for regression.

[0096] The cloud top optical flow field is used to quantify the translational velocity of convective clouds in the upper troposphere and is a key input for subsequent potential-optical flow coupling correction. The spatiotemporal transformer model proposed in this invention can simultaneously resolve the temporal evolution and spatial texture of infrared brightness temperature within a single network, and enhance the deep convection signal by using the discharge location provided by lightning cluster flow observations, thereby achieving accurate estimation of cloud top translation and local upwelling.

[0097] In principle, the infrared brightness temperature sequence of a geostationary satellite can be viewed as a two-dimensional field that varies over time. If cloud microphysical evolution is ignored, pixel grayscale satisfies the brightness temperature order preservation assumption. Let the brightness temperature sequence be... ,in Represents two-dimensional coordinates on the Earth's surface. Representing time; let the optical flow vector field within the short time window be... Applying a Taylor expansion to the brightness temperature of consecutive frames and discarding higher-order terms yields the brightness temperature conservation constraint:

[0098] Traditional optical flow algorithms solve directly based on the brightness temperature gradient. However, there is a grayscale mismatch in the rapid evolution of high cloud tops. This invention uses a spatiotemporal transformer model to learn the laws of cloud morphological change from a broader perspective, and introduces lightning information to improve the feature discrimination of deep convection regions.

[0099] The spatiotemporal transformer model structure consists of a stacked four-level encoder. First, six consecutive frames of brightness temperature images are divided into 16×16 pixel blocks, and each block is convolved and embedded to obtain a local feature vector. Within each encoder stage, time-position encoding is introduced. Spatial location coding Add to form input Taking a single attention head as an example, its output is:

[0100] in , , These are the query, key and value matrix, and dimension. The value is 64. In the four-level stacking, lower-level blocks focus on capturing fast textures at a four-kilometer scale, while higher-level blocks focus on the outer envelope of cirrus clouds at a hundred-kilometer scale. To fuse lightning information, the time-aligned lightning cluster stream is projected onto a grid with the same resolution as the brightness temperature and encoded as a binary pulse map. The pulse map is processed through a one-dimensional convolution layer to obtain the electrical activity guidance vector, which is then concatenated with the brightness temperature feature in the second-stage encoder to explicitly prompt the network to pay attention to deep convection towers.

[0101] The model output optical flow estimation uses a linear regression head to map multi-scale features to displacement components. The training loss consists of three parts: brightness-temperature conservation loss. Gradient smoothing loss Losses consistent with lightning Brightness temperature conservation loss is based on the forward-backward brightness temperature difference:

[0102] in Gradient smoothing loss constrains the second-order difference of optical flow in space. Lightning consistency loss calculates the overlap error of lightning pulse maps before and after motion to highlight deep convection dynamics.

[0103] During the inference phase, the network receives six brightness temperature frames and one lightning pulse image. After outputting optical flow, it uses bicubic interpolation to upscale to the assimilated grid resolution. The output is then normalized to between 0 and 1 and then scaled using a scaling factor. It is converted into the dimension of potential gradient.

[0104] Implementation details: The model has 43 million parameters and uses a single-precision A100 GPU for inference, with each iteration taking approximately 0.9 seconds. Training employs the Adam optimizer with a base learning rate of 0.0003, converging on 1.2 million batches of data. The lightning pulse map time alignment error is controlled within 1 second to ensure accurate positioning of the rapid ascent tower.

[0105] Optical flow validation on 30 tropical cyclones and 20 continental convection systems showed that the model reduced the mean absolute error of wind vectors by approximately 25% compared to the traditional pyramid Lucas-Kanade method; and reduced the error in the lightning-weighted region by approximately 35%. In the subsequent assimilation link, the mean square value of the angle between potential height and optical flow direction decreased from 37 degrees to 22 degrees, significantly improving phase consistency.

[0106] Example 8: A case study of a Bay of Bengal cyclone on August 15, 2024. The cloud top optical flow output by the spatiotemporal transformer network forms a circular radial divergence near the cyclone eye. After combining with a mutual information mask, the differential gradient is applied only in the high mutual information region. The results show that the geopotential line shrinks by about 12 kilometers towards the cyclone center, the symplectic form forecast of the 24-hour cyclone path deviation is reduced by 15 kilometers, and the prediction error of the maximum wind speed area radius is reduced by about 6 kilometers.

[0107] like Figure 2 As shown, the initial field pulses consistent with cloud motion are encoded and input into the symplectic decomposition spiking neural network deployed on neuromorphic hardware to generate a forecast field; the total energy of the system is monitored at preset intervals, and when the energy drift exceeds the threshold, the pulse firing threshold is adjusted and the topology preservation weight is increased simultaneously; the initial field is perturbed in parallel to generate an uncertainty field, and the forecast field, uncertainty field and energy drift index are output.

[0108] After completing the cloud motion uniformity correction, the initial field is encoded into a pulse sequence and fed into neuromorphic hardware to perform symplectic decomposition spiking neural network inference. This step aims to obtain a high-resolution forecast field with low power consumption and to close the total energy and uncertainty in real time during the inference process.

[0109] The neuromorphic inference principle utilizes a symplectic decomposition spiking neural network, which, by explicitly preserving the Hamiltonian structure within the pulse domain, approximates a reversible numerical integrator. In this invention, the weights of each convolutional layer are first decomposed into regularized and potential components using a symplectic mapping, then quantized to an 8-bit fixed-point format and programmed into the neuromorphic chip. The input tensor is converted into a pulse firing frequency via amplitude mapping and into a pulse symbol via polarity mapping, forming a dual-channel pulse code. During inference, the node membrane potential is updated according to the symplectic transition rule, ensuring that the total energy drift within a finite number of steps is only a higher-order time step term. Compared to traditional floating-point inference, this structure inherently supports event-driven computation in hardware, triggering calculations only when a pulse arrives, significantly reducing power consumption.

[0110] Energy monitoring and threshold adaptation: To further suppress energy drift caused by quantization and finite step size, this invention calculates the total system energy every 6 hours during the inference process. The meanings of the symbols are consistent with those described above. If an energy drift relative to the previous segment exceeds a threshold, then... It automatically performs two adjustments: first, it increases the neuromorphic chip pulse firing threshold. First, reduce the neuron firing rate; second, simultaneously increase the topology preservation penalty coefficient. This strengthens the constraints of observational screening on the spatial framework. Two complementary adjustment levels are used: the pulse threshold affects local energy injection, while the penalty coefficient affects the large-scale error structure, ensuring that energy returns to the conservation zone in subsequent periods.

[0111] For uncertainty estimation, to address initial value sensitivity, this invention adds 64 Gaussian perturbations to the assimilated initial field, with a variance of 2% of the mean of the basic variables. All perturbation fields are inferred in parallel with the master field on the neuromorphic board, and the ensemble mean and variance are output every 3 hours. The uncertainty field is obtained by square rooting the ensemble variance tensor and is released along with the master forecast. Due to the native multi-core parallel architecture of the neuromorphic hardware, the latency of 64 inference sets is only 0.5 seconds higher than that of a single set, which can meet the real-time requirements of operational applications.

[0112] The neuromorphic hardware employs four second-generation chips, totaling 80 neural cores. Each chip stores approximately 120 megabits of pulse weights, consuming 42 watts. Pulse encoding uses a linear amplitude-frequency mapping: a temperature increment of 0.1 Kelvin corresponds to a 1 Hz pulse frequency increment, with sign mapping distinguished by positive and negative membrane potentials. Pulse decoding is performed on the host GPU, with a pulse counter calculating the firing rate using a 5-millisecond sliding window, and then restoring the floating-point value through linear scaling.

[0113] Example 9, in a subtropical depression case in May 2024, shows that the master inference time was 48 seconds, energy monitoring triggered two threshold adjustments, and the total energy drift peak remained within 0.4%. The ensemble uncertainty field was used to generate probabilistic precipitation products, and its 24-hour cumulative precipitation critical probability curve showed a 0.06 decrease in the Brier score compared to observational verification, demonstrating the positive effect of parallel perturbation inference on uncertainty quantification. Compared to the GPU floating-point baseline scheme (600 watts power consumption, 45 seconds), the neuromorphic scheme reduced energy consumption by approximately 92% while maintaining the same computation time, achieving significant energy savings.

[0114] The terminal prediction closed loop of this invention is composed of parallel neuromorphic symplectic decomposition pulse inference, energy closed-loop feedback, and ensemble perturbation. Based on a hardware-level conservation structure, this closed loop constrains long-term drift through a threshold-penalty dual adaptive mechanism, while simultaneously outputting spatial explicit uncertainty, providing reliable input for downstream typhoon path probability cones, extreme precipitation risk maps, etc. In summary, this step further solidifies the advantages of the data correction model in multi-scale, low-energy consumption, and high-stability meteorological large-scale prediction scenarios.

[0115] Preferably, the convolution weights of the symplectic decomposition spiking neural network are first decomposed into regularization terms and potential energy terms according to the symplectic mapping rule, then fixed-point quantization is performed, and the weights are mapped to the pulse kernels of the neuromorphic hardware through pulse frequency and pulse polarity encoding.

[0116] In the neuromorphic inference stage, to ensure compatibility between convolutional computation and the atmospheric conservation equations, this invention transforms the weights of deep convolutional networks into symplectic decomposition pulse structures. The term "symplectic decomposition" treats the original floating-point convolution kernel as a generative mapping of a discrete Hamiltonian system, approximating a continuous Hamiltonian flow with a reversible combination of regularization and potential energy terms, thus naturally preserving higher-order conservation of numerical energy during inference. If the symbol... The floating-point weight tensor of a certain convolutional layer is first expanded into a two-dimensional matrix along the channel dimension. According to the Lie group exponential mapping theory, we can write:

[0117] in It is the discrete time step. It is a standard symplectic matrix. and These represent the gradient operators for the regularization term and the potential energy term, respectively. Decomposed into regularized weights under the meaning of exponential mapping. With potential energy weights Then, the approximation is completed by using fast Taylor series truncation. It primarily captures odd-symmetric convolutional quantities, corresponding to momentum exchange. The main focus is on capturing even-symmetric convolutional integrals, corresponding to positional coupling. This decomposition can be minimized. Solving for the Frobenius norm, This represents the mapping of the aforementioned three exponents. The solution satisfies... and Consistent within the second-order step size error range.

[0118] To adapt to neuromorphic chips, this invention... 8-bit fixed-point quantization is used. The quantization interval is limited by observing the distribution of the absolute values ​​of the weights. ,in To balance quantization error and power consumption, sign and amplitude are encoded separately: the sign is stored using 1 bit, and the amplitude is stored using 7 bits, with the amplitude using a logarithmic scale to compress the dynamic range. Quantization bias is compensated for by fine-tuning the weight details through offline retraining. Quantized data is packetized into the neuromorphic hardware spiking kernel register, with each convolutional kernel corresponding to a set of synaptic arrays in hardware.

[0119] The state tensor is input to the chip via pulse encoding: on the one hand, an amplitude-frequency mapping is used, i.e., a temperature of 0.05 Kelvin corresponds to a 1 Hz pulse increment; on the other hand, the sign is characterized by positive and negative membrane potentials corresponding to pulse polarity. The chip's internal synaptic voltage is accumulated, and when the accumulated voltage exceeds the firing threshold... A pulse is emitted and the membrane potential is reset. The convolution multiplication-addition process is naturally completed by the convergence of pulses from the synaptic array, eliminating the need for floating-point multiplication. Due to the canonical-potential alternation structure of the symplectic mapping weights, the chip updates the membrane potential in segments under event-driven conditions, which is equivalent to performing a symplectic integrator at the physical level, thereby keeping the system energy drift at only a certain level within the prediction time window. .

[0120] To avoid cumulative drift caused by quantization and finite step size, this invention monitors the total energy of the system in real time. ,in It is air density, It is the specific heat capacity at constant volume. It is temperature, It is the horizontal wind component, It is gravitational acceleration, It refers to the potential height. When the energy drift relative to the previous check period exceeds a threshold... At that time, the pulse firing threshold is automatically increased. Simultaneously increase the topology preservation penalty coefficient in the observed and screened targets to mitigate the growth of energy injection and background error.

[0121] In terms of parallel perturbation inference, this invention deploys 64 quantized weight instances on the same chip for uncertainty estimation, each receiving an initial field with random perturbations. Since the pulse core array can share the weight table, replicating instances only requires additional storage of the membrane potential, resulting in limited increase in GPU memory overhead. The uncertainty field is obtained by summing the aggregate variance at the end of the inference process. In a measured 15-day forecast, the main control and 64 perturbations in parallel consume approximately 2000 joules of energy, representing over 90% energy savings compared to the 29 kilojoules of traditional GPU floating-point inference.

[0122] Example 10: In a case study of a tropical storm in the North Atlantic in October 2024, the symplectic decomposition pulse network completed a 15-day integration in 48 seconds, with a maximum energy drift of 0.4% and a path deviation reduced by 17 km compared to the floating-point UNet model at the same resolution. The 70% confidence ellipse provided by the ensemble uncertainty showed a 0.08 improvement in the agreement between the ensemble uncertainty and the optimal path error bound. This demonstrates the advantages of the decomposition-quantization-pulse coding joint strategy in balancing energy conservation, real-time performance, and probabilistic forecasting capabilities.

[0123] Preferably, during the inference process of the spiking neural network, when the total energy drift of the system exceeds a preset threshold, the pulse firing threshold is automatically increased, and the weight of the observation topology importance weighted penalty term in the objective function is increased simultaneously to maintain energy conservation.

[0124] The core idea of ​​the system's total energy closed-loop control is to treat the pulse neuromorphic inference as a quasi-Hamiltonian dynamical system. When quantization errors or finite step size errors cause the conserved quantity to drift, a hardware-software joint adjustment method is used to quickly bring the energy balance back. First, a fixed detection interval is set at the inference end, for example, every 6 hours, to calculate the total energy of the three-dimensional atmospheric state in real time. ,in It is air density, It is the specific heat capacity at constant volume. It is temperature, It is the horizontal wind component, It is gravitational acceleration, This refers to the potential height. The relative drift between two consecutive detection results is calculated:

[0125] like Exceeding the preset threshold This triggers the dual-channel compensation mechanism.

[0126] The hardware channel increases the pulse firing threshold. Suppress the firing frequency of neurons in the future step. Increase amplitude. Depends on the drift direction: If This indicates that the energy level is too high. Incremental; if The threshold adjustment is equivalent to reducing the time step in the numerical integrator, directly decreasing the rate of energy injection in the next stage of simulation. Since the pulse core operates on an event-driven basis in hardware, the threshold modification is a register write operation with a delay in the microsecond range.

[0127] Software channel synchronization increases the topology penalty coefficient in the optimal transmission-topology preservation objective function:

[0128] in This is a proportionality constant. The increase... In the next observation selection cycle, it will favor retaining observations with high topological importance, indirectly increasing the energy constraint strength of the space framework and making the initial field energy generated by subsequent assimilation closer to the real atmosphere. In actual implementation, The update is written to the quantum annealing control word and automatically invoked during hardware solving.

[0129] To avoid oscillations caused by frequent adjustments, this invention requires that if all four detection intervals are satisfied... Then it will gradually recover in the opposite direction. and The value drops to 90% of its original value, forming a hysteresis zone.

[0130] Example 11: During a tropical cyclone, Set to 0.5%. Detected in the 12th hour. The system will Increase by 5 millivolts, while The value was increased by 0.02. The subsequent energy curve returned to 0.3% within two intervals, and the 24-hour geopotential height deviation decreased by approximately 9 meters compared to the uncompensated experiment. These results demonstrate that hardware-software dual compensation effectively suppressed the energy explosion caused by quantization errors and maintained the phase accuracy of the dynamic field.

[0131] Comprehensive evaluation shows that after introducing the energy drift criterion and dual-channel compensation, the maximum system energy drift over 15 days decreased from 1.3% to 0.4%, while the average typhoon path deviation decreased by approximately 12 kilometers. The power consumption increase due to hardware threshold adjustment is less than 3%, far lower than the re-assimilation-re-inference cost caused by excessive energy drift. Therefore, the closed-loop strategy of this invention maintains physical consistency while considering real-time performance and energy consumption constraints, ensuring the long-term stable operation of the data correction model on the neuromorphic platform.

[0132] Preferably, the uncertainty field is obtained by calculating the grid variance at each time step after performing parallel inference on a fixed number of perturbation initial fields.

[0133] The uncertainty field is used to quantify the magnitude of initial value uncertainty propagated to each time grid during the forecasting process, and serves as a direct basis for risk warning and probabilistic products. This invention obtains the uncertainty field through a three-step process: "fixed-quantity perturbation of the initial field - neuromorphic parallel inference - temporal variance statistics," balancing real-time performance and statistical sufficiency.

[0134] First, after generating an initial field consistent with cloud motion, independent and identically distributed random perturbations are superimposed on it. The perturbations use spatially correlated Gaussian noise, and the covariance is constructed using the local spectral method: the historical error power spectra of variables such as temperature, wind speed, and humidity are fitted to... ,in For wave number, and Take an empirical constant. Perform an inverse fast Fourier transform on the power spectrum to obtain the spatial correlation function, and then generate a random field for each variable in the lattice domain. Let the number of perturbation sets be... , record The initial value after the perturbation is The reason for choosing 64 is twofold: firstly, a single neuromorphic chip can handle 64 parallel pulse nuclei without significantly increasing latency; secondly, according to the central limit theorem, The time variance estimation error is now below 5%.

[0135] The second step is to The initial field from the master control is pulse-encoded and fed into the symplectic decomposition pulse neural network. The chip internally shares a weight table, only replicating the membrane potential matrix, resulting in an overall power consumption increase of no more than 10%. The inference results are displayed at the predicted time. Output set field To save bus transmission bandwidth, the aggregated data is first averaged cell-by-cell on the chip side. Sum of squares of deviations Only these two matrices are sent back to the host via a high-speed channel.

[0136] The third step is to divide by on the host side. The lattice variance is obtained, which is the uncertainty field. :

[0137] In the formula It contains multiple variable components such as temperature, wind speed, and geopotential height, and can be projected as the probability of cumulonimbus cloud top brightness temperature or cumulative precipitation probability as needed. This represents a state vector, with components varying depending on the variable, such as temperature. Potential height Horizontal wind component wait; Represents the corresponding variance vector; It is the average of the set.

[0138] In 30 cases across three categories—typhoons, fronts, and severe convection—the 64-way parallel inference reduced energy consumption by 92% compared to the traditional double-precision ensemble mode, with an overall latency increase of only 0.5 seconds. For products with a probability of 24-hour cumulative precipitation > 50 mm, the reliability of the inspected grid improved from 0.54 to 0.62; the average coverage of the path probability cone increased by 7 percentage points, demonstrating the operational value of uncertainty estimation using parallel perturbation inference.

[0139] Example 12, July 2024 South China Rainstorm: The ensemble variance shows a strip-shaped high-value area along the Lingnan coast, corresponding to the uncertainty grid points of the master forecast precipitation peak; observational verification shows that the actual rain area is indeed within the high-variance area, demonstrating the risk indication capability of the variance field. Combined with the system energy closed loop, this invention can simultaneously output energy indicators, uncertainty probability, and deterministic forecasts, providing "one-stop" multi-scale information for severe weather.

[0140] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects.

[0141] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A weather macro-model prediction method based on a data correction model, characterized in that, The method comprises the following steps: Collecting multi-source atmospheric observation data and completing time alignment to form a time-aligned observation set; and converting an initial field of a previous cycle into a coarse background field based on a graph structure neural weather model; Performing topological analysis on the time-aligned observation set, obtaining an observation selection vector by using an optimal transport-topology preserving binary optimization model and quantum annealing, constructing a weighted observation error covariance matrix, and obtaining an initial field after assimilation by diffusion implicit sampling under the joint action of mass conservation gradient, energy conservation gradient, geostrophic balance gradient and weighted observation gradient; Calculating mutual information between the initial field after assimilation and the geopotential height of the previous cycle to generate a mutual information mask, estimating a cloud top optical flow field, and injecting a gradient difference between the mutual information mask weighted geopotential height gradient and the cloud top optical flow field into diffusion fine tuning to obtain an initial field consistent in cloud motion; Pulse coding the initial field consistent in cloud motion and inputting the same into a symplectic decomposition pulse neural network deployed on neuromorphic hardware for inference to generate a prediction field; monitoring total system energy at a preset interval, adjusting a pulse firing threshold when the energy drift exceeds a threshold, and simultaneously increasing a topological preservation weight; and generating an uncertainty field by parallel inference of a perturbed initial field, and outputting the prediction field, the uncertainty field and the energy drift index.

2. The method of claim 1, wherein, Before time alignment, sequentially performing cloud pixel removal processing, neighborhood consistency verification processing based on spatial first-order difference and fixed time interval resampling processing on the multi-source atmospheric observation data.

3. The method of claim 1, wherein, The topological analysis uses a sparse kernel persistent homology algorithm to calculate the connected persistence of each observation at a zero-dimensional topological scale and the ring persistence at a one-dimensional topological scale, and the topological importance of the observation is obtained by weighted summation of the two.

4. The method of claim 3, wherein, The objective function of the optimal transport-topology preserving binary optimization model is composed of a second-order Wasserstein distance term between the observation empirical distribution and the coarse background field distribution and an observation topological importance weighted penalty term, and the quantum annealing solves the objective function according to a decreasing simulated temperature sequence and outputs the observation selection vector.

5. The method of claim 1, wherein, The diffusion implicit sampling adopts a fixed step reverse integral numerical solution method, simultaneously applies mass conservation gradient, energy conservation gradient and geostrophic balance gradient at each step, and corrects the state variable by the observation gradient determined by the weighted observation error covariance matrix.

6. The method of claim 1, wherein, The mutual information mask calculates the mutual information between the geopotential height of the initial field after assimilation and the geopotential height of the previous cycle by using a kernel density estimation method, and selects the grids with mutual information greater than a preset threshold to form the mask.

7. The method of claim 1, wherein, The cloud top optical flow field extracts multi-scale features of a satellite infrared brightness temperature sequence by a space-time transformer model, and obtains the multi-scale features and time-aligned lightning cluster flow observations by simultaneous regression.

8. The method of claim 1, wherein, The convolution weights of the symplectic decomposition pulse neural network are first decomposed into regular terms and potential terms according to the symplectic mapping rule, then fixed-point quantization is performed, and the pulse frequency and pulse polarity coding method are used to map the weights to the pulse core of the neuromorphic hardware.

9. The method of claim 4, wherein, During the inference of the pulse neural network, when the total system energy drift exceeds a preset threshold, the pulse firing threshold is automatically increased, and the weight of the observation topological importance weighted penalty term in the objective function is simultaneously increased to maintain energy conservation.

10. The method of claim 1, wherein, The uncertainty field is calculated by performing parallel reasoning on a fixed number of perturbed initial fields and then calculating the grid variance at each time. The uncertainty field is calculated by performing parallel reasoning on a fixed number of perturbed initial fields and then calculating the grid variance at each time.

Images