A neural computing-based atmospheric environment prediction method

Abstract

The application discloses a kind of atmospheric environment prediction methods based on neural computing, it is related to meteorology technical field, including the following steps: S1, constructs meteorological physics prior tensor;S2, meteorological physics prior tensor is input improved GMFlow network, utilizes turbulence cascade and physical tensor constraint to construct correlation volume, based on wind field guidance non-local aggregation and logarithm-polar coordinate transformation iterative refinement, output motion vector field;S3, constructs mass conservation physical deviation matrix;S4, outputs physical consistency optical flow field;S5, generates dynamic source increment concentration chart;S6, outputs double physical drive prediction frame;S7, output final atmospheric environment prediction result.The application overcomes the limitation of traditional method physical mechanism missing, high frequency detail fuzzy and transport simulation distortion, provides efficient solution for accurate prediction under complex atmospheric environment.

Description

Technical Field

[0001] This invention relates to the field of meteorological technology, and in particular to an atmospheric environment prediction method based on neural computation. Background Technology

[0002] With the continuous improvement of the spatiotemporal resolution of atmospheric environmental monitoring data, traditional numerical forecasting models and classical optical flow algorithms face severe computational challenges in nonlinear motion estimation and refined prediction of massive convection data. Existing deep learning-based motion estimation networks, such as the GMFlow and RAFT series models, while improving the ability to capture pixel-level displacements through hierarchical feature matching, mainly rely on the similarity measurement of image appearance features for correlation construction. This method, which is based solely on appearance similarity, ignores the complex physical dynamic mechanisms implicit in atmospheric fluids (such as continuity constraints, energy cascade transfer, etc.) and deep topographic-induced effects, leading to distortion of the physical structure when constructing the motion vector field, thus limiting the accuracy of atmospheric environmental prediction under severe weather evolution conditions. In addition, classical data-driven methods often fail to strictly follow the law of conservation of mass in fluid dynamics to constrain transport processes, resulting in non-physical material sources and sinks when dealing with pollution diffusion and turbulent mixing simulations, increasing the uncertainty of environmental prediction results.

[0003] Therefore, how to provide a neural computing-based method for atmospheric environment prediction is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0004] This invention proposes a neural computation-based atmospheric environment prediction method. It improves the GMFlow network by inputting meteorological and physical prior tensors, constructs a correlation volume using atmospheric turbulence cascade transfer mechanisms and divergence and curl physical tensor constraints, and iteratively refines the preliminary motion vector field based on the nonlocal adjacency matrix of the wind field UV components and logarithmic-polar coordinate transformation. The divergence is calculated using the fluid dynamics continuity equation, and a mapping is established between the divergence and mass conservation conditions to construct a mass conservation physical deviation matrix. The associated potential energy field is calculated based on this matrix, and the negative gradient direction of the associated potential energy field is used to backpropagate and correct the motion vector field, outputting a physically consistent optical flow field. Based on the physically consistent optical flow field, observation frames are reconstructed, incompressible residual feature maps are extracted, and mass conservation source terms are inverted to generate a dynamic source increment concentration map. This dynamic source increment concentration map is used as the pollution input source, and a grid deformation field is constructed in conjunction with the physically consistent optical flow field to perform Lagrange advection tracing. Furthermore, Eulerian space deformation modulation is performed using local divergence tensors to simulate diffusion and dilution effects, outputting a dual-physics-driven prediction frame. The dual-physics-driven prediction frames are transformed to the frequency domain, and the deviation between the energy spectral density of high-frequency components and the theoretical spectrum of atmospheric turbulence is calculated. Dynamic high-frequency compensation weights are then applied to the spatial domain to generate prediction frames, outputting the final atmospheric environment prediction results. This mechanism effectively eliminates physical biases in pure data-driven motion estimation by establishing a deep physical coupling path from "meteorological physical prior construction" to "motion vector physical correction" and then to "dual-physics-driven prediction." This ensures that the generated prediction results dynamically follow the laws of mass conservation in fluid mechanics and turbulent cascade transmission, achieving the technical effect of improving the physical consistency and accuracy of atmospheric environment prediction while maintaining high spatiotemporal resolution. This invention overcomes the limitations of traditional methods, such as the lack of physical mechanisms, fuzziness of high-frequency details, and distortion in transport simulation, providing an efficient solution for accurate prediction under complex atmospheric environments.

[0005] An atmospheric environment prediction method based on neural computing according to an embodiment of the present invention specifically includes: S1. The ground monitoring station data and satellite cloud image are fused to generate the wind field tensor. The fluid potential energy is calculated based on the terrain elevation and encoded as a position bias term. The meteorological and physical prior tensor is then constructed by superimposing the data. S2. Input the meteorological and physical prior tensor into the improved GMFlow network, construct the correlation volume using turbulence cascade and physical tensor constraints, and iteratively refine it based on wind field-guided nonlocal aggregation and log-polar coordinate transformation to output the motion vector field. S3. Introduce the fluid dynamics continuity equation, calculate the divergence of the motion vector field, establish a mapping between the divergence and the mass conservation condition, and construct the mass conservation physical deviation matrix. S4. Construct an adjoint potential energy field based on the mass conservation physical deviation matrix, and use the negative gradient direction of the adjoint potential energy field to perform backpropagation correction on the motion vector field, outputting a physically consistent optical flow field. S5. Based on the physical consistency optical flow field, reconstruct the observation frame, extract the incompressible residual feature map and invert the mass-conserving source term to generate a dynamic source increment concentration map. S6. Introduce a dynamic source incremental concentration map as a pollution input source, construct a grid deformation field based on a physically consistent optical flow field to perform Lagrange advection tracing; use the local divergence tensor to perform Eulerian space deformation modulation to simulate the diffusion and dilution effect, and output a dual-physical-driven prediction frame. S7. Convert the dual-physics-driven prediction frame to the frequency domain, calculate the deviation between the energy spectral density of the high-frequency components and the theoretical spectrum of atmospheric turbulence, generate a prediction frame with dynamic high-frequency compensation weights applied to the spatial domain, and output the final atmospheric environment prediction result.

[0006] Optionally, S1 specifically includes: S11. Spatial interpolation is performed on the wind field data from the ground monitoring station to extract the optical flow field information from the satellite cloud image. The basic wind field tensor containing the UV components of the wind field is generated by pixel-level fusion. S12. Analyze the geographic elevation data of the study area, construct a gridded elevation matrix and calculate the fluid gravitational potential energy at the grid points; establish the conversion relationship between potential energy and kinetic potential energy based on Bernoulli's equation, and solve for the fluid kinetic potential energy bias; perform high-dimensional feature mapping on the fluid kinetic potential energy bias to generate a location encoding vector. S13. The location encoding vector is superimposed with the basic wind field tensor, and the terrain elevation constraint is introduced into the feature representation of the wind field tensor to generate a meteorological physical prior tensor containing wind field UV components and terrain location information.

[0007] Optionally, the improved GMFlow network includes a feature pyramid extraction layer, a group correlation calculation layer, a Transformer feature extraction layer, an iterative update layer, and a motion vector field generation layer: The feature pyramid extraction layer is used to extract multi-scale coefficients by convolution of two-dimensional wavelet basis functions and meteorological physical prior tensors, and to perform physical mapping by combining atmospheric turbulence cascade mechanism. The feature pyramid is then output by fusion of attention weight modulation and feature decomposition. The group correlation calculation layer is used to decouple the feature pyramid into advection and deformation groups, generate advection correlation volume and physical constraint correlation volume based on divergence and curl respectively, introduce wind field UV components to calculate inner product weights and fuse them to generate multidimensional correlation volume. The Transformer feature extraction layer is used to convert the multidimensional correlated volume into a point cloud, construct an adjacency matrix and physical field superpixel clusters based on wind field UV components, decouple the background and turbulence components, and decode them after aggregation through an attention mechanism to output the motion context feature tensor. The iterative update layer is used to calculate the divergence and curl field of the motion context feature tensor to generate a physical constraint mask. It combines the logarithmic-polar coordinate transformation features to perform mask weighting and fusion, and uses residual iteration to refine the motion feature map. The motion vector field generation layer is used to regress and upsample the motion feature map to output the final motion vector field.

[0008] Optionally, S3 specifically includes: S31. Construct a fluid dynamics continuity equation operator, input the motion vector field into the fluid dynamics continuity equation operator, and calculate the velocity gradient components along the coordinate axis direction at each grid point in space. S32. Receive the velocity gradient components along the coordinate axes of each grid point in the space, perform tensor summation on the velocity gradient components along each coordinate axis, and solve for the divergence value of the motion vector field in three-dimensional space; map the divergence value according to the spatial grid arrangement rules to generate a full-field divergence feature map. S33. Compare the whole field divergence feature map with the preset mass conservation zero threshold element by element, and calculate the deviation between the divergence value and the zero threshold. S34. Map the deviation back to spatial grid coordinates to generate a two-dimensional tensor characterizing the degree of fluid mass conservation deviation, and construct the mass conservation physical deviation matrix.

[0009] Optionally, S4 specifically includes: S41. Calculate the accompanying potential energy field, including receiving the mass conservation physical deviation matrix, initializing the mass conservation physical deviation matrix as Lagrange multiplier variables; applying Poisson equation constraints to the Lagrange multiplier variables, solving for the potential energy values ​​at each grid point; arranging the solved potential energy values ​​according to spatial location to obtain the accompanying potential energy field. S42. Obtain the physical correction force field, including calculating the gradient vector of the accompanying potential energy field at each grid point in space, and negating the gradient vector to obtain the physical correction force field. S43. Update the motion vector field, including superimposing the physical correction force field with the motion vector field, and using the superposition result to update the vector value of each pixel in the motion vector field. S44. Iteratively execute the operations of calculating the accompanying potential energy field, obtaining the physical correction force field, and updating the motion vector field until the motion vector field meets the preset physical convergence condition, and output a physically consistent optical flow field.

[0010] Optionally, S5 specifically includes: S51. Extract the observation frame from the previous moment as a reference benchmark, and use the motion vector in the physically consistent optical flow field to perform pixel-level spatial distortion transformation on the observation frame from the previous moment to reconstruct the predicted observation frame at the current moment. S52. Calculate the pixel-by-pixel difference between the reconstructed predicted observation frame and the actual observation frame at the current time to generate the original residual map representing the reconstruction error. S53. Construct an incompressible fluid constraint operator, perform a two-dimensional convolution operation between the incompressible fluid constraint operator and the original residual map in the spatial domain, and calculate the weighted response value in the neighborhood of each pixel; substitute the weighted response value into the incompressible fluid divergence constraint equation to calculate the divergence deviation of the weighted response value; set a divergence deviation threshold, compare the calculated divergence deviation with the threshold, and filter out pixels with divergence deviation less than the threshold; remove the weighted response values ​​corresponding to pixels with divergence deviation greater than or equal to the threshold, and retain the response features that satisfy fluid continuity; map the retained response features back to the pixel space to generate an incompressible residual feature map. S54. Construct a mass conservation source term inversion network, input the incompressible residual feature map into the mass conservation source term inversion network, and calculate the source term intensity distribution caused by mass change. S55. Normalize the source term intensity distribution and map it to the color space to generate a dynamic source increment concentration map that characterizes the dynamic mass source concentration distribution.

[0011] Optionally, S6 specifically includes: S61. Extract the dynamic source incremental concentration map, analyze the grid index coordinates and pixel gray values, calculate the pollution source intensity value according to the gray-intensity quantization mapping relationship, obtain the distribution feature data containing spatial location coordinates and pollution source intensity values, and generate the initial concentration field based on the distribution feature data. S62. Read the motion vector data in the physically consistent optical flow field, use the motion vector data to construct a geometric transformation matrix that maps the grid coordinates to the grid coordinates at the next moment, and generate the grid deformation field; S63. Based on the grid deformation field and the initial concentration field, the semi-Lagrange inverse tracking algorithm is used to calculate the inverse position of the pixel at the current time in the previous time. The pixel value at the inverse position is obtained by bilinear interpolation. Lagrange advection tracking is performed to generate an advection prediction intermediate frame. S64. Calculate the first-order partial derivatives of the advection prediction intermediate frame in the horizontal and vertical directions respectively, and construct a local divergence tensor based on the first-order partial derivatives; calculate the trace of the local divergence tensor, input the value of the trace into the exponential decay function to calculate the diffusion coefficient, and construct a diagonalized diffusion modulation matrix using the diffusion coefficient; perform an element-wise multiplication of the diffusion modulation matrix and the advection prediction intermediate frame using the Hadamard product operation to update the pixel values ​​of the advection prediction intermediate frame; apply a Gaussian kernel to perform convolution smoothing on the updated advection prediction intermediate frame to simulate the diffusion dilution effect in Eulerian space and generate a diffusion modulation prediction frame. S65. Perform pixel-level weighted fusion of the advection prediction intermediate frame and the diffusion modulation prediction frame to output the dual physical drive prediction frame.

[0012] Optionally, S7 specifically includes: S71. Perform a two-dimensional discrete Fourier transform on the dual-physics-driven prediction frame and divide the frequency domain coefficient matrix. S72. Calculate the energy spectral density and theoretical spectral density values ​​based on the frequency domain coefficient matrix, and output the energy spectral density deviation. S73. Perform normalization and inverse proportional function nonlinear mapping on the energy spectral density deviation, and output dynamic high-frequency compensation weights; S74. Perform an inverse transformation on the dynamic high-frequency compensation weights to map them back to the spatial domain to generate a spatial domain compensation weight map. Multiply the spatial domain compensation weight map with the dual physical drive prediction frame to obtain the prediction frame with enhanced high-frequency details. S75, output the prediction frame with enhanced high-frequency details as the final atmospheric environment prediction result.

[0013] The beneficial effects of this invention are: (1) This invention achieves accurate representation and decoupling of atmospheric dynamic features and topographic flow field constraints by constructing a deep coupling mechanism between meteorological physics prior tensors and an improved GMFlow network. Low-frequency approximation coefficients and high-frequency detail coefficients are extracted using two-dimensional wavelet basis functions and meteorological physics prior tensors convolution. Energy attenuation weights are constructed based on the atmospheric turbulence cascade transfer mechanism, physically mapping high-frequency details and low-frequency approximations to small-scale turbulence features and large-scale background flow field features, respectively. Combining grouped correlation calculation and the Transformer feature extraction module, the features are decoupled into advection feature groups and deformation feature groups. A correlation volume containing physical constraint dimensions is constructed using local divergence tensors, curl tensors, and meteorological prior wind field vectors. This mechanism transforms discrete wind field data into a feature space conforming to fluid physics laws, deeply capturing the dynamic characteristics of multi-scale flows and providing a high-quality feature representation that simultaneously possesses the accuracy of large-scale background flow fields and small-scale turbulence details.

[0014] (2) This invention establishes a physical closed-loop system from motion vector consistency verification to enhanced transport prediction accuracy by employing fluid dynamics continuity equation correction, Lagrange-Euler dual drive, and frequency domain turbulence spectrum compensation mechanism. The fluid dynamics continuity equation correction layer calculates the divergence of the motion vector field, constructs a mass conservation physical deviation matrix, and uses the negative gradient direction of the accompanying potential energy field to perform backpropagation correction on the motion vector field, ensuring that the motion vector field strictly follows the mass conservation law. The dual physical drive layer constructs a grid deformation field based on the physically consistent optical flow field to perform Lagrange advection tracking, and uses the local divergence tensor to perform Euler space deformation modulation to simulate the diffusion and dilution effect, accurately characterizing the transport path of pollutants. The frequency domain turbulence spectrum compensation layer calculates the deviation between the energy spectral density of high-frequency components and the theoretical spectrum of atmospheric turbulence, and generates a prediction frame in the spatial domain with dynamic high-frequency compensation weights. This system achieves full-process physical enhancement from motion estimation to detailed reconstruction through accompanying potential field gradient correction, Lagrange-Euler cooperative transport, and frequency domain energy spectrum constraints, ensuring that the final atmospheric environment prediction results have extremely high physical consistency and spatiotemporal resolution. Attached Figure Description

[0015] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is an overall flowchart of an atmospheric environment prediction method based on neural computing proposed in this invention; Figure 2 This is a flowchart illustrating the working principle of the improved GMFlow network, a neural computation-based atmospheric environment prediction method proposed in this invention. Detailed Implementation

[0016] The invention will now be described in further detail with reference to the accompanying drawings. These drawings are simplified schematic diagrams, illustrating only the basic structure of the invention, and therefore only show the components relevant to the invention.

[0017] refer to Figure 1 and Figure 2 A neural computation-based method for predicting atmospheric environment conditions, specifically including: S1. The ground monitoring station data and satellite cloud image are fused to generate the wind field tensor. The fluid potential energy is calculated based on the terrain elevation and encoded as a position bias term. The meteorological and physical prior tensor is then constructed by superimposing the data. S2. Input the meteorological and physical prior tensor into the improved GMFlow network, construct the correlation volume using turbulence cascade and physical tensor constraints, and iteratively refine it based on wind field-guided nonlocal aggregation and log-polar coordinate transformation to output the motion vector field. S3. Introduce the fluid dynamics continuity equation, calculate the divergence of the motion vector field, establish a mapping between the divergence and the mass conservation condition, and construct the mass conservation physical deviation matrix. S4. Construct an adjoint potential energy field based on the mass conservation physical deviation matrix, and use the negative gradient direction of the adjoint potential energy field to perform backpropagation correction on the motion vector field, outputting a physically consistent optical flow field. S5. Based on the physical consistency optical flow field, reconstruct the observation frame, extract the incompressible residual feature map and invert the mass-conserving source term to generate a dynamic source increment concentration map. S6. Introduce a dynamic source incremental concentration map as a pollution input source, construct a grid deformation field based on a physically consistent optical flow field to perform Lagrange advection tracing; use the local divergence tensor to perform Eulerian space deformation modulation to simulate the diffusion and dilution effect, and output a dual-physical-driven prediction frame. S7. Convert the dual-physics-driven prediction frame to the frequency domain, calculate the deviation between the energy spectral density of the high-frequency components and the theoretical spectrum of atmospheric turbulence, generate a prediction frame with dynamic high-frequency compensation weights applied to the spatial domain, and output the final atmospheric environment prediction result.

[0018] In this embodiment, S1 specifically includes: S11. Read the wind speed and direction values ​​recorded by the ground monitoring stations, calculate the reciprocal of the distance from each grid point to the surrounding stations using the inverse distance weighting method, use the reciprocal of the distance as a weight to perform a weighted summation of the wind speed values ​​of the surrounding stations, divide the weighted summation result by the sum of all weights to obtain the interpolated wind speed value of the grid point, extract the brightness distribution information of the pixel blocks in the satellite cloud image, calculate the displacement vector of the pixel block in the previous and next frames using the block matching method, convert the displacement vector into horizontal and vertical velocity components, align and fuse the interpolated wind speed value and velocity components according to the grid coordinates to generate a basic wind field tensor containing the UV components of the wind field. S12. Read the terrain elevation values ​​from the geographic information system, construct an elevation matrix with the same resolution as the wind field grid, calculate the product of the elevation value of each grid point in the matrix and the gravitational acceleration constant 9.81 to obtain the fluid gravitational potential energy, select 1.225 as the air density constant, divide the fluid gravitational potential energy value by the air density constant to calculate the potential energy value per unit mass of fluid, divide the potential energy value by the constant 10000 for normalization, calculate the natural logarithm using the ratio of the constant 1 to the normalized value, use the value of the natural logarithm as a correction factor, and add the correction factor to the normalized value. The feature bias reflecting the influence of terrain height is calculated. The feature bias is then multiplied by the 128 column vectors of the first-layer transformation matrix and superimposed with the first-layer bias vector. The result is mapped to a 128-dimensional feature space through a hyperbolic tangent activation function. The activated value is then multiplied by the 64 column vectors of the second-layer transformation matrix and superimposed with the second-layer bias vector. This is then compressed to 64 dimensions again through a hyperbolic tangent activation function. Finally, the value is multiplied by the 32 column vectors of the third-layer transformation matrix and superimposed with the third-layer bias vector. The result is linearly projected to a 32-dimensional space and arranged to generate a position encoding vector. S13. Add the values ​​of the location encoding vector to the corresponding values ​​of the basic wind field tensor element by element, use the calculation results as new feature values, embed the physical constraints of terrain elevation into the data structure of the wind field tensor, and generate a meteorological physical prior tensor containing wind field UV components and terrain location information.

[0019] In this embodiment, the improved GMFlow network includes a feature pyramid extraction layer, a group correlation calculation layer, a Transformer feature extraction layer, an iterative update layer, and a motion vector field generation layer: The feature pyramid extraction layer reads the meteorological and physical prior tensor, initializes the first approximate convolution kernel with values ​​of 0.5 and 0.5, and initializes the second detail convolution kernel with values ​​of 0.5 and -0.5. The meteorological and physical prior tensor is then subjected to one-dimensional sliding convolution with both the first and second detail convolution kernels along the horizontal dimension to obtain horizontal approximate and horizontal detail feature maps. The horizontal approximate feature maps are then subjected to one-dimensional sliding convolution with both the first and second detail convolution kernels along the vertical dimension to obtain the final low-frequency approximation coefficients and vertical detail coefficients. The horizontal detail feature maps are then subjected to one-dimensional sliding convolution with the first approximate convolution kernel along the vertical dimension to obtain horizontal detail coefficients. Finally, the horizontal detail feature maps are subjected to one-dimensional sliding convolution with the second detail convolution kernel along the vertical dimension to obtain diagonal detail coefficients. The vertical detail coefficients and horizontal detail coefficients are then combined... Multiply the number and diagonal detail coefficients by a constant of 2.5, and multiply the low-frequency approximation coefficients by a constant of 0.5 to construct an attention matrix with the same size as the low-frequency approximation coefficients. Initialize all element values ​​in the matrix to 0.01. Add the weighted low-frequency approximation coefficients, vertical detail coefficients, horizontal detail coefficients, and diagonal detail coefficients along the channel dimension. Replace the values ​​in the attention matrix with the sum. Perform element-wise multiplication of the attention matrix with the weighted low-frequency approximation coefficients, vertical detail coefficients, horizontal detail coefficients, and diagonal detail coefficients to obtain spatially aligned dual-scale features. Perform singular value decomposition on the spatially aligned low-frequency component feature matrix along the length and width dimensions, retain the 32 largest singular values, set the remaining singular values ​​to zero, and reconstruct the low-frequency component feature matrix using the zeroed singular value matrix to output the feature pyramid. The group correlation calculation layer reads the current frame feature map and the reference frame feature map from the feature pyramid. It divides the total number of channels in the feature map by 2 and rounds down to obtain the group index. The channel data before the index value is extracted as the advection feature group, and the channel data after the index value is extracted as the deformation feature group. The pixel feature vector at coordinates 0, 0 in the current frame advection feature group is multiplied by the corresponding position value of the pixel feature vector at the same coordinates in the reference frame advection feature group, and the multipliers are summed. The dot product value is calculated by traversing all pixel coordinates to generate the advection correlation volume. In the deformation feature group, a square neighborhood window with a side length of 3 is constructed centered on the pixel at coordinates 0, 0. The difference between the right and left pixel feature values ​​within the window is calculated to obtain the horizontal difference, and the difference between the lower and upper pixel feature values ​​is calculated to obtain the vertical difference. The divergence tensor value is obtained by adding the horizontal difference value and the vertical difference value. The curl tensor value is obtained by subtracting the horizontal difference value from the vertical difference value. The divergence and curl tensor values ​​are calculated by traversing all pixel coordinates. A two-dimensional tensor composed of divergence and curl values ​​is constructed. The Euclidean distance between the tensor of the current frame and the tensor of the reference frame is calculated to generate the physical constraint correlation volume. The meteorological prior wind field vector is read. The dot product value of the wind field vector of the current frame and the wind field vector of the reference frame is calculated. The dot product value is divided by the product of the vector magnitude to obtain the inner product weight of the wind field tensor. All channel values ​​of the advection correlation volume are multiplied by the inner product weight. All channel values ​​of the physical constraint correlation volume are multiplied by the difference between the constant 1 and the inner product weight. The weighted advection correlation volume value and the physical constraint correlation volume value are added to each other and fused to generate a correlation volume containing the physical constraint dimension. The Transformer feature extraction layer reads data points from the correlation volume, retaining the row, column, and channel indices of each data point. It then combines these indices to construct an unstructured point cloud dataset. Based on the row and column indices, it determines the two-dimensional pixel coordinates of each data point in the original image. The wind field UV components are read, and the Euclidean distance between each data point and the center point is calculated. Points with a distance less than 5 are identified as neighbor nodes, and their indexes are set to 1, while the rest are set to 0. A nonlocal adjacency matrix is ​​constructed. Finally, the feature vector of the center node is read and compared with that of its neighbor nodes. The feature vectors of a point are multiplied at corresponding positions to obtain the Hadamard product. The Hadamard products of all neighboring nodes are summed along the neighbor dimension to obtain the aggregated feature vector. The absolute values ​​of the differences between the UV components of the wind field in the horizontal and vertical directions are calculated and summed to obtain the wind field gradient. Regions with gradient values ​​less than 0.1 are grouped into the same cluster. The arithmetic mean of all feature values ​​within the cluster is calculated to obtain the mean vector. The difference vector between each feature vector and the mean vector within the cluster is calculated. The outer product matrix is ​​calculated by multiplying the difference vector with its transpose. All features within the cluster are then grouped together. The values ​​at corresponding positions in the outer product matrix are summed to generate a covariance matrix. The mean vector and the covariance matrix are concatenated along the channel dimension to generate a block-level feature vector. The first half of the channel data in the block-level feature vector is set as the background component, and the second half of the channel data is set as the turbulence component. The background component feature vector is multiplied by the key transformation matrix to generate 64 key vectors. The background component feature vector is multiplied by the value transformation matrix to generate 64 value vectors. The turbulence component feature vector is multiplied by the query transformation matrix to generate the query vector. The dot product value of the query vector and each key vector is calculated, and the dot product value is divided by... The value is scaled by a constant 8. The scaled value is then input into an exponential function to calculate the attention weight. The attention weight is used to sum the 64 value vectors to obtain a weighted block-level feature vector. The two-dimensional pixel coordinates of the weighted block-level feature vector in the point cloud are read. The coordinate values ​​are used to determine the four neighboring pixels in the pixel-level grid. The distance difference between the coordinates and the coordinates of the four neighboring pixels is calculated. The distance difference is used as the bilinear interpolation weight. The weighted block-level feature vector values ​​are multiplied by the corresponding weights and summed to map back to the pixel-level grid coordinates. The motion context feature tensor is then output. The iterative update layer reads the motion context feature tensor and the upsampled features of the multi-scale feature pyramid, concatenates the two tensors along the channel dimension, reads the motion vector field of the current iteration step, subtracts the horizontal value of the right-hand neighbor pixel from the horizontal value of all pixels in the motion vector field to obtain the horizontal first-order difference, subtracts the vertical value of the lower neighbor pixel from the vertical value of all pixels in the motion vector field to obtain the vertical first-order difference, adds the horizontal and vertical first-order differences to obtain the local divergence field, subtracts the vertical first-order difference from the horizontal first-order difference to obtain the local curl field, sets elements with values ​​greater than 0.8 in the local divergence and local curl fields to 1, and sets elements with values ​​less than or equal to 0.8 to 0, generating a physical constraint mask. The hidden layer state values ​​of the Transformer decoder are read, and multiplied by the corresponding values ​​in the physical constraint mask to obtain a result containing... The physical correction feature components are read from the feature map of the Transformer decoder. The Euclidean distance from the origin of the feature map grid coordinates to each pixel coordinate is calculated. The distance value is input into the logarithmic function to calculate the logarithmic radius. The angle from the origin of the grid coordinates to each pixel coordinate is calculated to obtain the polar angle. The coordinate transformation index constructed by the logarithmic radius and the polar angle is applied to the decoder feature map to map the decoder feature map to obtain the decoder polar coordinate features. The previously obtained motion vector logarithmic polar coordinate domain features and decoder polar coordinate features are concatenated in the channel dimension. Bilinear interpolation is used to map the concatenated features from the logarithmic polar coordinate system back to the Cartesian coordinate system to obtain the fused features in the Cartesian coordinate system. The fused features are dimensionality reduced using a convolutional layer to obtain the motion vector residual. The motion vector residual value is added element-wise with the current motion vector estimate value to output the refined motion feature map. The motion vector field generation layer reads the refined motion feature map and performs convolution operations on the feature map using a convolutional layer with a kernel size of 3, a stride of 1, and a padding pixel of 1. It outputs a regression feature map with 2 feature channels. The width and height values ​​of the regression feature map are enlarged to match the width and height values ​​of the original image using the nearest neighbor interpolation algorithm. All pixel position values ​​of the first channel in the regression feature map are extracted and set as horizontal displacement values. All pixel position values ​​of the second channel in the regression feature map are extracted and set as vertical displacement values. The horizontal displacement values ​​and vertical displacement values ​​are combined to generate the motion vector field.

[0020] The improved GMFlow network architecture proposed in this step is similar to the traditional GMFlow network architecture in that it is based on a pyramid feature extraction and iterative refinement mechanism. That is, it extracts semantic information at different resolutions by constructing a multi-scale feature pyramid, searches for matching relationships between features using the correlation volume, and uses the Transformer architecture to aggregate and update global features.

[0021] The difference lies in that this invention breaks through the limitations of traditional methods that rely solely on image appearance features for matching or ignore prior constraints of atmospheric physics. In the feature extraction layer, this invention adds a meteorological physics prior tensor input step, using two-dimensional wavelet basis functions to extract low-frequency approximation coefficients and high-frequency detail coefficients, rather than simply relying on image features. In the correlation calculation layer, features are decoupled into advection feature groups and deformation feature groups, calculating the advection dot product and physical divergence and curl tensors respectively, instead of relying solely on feature similarity. In the feature extraction and update layer, a nonlocal adjacency matrix is ​​constructed based on the UV components of the wind field, and block-level feature vectors are generated using the mean and covariance for attention calculation. Furthermore, in the iterative update layer, the divergence and curl of the motion vectors are used to generate a physical constraint mask and a logarithmic-polar coordinate transformation, rather than simple vector aggregation and residual update.

[0022] The beneficial effects of this improvement are that, by introducing a meteorological physical prior tensor and decoupling it from dual-scale physical features, the background flow field and turbulence features in atmospheric dynamics can be integrated into the optical flow estimation process. This breaks through the limitation of traditional methods that ignore physical conservation laws under complex meteorological conditions, leading to prediction bias, and achieves precise constraints from visual feature similarity to physical field dynamics. This design significantly enhances the model's ability to perceive non-rigid deformation of the wind field and turbulent motion, and can more accurately capture the nonlinear evolution of meteorological elements. Based on the transformation of the logarithmic-polar coordinate system and the correction of the physical constraint mask, the translation invariance and physical consistency of motion estimation are effectively improved, enhancing the robustness and computational accuracy of the system in complex meteorological motion prediction tasks.

[0023] In this embodiment, S3 specifically includes: S31. Construct a horizontal difference operator containing the value [1,-1] and a vertical difference operator containing the value [1,-1]. Read the motion vector field, use the horizontal difference operator to perform a horizontal sliding window calculation on the motion vector field, subtract the horizontal motion vector value of the adjacent pixel on the right side of the same row from the horizontal motion vector value of each pixel position to obtain the horizontal velocity gradient component. Use the vertical difference operator to perform a vertical sliding window calculation on the motion vector field, subtract the vertical motion vector value of the adjacent pixel below the same column from the vertical motion vector value of each pixel position to obtain the vertical velocity gradient component. Output the velocity gradient components of each grid point in space along the coordinate axis direction. S32. Read the horizontal and vertical velocity gradient components along the coordinate axes of each grid point in space. Take any grid point position, obtain the value of the horizontal velocity gradient component at that grid point position, obtain the value of the vertical velocity gradient component at that grid point position, perform an addition operation on the horizontal and vertical velocity gradient components at that grid point position, obtain an addition result value, assign the addition result value to that grid point position as the divergence value, traverse all grid point positions and repeat the above steps of obtaining, adding and assigning values, solve the divergence value of the motion vector field at each grid point position, fill the divergence values ​​of all grid point positions into the numerical matrix according to the grid row and column order, and generate the whole field divergence feature map. S33. Read the full-field divergence feature map, extract the divergence value of each pixel position in the feature map, subtract the constant 0 from the divergence value to obtain the difference, calculate the absolute value of the difference, and obtain the deviation between the divergence value and the zero threshold. S34. Fill the calculated deviation into the corresponding spatial grid coordinate position, construct a two-dimensional numerical matrix with the same size as the whole field divergence feature map, and generate the mass conservation physical deviation matrix.

[0024] In this embodiment, S4 specifically includes: S41. Read the mass conservation physical deviation matrix, use it as the source term of the right-hand side of the Poisson equation, initialize the Lagrange multiplier matrix, construct a discrete Laplacian operator convolution kernel containing values ​​of [0,1,0], [1,-4,1], and [0,1,0], and use the discrete Laplacian operator convolution kernel to perform a sliding window traversal on the Lagrange multiplier matrix to calculate the weighted sum of the values ​​of the current center pixel and its four neighboring pixels (up, down, left, and right) to obtain the Laplacian operator output value. Read the mass conservation physical deviation. The deviation value at the corresponding position in the matrix is ​​used to calculate the difference between the Laplacian operator output value and the deviation value to obtain the residual error. The Jacobian iteration step size is set to 0.25. The value of the current center pixel is subtracted from the product of the residual error and the iteration step size to obtain the updated Lagrange multiplier value. The updated Lagrange multiplier value is filled into the corresponding grid position. All grid points are traversed and the above steps of calculating the difference and updating the value are repeated until the values ​​of all grid points are updated, generating an adjoint potential energy field containing the potential energy values ​​of each grid point. S42. Read the associated potential energy field, construct a horizontal difference operator containing values ​​of [1,-1], use the horizontal difference operator to perform sliding window calculation on the associated potential energy field, subtract the potential energy value of the adjacent grid point on the right side of the same row from the potential energy value of each grid point to obtain the horizontal gradient value, construct a vertical difference operator containing values ​​of [1,-1], use the vertical difference operator to perform sliding window calculation on the associated potential energy field, subtract the potential energy value of the adjacent grid point below the same column from the potential energy value of each grid point to obtain the vertical gradient value, take the opposite of the horizontal gradient value and the vertical gradient value respectively to obtain the negative horizontal gradient value and the negative vertical gradient value, fill the negative horizontal gradient value and the negative vertical gradient value into the corresponding grid point position to generate the physical correction force field; S43. Read the physical correction force field and motion vector field. Take the position of the first pixel in the motion vector field. Read the horizontal and vertical motion vector values ​​of the pixel position. Read the horizontal correction force value and vertical correction force value of the same pixel position in the physical correction force field. Calculate the sum of the horizontal motion vector value and the horizontal correction force value to obtain the updated horizontal value. Calculate the sum of the vertical motion vector value and the vertical correction force value to obtain the updated vertical value. Write the updated horizontal value and the updated vertical value into the horizontal and vertical channels of the pixel position in the motion vector field, respectively, overwriting the original values. Take each pixel position in the motion vector field in the order of pixel index and repeat the above reading, calculation, summation and writing overwrite steps to complete the update of the motion vector field. S44. Record the vector values ​​of all pixels in the currently updated motion vector field. Calculate the sum of the squares of the horizontal and vertical values ​​of each pixel, and then take the square root of the sum to obtain the current pixel's magnitude value. Summate the magnitude values ​​of all pixels to obtain the current total magnitude value. Read the motion vector field output from the previous iteration, calculate the previous total magnitude value, and subtract the previous total magnitude value from the current total magnitude value to obtain the magnitude difference. Determine whether the absolute value of the magnitude difference is less than 0.0001. If the absolute value of the magnitude difference is less than 0.0001, stop the operation and output the current numerical arrangement result as the physically consistent optical flow field. If the absolute value of the magnitude difference is greater than or equal to 0.0001, re-execute the steps of calculating the accompanying potential energy field, obtaining the physical correction force field, and updating the motion vector field.

[0025] The mass conservation physical deviation-driven optimization process proposed in this step is similar to the traditional variational optical flow algorithm in that it is based on the theory of minimizing fluid dynamics energy. That is, by introducing physical constraint terms to modify the motion vector field, the optimal solution is solved by numerical calculation method of partial differential equations, and the optical flow estimation accuracy is gradually optimized by iterative update strategy.

[0026] The difference lies in that this invention breaks through the limitations of traditional methods that rely solely on the assumption of brightness conservation or ignore the divergence caused by the incompressibility of fluids. While traditional models directly use a weighted average of data and smoothing terms, this invention adds a mass conservation deviation detection step, calculating the divergence value of the motion vector field as a physical constraint deviation, rather than a single data error. In the physical correction force calculation step, the negative gradient direction is calculated using the accompanying potential field obtained from solving the Poisson equation, instead of simple gradient descent or neighborhood averaging. Finally, in the vector field update step, based on the superposition of the physical correction force field and the original motion vector field and the magnitude convergence criterion, a physically consistent optical flow field satisfying the divergence-free condition is jointly output, rather than a conventional optical flow field that only satisfies the condition of unchanged grayscale.

[0027] The beneficial effects of this improvement are that, by constructing the mass conservation deviation matrix and solving the Poisson equation in reverse, the physical law of fluid incompressibility can be integrated into the optical flow estimation process. This breaks through the limitations of traditional methods that cause prediction distortion due to violations of physical conservation in occluded or discontinuous motion regions, achieving an accurate mapping from apparent grayscale matching to intrinsic physical consistency. This design significantly enhances the model's ability to perceive the intrinsic physical laws of complex fluid motion and rigid body motion, and can more accurately eliminate non-physical noise and outliers in the motion vector field. The closed-loop iterative mechanism based on divergence convergence effectively improves the physical reliability of motion estimation and enhances the robustness and accuracy of the system in low-texture regions, motion occlusion boundaries, and fluid dynamic tracking tasks.

[0028] In this embodiment, S5 specifically includes: S51. Read the image data of the observation frame from the previous moment and the physically consistent optical flow field. Take the coordinates of any pixel in the physically consistent optical flow field and extract the horizontal and vertical motion vector values ​​at that coordinate position. Using the pixel at that coordinate position in the previous observation frame as the center, determine the corresponding target coordinate position at the current moment based on the horizontal and vertical motion vector values. Copy the pixel gray value at that coordinate position in the previous observation frame to the target coordinate position at the current moment. Traverse all pixel coordinates in the previous observation frame and repeat the above steps of extracting values, determining target coordinates, and copying gray values. Use the bilinear interpolation algorithm to process the gray values ​​at non-integer coordinate positions and generate the predicted observation frame at the current moment. S52. Read the predicted observation frame and the actual observation frame at the current time. Take a pixel at the same coordinate position in the two observation frames. Read the predicted gray value of the pixel in the predicted observation frame. Read the actual gray value of the pixel in the actual observation frame. Subtract the predicted gray value from the actual gray value to get the gray value difference. Determine the absolute value of the gray value difference as the residual value at the coordinate position. Traverse all coordinate positions and repeat the above steps of reading values, subtracting and taking the absolute value to generate the original residual map containing the residual values ​​of all coordinate positions. S53. Construct a 3x3 mean convolution kernel containing the values ​​[1,1,1], [1,1,1], [1,1,1] as a spatial smoothing operator. Use the mean convolution kernel to perform a sliding window traversal on the original residual map, sum the residual values ​​of the 9 pixel positions covered by the sliding window, divide the sum by 9 to obtain the weighted response value of the center pixel, calculate the first-order difference in the horizontal direction and the first-order difference in the vertical direction of the weighted response value of the center pixel, and then convert the first-order difference in the horizontal direction and the first-order difference in the vertical direction respectively. The first-order difference is squared and the results are summed to obtain the sum of squares. The square root of the sum of squares is then taken to obtain the magnitude of the gradient. The gradient magnitude threshold is set to 0.05. The magnitude of the gradient is checked to see if it is less than 0.05. If it is less than 0.05, the weighted response value of the pixel is retained. If it is greater than or equal to 0.05, the weighted response value of the pixel is set to zero. The processed weighted response value is filled into the corresponding grid coordinate position to generate an incompressible residual feature map. S54. Construct a fully convolutional neural network model with five convolutional layers as a quality-conserving source term inversion network. Input the numerical matrix of the incompressible residual feature map into the fully convolutional neural network model. Use the first convolutional kernel to extract the residual spatial texture features, use the second convolutional kernel to downsample the texture features, use the third convolutional kernel to calculate the deep semantic mapping of the features, use the fourth convolutional kernel to upsample the features to restore the resolution, and use the fifth convolutional kernel to output the source term value of each pixel. Traverse all pixels and repeat the above convolutional calculation steps to obtain the source term intensity distribution matrix caused by quality changes. S55. Read the source term intensity distribution matrix, calculate the maximum value of all source term values ​​in the matrix, divide the source term value of each pixel in the matrix by the maximum value to obtain the normalized value, multiply the normalized value by 255, round the multiplication result, map the rounded value to a grayscale value or pseudocolor value in the RGB color space, fill the mapped color value into the corresponding pixel position, and generate a dynamic source increment concentration map representing the dynamic quality source concentration distribution.

[0029] In this embodiment, S6 specifically includes: S61. Read the dynamic source incremental concentration map, extract the grid row and column coordinate index and the corresponding pixel gray value of each pixel in the image, read the preset gray-intensity quantization lookup table, use the pixel gray value to find the corresponding pollution source intensity value in the lookup table, map the grid row and column coordinate index to two-dimensional spatial location coordinates, construct a blank two-dimensional matrix with the same size as the dynamic source incremental concentration map, combine the two-dimensional spatial location coordinates and the corresponding pollution source intensity value to form distribution feature data, fill the pollution source intensity value in the distribution feature data into the corresponding spatial location coordinates in the blank two-dimensional matrix, and generate a pollution source term matrix as the initial concentration field. S62. Read the physically consistent optical flow field, extract the horizontal and vertical motion vector values ​​of each grid point in the optical flow field, construct a blank grid matrix with the same size as the observed grid at the current time, add the horizontal and vertical motion vector values ​​to the coordinate index of each grid point at the current time, calculate the grid point coordinate index corresponding to the next time, establish the correspondence between the coordinate index at the current time and the coordinate index at the next time, fill the correspondence into the blank grid matrix, and generate a grid deformation field containing grid point displacement information. S63. Read the concentration field of the previous moment, take any grid point in the concentration field of the previous moment as the pixel point of the current moment, extract the coordinate position of the pixel point, use the negative number of the motion vector value stored in the grid deformation field to calculate the inverse position coordinate of the pixel point in the previous moment, determine whether the inverse position coordinate is an integer, if it is an integer, directly extract the value of the coordinate in the concentration field of the previous moment, if it is not an integer, read the values ​​of the four adjacent pixels around the inverse position coordinate, use the bilinear interpolation algorithm to calculate the weighted average of the four values, assign the calculated value to the pixel point of the current moment, traverse all pixels and repeat the above steps of calculating the inverse position and interpolation value to generate the advection prediction intermediate frame; S64. Read the intermediate frame of the advection prediction, construct a horizontal difference operator and a vertical difference operator containing the values ​​[1,-1], use the horizontal difference operator to calculate the first-order partial derivative of the intermediate frame of the advection prediction in the horizontal direction, use the vertical difference operator to calculate the first-order partial derivative of the intermediate frame of the advection prediction in the vertical direction, calculate the sum of the squares of the first-order partial derivatives in the horizontal direction and the squares of the first-order partial derivatives in the vertical direction to obtain the square value of the magnitude of the gradient, add the constant 1 to the square value of the magnitude of the gradient to obtain the denominator value, divide the constant 1 by the denominator value to obtain the diffusion coefficient of the current pixel, use the diffusion coefficient to construct a diagonalized diffusion modulation matrix with the diagonal elements being the diffusion coefficient and the off-diagonal elements being 0, read the pixel value of the intermediate frame of the advection prediction, multiply the pixel value with the value at the corresponding position in the diagonalized diffusion modulation matrix, update the pixel value, construct a 3x3 Gaussian smoothing kernel to perform sliding window convolution calculation on the intermediate frame of the advection prediction after updating the value, and generate a diffusion modulation prediction frame. S65. Read the pollution source term matrix, the advection prediction intermediate frame, and the diffusion modulation prediction frame. Take a pixel at the same coordinate position in the three matrices. Read the source term value of the pixel in the pollution source term matrix. Read the first prediction value of the pixel in the advection prediction intermediate frame. Read the second prediction value of the pixel in the diffusion modulation prediction frame. Multiply the first prediction value and the second prediction value by a weighting coefficient of 0.5 to obtain the advection-diffusion fusion value. Add the advection-diffusion fusion value to the source term value to obtain the fusion value. Fill the fusion value into the pixel position. Traverse all pixel positions and repeat the above reading, multiplication, and addition steps to generate a dual physical drive prediction frame.

[0030] The dual physics-driven prediction process proposed in this step is similar to the traditional fluid optical flow prediction framework in that it is based on the fluid dynamics mass conservation theory. That is, it uses the motion vector field to describe the spatiotemporal evolution of the fluid medium, uses numerical difference algorithms to calculate the changes in physical quantities of the flow field, and uses spatiotemporal interpolation and filtering techniques to process nonlinear flow data.

[0031] The difference lies in that this invention breaks away from the limitations of traditional methods that rely solely on a single advection transport assumption or ignore the diffusion mechanism of substance concentration. Instead of directly using the semi-Lagrange algorithm for pixel displacement in traditional models, this invention adds a dynamic source increment inversion step. It calculates the pollution source intensity mapped from the incompressible residual feature map as the initial driving term for concentration field evolution, rather than directly using the observation image from the previous moment. In the diffusion modulation step, an adaptive diffusion coefficient constructed based on the inverse of the gradient magnitude is used to nonlinearly modulate the advection results, instead of simple fixed-weight smoothing. Finally, in the prediction frame generation step, a dual-physical-driven prediction frame is output based on the joint weighted fusion of the advection prediction intermediate frame and the diffusion-modulated prediction frame, rather than a single advection inference result.

[0032] The beneficial effects of this improvement are that, through dynamic source incremental inversion and gradient adaptive diffusion modulation, the spatial heterogeneity of the injection intensity and local concentration gradient of the fluid pollution source can be integrated into the prediction process. This breaks through the limitation of traditional methods that ignore source term changes and edge details when dealing with complex fluid diffusion, resulting in prediction ambiguity. It achieves an accurate mapping from pure motion advection deduction to coupled simulation of mass transport and diffusion. This design significantly enhances the model's ability to perceive the evolution of mass concentration distribution within the fluid and can more accurately capture the nonlinear mixing law of the fluid in a non-uniform field. The fusion processing based on the dual physical constraints of advection and diffusion effectively improves the physical realism and detail clarity of the prediction results, and enhances the robustness and accuracy of the system in complex environment fluid simulation and pollutant diffusion prediction tasks.

[0033] In this embodiment, S7 specifically includes: S71. Read the dual physical drive prediction frame, obtain the total number of row pixels and column pixels of the prediction frame, construct a complex matrix with the same size as the prediction frame, perform a two-dimensional discrete Fourier transform operation on the gray value of each pixel in the prediction frame, calculate the frequency components of the gray value in the horizontal and vertical directions, fill the calculated complex results into the corresponding positions of the complex matrix, generate a frequency domain coefficient matrix, set a low frequency cutoff radius of 20 pixels at the center of the frequency domain coefficient matrix, classify the frequency point values ​​less than 20 pixels from the center into the low frequency component coefficient matrix, and classify the frequency point values ​​greater than or equal to 20 pixels from the center into the high frequency component coefficient matrix. S72. Read the high-frequency component coefficient matrix, traverse each coordinate position in the matrix, extract the real and imaginary parts of the complex coefficients at that position, calculate the square of the real part, and calculate the square of the imaginary part. Add the squares of the real and imaginary parts to obtain the high-frequency energy value at that frequency point. Fill the high-frequency energy values ​​of all frequency points into the corresponding matrix coordinates to generate the energy spectral density distribution map of the high-frequency components. Obtain the row and column indices of the current frequency point in the matrix. Subtract the row index of the frequency domain center point from the row index to obtain the row distance, and subtract the column index of the frequency domain center point from the column index to obtain the column distance. Calculate the sum of the squares of the row distance and the squares of the column distance, and perform an operation on the sum. The square root operation is performed to obtain the radial distance value at the current frequency point. This radial distance value is multiplied by a preset frequency domain resolution coefficient to obtain the radial spatial frequency value at the current frequency point. The radial spatial frequency value is then multiplied by a preset turbulence internal scale parameter value to obtain the first product term. The square of the first product term is calculated, and 1 is added to the squared result to obtain the first denominator term. The fourth power of the radial spatial frequency value is multiplied by the cube power of the preset turbulence external scale parameter value to obtain the second product term. The constant value 0.033 is multiplied by a preset turbulence refractive index structure constant to obtain the third product term. The third product term is divided by the first denominator term to obtain the first intermediate result. The first intermediate result is divided by the second product term to calculate the theoretical spectral density value at the corresponding frequency point. S73. Read the energy spectral density distribution map and the theoretical spectral density value. Take a frequency point at the same coordinate position, read the energy spectral density value and the theoretical spectral density value at that point, subtract the theoretical spectral density value from the energy spectral density value and take the absolute value to obtain the spectral density difference. Traverse all frequency points and repeat the above subtraction and absolute value taking steps. Calculate the average value of all spectral density differences as the mean deviation. Calculate the sum of squares of the differences between all spectral density differences and the mean deviation. Divide the sum of squares by the total number of frequency points to obtain the variance. Use the variance value to normalize the spectral density difference of the current pixel. Substitute the normalized value into an inverse proportional fraction with a constant of 1 for calculation, that is, divide 1 by 1 plus the sum of the normalized values. Output the result as the dynamic high-frequency compensation weight. S74. Construct a frequency domain weight matrix composed of dynamic high-frequency compensation weights for all frequency points. Perform a two-dimensional discrete Fourier inverse transform operation on the frequency domain weight matrix to map the frequency domain weight matrix back to the image pixel space, generate a spatial domain compensation weight map, read the weight value of each pixel in the spatial domain compensation weight map, read the gray value of the same pixel in the dual physical drive prediction frame, multiply the weight value and the gray value, update the gray value of the pixel with the multiplication result, traverse all pixels and repeat the above steps of reading weight, multiplication and updating to obtain the prediction frame after high-frequency detail enhancement. S75. Mark the high-frequency detail-enhanced prediction frames as the final output of the atmospheric environment prediction model, and generate the final atmospheric environment prediction results to characterize the atmospheric environment state and evolution trend.

[0034] Example 1: To verify the effectiveness of this invention in predicting complex atmospheric environments, the method of this invention was applied to the refined atmospheric environment prediction system of a provincial meteorological bureau (hereinafter referred to as "Meteorological Bureau M"). In traditional meteorological environment prediction systems, dynamic calculations based on numerical models (such as WRF) or visual inference based on classical optical flow methods are typically employed. Numerical models suffer from long computation times and low spatial resolution when dealing with drastic local convective changes, while purely data-driven optical flow methods often lack fluid physics constraints, resulting in non-physical distortions in cloud motion estimation, leading to significant deviations in the prediction results regarding mass conservation and detailed texture. To address these issues, Meteorological Bureau M decided to adopt the atmospheric environment prediction method based on neural computation proposed in this invention.

[0035] During implementation, Meteorological Bureau M first acquired multi-source wind field observation data using automatic weather stations deployed throughout the province, while simultaneously acquiring high-resolution meteorological satellite cloud images. After preprocessing operations such as spatiotemporal registration, missing measurement interpolation, and radiometric calibration, a basic dataset containing dynamic and textural features was generated. At the same time, meteorological experts at Meteorological Bureau M performed precise flow field labeling and source term intensity inversion on historical typical weather processes, serving as a benchmark for model training and physical correction.

[0036] Meteorological Bureau M constructs a meteorological physics prior tensor, extracts low-frequency approximation coefficients and high-frequency detail coefficients using two-dimensional wavelet basis functions, and maps high-frequency details to small-scale turbulence features using an atmospheric turbulence cascade mechanism, while mapping low-frequency approximations to large-scale background flow field features. A feature pyramid is then output through attention weight modulation. Next, the features are decoupled into advection feature groups and deformation feature groups. Inner product weights are calculated using wind field UV components, and the advection correlation volume is fused with the physical constraint correlation volume based on divergence and curl to generate a multidimensional correlation volume. Subsequently, the correlation volume is converted into a point cloud, and a nonlocal adjacency matrix and physical field superpixel clustering are constructed based on wind field UV components. The background and turbulence components are decoupled, aggregated through an attention mechanism, and then decoded to output a motion context feature tensor.

[0037] In the core prediction and correction stage, this invention generates a physical constraint mask by iteratively updating the layer to calculate the divergence and curl field of the motion context feature tensor. It then performs residual iterative refinement using logarithmic-polar coordinate transformation features, outputting a preliminary motion vector field. Subsequently, the fluid dynamics continuity equation is introduced to calculate the divergence, constructing a mass-conserving physical deviation matrix. The adjoint potential field is obtained by solving the Poisson equation, and the motion vector field is corrected by backpropagation using the negative gradient direction, outputting a physically consistent optical flow field. Based on this optical flow field, the observation frame is reconstructed, the incompressible residual feature map is extracted, and the mass-conserving source term is inverted to generate a dynamic source increment concentration map. Finally, a dual-physical-driven prediction frame is output by combining Lagrange advection tracing and Eulerian space diffusion modulation. This frame is then transformed to the frequency domain, and dynamic high-frequency compensation weights are generated using the deviation between the high-frequency energy spectral density and the atmospheric turbulence theoretical spectrum, outputting the final atmospheric environment prediction result.

[0038] During implementation, the technical team at the Meteorological Bureau M discovered that, compared to traditional numerical model solving methods and conventional optical flow prediction methods, the method of this invention significantly improves the physical consistency and spatiotemporal resolution of atmospheric environment predictions. Traditional methods cannot effectively balance computational efficiency and physical conservation, and their predictions of small-scale turbulence details are often ambiguous. In contrast, the method of this invention, through coupling of meteorological and physical prior tensors, correction of the accompanying potential field, and compensation of the frequency domain turbulence spectrum, effectively achieves high-precision dynamic predictions that conform to the laws of fluid mechanics.

[0039] To further verify the actual performance of the method of the present invention, the meteorological bureau M conducted a detailed comparative test between the method of the present invention and the traditional method. The specific performance data is shown in Table 1: Table 1. Performance Comparison of Atmospheric Environment Prediction Method M by the Meteorological Bureau

[0040] As shown in Table 1, the performance of the atmospheric environment prediction system was comprehensively improved after applying the method of this invention. The estimation error of the motion vector field decreased from 4.5 pixels / frame in the traditional method to 0.8 pixels / frame, and the mass conservation divergence bias decreased from 15.2% to 2.1%, significantly improving the physical accuracy and reliability of the flow field estimation. The effective duration of short-term predictions was extended from 30 minutes to 65 minutes, and the cloud edge ambiguity decreased from 6.3 pixels to 1.5 pixels, significantly enhancing the timeliness and detail clarity of the system's predictions. In addition, the goodness of fit of pollution diffusion concentration increased from 81.5% to 95.4%, and the underreporting rate of extreme weather decreased from 18.5% to 3.2%, effectively improving disaster prevention and mitigation capabilities. The system operation and maintenance cost decreased from 2 million yuan / year to 1.2 million yuan / year, and the satisfaction rate of meteorological services also significantly improved, from 88.0% to 98.0%.

[0041] Through the method of this invention, the meteorological bureau successfully achieved high-precision physical perception and refined prediction of atmospheric environmental motion evolution. It effectively solved the problems of missing physical mechanisms and loss of high-frequency details in traditional methods, ensuring the accuracy of weather forecasts, greatly improving the intelligence and digitalization level of atmospheric environmental monitoring, significantly reducing the burden of manual correction for forecasters, enhancing the stability and robustness of the prediction system, and providing strong technical support for the construction of smart meteorology.

[0042] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for predicting atmospheric environment based on neural computing, characterized in that, Includes the following steps: S1. The ground monitoring station data and satellite cloud image are fused to generate the wind field tensor. The fluid potential energy is calculated based on the terrain elevation and encoded as a position bias term. The meteorological and physical prior tensor is then constructed by superimposing the data. S2. Input the meteorological and physical prior tensor into the improved GMFlow network, construct the correlation volume using turbulence cascade and physical tensor constraints, and iteratively refine it based on wind field-guided nonlocal aggregation and log-polar coordinate transformation to output the motion vector field. S3. Introduce the fluid dynamics continuity equation, calculate the divergence of the motion vector field, establish a mapping between the divergence and the mass conservation condition, and construct the mass conservation physical deviation matrix. S4. Construct an adjoint potential energy field based on the mass conservation physical deviation matrix, and use the negative gradient direction of the adjoint potential energy field to perform backpropagation correction on the motion vector field, outputting a physically consistent optical flow field. S5. Based on the physical consistency optical flow field, reconstruct the observation frame, extract the incompressible residual feature map and invert the mass-conserving source term to generate a dynamic source increment concentration map. S6. Introduce a dynamic source incremental concentration map as a pollution input source, construct a grid deformation field based on a physically consistent optical flow field to perform Lagrange advection tracing; use the local divergence tensor to perform Eulerian space deformation modulation to simulate the diffusion and dilution effect, and output a dual-physical-driven prediction frame. S7. Convert the dual-physics-driven prediction frame to the frequency domain, calculate the deviation between the energy spectral density of the high-frequency components and the theoretical spectrum of atmospheric turbulence, generate a prediction frame with dynamic high-frequency compensation weights applied to the spatial domain, and output the final atmospheric environment prediction result.

2. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, S1 specifically includes: S11. Spatial interpolation is performed on the wind field data from the ground monitoring station to extract the optical flow field information from the satellite cloud image. The basic wind field tensor containing the UV components of the wind field is generated by pixel-level fusion. S12. Analyze the geographic elevation data of the study area, construct a gridded elevation matrix and calculate the fluid gravitational potential energy at the grid points; establish the conversion relationship between potential energy and kinetic potential energy based on Bernoulli's equation, and solve for the fluid kinetic potential energy bias; perform high-dimensional feature mapping on the fluid kinetic potential energy bias to generate a location encoding vector. S13. The location encoding vector is superimposed with the basic wind field tensor, and the terrain elevation constraint is introduced into the feature representation of the wind field tensor to generate a meteorological physical prior tensor containing wind field UV components and terrain location information.

3. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, The improved GMFlow network includes a feature pyramid extraction layer, a group correlation calculation layer, a Transformer feature extraction layer, an iterative update layer, and a motion vector field generation layer. The feature pyramid extraction layer is used to extract multi-scale coefficients by convolution of two-dimensional wavelet basis functions and meteorological physical prior tensors, and to perform physical mapping by combining atmospheric turbulence cascade mechanism. The feature pyramid is then output by fusion of attention weight modulation and feature decomposition. The group correlation calculation layer is used to decouple the feature pyramid into advection and deformation groups, generate advection correlation volume and physical constraint correlation volume based on divergence and curl respectively, introduce wind field UV components to calculate inner product weights and fuse them to generate multidimensional correlation volume. The Transformer feature extraction layer is used to convert the multidimensional correlated volume into a point cloud, construct an adjacency matrix and physical field superpixel clusters based on wind field UV components, decouple the background and turbulence components, and decode them after aggregation through an attention mechanism to output the motion context feature tensor. The iterative update layer is used to calculate the divergence and curl field of the motion context feature tensor to generate a physical constraint mask. It combines the logarithmic-polar coordinate transformation features to perform mask weighting and fusion, and uses residual iteration to refine the motion feature map. The motion vector field generation layer is used to regress and upsample the motion feature map to output the final motion vector field.

4. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, S3 specifically includes: S31. Construct a fluid dynamics continuity equation operator, input the motion vector field into the fluid dynamics continuity equation operator, and calculate the velocity gradient components along the coordinate axis direction at each grid point in space. S32. Receive the velocity gradient components along the coordinate axes of each grid point in the space, perform tensor summation on the velocity gradient components along each coordinate axis, and solve for the divergence value of the motion vector field in three-dimensional space; map the divergence value according to the spatial grid arrangement rules to generate a full-field divergence feature map. S33. Compare the whole field divergence feature map with the preset mass conservation zero threshold element by element, and calculate the deviation between the divergence value and the zero threshold. S34. Map the deviation back to spatial grid coordinates to generate a two-dimensional tensor characterizing the degree of fluid mass conservation deviation, and construct the mass conservation physical deviation matrix.

5. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, S4 specifically includes: S41. Calculate the accompanying potential energy field, including receiving the mass conservation physical deviation matrix, initializing the mass conservation physical deviation matrix as Lagrange multiplier variables; applying Poisson equation constraints to the Lagrange multiplier variables, solving for the potential energy values ​​at each grid point; arranging the solved potential energy values ​​according to spatial location to obtain the accompanying potential energy field. S42. Obtain the physical correction force field, including calculating the gradient vector of the accompanying potential energy field at each grid point in space, and negating the gradient vector to obtain the physical correction force field. S43. Update the motion vector field, including superimposing the physical correction force field with the motion vector field, and using the superposition result to update the vector value of each pixel in the motion vector field. S44. Iteratively execute the operations of calculating the accompanying potential energy field, obtaining the physical correction force field, and updating the motion vector field until the motion vector field meets the preset physical convergence condition, and output a physically consistent optical flow field.

6. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, S5 specifically includes: S51. Extract the observation frame from the previous moment as a reference benchmark, and use the motion vector in the physically consistent optical flow field to perform pixel-level spatial distortion transformation on the observation frame from the previous moment to reconstruct the predicted observation frame at the current moment. S52. Calculate the pixel-by-pixel difference between the reconstructed predicted observation frame and the actual observation frame at the current time to generate the original residual map representing the reconstruction error. S53. Construct an incompressible fluid constraint operator, perform a two-dimensional convolution operation between the incompressible fluid constraint operator and the original residual map in the spatial domain, and calculate the weighted response value in the neighborhood of each pixel; substitute the weighted response value into the incompressible fluid divergence constraint equation to calculate the divergence deviation of the weighted response value; set a divergence deviation threshold, compare the calculated divergence deviation with the threshold, and filter out pixels with divergence deviation less than the threshold; remove the weighted response values ​​corresponding to pixels with divergence deviation greater than or equal to the threshold, and retain the response features that satisfy fluid continuity; map the retained response features back to the pixel space to generate an incompressible residual feature map. S54. Construct a mass conservation source term inversion network, input the incompressible residual feature map into the mass conservation source term inversion network, and calculate the source term intensity distribution caused by mass change. S55. Normalize the source term intensity distribution and map it to the color space to generate a dynamic source increment concentration map that characterizes the dynamic mass source concentration distribution.

7. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, S6 specifically includes: S61. Extract the dynamic source incremental concentration map, analyze the grid index coordinates and pixel gray values, calculate the pollution source intensity value according to the gray-intensity quantization mapping relationship, obtain the distribution feature data containing spatial location coordinates and pollution source intensity values, and generate the initial concentration field based on the distribution feature data. S62. Read the motion vector data in the physically consistent optical flow field, use the motion vector data to construct a geometric transformation matrix that maps the grid coordinates to the grid coordinates at the next moment, and generate the grid deformation field; S63. Based on the grid deformation field and the initial concentration field, the semi-Lagrange inverse tracking algorithm is used to calculate the inverse position of the pixel at the current time in the previous time. The pixel value at the inverse position is obtained by bilinear interpolation. Lagrange advection tracking is performed to generate an advection prediction intermediate frame. S64. Calculate the first-order partial derivatives of the advection prediction intermediate frame in the horizontal and vertical directions respectively, and construct a local divergence tensor based on the first-order partial derivatives; calculate the trace of the local divergence tensor, input the value of the trace into the exponential decay function to calculate the diffusion coefficient, and construct a diagonalized diffusion modulation matrix using the diffusion coefficient; perform an element-wise multiplication of the diffusion modulation matrix and the advection prediction intermediate frame using the Hadamard product operation to update the pixel values ​​of the advection prediction intermediate frame; apply a Gaussian kernel to perform convolution smoothing on the updated advection prediction intermediate frame to simulate the diffusion dilution effect in Eulerian space and generate a diffusion modulation prediction frame. S65. Perform pixel-level weighted fusion of the advection prediction intermediate frame and the diffusion modulation prediction frame to output the dual physical drive prediction frame.

8. The atmospheric environment prediction method based on neural computing according to claim 1, characterized in that, Specifically, S7 includes: S71. Perform a two-dimensional discrete Fourier transform on the dual-physics-driven prediction frame and divide the frequency domain coefficient matrix. S72. Calculate the energy spectral density and theoretical spectral density values ​​based on the frequency domain coefficient matrix, and output the energy spectral density deviation. S73. Perform normalization and inverse proportional function nonlinear mapping on the energy spectral density deviation, and output dynamic high-frequency compensation weights; S74. Perform an inverse transformation on the dynamic high-frequency compensation weights to map them back to the spatial domain to generate a spatial domain compensation weight map. Multiply the spatial domain compensation weight map with the dual physical drive prediction frame to obtain the prediction frame with enhanced high-frequency details. S75, output the prediction frame with enhanced high-frequency details as the final atmospheric environment prediction result.

Images